Global multiquadric collocation method for groundwater contaminant source identification

Environmental Modelling & Software 26 (2011) 1611e1621

Contents lists available at ScienceDirect

Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft

Global multiquadric collocation method for groundwater contaminant source identification Zi Li, Xian-zhong Mao* Research Center for Environmental Engineering and Management, Graduate school at Shenzhen, Tsinghua University, Shenzhen 518055, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 April 2010 Received in revised form 14 July 2011 Accepted 16 July 2011 Available online 23 August 2011

In this paper, a radial basis collocation method (RBCM) based on the global space-time multiquadric (GST-MQ) is developed to solve the inverse problem of groundwater contaminant source identification. This deterministic method directly induces the problem to a single-step solution of a system of linear algebraic equations in the entire space-time domain. The least-square-based radial basis collocation method (LS-RBCM) is introduced to overcome the ill-posedness of the linear system. To increase confidence in the GST-MQ solutions, the sensitivity analysis with respect to calculation parameters, observation data and model uncertainty is investigated. The performance of the application examples of contaminant source identification in one- and two-dimensional porous media demonstrates that the proposed method poses the meshfree advantage and direct identification of contaminant source with efficient methodology. The GST-MQ method provides a robust tool for estimating spatial plume distribution as well as the release history of groundwater contaminant source from concentration measurements. Ó 2011 Elsevier Ltd. All rights reserved.

Keywords: Groundwater Contaminant source identification Inverse problem Radial basis collocation method Global space-time multiquadric method Sensitivity analysis

1. Introduction Groundwater, which accounts for one third of the world population’s source of drinking water, is extremely susceptible to contamination. Major contamination sources are waste disposal sites, landfills and underground storage tanks and other contaminants such as pesticides, fertilizers and salt water intrusion (Commission on Geosciences, Environment and Resources, 1993). In order to decrease the detriment to human health and ecosystems arising from groundwater contamination, effective efforts should be immediately made to clean up groundwater systems once a contaminant plume is detected. The cost of mitigating groundwater contamination problems is staggering; in 1994 the National Academy of Sciences estimated that over a trillion dollars will be spent during the next 30 years on the remediation of contaminated soil and groundwater (USEPA, 1999). It is essential to develop inverse models for identifying groundwater contaminant source in order to reduce the cost requirements associated with the complex process of remediation. The inverse modeling approaches are commonly used as a helpful tool for the identification of contaminant sources and

* Corresponding author. E-mail address: [email protected] (X.-z. Mao). 1364-8152/$ e see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2011.07.010

evaluation of parameters in many sectors of environmental science, especially in air environment (Peña et al., 2010; Keats et al., 2010), and water environment (Golosov and Kirillin, 2010; Rodríguez and Martos, 2010; Shrestha et al., 2008). Inverse problems are ill-posed since contaminant fate and transport processes are irreversible, therefore, the inversion solutions are significantly sensitive to the errors in the observation data. In the last 20 years many approaches, including optimization methods, geostatistical and probabilistic modeling and deterministic direct methods, have been developed to address the groundwater contaminant source identification problem. Wagner (1992) used nonlinear maximum likelihood estimation to perform simultaneous model parameter estimation and source characterization. Mahar and Datta (1997) extended nonlinear optimization with embedding techniques to source identification as well as design of groundwater quality monitoring networks. Aral et al. (2001) proposed a progressive genetic algorithm for the solution of the nonlinear optimization model for identifying contaminant source location and release history in aquifers. Singh et al. (2004) developed artificial neural networks based on a back propagation algorithm to estimate unknown pollutant sources. Huang et al. (2008) applied the conjugate gradient method to determine the unknown space and time dependent contaminant source. Bagtzoglou et al. (1992) modeled the reversed time transport equation using the random walk particle method, which

1612

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

successfully assessed the relative importance of each potential source. Snodgrass and Kitanidis (1997) combined Bayesian theory with geostatistical techniques to derive the best estimate of the release history. Michalak and Kitanidis (2004) extended the geostatistical approach with the adjoint state method to allow for the recovery of the antecedent contaminant distribution. Shlomi and Michalak (2007) improved the geostatistical tools by combining geostatistical kriging and inverse modeling approaches to estimate the spatial distribution of contaminant plumes; the minimum relative entropy inversion technique was introduced by Woodbury and Ulrych (1996) to recover the release and evolution history of a plume. The deterministic method has been used to directly solve the governing equations backward for the reconstruction of the plume’s spatial and temporal history. For example, Skaggs and Kabala (1994) attempted to find the history of the contaminant plume using the Tikhonov regularization method; Skaggs and Kabala (1995) employed the quasi-reversibility to solve the equations with a negative time step; Atmadja and Bagtzoglou (2001) enhanced the backward beam equation method to solve the advection-dispersion equation within the contaminant source identification context. Recently the radial basis collocation method has emerged as a competitive alternative to discretization methods for the inverse problems. As a domain-type numerical method, the RBCM is applicable to practical applications. Hon and Wu (2000) utilized harmonic functions with shift invariability instead of the radial basis functions (RBFs) in the Hermite-Birkhoff collocation method to determine an unknown boundary. Cheng and Cabral (2005) chose the inverse multiquadric (MQ) as the RBFs to solve several types of ill-posed inverse boundary problems for the Laplace equation. Li and Mao (2011) proposed a RBCM method based on the GST-MQ function to solve the inverse heat conduction problems. The present work extends the GST-MQ collocation method to the identification of groundwater contaminant source. This paper is organized as follows. Section 2 describes the contaminant transport processes as well as the source identification problems. Section 3 presents the GST-MQ function approximation and the formulation of the LS-RBCM scheme. Section 4 investigates the sensitivity of the GST-MQ solution with respect to the calculation parameters, observation data and model uncertainty involved in contaminant source identification problems. Section 5 includes application examples in one- and two-dimensional porous media to demonstrate the efficiency of the proposed method. Finally, conclusions are given in Section 6.

vCðx; tÞ v2 Cðx; tÞ vCðx; tÞ þU D þ KCðx; tÞ vt vx vx2 ¼ 0; 0 < x < L; 0 < t < T;

F½Cðx; tÞ ¼

