Coupled groundwater flow and contaminant transport simulation in a confined aquifer using meshfree radial point collocation method (RPCM)

Engineering Analysis with Boundary Elements 66 (2016) 20–33

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Coupled groundwater flow and contaminant transport simulation in a confined aquifer using meshfree radial point collocation method (RPCM) L. Guneshwor Singh a, T.I. Eldho b,n, A. Vinod Kumar c a

Health Physics Division, Bhabha Atomic Research Centre, Trombay, Mumbai, India Civil Engineering Department, Indian Institute of Technology Bombay, Powai, India c Radiation Safety Systems Division, Bhabha Atomic Research Centre, Mumbai, India b

art ic l e i nf o

a b s t r a c t

Article history: Received 14 February 2015 Received in revised form 22 November 2015 Accepted 2 February 2016

In this study, a meshfree radial point collocation method is used to model the contaminant transport through confined aquifer. The discretization of the governing equations is done by a point collocation method and radial basis functions (RBF) are used as the interpolation function. For comparative study, two widely used radial basis functions namely multi-quadrics and exponential RBF are used. A local circular support domain is employed to construct the shape functions. In the model, no information on nodal inter-relationship is required for shape function construction except the nodal coordinates, unlike in finite-difference (FDM) or finite-element (FEM) based methods. The developed model is validated through benchmark problems in one and two dimensions. Further, application of the model for advective transport with high Peclet number has been studied and the model has been found to be effective in handling the instability of high Peclet problems. For the field problem considered, the results obtained from the model have been compared with the FEM solution and was found to be satisfactory. This method is relatively easy to implement and offers better accuracy with acceptable computational time. Considering the significant advantages offered by this method, it can serve as a good alternative to the conventional methods. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Groundwater flow Contaminant transport simulation Confined aquifer Radial basis functions Local support domain Shape parameter Point collocation method

1. Introduction Groundwater modeling has nowadays become a major part of most projects dealing with groundwater development, protection and remediation. With the improvements in computer hardware and software, the role of groundwater models will continue to increase. Most of the currently available codes or software for groundwater modeling is based on finite difference (FDM) and finite element (FEM) methods. In these models, one has to construct a grid or mesh before simulation which is tedious and time consuming. For many problems, the creation of mesh becomes the major component, being computationally expensive and also difficult to create specially for three-dimensional cases. Moreover, because of the need to create mesh, adaptive analysis of a problem is also difficult to carry out by mesh based methods as several remeshing or re-zoning are required during the analysis. Mesh based methods are also not suitable or lead to considerable difficulty in certain problems such as large deformation problems, crack n

Corresponding author. E-mail address: [email protected] (T.I. Eldho).

http://dx.doi.org/10.1016/j.enganabound.2016.02.001 0955-7997/& 2016 Elsevier Ltd. All rights reserved.

propagation, breakage of materials etc. Meshfree methods, on the other hand, do not suffer from the above shortcomings and offers itself as an alternative approach to the grid or mesh based techniques such as FDM and FEM. Because of the advantages it offers, it is increasingly becoming popular. Meshfree methods come in many varieties. Several of these approaches or types of meshfree methods have been discussed by Liu and Gu [15]. Radial basis functions (RBF) were originally proposed for scattered data interpolation [10] and was found to have exponential convergence [16]. It was Kansa [12] who first introduced a meshfree scheme for solving partial differential equations (PDE) using a globally supported multi-quadrics (MQ) RBF interpolation. In this method, a point collocation discretization of the governing PDEs, known as strong form approach was employed. Since then the technique was successfully applied for solving various two and three dimensional PDEs. But Kansa's globally supported RBF collocation method suffers from ill-conditioning and instability as the number of nodes increases. Kansa and Hon [13] examined several approaches to overcome these difficulties and improve accuracy such as block-partitioning, using matrix pre-conditioner, variable MQ shape parameters, multi-zone method, knot adaptivity, multi-level approximation etc. They concluded that each of the

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

above techniques improves the conditioning of the system matrix and the solution accuracy, and recommended a hybrid combination of the techniques to obtain stable and very accurate solutions. Another issue with the collocation technique is the instability and drastic loss of accuracy at the boundary when it involves singularity or derivative (Neumann) boundary conditions [3,33]. Fedoseyev et al. [5] proposed an improved Kansa–MQ method using PDE collocation on the boundaries which deal with this issue. It was shown that their approach leads to considerable improvement in the accuracy and substantially better conditioning of the resulting system matrix over the Kansa approach [12]. Gutierrez and Florez [9] forced the PDE and the normal derivative to be satisfied at the Neumann boundary points and observed that it yields very accurate solutions.This formulation is named double collocation approach. Liu and Gu [15] discussed different strategies of handling the derivative boundary conditions namely direct collocation, using fictitious points, Hermite-type collocation, method of regular grids, dense nodes and weak–strong form method. They compared the relative performances of each strategy. Based on the sample problems they studied, they found that though the accuracies of the methods are different, their convergence rates are nearly the same. Though no particular method was specifically recommended, they observed that the method of regular grids converged poorly while the Hermite collocation method worked well for both regular and irregular nodes. Another recently proposed collocation approach is the local RBF approach [14]. It uses compact sub-domains, called local support domains, for function interpolation. This is in contrast to the global approach where the entire domain is used for interpolation. This approach can circumvent the ill-conditioning problem and also produces a very sparse discretized system equation which can be efficiently solved. For large systems, it is very important and desirable that the discretized equation system be sparse. Praveen Kumar and Dodagoudar [23] used the meshfree radial point interpolation method to model two dimensional contaminant transport through saturated media. The shape function was constructed by employing thin plate spline radial basis function (TPS-RBF) augmented with polynomials using local support domains. Galerkin weak formulation was then applied to discretise the governing equation. Numerical results from the method were compared with analytical, FEM and experimental results, and was found to be satisfactory. Additionally it was observed that the meshfree method results were free from numerical oscillations and insensitive to Peclet constraints. Meenal and Eldho [17] developed a meshfree point collocation method (PCM) based model using multi-quadrics RBF to solve the groundwater flow equation in an unconfined aquifer and applied it to a case study using field data. The authors have extended the method to solve the groundwater contaminant transport equation [18]. As further application of their method, the meshfree point collocation method based groundwater flow and transport model was coupled to an optimization model [19] and used to solve a groundwater contamination remediation problem. The authors further developed simulation-optimization (SO) model for an unconfined aquifer in Meenal and Eldho [20] and a multi-objective SO model in Meenal and Eldho [21] for groundwater remediation. Rectangular local support domains were employed in their work. However as noted by the authors, their method requires the creation of fictitious or dummy node to keep uniform shape and size of the local support domains, and for shape function creation for nodes in the periphery or boundary. In this study, we further explore the application of radial point collocation method (RPCM) to solve groundwater flow problems and address a few of the issues noted above. It is used to model a confined aquifer with given field measured parameters. Circular

