NUMERICAL SIMULATION OF GROUNDWATER FLOW AND POLLUTION TRANSPORT USING THE DUAL RECIPROCITY AND RBF METHOD NUMERICAL SIMULATION OF GROUNDWATER FLOW AND POLLUTION TRANSPORT USING THE DUAL RECIPROCITY AND RBF METHOD

Transport of pollution is strongly influenced by groundwater flow. All numerical models of transport equation must be based on the groundwater flow models. The most of them suffer from numerical diffusion and/or oscillations of solution. Both phenomena lower the of results of these models. Plenty of authors are working on improving the quality of the numerical simulation. This paper presents the connection between a dual reciprocity method used to model groundwater flow and the meshless radial basis function methods used for a transport model. Both methods are one of the latest tools for modeling these phenomena.


Introduction
At the dawn of the use of numerical models, they were mostly used to study a groundwater flow in porous media. Today, our environment suffers more and more from by-products of man's industrial activities. Groundwater is one of the things that have sustained extensive damage in the past decade. Even ordinary events in our everyday life can cause pollution of soil or surface water. Numerical models have become a great tool that can be used to predict spreading of pollutants and assess different remediation projects.
Before we start to prepare the model of pollutant transport we need to know groundwater velocity field because the velocity is one of the main phenomena which influences the spreading of pollution in porous media. Therefore, we focused also on the proper numerical method for simulating the groundwater flow. The boundary element method (BEM) is one of the most suitable methods for this purpose because it allows computing the velocity in any point inside the domain. Unfortunately BEM is not very suitable to solve time-dependent groundwater flow and, therefore, we decided to use the dual reciprocity method (DRM) which was developed like an enhancement of BEM (see [6], [4]).
In this paper we used DRM for computing the velocity field and then the computed field is used as an input to the transport model. The groundwater velocity caused an advection and also increased the dispersion in porous media (see [1]). The advection can cause severe problems in numerical models of pollution transport. These problems were solved by different methods such as Galerkin characteristic method (see [5]) or random walk method [4]. In the last few years, there has been tremendous attention focused on the development of meshless methods to alleviate meshing problems associated with finite element and finite volume methods. The radial basis function method (RBF) was used for simulation of pollutant transport and this method seems to be quite successful.