where F is the advection-dispersion-reaction operator; C(x, t) is the contaminant concentration at the position x at time t; D is the hydrodynamic diffusion coefficient; U is the Darcy flow velocity; K is the decay parameter. As shown in Fig. 1, if the source is known at the boundary x ¼ 0, the plume, which travels downstream in the porous media due to the advection and dispersion in the field (0, L), can be predicted subject to the initial and boundary conditions

Cðx; 0Þ ¼ 0; 0 < x < L;

(2)

Cð0; tÞ ¼ sðtÞ; 0 < t < T;

(3)

CðL; tÞ ¼ 0; 0 < t < T;

(4)

where s(t) is the release history at the location x ¼ 0 from which the plume originated, and the concentration of the cross-section at the sufficiently far field x ¼ L is assumed to be zero. If the transport parameters are all prescribed, it is a direct problem that can be easily resolved using conventional procedures such as the finite difference method or the finite element method. When it is defined in the semi-infinite domain, L/N, the analytical solution (Jury and Roth, 1990) of Eqs. (1)e(4) is expressed in the integral form as

Zt Cðx; tÞ ¼

sðsÞvðx; t  sÞds;

0

in which, vðx; tÞ ¼

x pffiffiffiffiffiffiffiffiffiffiffi exp 2 pDt 3

2.1. Advection-dispersion-reaction equation The processes affecting the movement of chemical contaminants should be understood for predicting contaminant concentrations in groundwater. The significant transport process is advection, which is the transport of the contaminant caused by the flow of subsurface water. Mechanical dispersion causes the contaminant to spread out, while molecular diffusion causes the contaminant to move from high concentration to low concentration. The contaminant may also react with other compounds in the soil grains. The processes of advection, hydrodynamic dispersion, and chemical reaction may be described by the advection-dispersion-reaction equation. Taking a one-dimensional homogenous transport for example, the governing equation is expressed as

(5) " 

ðx  UtÞ2  Kt 4Dt

#

2.2. Contaminant source identification problem When a contamination accident happens in groundwater, the source information s(t) is unknown, and the contaminant plume has evolved a distance before being discovered. It is supposed that the measurement of time series concentration at observation sites or/and spatial distribution of the plume at a certain time are given, a contaminant source identification problem arises to reconstruct the missing source information s(t). For example, as illustrated in Fig. 1, the measured concentrations q(t) at observation site x ¼ l

Cðl; tÞ ¼ qðtÞ; 0 < t < T 2. Mathematical modeling for the contaminant source identification

ð1Þ

(6)

is given, Eq. (1), together with the initial boundary (2), far field condition (4) and the measurement (6), produces a solvable inverse boundary value problem. The study will focus on this type of inverse problems. In cases where neither the source location nor the release function are available, it is necessary to provide following additional measurement information at another observation site x ¼ l* to identify the contaminant source.

Fig. 1. Definition sketch of source identification in one-dimensional porous media.

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

  C l* ; t ¼ q1 ðtÞ; 0 < t < T:

(7)

If the initial condition (2) is unknown, the problem can be categorized as time inversion problem. It is necessary to provide following additional spatial distribution of the plume at a certain time t ¼ t* to recover the initial condition (2).

  C x; t * ¼ q2 ðxÞ; 0 < x < L:

(8)

3. Formulation of GST-MQ collocation method 3.1. Global space-time multiquadric function Since Kansa (1990) successfully used the RBCM method based on Hardy (1971) to solve partial differential equations, the approach has been applied to solve many practical problems (Zerroukat et al.,1998; Hon et al., 1999). Lee et al. (2009) gave the theoretical solvability for Kansa’s method. And if the method is modified, it becomes possible to give error bounds of the numerical solution (Ling and Schaback, 2008). Due to its truly meshless property and independence of spatial dimension, the RBCM has been directly applied to the solution of ill-posed problems in recent years. Cheng and Cabral (2005) found that the RBCM can be applied to the inverse problems using the same numerical procedure as the direct problems, which involves a singlestep solution of a system of linear algebraic equations. The method is much more efficient than the iterative optimization techniques for solving ill-posed problems using conventional finite difference method or finite element method. The RBCM utilizes a RBF to interpolate the numerical solutions at scattered collocation points and does not require the generation of a grid, which is a non-trivial task for complicated domains for the grid-based methods. In comparison with other RBFs, such as thin plate spline, the MQ interpolant has exponential convergence with minimal semi-norm errors (Madych and Nelson, 1990; Wu and Schaback, 1993). Therefore, Kansa’s method can provide a highly accurate numerical solution with a relatively small number of collocation points (Hon et al., 1999). Li et al. (2003) investigated the elliptic problems with Dirichlet boundary conditions, and reported that the accuracy of the MQ solution outperforms the finite element method by about 1000 folds. In this paper, we follow the idea of global MQ collocation method in Li and Mao (2011) to identify the time dependent contaminant source release in groundwater systems. Taking a onedimensional problem for example, we consider N and M nodes uniformly distributed in the spatial domain (0, L) and time interval (0, T), respectively. The GST-MQ function is given by

4j ðx;tÞ ¼ 1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2  2 c2 þ x  xj þl t  tj ; j ¼ 1;2;.; N  M;

(9)

where xj, tj are the position and time of collocation points, respectively; N  M is the total number of the collocation points distributed in the global space-time domain; l is the scaling factor that normalizes time to be a pseudo-space; c is the shape parameter. In RBCM method, it is assumed that the concentration C(x, t) is approximated by the following linear combination of the GST-MQ function

Cðx; tÞ ¼

NM X

4j ðx; tÞaj ;

(10)

j¼1

where {aj} represents the coefficients to be determined by the collocation method. If we denote the vectors as

h

Fðx; tÞ ¼ 4j ðx; tÞ

i

1613

 ; a ¼ aj NM;

;

(11)

Cðx; tÞ ¼ Fðx; tÞa; G½Cðx; tÞ ¼ G½Fðx; tÞa;

(12)

1; NM

1

we have