21

Local Support Domains

: Field Node

: Point of interest

Fig. 1. Local support domains used in mesh free methods.

local support domain is employed here, giving the advantage to change only the radius in order to include or exclude more nodes in the construction of shape functions without affecting other parts of the program. In the present work, there is no need for creation of dummy nodes as reported by Meenal and Eldho [17] to deal with the irregular boundary problem. Another major issue in the previous studies is in the treatment of the derivative boundary condition which causes poor accuracy and instability problems in point collocation methods. In this study, the derivative boundary condition has also been discretized using a direct collocation method which gives good accuracy without the need of complicated treatments or the creation of fictitious or dummy nodes. Both the multi-quadrics and Gaussian RBFs are used in this study for interpolation and their performances are compared. We have also demonstrated how the numerical instability associated with high Peclet transport problems can be handled easily in the meshfree RPCM method. The formulation and implementation of the method presented in this work are organized as follows. Section 2 gives the numerical formulation of the meshfree RPCM method. The discretization of the transient governing equations for groundwater flow and transport equation by the meshfree RPCM method is presented. In Section 3, the proposed method is verified for one and two dimensional problems and its capability to handle highly advective transport problem is demonstrated. In Section 4, the proposed method is applied to a field case study. Section 5 discusses the results from the case study. Conclusions are given in Section 6.

2. Governing equations and meshfree formulation 2.1. Groundwater flow The governing equation describing the groundwater flow in a two dimensional inhomogeneous confined aquifer is given as [2]:

   ∂ ∂h ∂ ∂h ∂h Tx þ Ty ¼ S þ Q w δðx  xi Þ y  yi q ð1Þ ∂x ∂x ∂y ∂y ∂t Following initial conditions are used for transient analysis, hðx; y; 0Þ ¼ h0 ðx; yÞ; x; y A Ω

ð2Þ

Generally, the boundary conditions (BC) can be of two types, the prescribed head or flux. It can be written as: hðx; y; t Þ ¼ h1 ðx; y; t Þ; x; y A ∂Ω1 ðDirichlet BCÞ T

∂h ¼ q1 ðx; y; t Þ; x; yA ∂Ω2 ðNeumann BCÞ ∂n

ð3Þ ð4Þ

where hðx; y; t Þ is the Piezometric head (m) which is the state variable, Tx, Ty are the transmissivities (m2/day) in x and y

22

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

directions, S is the storage coefficient, Qw is the source or sink term (m3/day/m2). The flow region is represented by Ω while the ∂ denotes the normal boundary is denoted by ∂Ω. In Eq. (4), ∂n derivative to the boundary, h0 ðx; yÞ is the initial head in the flow domain (m), h1 ðx; y; t Þ is the known head value of the boundary head (m) and qðx; y; t Þ is the recharge rate (m/day) while q1 ðx; y; tÞ is the known inflow rate (m3/day/m). 2.2. Contaminant transport The governing equation for contaminant transport of a single chemical constituent in groundwater is given by [7]:      ∂C ∂ ∂C ∂ ∂C ∂ ∂ Dxx þ Dyy  ðV x C Þ  V yC R ¼ ∂t ∂x ∂x ∂y ∂y ∂x ∂y 

c0 w q C  RλC þ w θb θ

ð5Þ

where V x ; V y are respectively the seepage velocity in x and y direction (m/d); Dxx ; Dyy are the components of dispersion coefficient [m2/day] tensor in x and y directions respectively, C is the concentration of the dissolved species (kg/m3), λ is the reaction rate constant (day  1), w is the elemental recharge rate with solute concentration c 0 , θ is the porosity , b is the aquifer thickness while R ¼ 1 þ ρb K d =θ is the retardation factor, in which ρb is the media bulk density, K d is the sorption coefficient and qw is the volumetric pumping rate from a source. For transient analysis, initial condition is of the form, C ðx; y; 0Þ ¼ f ðx; y; 0Þ; ðx; yÞ A Ω

ð6Þ

while the general boundary conditions are of the form, C ðx; y; t Þ ¼ g 1 ðx; y; t Þ;ðx; yÞ A ∂Ω1 ðDirichlet boundary conditionÞ  Dxx

ð7Þ

   ∂C ∂C nx þ Dyy ny ¼ g 2 ðx; y; t Þ; ðx; yÞ A ∂Ω2 ∂x ∂y

ðNeumann boundary conditionÞ

ð8Þ

Here f is a given function in Ω; g 1 ; g 2 are functions along boundaries; and nx ; ny are the components of the unit outward normal vector to the given boundary. The seepage velocities V x ; V y are evaluated from the solutions of the flow equations using the following relations [2], Vx ¼ 

K x ∂h K y ∂h ;V ¼  θ ∂x y θ ∂y

ð9Þ

where θ is the porosity of the medium while K x and K y are the hydraulic conductivities. The elements of the dispersion coefficient tensor are evaluated from the longitudinal dispersivity ðαL Þ and

transverse dispersivity ðαT Þ using the relations [2]: Dxx ¼

αL V 2x þ αT V 2y V

; Dyy ¼

αL V 2y þ αT V 2x V

; Dxy ¼ Dyx ¼

ðαL  αT ÞV x V y V ð10Þ