Basic equations
The planar unsteady groundwater flow with a confined surface is governed by the following differential equation (1) where Φ is the potential of groundwater flow (piezometric head) [m], T x , T y are the transmissivity coefficients [m 2 s Ϫ1 ], S is the storativity coefficient [m Ϫ1 ], Q m are specific discharges of sources or sinks [s Ϫ1 ], and δ is the Dirac function (see e.g. [1]).
The boundary conditions of equation (1) are divided into three basic groups -boundary conditions of the 1st kind (the Dirichlet conditions), where the potential's value on the part of boundary Γ 1 is given, i.e. Φ ϭ Φ 0 , -boundary conditions of the 2nd kind (the Neumann conditions), where the value of the potential's normal derivative (or flux) on the boundary Γ 2 is given, i.e. ( where n x , n y are the components of a outer normal vector which is perpendicular to the boundary, where the given potential's normal derivative (flux) is a function of the potential, i.e.
The initial condition is created by the prescribed value of potential at time t 0 ϭ 0.
The transport of pollution in porous media is an irreversible, non-stationary process that creates a transition zone where the concentration of pollutant continuously changes from the minimal to the maximal value. This phenomenon is called dispersion and it is a macroscopic reflection of a real movement of particles in pores and a reflection of different physical and chemical phenomena.
The simplest governing equation of transport in porous media is (4) where C is concentration of pollutant, D ij is the tensor of dispersion coefficients and v P is the porous velocity of groundwater flow.
Eq. (4) is the basic differential equation of a conservative transport of pollutants in a porous medium that is a transport where the pollutant does not react with the environment and it does not spontaneously decay.
The boundary conditions can be again divided into three groups (see [2]) -Dirichlet boundary condition when the known value of concentration is prescribed on the part of boundary Γ 1 , i.e. C ϭ C 0 , -Neumann boundary condition (impervious boundary) when the flux of pollutant is zero across the boundary Γ 2 , i.e.
-Cauchy boundary condition when the known flux of pollutant is prescribed on the boundary Γ 3 , i.e.
We present the connection of two different methods of solution in this paper. Equation (1) of groundwater flow is solved by a dual reciprocity method; the transport of pollution (Eq. 4) is solved by the local RBF collocation method.

Dual Reciprocity Method
The boundary element method (BEM) is very suitable for solving the groundwater flow equation because it can solve exactly the influence of sources and sinks (see [4]). The dual reciprocity method (DRM) developed from BEM in order to remove major disadvantages connected with the non-zero right side of Eq. (1) (see [6]). Equation (1) should be transformed by the following transformation of coordinates to the Poisson's equation . (8) To solve this equation we use the weighted residuals method (see [8]).
The solution of Eq. (8) can be expressed as the sum of the solution of a homogenous equation and a particular solution Φ . The solution of the homogeneous part of this equation can be expressed using the BEM as (9) where w k are the weighting functions defined for a planar problem as (10) We can now divide the boundary Γ to the set of elements like in the BEM and try to find the approximate solution of potential and flux as , where N i is a set of base functions. We use the simplest form of the base function which is constant in the whole element, i. e. N i ϭ 1 in this paper because these functions are very suitable for solving non-homogeneous areas. Then Eq. (9) can be written as (12) In matrix notation we can write where N is the number of boundary nodes and L is the number of inner points, α i is a set of unknown coefficients and f i are approximation functions. The particular solutions from the sequence must fulfil these equations .
The approximation functions fi can be chosen as The sequence of particular solutions Φ i has the form of a series , and the exterior normal derivative is .
If we substitute relations (17) and (18) to Eq. (15), we obtain the governing equation . ( Then we apply approaches similar to BEM and afterwards we acquire where u n , q n and u nϩ1 , q nϩ1 are the values of variable u, q at time intervals n and nϩ1, respectively. The values u and q are parameters which position the value of u and q, respectively between the time intervals n and nϩ1. If we choose u ϭ q ϭ 0.5 we get the well-known Crank-Nicholson scheme.
The DR method requires a certain number of points to be given inside the domain of solution. The inclusion of boundary conditions is the same as in the BEM because the DR method is a variant of BEM.

Local RBF Method
Radial basis functions (RBF) are initially known as a powerful tool for approximating multivariate functions on a scattered data. Due to their mesh-free nature RBF have received an increasing attention for solving partial differential equations (PDE) of different kinds. The first trial of such exploration was made by Kansa [3].
Full exploitation of the RBF method was constrained by the progressive ill-conditioned coefficient matrix as the number of nodes increases. To remove this difficulty, Shu at al. (see [7]) suggested using the local RBF method in which the approximation is formed by using only several local supporting points.
The unknown concentration function C(x) is approximated in a sub-domain Ω i which forms the neighbourhood or support of a reference point x i by weighted sum of multiquadric functions and polynomials , where λ j and γ j are weights, ϕ(r j ) are the RBF basis functions, and p j is a basis for polynomial space with degree mϪ1, m is the order of ϕ.
Multiquadric functions are one of the most popular radial functions used for this purpose and they are defined as (25) where r j is a distance from a point and ε is a so-called shape factor of multiquadric function. The order m of multiquadric functions is equal to one and, therefore, we need one additional condition to make the interpolation problem well-posed. This condition is (26) The interpolation problem can be written in a matrix form as (27) The differential quadrature (DQ) method was used to approximate derivatives. The derivative at a point x i can be approximated by DQ as and we can compute the weighting coefficient as The same principle can then be used for derivatives of higher order.
The transport equation (4) can be solved in time using the similar time approximation like in Eq. (23) We can now write Eq. (4) as , where we used two operators The solution of PDE is based on a collocation method, i. e., we try to find the solution in the set of nodes placed inside the domain Ω and on boundary Γ. The domain is discretized by N ϭ ϭ N I ϩ N B nodes (N I is the number of interior and N B is the number of boundary nodes). For each node x i , i ϭ 1 … N we assume the N i nodes are used in the support domain to interpolate the unknown variable C. Now we can obtain the weights w 1 and w 2 for operators Y 1 and Y 2 as

Numerical examples
For the validation of the proposed numerical approach, several standard test problems are considered.

Pumped well in circular aquifer
First of all we validate the groundwater flow model. It is tested on the interesting example of a pumped well in a homogeneous aquifer. This example can be solved analytically only for an infinite aquifer. The drawdown can be computed as where T is the transmissivity coefficient, Q is the pumped discharge and u is defined as Here r is the distance from the pumped well, S is coefficient of storativity, and t is time.
The numerical model can not model the infinite domain. We prepared a circular domain with a diameter of 20 m and the Dirichlet boundary condition around. The boundary was divided into 20 elements and 32 internal points (see Fig. 1 Table 1. These results are quite satisfactory and we can see that the DRM method gives very good values also in the vicinity of the pumped well on the contrary to other numerical methods (e. g. finite elements).

Convective-diffusion transport equation
In this section, the solution of the convection-diffusion equation with a constant diffusion coefficient and convection velocity is presented. Only a longitudinal diffusion and convection are considered. The governing equation is (37) The analytical solution of the above transport problem in a onedimensional semi-infinity domain with initial and boundary conditions is given by the following expression (see Bear (1972)) where erfc is an error function.
We prepared a numerical model using a local RBF method. The semi-infinity domain was treated in a numerical model as a rectangular one of 10x4 m. The network had 84 boundary points and 585 internal points. We can compare only short time intervals because of the influence of the rectangular domain boundary at x ϭ 10. The boundary condition C(0,t) ϭ 100, the coefficient of diffusion is 0.1 and velocity v ϭ 1.
Comparison of results of numerical and analytical solution for time t ϭ 2 and t ϭ 4 is presented at Fig. 2.
The influence of the shape factor ε (Eq.25) was also studied and values of relative errors for different values of the shape factors are presented in Fig. 3. The finding of the optimal value is a very interesting and complicated problem and has not been satisfactorily solved yet.

Conclusions
Numerical models became a great tool that can be used to predict spreading of pollutants and assess different remediation projects. There exist plenty of commercial software systems to solve the transport equation but the new numerical methods seems to be more effective for simulation this phenomena. The paper try to present two non-usual methods for solution of coupled problem of groundwater flow and pollution transport in porous media.
The dual reciprocity method and truly meshless local RBF method is used for contaminant transport modelling. The presented numerical results show the simplicity and quite good accuracy of both methods.
The advantage of the dual reciprocity method is the introduction of point sources (wells) into the domain. The sources can be placed arbitrarily inside the domain without having the change the network of elements. The computed values of potential in the point source and its neighbourhood are very precise.
The truly meshless method using the local RBF interpolation connected with collocation solution of the governing dispersion equation is used to solve the transport equation. Numerical results show a good accuracy of this method. The study of influence of the shape factor will continue to improve the usability of the local RBF method in modelling of contaminant transport.