where G is a partial differential operator that acts on the GST-MQ function vector F with respect to x or t. Apparently, the behavior of the solution C(x,t) is insensitive to the direction of time. As a result, the time dependent ill-posed problem can be solved as a direct problem (Li and Mao, 2011). The convergence property of inverse MQ function is dependent on the magnitude of shape parameter c. Early works recommended the choice of the value of c to be proportional to the mean distance dave from each data point to its nearest neighbor, and Hardy (1971) used the formula c ¼ 0.815dave. Carlson and Foley (1991) suggested that a small value of c should be used if the function to be interpolated varies rapidly, but a large value be used if it is a smooth function. Following this idea, Kansa and Carlson (1992) observed that the computational accuracy of the interpolant can be improved by varying the shape parameter given by c2j ¼ c2min (c2max/c2min)(j1/N1), where cmin and cmax are preset. Golberg et al. (1996) used statistical method of cross validation to determine the optimal shape parameter. Unfortunately, it is still an open question to choose the optimal shape parameter. In GST-MQ function, the shape parameter c is adjusted as follows (Li and Mao, 2011)

c ¼ c0

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðDxÞ2 þl ðDtÞ2 ;

(13)

where Δx, Δt are the space and time step; c0 is a constant.

3.2. Least-square-based radial basis collocation method The conventional RBCM requires the variable to satisfy the governing equation and boundary conditions exactly at all nodes, which may cause the instability for the solution (Kee et al., 2008). Li and Mao (2011) used the LS-RBCM method (Mao and Li, 2010) to provide a kind of ‘relaxation’ effect to search for a best-fit solution of the system based on two sets of collocation points: One is satisfied with the governing equation and another is for the given conditions. We follow the numerical procedures in Li and Mao (2011) to identify the contaminant source problem. Firstly, the concentration C(x, t) approximated by the GST-MQ function (10) is required to satisfy the governing Eq. (1), and collocating on all collocation points in the entire space-time domain [0, L]  [0, T], then we can obtain a square linear system with N  M algebraic equations

F½Fðxi ; ti Þa ¼ 0; i ¼ 1; .; N  M:

(14)

Secondly, the concentration (10) is substituted into the definite conditions (2) and (4) and measured conditions (6), and collocating on the corresponding collocation points, then we have

Fðxi ; 0Þa ¼ 0; i ¼ N  M þ 1; .; N  M þ N;

(15)

FðL; ti Þa ¼ 0; i ¼ N  M þ N þ 1; .; N  M þ N þ M;

(16)

Fðl;ti Þa ¼ qðti Þ; i ¼ N  M þ N þ M þ 1;.; N  M þ N þ M þ W; (17) where W is the number of measurement data.

1614

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

Then Eqs. (14)e(17) constitute an over-determined system of linear algebraic equations with N  M unknowns and (N  M þ N þ M þ W) equations, which can be expressed in the matrix form as

Aa ¼ b

(18)

where A is coefficient matrix, which includes the governing equation and definite conditions; b is the vector containing the prescribed boundary conditions and observation data; a is the vector of unknowns aj. If the additional observed information (7) or (8) is known for input, by collocating on the corresponding collocation points, we have





¼ N  M þ N þ M þ W þ 1;.; N  M þ N þ M þ W þ W1 ; (19) or



F xi ;t * a ¼ q2 ðxi Þ; i ¼ N M þN þM þW þ1;.; N M þN þM þW þW2 ; (20) where W1, or W2 is the number of supplementary measurement. They can also be included into the linear algebraic equations system (18). Finally, the unknown coefficients a in Eq. (18) can be determined by minimizing L2 norm using the least-square method and the singular value decomposition with rank reduction is introduced to stabilize the solution of the following over-determined system (Mao and Li, 2010).

n

o

min Aa  bk2

(21)

Then it can provide a best-fit approximation to the solution of the system. Once the solution of coefficient vector a is obtained, the contaminant spatial plume distributions C(x, t) at any time within the porous media can be determined by substituting a into Eq. (10), and the contaminant release history can be estimated as

sðtÞ ¼ Fð0; tÞa:

4. Sensitivity analysis To increase confidence in the GST-MQ collocation method for contaminant source identification problems, a sensitivity and convergence analysis is performed on the calculation parameters, observation data and model uncertainty. We suppose a contamination accident in the porous media with D ¼ 10, U ¼ 1 and K ¼ 0.001, as illustrated in Fig. 1. The source release s(t) is assumed to be a simple smooth function

" sðtÞ ¼ CS exp

F l* ;ti a ¼ q1 ðti Þ; i



identification problems in the multi-dimensional porous media using the same procedure.

(22)

Since the solution of the inverse problem is calculated simultaneously in global space and time domain [0, L]  [0, T], this method can be generally extended to solve other contaminant source



# ðt  0:4TÞ2 ; 0
(23)

where CS is the peak value of concentration taken to be CS ¼ 100. The contaminant plume has an evolution history lasting for T ¼ 200 but does not exceed the downstream L ¼ 400. The observation data (6) is based on the analytical solution (5). The deviation between the true and estimated release function s(t) is measured by the relative root mean square error

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi,vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u PS PS u1 X u1 X ES ¼ t ½sest ðtk Þ  strue ðtk Þ2 t ½strue ðtk Þ2 PS PS k¼1

(24)

k¼1

where {tk} is the set of test points; PS is total number of test points; sest is the estimated solution obtained from Eq. (22); strue is the true solution given by Eq.(23). We take PS ¼ 200 and {tk} is uniformly distributed in [0, T]. 4.1. Calculation parameters In this section, the sensitivity of the inverse solutions with respect to the parameters, such as the shape parameter c0, scaling factor l, time step Δt and collocation point number N, is investigated, and the accuracy and stability of the method are also discussed. It is assumed that W ¼ 41 error-free measurement are given at the observed site l ¼ 50. Li and Mao (2011) investigated the convergence of the GST-MQ solution with respect to the shape parameter c0 and the scaling factor l. The errors change with c0 and l in V or U-shaped distribution in the inverse heat conduction problems. Fig. 2a shows that the errors ES change with c0 for l ¼ 2.4, 6t ¼ 5, N ¼ 41 and l for c0 ¼ 5.8, 6t ¼ 5, N ¼ 41 in V-shaped distribution. The optimal value c0 is about 5.0 w 6.0, and the optimal l is about 2.0. The above optimum values c0 and l are the local optimum ones since they are obtained by taking the other constant. Fig. 2b illustrates that the