where V 2 ¼ V 2x þ V 2y . Eqs. (9) and (10) together provide the linkage or coupling between the groundwater flow and transport equations i.e. the velocities computed from the flow equation are used as inputs to the transport equation through these equations. 2.3. Radial basis function interpolation To solve any PDE, first the unknown field variable is approximated using trial or shape functions. A local support domain of a point x determines the number of nodes to be used to approximate the function value at x. It is a set of nodes in the neighborhood of the point. The constructed shape function will not be used or regarded as zero outside this local support domain; hence it is called locally supported. The support domain can have different shapes and also its dimension and shape can be different for different points of interest x, as shown in Fig. 1. Circular or rectangular support domains are the most commonly used shapes. The approximation of a function h(x) within a local support domain can be constructed as a linear combination of n radial basis functions and m polynomial basis functions as given below [15]: hðxÞ ¼

n X

ai Ri ðxÞ þ

i¼1

m X

P j ðxÞbj ¼ RT ðxÞa þ P T ðxÞb

ð11Þ

j¼1

where Ri ðxÞ is a radial basis function (RBF) such as a multi-quadrics or Gaussian, n is the number of points in the support domain, ai and bj are unknown coefficients to be determined. P j ðxÞ are polynomial basis functions. In the RBF Ri ðxÞ, the only variable is the distance between the point of interest x, called the data site and a node at xi , known as the center point, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 ð12Þ r i ¼ ðx xi Þ2 þ y yi ðin case of 2 D problemsÞ There are several types of radial basis functions. Kansa [12], Franke and Schaback [6], Sharan et al. [27] etc. have extensively investigated the characteristics of RBFs. In this paper, two of the widely used RBF functions viz., the multi-quadrics RBF (MQ-RBF) and the exponential or Gaussian RBF (EXP-RBF) have been used for function interpolation. These RBFs are defined as below,

q Ri ðx; yÞ ¼ r 2i þC 2s ð13Þ

= 500 = 500 / /

= 0

/

x = 1000 m

x=0

Local Support Domain (size = 2*Nodal distance)

Point of interest Fig. 2. One dimensional transport problem.

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

"

 2 # r Ri ðx; yÞ ¼ exp  αc i dc

ð14Þ

where q; C s ; αc are known as the shape parameters of the RBF and dc is the characteristic length that relates to the nodal spacing in the local support domain. The standard MQ-RBF has q ¼ 70:5, but this parameter can also be left as an additional open parameter [31] . In this study, the parameter q has been kept at 0.98 as in Liu and Gu [15]. It was observed that changing this parameter to the standard value of 0.5 did not produce any significant effect on the accuracy and stability of the solutions in this study. The shape parameter of the MQ-RBF, C s is usually defined in terms of the characteristic length (dc ) i.e. C s ¼ αc dc . The value of the shape parameter of the RBFs has a very profound impact on the quality of the interpolation within a support domain and hence on the overall solution accuracy. Choosing the optimum value of this parameter is very crucial. Currently there is no unique criterion or established method for finding the optimum value of this parameter and is generally found by sensitivity or parametric studies. Some strategies of choosing the shape parameter can be found in Schaback and Wendland [26], Rippa [25] and Wright [30]. The unknown coefficients ai and bj in Eq. (11) can be determined by enforcing the interpolation function to pass through all n nodes within the support domain [15]. The interpolation of the function at the kth point has the following form: n m X     X ai Ri xk ; yk þ bj P j ðxk ; yk Þ; k ¼ 1; :::; n h xk ; yk ¼ hk ¼ i¼1

ð15Þ

which yields n simultaneous linear algebraic equations with n þ m unknowns. The additional m equations can be obtained from the following m constraint equations imposed on the polynomial term to guarantee unique approximation [8]: n X

  P j xi ; yi ai ¼ P Tm a ¼ 0 ; j ¼ 1; 2; ……m

ð16Þ

i¼1

The system of Eqs. (15) and (16) can be solved to obtain the unknown coefficients ai and bj . Details of solving the equations can be found in Liu and Gu [15]. Substituting these coefficients back into Eq. (11), the interpolation can be expressed as, T

hðxÞ ¼ Φ ðxÞhs

ð17Þ

where ϕðxÞ, are known as the shape functions and expressed as: 

ΦT ðxÞ ¼ ϕ1 ðx; yÞ ϕ2 ðx; yÞ ::: ϕn ðx; yÞ



ð18Þ

 T and, hs ¼ h1 h2 ⋯hn is the vector of nodal head values at the support domain nodes. The shape functions depend only upon the position of the nodal points. The derivatives of hðx; yÞ at any point xI ðxI ; yI Þ can be easily determined as below: hðxI Þ ¼ Φ hs ¼ T

n X

ϕi hi

2.4. Discretization of governing equations using point collocation method Using the shape functions derived in the previous section i.e. Eq. (17), the piezometric head can be approximated as, hðx; t Þ ¼

n X

ð19Þ

ϕðxÞhi ðt Þ ¼ ΦT ðxÞhs ðtÞ

ð21Þ

i¼1

in which n is the number of nodes used in the local support domain. The discretization of the governing equation is achieved by simple collocation at all the internal nodes and using the above approximation for the piezometric head [15]. Thus by collocation at the point xr ðxr ; yr Þ, the governing Eq. (1), for a homogeneous and isotropic aquifer, is discretized as: " #   T T ∂2 Φ ∂2 Φ ∂h T ðxr Þ þ ð t Þ ¼ S ð x Þ þ Q w δðxr  xi Þ  qðxr Þ ð22Þ h s r ∂t xr ∂x2 ∂y2 Using nodal indices, the above equation can be re-written as, " #   T T ∂2 Φ ∂2 Φ ∂h þ þ Q w δ ðxr  xi Þ  qr ð23Þ hs ðt Þ ¼ Sr Tr 2 2 ∂t r ∂x ∂y The time discretization can be done by using the θ-method which is the most commonly used algorithm for time integration. In this method, the time derivative is replaced by a simple forward difference while the solution is replaced by a weighted value of the previous time-step solution and current solution as below [24]: t þ Δt

j¼1

23

∂h h ¼ ∂t

t

h Δt

ð24aÞ

and, h ¼ θh

t þ Δt

  t þ 1θ h

ð24bÞ