Fig. 2. (a) The error ES with respect to c0 (dCd) and l (d:d); (b)The ES distribution with respect to c0 and l.

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

1615

Fig. 3. (a) The error ES with respect to Δt (dCd) and N (d:d); (b) The ES distribution with respect to Δt and N.

error ES varies in a cone-shaped distribution with respect to c0 and l. The global optimal value ranges c0 ¼ 3.3 w 5.0 and l ¼ 3.0 w 6.0, as shown in the core of the cone in Fig. 2b. When either of c0 and l is very small or large, the ES goes up very quickly. When c0 is very small, the interpolation between two neighboring collocation points is nearly linear and so the surface fitted to the data contains sharp corners. As it increases, more collocation points are effectively involved in the interpolation and then the sharp corners spread out to form a smooth surface. When it is too large and reaches a critical value, the basis function tends to be a constant, and the matrix A becomes very highly ill-conditioned (Hon et al., 1999). When l is very small or large, the solution is almost insensitive to the temporal or spatial collocation points. Therefore the GST-MQ function fails to give a reasonable approximation. It is not easy to find the optimum values c0 and l for the numerical solution in practice since we never know the true solution. Fortunately, the reasonable numerical solution exists within a wide range of the parameters c0 and l, for example, c0, l ¼ 0.1 w 100, although it may not be the “best” one. In practice, we can use the statistical method of cross validation to obtain the optimal parameter for the “best” solution as shown in Golberg et al. (1996). The time step in the GST-MQ method has lost the implications in time marching numerical procedure, due to the behavior of the solution insensitive to the direction of time. Li and Mao (2011) showed that there is a maximum time step to guarantee the convergence of the solution. Fig. 3a investigates the error ES with

a 75

4.2. Observation data Since the contaminant source identification problem is an illposed inverse problem, the accuracy of the estimated release history will be influenced by the measurement quality and quantity. In practical situations, the measurement inevitably contains noise due to the technical difficulties associated with measurement. In this section, the sensitivity of factors associated with observations including noise level e, observation number W, observation site l and additional observation is investigated to show the stability and convergence of the method.

b

True e=0.001 e=0.01 e=0.04

100

respect to the time step Δt for c0 ¼ 5.8, l ¼ 2.4, N ¼ 41. When Δt increases from 1.0 to 4.0, the error ES goes down from 0.0607 to 0.0120; When Δt > 4.0, the error goes up almost proportionally. The computational accuracy and efficiency of the GST-MQ method also depends on the mesh size. Fig. 3a illustrates the error ES with respect to the collocation point number N for c0 ¼ 5.8, l ¼ 2.4, Δt ¼ 5. The estimated solution converges rapidly when N increases from 11 to 61 and the error ES keeps stable when N > 61. The pair-wise sensitivity of time step 6t and collocation point number N illustrates that the more collocation points and smaller time step used, the more reasonable estimations of release history can be obtained, as shown in Fig. 3b. When N > 51 and Δt < 6, the error ES is less than 0.01. The result demonstrates that the method can achieve a good numerical solution at a relatively small number of collocation points because of the exponential convergence property of MQ approximation.

s(t)

Es

50

10

10

0

-1

25

0 -2

-25 0

40

80

t

120

160

200

10 -5 10

10

-4

10 e

Fig. 4. (a) The estimated s(t) with respect to e; (b) The error ES with respect to e.

-3

10

-2

10

-1

1616

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

a

b 10

100

0

True W=10 W=20 W=40

75

Es

s(t)

50 10

-1

25

0

-25 0

40

80

t

120

160

200

10

-2

10

1

10

W

2

Fig. 5. (a) The estimated s(t) with respect to W; (b) The error ES with respect to W.

We take the calculation parameters as c0 ¼ 5.8, l ¼ 2.4, Δt ¼ 5, N¼41. The noisy observations are simulated by

qj;noisy ¼ qj;true þ e$uj $qj;true ; j ¼ 1; .; W

(25)

where qj,noisy is the measurement with noise; qj,true is true data; e is the noise level; uj is the random number according to standard normal distribution. Fig. 4 investigates the effect of measurement error level e on the estimated release history s(t) for W ¼ 40, l ¼ 50. When the measurement error is very small, say e¼0.001, 0.01, the estimated solutions match the true ones very well, including the shape and peak. The error ES is as small as 0.0181, 0.0875, respectively. When the error level e ¼ 0.04, the oscillations appear after t ¼ 100 in the estimated release history. The error ES goes to 0.3395. There exists a critical point e0 (¼0.001) as shown in Fig. 4b. The same phenomena can be found in Li and Mao (2011). When e is smaller than e0, the solutions are not sensitive to e and the accuracy depends on the truncated error of the method; however, when e is greater than e0, the accuracy depends linearly on the input error level. It shows the LS-RBCM possesses good stability of the solution to resist the large noise from the input. To guarantee the convergence of the source identification problem, the observation data should be sufficiently provided. Li and Mao (2011) showed that there exists a maximum observation number W to guarantee the full convergence of the solution, which means that too few measurement data are difficult to pin the

b

a True l=50 l=75 l=100

100

10

10

Es

75

s(t)

solution down, but more information do not improve the solution significantly. The same phenomena can be found in Fig. 5. The error ES is quite large when W < 30, but tends to be convergent as W increases. When W > 30, the solution is fully convergent. The error ES of the solutions for W ¼ 10, 20, 40 are 0.3395, 0.0562, 0.0145 for e ¼ 0.001, l ¼ 50, respectively. Mathematically, the suitable observation condition should be provided to make the inverse problem identifiable. If the condition is too limited, it will be an under-determined problem; and if the condition is provided sufficiently, it will be an over-determined problem. Fig. 5b shows that, although the observation condition is too limited to make the inverse problem identifiable when W < 30, the RBCM can obtain a reasonable solution. This is more important in real world applications, because available data are always limited in practice. The optimal observation site may help to reconstruct the problem. Li and Mao (2011) showed that when the observation site moves closer to the source location, the better agreement can be obtained. As shown in Fig. 6, when the observation site l < 45, the error ES keeps a level below 0.0076 and the error is less sensitive to the observation site; When l > 45, the error ES becomes extremely sensitive and rises quickly. The errors ES of the estimated history are 0.0181, 0.0619 and 0.5336 for observation site l ¼ 50, 75, 100 when e ¼ 0.001, W ¼ 40, respectively. The estimated history recovers the true one reasonably except t > 160 for l ¼ 100. Under some circumstances, it is necessary to provide additional observation data to improve the recovery of the plume’s release