Here, θ is a relaxation parameter which lies in the interval [0, 1] and is used to control the accuracy and stability of the algorithm. This method is unconditionally stable for all θ Z 1=2 [24]. Commonly used values of θ are 0, 1/2 2/3 and 1. When θ o1/2, the algorithm is only conditionally stable. When θ ¼ 1=2, it is known as the Crank–Nicholson time stepping method and the same is used in this study. Substituting Eq. (24a) and (24b), the discretized Eq. (23) can be written as, " # T T ϕT hs t þ Δt  htr ∂2 ϕ ∂2 ϕ t þ Δt Sr þ Q w δðxr xi Þ  qr ¼ θT r þ hs Δt ∂x2 ∂y2 " # T T   ∂2 ϕ ∂2 ϕ t þ 1θ Tr þ ð25Þ hs ∂x2 ∂y2 t

It is assumed that the nodal heads, hr at time t are known and is used as the initial condition to advance the solution to the next time level, t ¼ t þ Δt. Re-arranging the terms, " !# T T ∂2 ϕ ∂2 ϕ T t þ Δt Sr ϕ  θΔtT r þ hs ∂x2 ∂y2 " #   ∂ 2 ϕT ∂ 2 ϕT t t ¼ Sr hr þ Δt 1  θ þ ð26Þ hs  Q w δðxr  xi Þ þ qr ∂x2 ∂y2

i¼1

∂hI ∂Φ ¼ hs ¼ ∂x ∂x T

n X ∂ϕ

i

∂ hI ∂ Φ hi ; 2 ¼ hs ¼ ∂x2 ∂x 2

2

T

n X ∂2 ϕ

i

hi

ð20aÞ

T T n n X X ∂hI ∂Φ ∂ϕi ∂2 hI ∂2 Φ ∂2 ϕi ¼ hs ¼ hi ; 2 ¼ hs ¼ hi 2 ∂y ∂y ∂y ∂y ∂y ∂y2 i¼1 i¼1

ð20bÞ

i¼1

∂x

i¼1

∂x2

The above equation is established for all the internal nodes. A similar discretization can be performed for the transport equation using the RBF point interpolation approximation of the concentration, Cðx; tÞ as below: C ðx; t Þ ¼

n X i¼1

ϕi ðxÞC i ðtÞ ¼ ΦT C s

ð27Þ

24

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

Fig. 3. (a) Concentration profile at 10 m downstream as compared to the analytical result for various values of the shape parameter of the MQ-RBF (number of nodes: 301 size of support domain: 4*dc). (b) Concentration profile at 10 m downstream as compared to the analytical result for various values of the shape parameter of the EXP-RBF (number of nodes: 301; size of support domain: 4*dc).

This discretization yields the following equation at the node r, assuming homogeneous dispersion and uniform flow field, and neglecting the source and reaction terms in Eq. (5): " !# T T T T ∂2 Φ ∂2 Φ ∂Φ ∂Φ T  V yr þ Dyyr  V xr C st þ Δt RΦ  θΔt Dxxr ∂x ∂y ∂x2 ∂y2 " # T T T T   ∂2 Φ ∂2 Φ ∂Φ ∂Φ t  V yr ¼ RC r þ 1  θ Δt Dxxr þDyyr  V xr C ts ∂x ∂y ∂x2 ∂y2 ð28Þ For heterogeneous problems, the model domain is divided into zones of different transmissivities (for flow equation) and dispersivities (for transport equation) depending on the variation of the transmissivity and dispersivity respectively. The values of transmissivity and dispersivity of a zone are considered as same for all the nodes lying in that particular zone. The flow and transport equations are formulated by considering the variations of transmissivity and dispersivity. The seepage velocities V x ; V y can be evaluated from the groundwater flow equation solutions using Eq. (9) (Darcy's law) as, K r ∂Φ K r ∂Φ h ;V ¼  h θ ∂x s yr θ ∂y s T

V xr ¼ 

T

ð29Þ

The dispersion coefficients are evaluated using Eq. (10). Eqs. (26) and (28) are formed for all the internal nodes. For the

boundary nodes, the proper boundary conditions have to be imposed. This can be done as illustrated below. The Dirichlet boundary condition can be easily discretized by directly substituting the approximation equations. For example for the flow model, Eq. (17) can be directly substituted into the boundary Eq. (2) i.e. hðx; yÞ ¼ h0

ð30Þ

or, on using Eq. (17)

ϕT hs ¼ h0

ð31Þ

Implementing the derivative boundary conditions requires some special considerations as the presence of derivative (Neumann) boundary conditions can drastically reduce the accuracy and stability of the solution of collocation or strong form methods [15]. A number of strategies have been proposed by Liu and Gu [15] to deal with derivative boundary conditions such as direct collocation, method of fictitious points, Hermite-type collocation, method of regular grids, dense nodes method and weak–strong form method. In this study, the direct collocation approach has been implemented. In this approach a discretization scheme, similar to that for the internal nodes, is constructed at the derivative boundary nodes. For example, the derivative boundary

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

25

Fig. 4. Concentration profile at downstream distance (a) x¼ 10 m and (b) x ¼ 25 m. Table 1 Error measures of the concentration profile plots of Fig. 4. RBF

x ¼10 m RMSE

x ¼25 m RMSE

Multi-quadrics Exponential

1.4 1.6

0.4 0.6

Size of local support domain (multiples of nodal separation)

2 3 4 5

condition, ∂h AðxÞ þ BðxÞh þ qb ðxÞ ¼ 0 ∂x can be discretized by direct collocation as, " # T ∂Φ T þ BðxÞΦ hs ¼  qb ðxÞ AðxÞ ∂x

Table 2 Effect of different nodal separations and sizes of the local support domain on the RMS error of meshfree RPCM method (concentrations taken at the end of simulation time). Number of nodes (N ¼101) Pe¼ 72

Number of nodes (N ¼ 201); Pe¼36

Number of nodes (N ¼ 401) ; Pe¼18

RMSE

e0 (%) RMSE

e0 (%)

RMSE

e0 (%)

12.54 7.7 5.7 5.2

5.8 3.5 2.6 2.4

1.5 0.5 0.5 0.46

1.2 0.5 0.2 0.2

0.5 0.2 0.09 0.1

3.3 1.2 1.1 1.0

ð32Þ the groundwater flow equation can be assembled in matrix form as, ð33Þ

Let N be the total number of nodes including internal and boundary nodes. Then nodal discretized equations derived above for

 t þ Δt  t ½K NN h N1 ¼ ½F NN h N1 þ fQ gN1

ð34Þ

A similar system of equation follows for the contaminant transport equation. These equations can be solved using any of the matrix inversion methods or by using iterative methods.

26

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

Fig. 5. (a) Stabilizing effect of enlargement of local support domain for high Peclet transport problem (number of nodes: 101 ; concentrations obtained after simulation of 2 years). (b) A combination of enlargement of local support domain and increasing the nodal density has nearly removed the numerical instability observed in the previous figure (number of nodes: 401; concentrations obtained after simulation of 2 years).

(0, 0)

Point source: contaminant injection

(150, 150)

(450, 300)

Constant-head boundary

Constant-head boundary

No-flow boundary

3. Model development and verification Based on the above formulations, a coupled flow and transport RPCM model (CFTM-RPCM) has been developed in MATLAB. The meshfree collocation formulation is verified for one and two dimensional cases through standard benchmark problems as described in Sections 3.1 and 3.3. In Section 3.2, the capability of the meshfree model to handle high Peclet transport problem is demonstrated. In case of the two dimensional problem in Section 3.3, the local support domain used is circular in shape and its size can be easily controlled by varying the radius, thereby controlling the number of nodes within the local support domain.

No-flow boundary

Fig. 6. Confined aquifer considered as benchmark for coupled flow and transport simulation.

The presence of singular sources or sinks such as pumps leads to instability of the solution and loss of accuracy in the whole domain, known as pollution effect. This has been studied by Tornberg and Engquist [29], Ashyraliyev et al. [1], Jung [11] etc. and several techniques have been proposed to overcome this problem. In the present study, the Gaussian function regularization as suggested in Jung [11] has been implemented which is a simple yet very effective method. It yields a first order convergence near the jump discontinuity.

3.1. One dimensional transport problem For verification of one dimensional transport problem, the following combined advection–dispersion equation is considered, ∂C ∂C ∂2 C ¼ v þ DL 2 ∂t ∂x ∂x

ð35Þ

The boundary and initial conditions are, Initial condition: C ðx; 0Þ ¼ 0; x Z 0; boundary conditions: C ð0; t Þ ¼ C 0 ; t Z 0 and C ð1; t Þ ¼ 0; t Z0

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

27

Fig. 7. Comparison of the meshfree RPCM and FEM solutions for two-dimensional transport from a continuous point source in a confined aquifer.

Fig. 8. Location map of the study area.

The analytical solution to this problem has been given by Ogata and Banks [22] as,
      C0 x vt vx x þ vt C ðx; t Þ ¼ erfc pffiffiffiffiffiffiffiffi þ exp ð36Þ erfc pffiffiffiffiffiffiffiffi DL 2 2 DL t 2 DL t where x is the distance from the injection point. To verify the RPCM method for this 1D transport problem, a one dimensional domain as depicted in Fig. 2(a) is considered. Solute (contaminant) of concentration C0 ¼500 mg/L is injected at x¼ 0 and it is required to find the concentration at different downstream distances (x 40) after 2 years. The parameters used for this problem are [4]: pore flow velocity (v)¼0.0259 m/day and longitudinal dispersion coefficient (DL )¼0.0360 m2/day. The modeling domain is considered to be 1000 m in length, sufficient enough for the boundary condition C ð1; t Þ ¼ 0; t Z 0 to

be satisfied i.e. far away from the points where concentration is required. Considering the very small pore flow velocity of 0.0259 m/day, the above assumption regarding zero concentration boundary condition is justified since the simulation time considered is only 2 years. Fig. 2b illustrates the construction of local support domain for a support domain size equal to twice the nodal separation (dc ). For finding the optimum value of the shape parameters for the RBF function used in the interpolation, a sensitivity analysis is performed. The RBF functions used in this study are the multiquadrics and exponential or Gaussian RBF. Fig. 3(a) and (b) shows respectively the results of the sensitivity analysis for each of the RBFs. It is found that for the MQ-RBF, the values of the shape parameter (C s ¼ αc dc ) in the range of 3–7 times the nodal separation gave very good accuracy whereas for the exponential

28

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

Fig. 9. Nodal distribution map and boundaries of the field case study (numbers in the model area represents the node numbering used in the program).

Fig. 10. Transmissivity distribution map of the study area.

RBF very good accuracy was observed when the shape parameter (αc ) was in the range from 0.1 to 0.08. It was observed that the exponential RBF was more sensitive to changes in the shape parameter values as compared to MQ-RBF. The values of the shape parameter can further be refined using some of the recent approaches suggested in Schaback and Wendland [26], Rippa [25], Wright [31] etc. to obtain the optimum value. With the shape parameters determined, the 1D transport equation is solved and Fig. 4 shows the concentration profiles at two downstream distances as compared to the exact solution of Ogata and Bank [22]. Results of using MQ-RBF and EXP-RBF as interpolation functions are also shown for comparison.

In the above figures, number of nodes used for the simulation is 401 and the size of the support domain is 4 times the nodal separation (dc ). Since the analytical solution is known, the root mean square error (RMSE) as defined below is used as the measure of error: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  num 2 C i  C exact i RMSE ¼ N i¼1 N X

ð37Þ

Table 1 shows the RMS error measures of the above plots: From Fig. 4a and b and Table 1, it can be seen that there is very good agreement between the mesh free calculated concentration values and the analytical results. It is also seen that both the MQRBF and EXP-RBF gave almost equally good results.

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

3.2. Model application to high Peclet number problems Since in case of transport problems, higher Peclet numbers (Pe ¼ DvL  Δx) causes numerical instability, it is important to analyze the capability of the meshfree RPCM method to handle such instability issues. Two simple approaches to handle the numerical instability are considered in this paper [15]. One obvious way is to reduce the nodal separation (Δx), thereby reducing Pe itself. The other is to increase the size of the local support domain so as to capture the upstream information which is the source of the numerical instability. This later approach is not easy to implement in FDM or FEM where the mesh is predefined but in meshfree RPCM it can be easily implemented. The one-dimensional transport problem [22] from the previous section is analyzed here to examine the capability of the meshfree Table 3 Comparison of field measured heads at selected monitoring bore wells to FEM and Meshfree RPCM models. Borewell no. (n)

Field measured head (m)

FEM model prediction (m)

Meshfree RPCM model prediction (m)

1 2 3 4 5 6 7 8 9 10 11 14 15 16 17 18 19 20 21