10

50

25

10

1

0

-1

-2

0 0

40

80

t

120

160

200

10

-3

10

1

Fig. 6. (a) The estimated s(t) with respect to l; (b) The error ES with respect to l.

l

10

2

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

t

1

10

*

where P represents the model parameters D, U and K; Pj,err, Pj,exact is the erroneous and exact value at each collocation point on the global domain, respectively; r is the error level. The parameters related to calculation procedures and observation data are specified as c0 ¼ 5.8, l ¼ 2.4, Δt ¼ 5, N ¼ 41, e ¼ 0.001, W ¼ 40 and l ¼ 50. Fig. 8 shows the comparison of the estimated release function s(t) with respect to the error level r in D, U, and K. For parameter D, when the error level r ¼ 0.01, 0.1, the estimated histories coincide with the true history very well and the error ES is as low as 0.0219, 0.0712; When r ¼ 0.4, the deviations of estimated history occur near t ¼ 60, 100 and 200, and the ES is 0.1571. For parameter U, when r ¼ 0.04, 0.1, the numerical history under-estimates the peak, and the deviation increases very quickly after t ¼ 180. For parameter K, at larger error level, say r ¼ 1.0, the numerical history can get the peak very well, but starts to deviate near t ¼ 200. The performance of different errors in D, U, and K is presented in Fig. 9. It is observed that the same error level of the parameters D, U, and K may have different influence to the estimated release history. For parameters K, the ES errors are not sensitive to the error level when r < 0.06; whereas the ES error increase with the error level increase when r > 0.06, and the gradient is still very small. For parameter D and U, the ES errors are not sensitive to the error level when r < 0.004; whereas the ES errors increase quickly with the parameter error levels increase when r > 0.004, the gradient is large. As show in Fig. 9, obviously, the source release history is very sensitive to the uncertainty of the parameter U (advection), but less sensitive to the parameter K (reaction).

2

10

ES0 -1

Es

10

-2

10

1

10

2

10

*

l

Fig. 7. The error ES with respect to l* (dCd) and t* (d:d).

history. If the additional measurement (7) is included for input at another observation site l*, the error ES with respect to observation location l* is shown in Fig. 7. For comparison, the error ES0¼0.1712 of the estimated solutions without the assistance of additional observations is plotted for e ¼ 0.02, W ¼ 40, l ¼ 50. It is observed that, when l* < 50, the estimated release history is improved greatly; however, the additional measurement from the site l* > 50 is trivial to improve the estimated release history. If the measured spatial distribution of the plume (8) at time t* is given for the additional information, the error ES with respect to observation time t* is shown in Fig. 7. It is found that, the error is almost constant and close to ES0. The additional observation of spatial distribution improves the estimated solutions slightly.

5. Application examples In this section, one- and two-dimensional source identification problems are chosen as application examples to illustrate the efficiency of the proposed GST-MQ method. Following relative root mean square error between the true and estimated contaminant distributions C(x, T0) is used to evaluate the accuracy of the method

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi,vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u PC PC u1 X u1 X 2 t t EC ¼ ½Cest ðxk ;T0 Þ  Ctrue ðxk ;T0 Þ ½Ctrue ðxk ;T0 Þ2 PC PC

4.3. Model uncertainty The behavior of contaminant in groundwater is mainly described by the processes including advection, hydrodynamic dispersion and chemical reaction, so the reliability of source identification should depend on the accuracy of these parameters. This section examines the sensitivity of the model uncertainty of the diffusion coefficient D, Darcy velocity U, and decay parameter K. The error in the model parameters is simulated by

Pj;err ¼ Pj;exact þ r$uj $Pj;exact ; j ¼ 1; .; N  M

b

100 True

(27) where PC is total number of the test points; Cest is the estimated concentrations; Ctrue is the true concentrations. We take Pc¼121, 101  101 for the one- and two-dimensional examples, respectively. The test points are uniformly distributed in the spatial domain.

(26)

s(t)

s(t)

c

True r=0.01 r=0.04 r=0.1

75

50

100 True r=0.1 r=0.4 r=1

75

50

50

25

25

25 0

0 0

k¼1

100

r=0.01 r=0.1 r=0.4

75

k¼1

s(t)

a

1617

0 40

80

120 t

160

200

0

40

80

120

160

200

0

t

Fig. 8. The estimated release history s(t) with respect to the error r in D (a), U (b) and K (c).

40

80

120 t

160

200

1618

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

10

amount of the released contaminant assigned by AS¼1000. We simulate the plume history by T ¼ 120, and the plume fronts do not reach the end cross-section at L ¼ 240. By substituting Eqs. (28)e(30) into Eq.(5), the exact solutions subject to three impulse boundary conditions are

0

Es

D U K

10

" pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! CS x  t U 2 þ 4KD pffiffiffiffiffiffi C1 ðx; tÞ ¼ erfc 2 2 Dt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!#  Ux x þ t U 2 þ 4KD pffiffiffiffiffiffi erfc ; þ exp D 2 Dt  C Ux C2 ðx; tÞ ¼ S exp f½erfcðB1 þ B3 Þ 2D 2

-1

-2

10 -4 10

10

-3

10 r

-2

10

-1

10

(31)

 erfcðB2 þ B4 ÞSðt  TS ÞexpðB0 Þ þ ½erfcðB1  B3 Þ

0

 erfcðB2  B4 ÞSðt  TS ÞexpðB0 Þg;

ð32Þ

in which,

Fig. 9. The error ES with respect to the error r in D, U, and K.

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U2 x x U2t U 2 ðt  TS Þ þ K ; B1 ¼ pffiffiffiffiffiffi; B2 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi; B3 ¼ þ Kt ; B4 ¼ þ Kðt  TS Þ; 4D 4D 4D 2 Dðt  TS Þ 2 Dt # " AS ðx  UtÞ2 C3 ðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffi exp   Kt 4Dt 4pDt x B0 ¼ pffiffiffiffi D