47.995 48.383 48.147 48.112 46.004 45.799 45.129 46.194 45.82 45.719 45.704 45.619 45.469 45.439 45.387 45.699 45.319 45.139 44.934

48.0 48.0 47.8 47.2 46.7 46.7 46.9 47.0 46.7 46.6 46.6 46.6 46.6 46.5 46.6 46.0 46.0 46.1 46.0

48.3 48.4 48.4 47. 4 46.7 46.8 46.9 47.1 46.8 46.6 46.6 46.7 46.6 46.6 46.6 45.9 46.0 46.1 46.0

n

29

RPCM method to handle high Peclet transport simulation. For the purpose of analysis here, the pore flow velocity is increased by 10 times to 0.259 m/day resulting in a Peclet number value of 72, for a nodal separation of 10 m (corresponding to 100 equal subdivision of the domain). In addition to the RMSE defined in Eq. (37), an additional measure of error [15] as defined below is also used for the analysis viz., vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPN  num 2 u C exact i 1 Ci e0 ¼ t i ¼P N exact 2 i ¼ 1 ðC i ¼ 1 Þ

ð38Þ

Table 2 shows the results of the two approaches of handling numerical instability where the meshfree RPCM results are compared to the exact solution. Fig. 5a and b show the stabilization of the numerical instability by using a combination of enlargement of local support domain and using denser node distribution. MQ-RBF interpolation has been used in this analysis. But the results will apply to EXP-RBF interpolation also. It can be seen that the combination of enlargement of the local support domain and reducing the nodal separation is very effective in handling numerical instability present in high Peclet number transport problems. The advantage of using meshfree RPCM method is that this approach can be implemented very easily. Other approaches for stabilization can be found in Liu and Gu [15].

Table 4 Comparison of head output of Meshfree RPCM and FEM models. Nodes Head from mesh free RPCM model (m)

Head from FEM model (m)

Percentage difference

448 558 662 552 670 730 545

46.8 46.8 46.5 46.5 46.8 46.5 46.1

0.21 0.21 0.22 0.22 0.21 0 0

46.9 46.9 46.6 46.6 46.9 46.5 46.1

Ref. Singh and Sarma [28].

Fig. 11. Comparison of head contour predicted by Mesh free RPCM and FEM.

30

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

Fig. 12. Comparison of concentration plume predicted by meshfree CFTM-RPCM and FEM. Table 5 Comparison of concentration predicted by CFTM-RPCM and FEM models. Node Concentration (ppm) after 5 years

448 558 662 552 670 730 545

Concentration (ppm) after 10 years

RPCM

FEM

Percentage difference

RPCM

FEM

Percentage difference

442 985 905 732 781 785 502

451 976 927 720 757 804 538

2.0 1.0 2.4 1.6 3.1 2.4 7.2

442 986 959 733 792 964 618

451 975 944 722 767 930 662

2.0 1.1 1.6 1.5 3.2 3.5 7.1

3.3. Two dimensional transport problem For verification of the 2D coupled groundwater flow and solute transport model, a two-dimensional transport of solute injected continuously from a point source in a steady-state uniform flow field is considered [32]. In this problem, it is assumed that a contaminant enters at a point, such as an injection well or toxic spill, and spreads through the aquifer with time. The study area is shown in Fig. 6. The flow model is surrounded by constant-head boundaries on the east and west borders and no-flow boundaries on the north and south borders. The head values at the constant-head boundaries are arbitrarily chosen to establish the required hydraulic gradient so as to give a uniform flow field with groundwater seepage velocity (v)¼1/3 m/day. The simulation period is chosen so that the plume developed from the point source does not reach the boundaries. The model parameters used in the simulation are listed below [32]: Cell width along rows (Δx) ¼10 m, cell width along columns (Δy)¼10 m, layer thickness (Δz)¼ 10 m, porosity (θ) ¼0.3, longitudinal dispersivity¼10 m, ratio of transverse to longitudinal dispersivity¼0.3, volumetric injection rate ¼1 m3/day, concentration of the injected water ¼1000 ppm, simulation time (t)¼ 365 days The study area has been divided into a 46  31 uniformly spaced nodes, along the x and y axes respectively and meshfree

RPCM formulation as formulated in Section 2 was applied. The Crank–Nicholson scheme was implemented for time discretization with a time step of 5 days. The same problem was solved with FEM using 1092 triangular elements and 585 nodes in COMSOL Multiphysics (www.comsol.com) for the purpose of comparison. As in the 1D case, a numerical sensitivity analysis with respect to the size of support domain and the shape parameter value was conducted to arrive at an optimum value of these parameters taking into account the balance between accuracy, stability and computational time. While the support domain size mostly affects the computational time, the shape parameter has a profound impact on the accuracy and stability of the method. A support domain size (radius) of 3–4 times the nodal separation yielded sufficient accuracy within acceptable computational time. For finding the optimum value of the shape parameter, a sensitivity analysis is performed for both the RBF functions viz. MQRBF and EXP-RBF. It is found that for the MQ-RBF, the values of the shape parameter (C s ¼ αc dc ) in the range of 3–7 times the nodal separation gave very good accuracy whereas for the EXP-RBF the optimum values of the shape parameter (αc ) was in the range from 0.1 to 0.08. It can be noted that the ranges of the shape parameters are same as for the 1D transport case. Fig. 7 shows the contour plot of the contaminant distribution at the end of the simulation time i.e. 365 days obtained by meshfree RPCM method using both the MQ-RBF and EXP-RBF as interpolation functions. It is observed that both MQ-RBF and EXP-RBF produces almost identical results and there is no way to conclude which RBF function is superior to the other. This means that we can freely choose either of the RBF functions provided the proper shape parameter values are used. The result from the FEM (COMSOL) simulation is also shown alongside for comparison. From the figure, it can be seen that there is good agreement between the meshfree RPCM and FEM results.

4. Field case study To test the effectiveness of the developed model (CFTM-RPCM) for a field problem, a case study of a confined regional aquifer is considered for coupled flow and contaminant transport modeling.

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

Fig. 13. Concentration profile at selected nodes along the path of the contaminant plume.

Fig. 14. Propagation of the contaminant plume in the aquifer (number in the figure shows concentration in ppm).

31

32

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