5.1. One-dimensional examples

The true data are calculated from Eqs. (31)e(33). The noisy measurement data are simulated by the random errors in Eq. (25). The other parameters related to the calculation procedures and input observations are defined as Δt ¼ 5, N ¼ 49, e ¼ 0.02, W ¼ 24 and l¼30. Fig. 10 presents the comparison between the estimated and true release histories s(t) in three release history cases. Fig. 11 shows the comparison between the estimated and true spatial distributions C(x, T0) in the region [0, L/2] at time T0 ¼ 30, 60, 90, 120 for three different release sources. The contaminant release s1(t) is a constant release, so it is relatively easy to estimate. The estimated history converges to the true solution quickly except that there is a big discrepancy at t ¼ 0 as shown in Fig. 10a, which is caused by the discontinuities from the zero initial conditions to the constant release s1(0) ¼ 100. The GSTMQ formulation is continuously differentiable, so the numerical solution keeps the continuity to approximate the discontinuities. The error ES of the release history s1(t) is 0.1049. The estimated spatial distributions C1(x, T0) at T0 ¼ 30, 60, 90, 120 converge to the true solutions very well except the area near the source at T0 ¼ 30, as shown in Fig. 11a. The EC errors for C1(x, T0) are 0.0536, 0.0133, 0.0177 and 0.0185 at T0 ¼ 30, 60, 90, 120, respectively.

Three release history cases are considered in one-dimensional contaminant source identification problems, including constant, intermittent and transient release history, representing different typical types of source release history in practical contamination incidents. It is assumed that the contaminant plume travels in the porous media with D ¼ 1, U ¼ 1, and K ¼ 0.0001, as shown in Fig. 1. Three synthetic source functions are selected as the true contaminant release history

s1 ðtÞ ¼ CS ; 0 < t < T;

(28)

s2 ðtÞ ¼ CS ½SðtÞ  Sðt  TS Þ; 0 < t < T;

(29)

s3 ðtÞ ¼ AS dðtÞ; 0 < t < T;

(30)

where s1(t) is a continuous source with a constant release rate CS ¼ 100; s2(t) is an intermittent source that ceases at time TS ¼ T/3, and S() is the step function; s3(t) is a transient source that releases at t¼0 instantly, d() is the Dirac-delta function, and AS is the total

b

120

80 60

60

40

40

30 Estimated True

20

s 2(t)

80

s 1(t)

Estimated True

100

100

c

120

s 3(t)

a

(33)

10

20 20 0 0

Estimated True 20

40

60

t

80

100

0 120

0

0 20

40

60

80

100

120

t Fig. 10. The estimated release history s1(t) (a), s2(t) (b) and s3(t) (c).

0

20

40

60

t

80

100

120

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

a100

b

80

c

100

C3(x,T0)

C2(x,T0)

C1(x,T0)

60

40

50 40

80

60

1619

30 20

40 10

20 20

0

0 0

20

40

60 x

80

100

120

0

20

40

60 x

80

100

120

0

20

40

60 x

80

100

120

Fig. 11. The estimated concentration distribution C1(x, T0) (a), C2(x, T0) (b) andC3(x, T0) (c) at time T0 ¼ 30(dCd), 60 (d-d), 90 (d:d), 120 (dAd).

Fig. 12. Definition sketch of source identification in two-dimensional porous media.

The contaminant release s2(t) is a square-wave step function. The continuously differentiable GST-MQ solution approximates the release discontinuity at time t¼0 and t¼TS. The estimated release history has reasonable agreement with the true one as shown in Fig. 10b. The error ES of the release history s2(t) is 0.3796. The estimated spatial distributions C2(x, T0) at T0 ¼ 30, 60, 90, 120 match the true solutions very well except the area near the source at t ¼ 30, as shown in Fig. 11b. The EC errors for C2 (x, T0) are 0.0550, 0.0437, 0.0138 and 0.0278 at T0 ¼ 30, 60, 90, 120, respectively. The contaminant release s3(t) is the most difficult case to deal with because the initial concentration goes to infinity. It is a highly ill-posed problem. The release history is estimated with a finite peak value, greatly under-estimating the true one as shown in Fig. 10c. The error ES expands dramatically to 11.9020. Fortunately, the estimated spatial distributions C3(x, T0) have a good agreement with the true solutions, as shown in Fig. 11c. The EC errors for spatial distributions C3 (x, T0) are 0.0978, 0.0927, 0.0916, 0.1109 at T0¼30, 60, 90, 120, respectively.

homogeneous media with thickness H ¼ 1 and length L ¼ 1 and it brings with a slight background contaminant C0 ¼ 2. As shown in Fig. 12, it involves the identification of the unknown space and time dependent contaminant source s(x, t) in a uniform width (L2eL1) located at the upper boundary. It is a non-point contaminant source identification problem. The total simulation time is chosen as T ¼ 0.2. The dimensionless formulation of the mass transfer problem can be expressed as

Rd

vCðx; y; tÞ v2 Cðx; y; tÞ v2 Cðx; y; tÞ vCðx; y; tÞ 2 þ A D  Vaq ¼ Dx y r vt vx vx2 vy2  Rd zCðx; y; tÞ; 0 < x < L; 0 < y < H; 0 < t < T; ð34Þ

Cðx; y; 0Þ ¼ 0; 0 < x < L; 0 < y < H;

(35)

Cð0; y; tÞ ¼ C0 ; 0 < y < H; 0 < t < T;

(36)

vCðL; y; tÞ ¼ 0; 0 < y < H; 0 < t < T; vx

(37)

vCðx; 0; tÞ ¼ 0; 0 < x < L; 0 < t < T; vy

(38)

vCðx;H;tÞ sðx;tÞ ½Sðx  L1 Þ  Sðx  L2 Þ; 0 < x < L; 0 < t < T; ¼ vy Ar Dy

(39)

5.2. Two-dimensional example

where Ar is the aspect ratio; z is the decay rate; the retardation factor Rd and the dispersion coefficients in x and y direction, Dx, Dy, are defined, respectively, as

This section presents a 2D dimensionless application example, which is similar to Huang et al. (2008). The flow is horizontal and steady with an average Darcy velocity Vaq ¼ 0.1 in the isotropic,

Rd ¼ 1 þ

rs Kd ; Dx ¼ aL Vaq ; Dy ¼ at Vaq q

Fig. 13. (a) The true release history s(x,t); (b) the estimated release history s(x,t); (c) The residual of estimated history s(x,t).

(40)

1620

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

Fig. 14. The true concentration distribution C(x, y, T0) (a) and estimated distribution C(x, y, T0) (b) at time T0 ¼ 0.10.

where rs is the soil bulk density; Kd is the distribution coefficient; q is the effective porosity; aL, at are the longitudinal and transversal dispersivity, respectively. These media parameter values are chosen as Ar ¼ 20, z ¼ 1.0, rs ¼ 1.8, Kd ¼ 0, q ¼ 0.3, aL ¼ 1 and at ¼ 0.2. The following release history s(x, t) is imposed as contaminant source

sðx; tÞ ¼ 30exp ð6tÞexp

n

o  180½x  ðL1 þ L2 Þ=22 ;

L1 < x < L2 ; 0 < t < T

ð41Þ

If we know the above release history, it is a direct problem of the mass transfer. The numerical solutions, which are obtained from the direct problem using alternating directional implicit method, are assumed to be the true measured concentrations. If the release history (41) is unknown, it arises to be an inverse problem of source release identification. It is supposed that 17 sensors are placed between L1 ¼ 0.2 and L2 ¼ 0.6 and at H0 ¼ 0.2, close to the bottom of the media. The noisy measurement data are simulated by the random error (25). The noise level is chosen as e¼0.10. The two-dimensional application example of groundwater contaminant source identification is solved using the same GSTMQ solution procedure as in one-dimensional examples. Nx  Ny collocation points are uniformly distributed in the spatial domain. The parameters related to calculation procedures and input observations are prescribed by Nx ¼ 21, Ny ¼ 21, Δt ¼ 0.02, c0 ¼ 0.35, l¼2.75, e¼0.10, W ¼ 10. As shown in Fig. 13a, the true contaminant release s(x, t) is a function of space and time, representing a non-point contaminant source. The estimated release history from GST-MQ collocation

method is illustrated in Fig. 13b. The estimated source release s(x, t) precisely recovers the source location in the region (0.2, 0.6), but under-estimates the peak of the concentration history by 67% as shown in Fig. 13c, which is similar as in one-dimensional cases. The error ES of the estimated release history is 0.4547. The true and estimated plume evolution distribution C(x, y, T0) at T0 ¼ 0.10, 0.20 are plotted in Figs. 14 and 15, respectively. It showed that the estimated plumes have good agreements with the true ones, except slight discrepancy in the estimated peak for every snapshot time. The EC errors are 0.1092, 0.0801 for T0¼0.10, 0.20, respectively. The numerical results in one- and two-dimensional cases for groundwater contaminant source identification show that, although the errors ES of the release history may be relatively large, especially in the cases of intermittent and transient release history, the estimated spatial distributions match the true ones very well. The error from the release history does not spread out along the downstream. The GST-MQ method possesses the advantage to solve the inverse problem, which gives the numerical plume distribution and the source release history simultaneously. All the computational codes are written in Matlab, and the singular value decomposition is adopted based on the Matlab code. Since the RBCM involves a single-step solution of a linear algebraic equations system, most computational CPU time is consumed for the system. The scale of the matrix system for the 1D problems is 1324  1225 and the CPU time is about 40 s on a PC platform with 3.16 GHz CPU and 4.00 GB RAM; The scale of the matrix for the 2D problem is 6232  4851 and the CPU time is about 3200 s on the same machine.

Fig. 15. The true concentration distribution C(x, y, T0) (a) and estimated distribution C(x, y, T0) (b) at time T0 ¼ 0.20.

Z. Li, X.-z. Mao / Environmental Modelling & Software 26 (2011) 1611e1621

6. Conclusions In this paper, we present an inverse model based on the global MQ collocation method for the identification of groundwater contaminant source. This GST-MQ is constructed by incorporating the time dimension into the inverse MQ function to resolve the problem in the entire space-time domain. Groundwater contaminant source identification problem is approximated as an overdetermined linear system using the LS-RBCM method and the SVD with rank reduction is introduced to stabilize the solution of the linear system. The behavior of the solution of the proposed MQ-RBF is insensitive to the direction of time, and, as a result, the time dependent ill-posed problem can be solved as a direct one. A sensitivity analysis with respect to the calculation parameters, observation data and model uncertainty is investigated to increase confidence in the GST-MQ solutions. The computational results show that the method possesses good performance to resist the noisy measurement data and model uncertainty. The error ES varies in a cone-shaped distribution with respect to the shape parameter c0 and scaling factor l. The estimated release history is very sensitive to the uncertainty of the Darcy velocity U when the error level r > 0.004, but slightly sensitive to the decay parameter K. One- and two-dimensional application examples demonstrate that the proposed method poses the meshfree property and allows for the spatial plume distribution and temporal history as well as direct identification of contaminant source with efficient methodology. Although the errors ES of the estimated release history may be relatively large, the estimated spatial distributions match the true ones very well. The error of the estimated release history does not spread out along the downstream. It provides a robust deterministic direct tool for identification of groundwater contaminant source from observed concentration. Acknowledgements The project is sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, China (The Project-sponsored by SRF for ROCS, SEM), and National Water Grant (No. 2008ZX07423-001). References Aral, M.M., Guan, J., Maslia, M.L., 2001. Identification of contaminant source location and release history in aquifers. Journal of Hydrological Engineering 6 (3), 225e234. Atmadja, J., Bagtzoglou, A.C., 2001. Pollution source identification in heterogeneous porous media. Water Resources Research 37 (8), 2113e2125. Bagtzoglou, A.C., Dougherty, D.E., Tompson, A.F.B., 1992. Application of particle methods to reliable identification of groundwater pollution sources. Water Resources Management 6 (1), 15e23. Carlson, R.E., Foley, T.A., 1991. The parameter R2 in multiquadric interpolation. Computers & Mathematics with Applications 21 (9), 29e42. Cheng, A.H.-D., Cabral, J.J.S.P., 2005. Direct solution of ill-posed boundary value problems by radial basis function collocation method. International Journal for Numerical Methods in Engineering 64 (1), 45e64. Commission on Geosciences, Environment and Resources, 1993. Ground water vulnerability assessment: predicting relative contamination potential under conditions of uncertainty. National Academies Press, Washington, D.C. Golberg, M.A., Chen, C.S., Karur, S.R., 1996. Improved multiquadric approximation for partial differential equations. Engineering Analysis with Boundary Elements 18 (1), 9e17. Golosov, S., Kirillin, G., 2010. A parameterized model of heat storage by lake sediments. Environmental Modelling & Software 25 (6), 793e801. Hardy, R.L., 1971. Multiquadric equations of topography and other irregular surfaces. Journal of Geophysical Research 76 (8), 1905e1915.