For this case study, the MQ-RBF has been used as the interpolation function. This choice is arbitrary since as demonstrated in previous examples, both the MQ-RBF and EXP-RBF produces nearly the same results. Fig. 8 shows the location map of the study area. The field measured hydrogeological parameters of the site are obtained from Singh and Sarma [28]. The site is bounded by a lake on the north, north-east, west and south-west boundaries. There are no water bodies on the rest of the boundary. The area to be modeled is approximately 4.5 km2. The aim of the study is to track the movement of any contaminant released inadvertently or by leakage from the underground tanks. Fig. 9 shows the case study site and the nodal distribution used for discretization of the governing equations. A total of 1008 nodes were used corresponding to a separation between the nodes of 49.6 m along x-direction (Δx) and 42.8 m along the y-direction (Δy). As shown in the figure, there is a recharge zone within the model domain which corresponds to a water pond and is known to be leaking water to the aquifer. The presence of this recharge zone to the aquifer is observed as a mound in the water table map of the area [28]. Also shown is an area within the complex with underground containers which is assumed to be the contaminant source location. The boundaries adjacent to the lake are treated as constant head boundaries while the rest of the boundary is treated as flux boundaries. The level of water in the lake is maintained by the irrigation department of the local government. This water level, adjusted to mean sea level, is used as the constant head value at the boundaries bounded by the lake. The flux boundary value is estimated and adjusted during calibration of the model. The best estimated value of this flux obtained from model calibration is 6.1  10  5 m2/d. The distribution of transmissivity is obtained from pumping tests carried out at five different sites within the modeling zone [28]. Fig. 10 shows the different transmissivity zones inferred from the field experiments. It is observed that transmissivity of the area varies from a minimum of 30.0 m2/day to a maximum of 170.0 m2/day. Other required aquifer parameters are also obtained from the same reference. The recharge through the pond is estimated from the make up level of the pond and is adjusted during the calibration. The final recharge value used was 0.045 m/day. During calibration of the model, the hydrogeologic parameters are adjusted within a range of 20% from their field estimated values.

5. Results and discussion A mesh free CFTM-RPCM groundwater flow model of the site is constructed using the hydrogeological parameters described above [28] and the node distribution as in Fig. 9. The governing groundwater flow Eq. (1) is solved in steady state for obtaining hydraulic heads at the nodes. The level of water in the lake is maintained as also the water level of the recharging pond. Being a confined aquifer, rainfall recharge is negligible. Therefore the groundwater flow is basically steady state. Calibration of the RPCM model was done against the field measured heads. Table 3 shows the comparision of actual field measured heads in selected monitoring borewells at the site to the meshfree RPCM model predicted heads, after calibration of the model. For comparision, using the same calibrated values of the hydro-geological parameters, an FEM model was constructed in COMSOL Multiphysics. The heads predicted by this FEM model is also shown in the same table. Fig. 11 compares the contour map of the hydraulic heads as predicted by mesh free CFTM-RPCM and FEM (COMSOL) models. It can be seen that both the models predicts more or less the same hydraulic head distribution in the modeling domain. This confirms the good agreement between RPCM model and FEM models.

Table 4 shows the comparison of meshfree RPCM flow model results and those from FEM at 7 selected nodes (shown in Fig. 9) viz. nodes 448, 545, 552, 558 ,662, 670 and 730 of the nodal distribution. These nodes are also used for validating the transport model. It is seen that there is good agreement between the heads predicted by the meshfree RPCM and FEM models. With the flow model validated, a non-reactive contaminant transport model is constructed and coupled to the flow model. Steady state head values from the flow model are considered as the initial head distribution. The source of contaminant is an areal source as shown in Fig. 9 and corresponds to a particular area in the domain with underground containers which may act as potential sources of contaminants. This model aims to track the spreading of the contaminants in a postulated structural failure scenario where the containers are leaking the contaminants to the groundwater continuously for the period of simulation. The contaminant transport is governed by Eq. (5). The longitudinal dispersivity (αL ) for this problem is considered as 20 m and the transverse dispersivity (αT ) is taken as 10% of the longitudinal dispersivity. The areal contaminant source is assumed to be leaking contaminant of concentration 1000 ppm. The boundaries of the model are set to a concentration value of zero i.e. Dirichlet boundary condition. The same nodes (Fig. 9) used for the flow model are also used for the transport model. The time step size used is Δt ¼ 5 days and it is assumed that the whole area was pristine when simulation starts i.e. have zero concentration as the initial condition. For time-stepping, the Crank–Nicholson method (θ ¼ 1=2) was implemented. The total simulation time was 10 years corresponding to 730 time steps. The size of the local support domain was 3 times the average of the nodal distances along x and y-directions while the shape parameter of the multiquadrics RBF was kept at 5 times the average nodal distance. These parameters were determined also in the verification problems discussed in the previous sections. Fig. 12 shows the contaminant plume after 5 years of continuous discharge. The output from FEM is also shown alongside for comparison. The path of the contaminant plume is in the south-eastern direction. Table 5 shows a comparison of the concentrations predicted by meshfree CFTM-RPCM and FEM models at the 7 selected nodes viz. nodes 448, 545, 552,558, 662, 670 and 730 in the direction of the contaminant plume. The locations of these nodes have been shown in Fig. 9. From Table 5 and Fig. 12, it can be seen that the meshfree CFTM-RPCM results agrees very well with the FEM predictions. Fig. 13 shows the concentration profiles for the above selected nodes. It is seen that the nodes 558, 662 and 730 receives the highest concentrations being directly along the centerline of the plume while node 448 receives the lowest concentration value being in the flank. Fig. 14 shows the snapshot of the contaminant plume at the end of 1, 3, 5 and 10 years of simulation showing the propagation of the contaminant. The movement of the contaminant plume is towards the south-eastern boundary of the model domain. This study has demonstrated that meshfree CFTM-RPCM model can be effectively used to model coupled groundwater flow and transport simulations as an alternative to the grid based methods such as FDM or FEM. This method requires only a set of scattered nodes, instead of a mesh with all the nodal interconnection information. The implementation of the model is very simple regardless of the dimension of the problem and complexity of the model geometry. Boundary conditions are also easily incorporated. 6. Conclusions In this study, meshfree models based on the point collocation method have been developed to simulate coupled flow and

L. Guneshwor Singh et al. / Engineering Analysis with Boundary Elements 66 (2016) 20–33

transport problems in groundwater using multi-quadrics and exponential RBFs as the interpolation function. The developed model has been applied to one and two dimensional benchmark problems for verification. The output from the model is compared to the output obtained by analytical and FEM models and found to be in very good agreement. A field case study of groundwater flow and transport problem was simulated using the developed model and the results are compared with those obtained from an established FEM software. It was found that the meshfree method results are in very good agreement with those obtained with the FEM. It was also observed that hydraulic heads predicted by the meshfree method shows good agreement with the available field measured values of the hydraulic heads. A few issues in the application of the meshfree methods though are yet to be addressed such as choosing the optimal value of shape parameters, size and shape of the local support domain etc. In particular, the value of shape parameter has a profound impact on the accuracy. But it was found that these issues do not pose any significant hurdle in application of the method. With a numerical sensitivity or parametric analysis, these parameters can be determined. On the other hand, this method does away with the need to construct a grid or mesh, is relatively easy to implement and offers better accuracy. As shown in this study, both the MQ-RBF and EXP-RBF gave almost equal accuracy and either of them can be used for the interpolation without concerning about losing accuracy, provided the proper shape parameter values are used. It has also been demonstrated that the proposed method can easily handle the numerical instability associated with high Peclet transport problems. Considering the significant advantages offered by this method, it can serve as a good alternative to the finite difference and finite element methods.

References [1] Ashyraliyev M, Blom JG, Verwer JG. On the numerical solution of diffusionreaction equations with singular source terms. J Comput Appl Math 2008;216 (1):20–38. [2] Bear J. Hydraulics of groundwater. New York: McGraw Hill Publishing; 1979. [3] Bernal F, Guitierrez G, Kindelan M. Use of singularity capturing functions in the solution of problems with discontinuous boundary conditions. Eng Anal Bound Elem 2009;33:200–8. [4] Delleur JW. The handbook of groundwater engineering. 2nd ed.. Boca Raton: CRC Press; 2007. [5] Fedoseyev A, Friedman M, Kansa E. Improved multiquadrics method for elliptic partial differential equation via PDE collocation on the boundary. Comput Math Appl 2002;43:439–55. [6] Franke C, Schaback R. Solving partial differential equations by collocation using radial basis functions. Appl Math Comput 1997;93:73–82. [7] Freeze RA, Cherry JA. Groundwater. Englewood Cliffs,New York: Prentice Hall; 1979. [8] Goldberg MA, Chen CS, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM. Eng Anal Bound Elem 1999;23:285–96.

33

[9] Gutierrez G, Florez W. Issues of the local radial basis collocation method implementation for solving second order partial differential equation. Mec Comput I 2008;XXVI:2241–52. [10] Hardy RL. Multi-quadric equations of topography and other irregular surfaces. J Geophys Res 1971;76:1905–15. [11] Jung Jae-Hun. A note on the spectral collocation approximation of some differential equations with singular source terms. J Sci Comput 2009;39:49–66. [12] Kansa E. Multiquadrics – a scattered data approximation scheme with applications to computational fluid-dynamics-II: solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19:147–61. [13] Kansa R, Hon YC. Circumventing the Ill-conditioning problem with multiquadric radial basis functions. Comput Math Appl 2000;39:123–37. [14] Lee C, Liu X, Fan S. Local multi-quadric approximation for solving boundary value problems. Comput Mech 2003;30:396–409. [15] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Netherlands: Springer; 2005. [16] Madych WR, Nelson SA. Multivariable interpolation and conditionally positive definite functions II. Math Comput 1990;54:211–30. [17] Meenal M, Eldho TI. Simulation of groundwater flow in an unconfined aquifer using meshfree point collocation method. Eng Anal Bound Elem 2011;35:700–7. [18] Meenal M, Eldho TI. Two-dimensional contaminant transport modeling using meshfree point collocation method (PCM). Eng Anal Bound Elem 2012;36:551–61. [19] Meenal M, Eldho TI. Simulation-optimization model for groundwater contamination remediation using meshfree point collocation method and particle swarm optimization. In: Sadhana J. Academy Proceedings in Engineering Sciences; 2012, 37(3). p. 351–69. [20] Meenal M, Eldho TI. Groundwater remediation optimization using a point collocation method and particle swarm optimization. Environ Model Softw 2012;32:37–48. [21] Meenal M, Eldho TI. Multi-objective groundwater remediation design using a coupled mesh free point collocation method and particle swarm optimization. J. Hydrol. Eng., ASCE 2014;19:1259–63. [22] Ogata A, Banks RB. A solution of the differential equation of longitudinal dispersion in porous media. Washington, D.C.: U.S. Geological Survey, U.S. Government Printing Office; 1961. [23] Praveen Kumar R, Dodagoudar GR. Two-dimensional modeling of contaminant transport through saturated porous media using the radial point interpolation method (RPIM). Hydrogeol J 2008;16:1497–505. [24] Richtmyer RD, Morton KW. Difference methods for initial value problems. New York: John Wiley & Sons; 1967. [25] Rippa S. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Math 1999;11:193–210. [26] Schaback R, Wendland H. Adaptive greedy techniques for approximate solution of large RBF systems. Numer. Algorithms 2000;24:239–54. [27] Sharan M, Kansa EJ, Gupta S. Application of the multiquadric method for numerical solution of elliptic partial differential equations. Appl Math Comput 1997;84:275–302. [28] Singh VS, Sarma MRK. Groundwater data generation at Kakrapar, Gujarat (Technical Report no. NGRI-2009-GW-673). Hyderabad (India): National Geophysical Research Institute; 2009. [29] Tornberg AK, Engquist B. Numerical approximation of singular source terms in differential equations. J Comput Phys 2004;200:462–88. [30] Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int J Numer Methods Eng 2002;54:1623–48. [31] Wright GB. Radial basis function interpolation: numerical and analytical developments. Boulder: University of Colorado; 2003. [32] Zheng C, Wang P. MT3DMS: a modular three-dimensional multispecies transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems. Alabama: University of Alabama; 1998. [33] Zuppa C, Cardona A. A collocation meshless method based on local optimal point interpolation. Int J Numer Methods Eng 2003;57:509–36.