1621

Hon, Y.C., Cheung, K.F., Mao, X.Z., Kansa, E.J., 1999. Multiquadric solution for shallow water equations. Journal of Hydraulic Engineering 125 (5), 524e533. Hon, Y.C., Wu, Z., 2000. A numerical computation for inverse boundary determination problem. Engineering Analysis with Boundary Elements 24, 599e606. Huang, C.H., Li, J.X., Kim, S., 2008. An inverse problem in estimating the strength of contaminant source for groundwater systems. Applied Mathematical Modeling 32 (4), 417e431. Jury, W.A., Roth, K., 1990. Transfer functions and solute movement through soil. Birkhauser Verlag, Boston. Kansa, E.J., 1990. Multiquadrics-a scattered data approximation scheme with applications to computational fluid dynamics-II. Computers and Mathematics with Applications 19, 147e161. Kansa, E.J., Carlson, R.E., 1992. Improved accuracy of multiquadric interpolation using variable shape parameters. Computers and Mathematics with Applications 24 (12), 99e120. Keats, A., Yee, E., Lien, F.-S., 2010. Information-driven receptor placement for contaminant source determination. Environmental Modelling & Software 25 (9), 1000e1013. Kee, B.B.T., Liu, G.R., Lu, C., 2008. A least-square radial point collocation method for adaptive analysis in linear elasticity. Engineering Analysis with Boundary Elements 32 (6), 440e460. Lee, C.F., Ling, L., Schaback, R., 2009. On convergent numerical algorithms for unsymmetric collocation. Advances in Computational Mathematics 30 (4), 339e354. Li, J., Cheng, A.H.-D., Chen, C.S., 2003. A comparison of efficiency and error convergence of multiquadric collocation method and finite element method. Engineering Analysis with Boundary Elements 27, 251e257. Li, Z., Mao, X.Z., 2011. Global space-time multiquadric method for inverse heat conduction problems. International Journal for Numerical Methods in Engineering 85 (3), 355e379. doi:10.1002/nme.2975. Ling, L., Schaback, R., 2008. Stable and convergent unsymmetric meshless collocation methods. SIMA Journal of Numerical Analysis 46 (3), 1097e1115. Madych, W.R., Nelson, S.A., 1990. Multivariate interpolation and conditionally positive definite functions II. Mathematics of Computation 54 (189), 211e230. Mahar, P.S., Datta, B., 1997. Optimal monitoring network and ground-water pollution source identification. Journal of Water Resources Planning and Management 123 (4), 199e207. Mao, X.Z., Li, Z., 2010. Least square-based radial basis collocation method for solving inverse problems of laplace equation from noisy data. International Journal for Numerical Methods in Engineering 84 (1), 1e26. doi:10.1002/nme.2880. Michalak, A.M., Kitanidis, P.K., 2004. Estimation of historical groundwater contaminant distribution using the adjoint state method applied to geostatistical inverse modeling. Water Resources Research 40 (8), W08302. Peña, P.A., Hernández, J.A., Toro, M.V., 2010. Computational evolutionary inverse lagrangian puff model. Environmental Modelling & Software 25 (12), 1890e1893. Rodríguez, J.A., Martos, J.C., 2010. SIPAR_ID: freeware for surface irrigation parameter identification. Environmental Modelling & Software 25 (11), 1487e1488. Shlomi, S., Michalak, A.M., 2007. A geostatistical framework for incorporating transport information in estimating the distribution of a groundwater contaminant plume. Water Resources Research 43 (3), W03412. Shrestha, S., Kazama, F., Newham, L.T.H., 2008. A framework for estimating pollutant export coefficients from long-term in-stream water quality monitoring data. Environmental Modelling & Software 23 (2), 182e194. Singh, R.M., Datta, B., Jain, A., 2004. Identification of unknown groundwater pollution sources using artificial neural networks. Journal of Water Resources Planning and Management 130 (6), 506e514. Skaggs, T.H., Kabala, Z.J., 1994. Recovering the release history of a groundwater contaminant. Water Resource Research 30 (1), 71e79. Skaggs, T.H., Kabala, Z.J., 1995. Recovering the history of a groundwater contaminant plume: method of quasi-reversibility. Water Resources Research 31 (11), 2669e2673. Snodgrass, M.F., Kitanidis, P.K., 1997. A geostatistical approach to contaminant source identification. Water Resources Research 33 (4), 537e546. U.S. Environmental Protection Agency (EPA), 1999. Safe Drinking Water Act Section 1429: Ground Water Report to Congress. Office of Water, Washington, D. C. Wagner, B.J., 1992. Simultaneously parameter estimation and contaminant source characterization for coupled groundwater flow and contaminant transport modeling. Journal of Hydrology 135, 275e303. Woodbury, A.D., Ulrych, T.J., 1996. Minimum relative entropy inversion: theory and application to recovering the release history of groundwater contaminant. Water Resources Research 32 (9), 2671e2681. Wu, Z., Schaback, R., 1993. Local error estimates for radial basis function interpolation of scattered data. IMA Journal of Numerical Analysis 13 (1), 13e27. Zerroukat, M., Power, H., Chen, C.S., 1998. A numerical method for heat transfer problems using collocation and radial basis functions. International Journal for Numerical Methods in Engineering 42 (7), 1263e1278.