Applications of numerical methods to simulate the movement of contaminants in groundwater.

This paper reviews mathematical models and numerical methods that have been extensively used to simulate the movement of contaminants through the subsurface. The major emphasis is placed on the numerical methods of advection-dominated transport problems and inverse problems. Several mathematical models that are commonly used in field problems are listed. A variety of numerical solutions for three-dimensional models are introduced, including the multiple cell balance method that can be considered a variation of the finite element method. The multiple cell balance method is easy to understand and convenient for solving field problems. When the advection transport dominates the dispersion transport, two kinds of numerical difficulties, overshoot and numerical dispersion, are always involved in solving standard, finite difference methods and finite element methods. To overcome these numerical difficulties, various numerical techniques are developed, such as upstream weighting methods and moving point methods. A complete review of these methods is given and we also mention the problems of parameter identification, reliability analysis, and optimal-experiment design that are absolutely necessary for constructing a practical model.


Introduction
In recent years, mathematical models and numerical methods have been used extensively to simulate the movement of contaminants to groundwater. Many books and review papers have contributed to this topic (1)(2)(3)(4)(5)(6)(7)(8). The distributed parameter model describing solute transport in porous media generally includes an advection-dispersion equation, a groundwater flow equation, and state equations (1). Analytical solutions can be derived only for a few classical models that are not suitable for complex situations normally encountered in the field. Consequently, the development of numerical solutions is required.
During the last two decades, a variety of numerical approaches were presented, such as the finite difference methods (FDM), the finite element methods (FEM), and other alternative methods (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14). However, at least three important problems have not been solved perfectly: * The movement of contaminants in groundwater is an inherently three-dimensional (3-D) process. Unfortunately, most simulations of field problems are now in 2-D because of the high expense of 3-D models. Therefore, 3-D numerical models that can decrease the de-*Environmental Science Center, Shandong University, The People's Republic of China.
Present address: Civil Engineering Department, UCLA, Los Angeles, CA 90024. mands both on human and machine resources need to be developed. * When using the regular numerical methods mentioned above to solve advection-dominated transport problems, two types of numerical errors always occur: the oscillation around the computed profile (overshoot) and the smearing of the concentration front (numerical dispersion). To decrease these numerical errors, many special techniques have been presented. Each has its own advantages and disadvantages in comparison with the others. Therefore, further work in this area is required. * Another difficulty in using mathematical models to simulate mass transport in porous media is the determination of model parameters. Because it is impossible to measure directly the pore velocity and dispersivities in the whole flow region, one has to use an indirect method to determine these input parameters according to observation data of heads and concentrations, i.e., to solve the inverse problem of the simulation of groundwater quality. The inverse problem of groundwater flow simulation has been extensively studied (15)(16)(17)(18)(19)(20)(21)(22)(23)(24). A review paper on this topic was presented recently by Yeh (25). In comparison, the body of published research for parameter estimation for mass transport problems is still small (26)(27)(28)(29)(30)(31)(32)(33)(34). In the solution of inverse problems, the primary obstacle is their ill-posed nature, i.e., the existence, uniqueness, and stability of the solution may not be satisfied. However, if the observation data are collected through a reasonable experiment design, the ill-posedness may be controlled under certain restricted conditions. In this paper, we focus our attention on the numerical methods with respect to the previously mentioned three problems. The distributed parameter models of mass transport in porous media are given in the second section. This paper does not deal with stochastic models (35)(36)(37)(38).
The third section is devoted to 3-D numerical models developed in recent years and mentions the mixed finite difference-finite element method (FD-FE), the multiple cell balance method (MCBM), and the alternating direction Galerkin method (ADG). The MCBM (13,(39)(40)(41) combines the idea of local mass balance coming from the integral FDM (42), and the idea of basic functions coming from FEM (43,44). It is very suitable for the solution of field problems.
In the fourth section, various numerical methods are reviewed. These have been developed to deal with the advection-dominated transport problems. We mention the upstream weighting FEM, the upstream weighting MCBM, and the method of characteristics (MOC) and its modified forms. A new method, the advection-control MCBM, recently presented by Sun and Liang (45), will be also introduced. It displays the ability of eliminating the overshoot and decreasing the numerical dispersion efficiently.
In the last section of the paper, the concepts and numerical methods of inverse problems are introduced. Here we discuss extended identifiability and its applications to experimental design and reliability estimation.

Mathematical Models
It is well known that mass transport in porous media is governed by the following advection-dispersion equation (1,2): The distribution of average velocities can be computed according to Darcy's law. For instance, in the saturated zone we have Vji = -!i4 0 axj (3) where Ku = components of hydraulic conductivity tensor (L/T); h = head (L); 0 = porosity (dimensionless). The distribution of head in Eq. (3) can be obtained by solving the groundwater flow problem. The advectiondispersion equation [Eq. (1)] must be completed by boundary and initial conditions depending on the considered problem. The initial condition is the concentration distribution at t = 0: where CO (x) is a known function, point x belongs to the flow region. The boundary conditions are generally divided into three types: a) Prescribed concentration boundary conditions C(x,t) = g1 (x,t), x E rI i,j = 1,2,3 where C = concentration of contaminants (M/L); Du = components of hydrodynamic dispersion coefficient tensor (L2/T); Vi = components of average velocity vector of fluid in porous medium (L/T); M = source or sink (M/ L2T); 0 = volumetric fraction of the fluid phase, in saturated zone; 0 = porosity; in unsaturated zone, 0 = water content; xi = space variables (L); t = time (T).
The repeated subscripts indicate summation.
For an isotropic porous medium, Dii may be expressed as: where aL, aT = longitudinal and transverse dispersivities of the isotropic medium, respectively (L); V = magnitude of the velocity (L/T); D* = molecular diffusivity of porous media (L/T); and bu = Kronecker delta, b& = O, when i + j; u = 1, when i = j.

Classical Models
Analytical solutions can be derived only for a few classical models that are often used to simulate controlled experiments in the laboratory or in the field to obtain dispersion parameters. Another application of analytical solutions is testing and contrasting numerical solutions. For this reason, two solutions are quoted here.
Consider the movement of a tracer in a semi-infinite column. The mathematical statement of the problem is ac 2C ac The analytical solution of the classical model was given by Fried and Combernous (46): Another example is the dispersion in radial flow when a fully penetrating recharge well is located at the origin with a constant rate Q (the molecular diffusion is neglected). The mathematical statement of the problem is ac a2c ac a t ar2 a r C(r,O) = 0, where r, = radius of the injection well; the radial velocity V = Q/2-rBOr, where B = thickness of the aquifer; and CO = concentration of the recharge water. The analytical solution of problem (10) is given by Tang and Babu (47). A new formula that is simpler than the former was presented recently by Hsieh (48). More analytical solutions can be found in Li (49), Bear (2) and Sun (8). Some new results are given by Van Genuchten (50), Mironenko and Pachepsky (51), Lindstrom and Boersma (52), and Mornch (53) (58), Chen (59,60), Lowell (61), and Dillon (62) for transport in fractured porous media.

Field Models
For a field model, the groundwater flow equation (GFE) and the advection-dispersion equation (ADE) must be solved individually or simultaneously. Several combinations of GFE and ADE are given below. These are often seen in practical problems.  (14) where Sy = specific yield (dimensionless); m = saturated thickness of the aquifer (L); W2 (L/T) and M2 (M/ L2 T) = sink or source; i,j = 1,2.
The GEF and ADE of a field model are coupled through Darcy's law. In tracer cases (the density p and viscosity ,u of the fluid are constant); the two equations can be solved individually. Otherwise (p and ,u depend upon the solute concentration) they must be solved simultaneously with the state equations (1,79).
The effects of ion exchange, chemical reactions, radioactive decay, and adsorption, as well as pumping and recharging of wastewater, are all included in the source or sink term of the ADE (2). For equilibrium transport, the linear adsorption isotherm F = kC is often used, where F is the solute concentration on the solid phase and k is the slope of the isotherm. The ADE [Eq. (1) The two above equations represent the mass transport and exchange in fluid phase and solid phase, respectively, where k and a are constants. [For review of adsorption models, see (7,(80)(81)(82)(83)(84)(85)(86).]

Numerical Methods
The previously mentioned field models are generally intractable analytically because of irregular geometry and a nonhomogeneous structure of the flow region. Numerical methods are always required. In the past, the practical application of numerical solutions was almost limited in 1-D and 2-D cases to avoid high expenses of 3-D models on both human and machine resources. However, dispersion of contaminants in groundwater is an inherently 3-D process. It is now recognized that even a low-order transverse dispersion acting over a long distance can substantially affect the configuration of a contaminant plume (14,87,88). Burnett and Frind (89) further examined the effects of dimension through a field-scale system that is simulated with a 3-D model, a 2-D vertical section model, and a 2-D vertically integrated, horizontal plane model. In fact, there are many instances of field problems where 2-D approximations of contaminant plumes are not adequate. Therefore, numerous researchers have turned their attention to 3-D numerical solutions in recent years (40,41,63,(90)(91)(92). The major problem in investigating 3-D numerical models is knowing how to decrease the required input data, storage space, and computation time.
second derivatives in Eq. (20) by Green's theorem, one can translate Eq. (20) into a system of ordinary differential equations.
In the Galerkin FEM, let {¢j} = {wi}. In eliminating the In this case, the input geometry information is almost

Mixed Finite Difference-Finite Element Method (FD-FE)
To decrease the high expenses of 3-D models both on human and machine resources, Sun (17) presented a mixed FD-FE method for 3-D GFE that can be easily extended to solving the 3-D ADE. A similar method was given by Babu and Pinder (90,94). The flow region is divided into triangular prism elements as mentioned above. All 3-D nodes can be numbered by a 2-D array (i,m) where m is the number of a layer that the node locates (m = 1,2, .. ,M) and i is the number of the node in the layer (i = 1,2, . .. J,). Then a double division is introduced to form an exclusive subdomain for each node that is a multiangular prism surrounding the node as shown in Figure 2. The local mass balance over the exclusive subdomain (di,m) of node (i,m) can be expressed as the integral of ADE [Eq. (1)]. Applying Green's theorem, it can be rewritten as boring nodes of node i on the same layer m. Aj,j and A can be computed directly and will be given later in Eq. (42) (32) [LZZaz- In Babu and Pinder (90), assume that where Vi(x,y) are 2-D basis functions (36) Using the Galerkin process in the x-y space and the finite difference approximation in the z space, they ob-

Multiple Cell Balance Method (MCBM)
The 2-D MCBM (13) is a variant of the linear FEM. It combines finite-element approximation with the concept of the local mass balance used in the integral finite difference method (IFDM), presented by Narasimhan (12), Narasimhan and Witherspoon (42), and Rasmuson et al. (66). In IFDM, the flow region is divided into polygon elements. The local mass balance of each element is represented by an integral form over the element and the derivatives arising in these integrals are replaced by simple finite differencing. Consequently, the balance equation of each element becomes an algebraic equation that includes the unknown concentrations of the element and its neighboring elements. IFDM is easy to understand, but it is difficult to deal with the tensor nature of dispersion coefficients and arbitrary orientation of flow velocities.
In MCBM, the flow region is divided into triangle elements, a twice division is then used to form the exclusive subdomain for each node as shown in Figure 3. The local mass balance of a node i can be represented by the integral of Eq. (14) over its exclusive subdomain (di): Assume that in each element (Aijk) the unknown function C can be represented as C(x,y,z,t) = Ci(t) *g(x,y) + Cj(t) *j(x,y) The MCBM has been extended to 3-D models (40,41). The division of the 3-D flow region is the same as that mentioned before. Assume that in each triangular prism element ( Fig. 1) we have 6 C (x,y,z,t) = , Cj(t)4I(x,y,z) (45) 1=1 are basis functions given in Eq. (27). Substituting Eq. (45) into the 3-D mass balance equation Eq. (29) all integrals in the equation can be calculated explicitly. Consequently, we obtain a set of equations in as the 2-D case. A layer-by-layer iteration process can be used to solve the algebraic system.   Table 1. The MCBM has been used to simulate field problems in China. For instance, a simulation model of groundwater quality at Xian city, Shaanxi province, was created in 1985. The MCBM almost keeps the advantages of the mixed FD-FE method on decreasing input data and computational effort. In addition, it is easy to modify in order to deal with field problems that have sharp fronts. A man-machine talking software package based is used in z space, where repeated subscripts indicate summation as before. This method takes full advantage of the nature of flow and mass transport in multilayer systems of several aquifers and aquitards.  (97) and Frind and Germain (98). In the 3-D ADG, the advection-dispersion equation is split into three symmetric steps corresponding to directions x, y, and z, respectively; the final matrix equations can be solved in each direction by means of the Thomas algorithm for tri-diagonal matrices. However, the division of the flow region in the ADG method must be limited in a special way.

Advection-Dominated Transport Problems
The numerical methods mentioned previously are all successful for dispersion-dominated transport problems. For instance, when the Peclet number Pe = 0.5 in the classical model [Eq. (8)], the Galerkin FEM, the lumped mass FEM, and the MCBM all gave numerical solutions that nearly coincide with the analytical solution of the model (Fig. 4). However, for advection-dominated problems, they suffer from two kinds of numerical difficulties: overshoot and numerical dispersion. The former is characterized by oscillations around the computed concentration profile, the latter, by the smearing of the concentration front (Fig. 5, where Pe = 100).
Gray and Pinder (99) assumed that the spatial distribution of concentration can be expressed as a Fourier series with components of different wavelengths. The behavior of a numerical scheme is then evaluated in terms of its ability to propagate each component of the series with the correct speed and amplitude. They discussed several FD and FE schemes for the classical model [Eq. (8)]. When advection dominates dispersion, the solutions generated from both FEM and FDM are characterized by overshoot, although the oscillations are less pronounced for the finite element solution. Gray and Pinder (99) also pointed out that overshoot and numerical dispersion are closely related in FD and FE schemes. Generally, one can select a numerical scheme that decreases the magnitude of one problem at the expense of increasing the effect of the other. Various improved numerical methods have been presented in recent years to overcome these numerical difficulties, such as upstream weighting methods (UWM) (100)(101)(102)(103), moving point methods (MPM) (104)(105)(106)(107)(108), and other alternative methods (109)(110)(111). Perhaps, the most simple and efficient approach is the mesh refinement technique that improves numerical solutions by decreasing the local Peclet number. Some adaptive forms of this technique were given by Shepherd and Gallagher (111) and Yeh (112). However, when Peclet number tends to infinity, the mesh refinement technique may not be economically feasible. A perfect numerical method for advection-dominated transport should eliminate overshoot and simultaneously improve the accuracy of numerical solutions with only a little increase of computational effort. Beside that, it should deal with complex 3-D field problems.

Upstream Weighting Method (UWM)
If one uses the FDM to solve 1-D advection-dispersion equations, the first and second derivatives, with respect to the space variable x at a grid i, are usually replaced by central difference approximations as follows: 10   Ci+1 -2Ci + Ci-1 (Ax)2 (47) are used as basis functions in the isoparametric element, -1 -e -1, then the weighting functions are given by wi (t) = +j (4) -aF(k) W2(V) = +2(t -aF(t) (53) where a is a parameter to be determined and then the oscillation of the numerical solution is eliminated at the expense of increasing the numerical dispersion significantly as shown by Gray and Pinder (99).
In order to obtain an optimal compromise solution, the weighting average of Eq. (47) and Eq. (49) is used: where 0.5 < a. < 1 is the upstream weighting factor. Substituting Eqs. (48) and (50) into Eq. (8), we obtain a numerical scheme that is called the first-order UWFDM. However, it is possible to present high-order upstream FDM, for example, Leonard (113) presented a third-order scheme. Its accuracy is higher then that of first order UWFDM, but at the expense of increased computational effort, because the coefficient matrix of finite difference equations in this scheme is no longer tri-diagonal. In 2-D models, an upstream weighting factor a is introduced into FD equations by changing Ci = 1/2 (C1 + Cj) to Cij = aCj + (1-a)Cj (51) where CU is the concentration of water flowing between nodes i and j; Ci and Cj are the concentrations at the nodes. If fluid is moving from node i to node j, then 0.5 < a < 1. Generally, let a = 0.75. Eq. (51) is also used in IFDM [Rasmuson et al. (66)] to damp the oscillations of numerical solutions with a = 0.65. Note that the tensor nature of the dispersion coefficient is not considered in the UWFDM. Therefore, it is difficult to deal with field problems. In the late 1970s, the idea of upstream weighting was extended to FEM. The terms "upwind FEM," "asymmetric weighting function FEM" and "Petrov-Galerkin FEM" are often It is interesting to note that for 1-D ADE with constant coefficients, the discretization equations of UWFEM using Eqs. (52) and (53) are exactly the same as that of the first-order UWFDM. The parameter a in Eq. (54) corresponds to the upstream weighting factor a in Eq. (50). Because the weighting functions [(Eq. (53)] are quadratic functions, the UWFEM needs more computational effort then that of the linear FEM. Some authors suggested that the upstream weighting functions should only be used for the advection term of Eq. (20). For the other terms of the equation, the basis functions are still used as weighting functions (114).
Extensions of UWFEM with linear basis and quadratic weighting functions to 2-D or 3-D cases are straightforward (115). The UWFEM using triangular elements was presented by Huyakorn (116). The weighting functions are wl(x,y) = Oi(x,y) + F(x,y) along the sides of elements that are not the flow di tions. Consequently, the effect of so-called "crossv dispersion" may decrease the accuracy of numerica lutions. In order to lessen the crosswind dispersio modified form named the streamline upwind/Pet Galerkin (SU/PG) method was presented (102,117,1 Generally, the accuracy of SU/PG method is higher t the PG method. If the set of weighting functions is required t( orthogonal to the set of basis functions, the corresp ing method is named as the orthogonal upstr weighting (OUW) finite element scheme. Yeh ( noted that the OUW scheme provides an alternativ4 large-scale problems because it gives a convergent p successive overrelaxation (SOR) computation for all clet numbers.
Although UWFEMs have been extended to 2-D even 3-D cases, it is necessary to remark that t practical applications are limited because of comp tional complexity.
Sun and Yeh (39) presented a new upstream wei ing method for 2-D models based on the MCBM n tioned above. In this method, the center of each angular element is divided into three subtriang elements by linking the center and vertices of the ment. In Figure 7, m is an invented node. The con tration at the invented node is defined as the weigl average of the three vertices of the element: numerical examples including a field problem of elecontaminant transport in a stream-aquifer system. Beicen-cause upstream weighting coefficients in the method are hted introduced into basis functions directly, the representation of unknown functions in each element remains linear. The coefficient matrix of discretization equations can be computed without numerical quadrature. Addi- (57) tionally, the method is specially available for the solution of field problems. In fact, the method has been used to predict groundwater pollution at several aquifers in ae China. are The 3-D form of UWMCBM was given by Sun et al.
Odes (40). Several numerical examples were given by Wang local et al. (41). A multilayer transport problem explains that the the elimination of oscillations of numerical solutions is For also important for field problems. However, all UWMs previously mentioned have two disadvantages: overshoot is controlled, but at the expense of increasing numerical dispersion, and the upstream weighting coefficients arising in each UWM have to be designated artificially.

Moving Point Methods (MPM)
All UWMs can be seen as Eulerian methods. The advantage of Eulerian methods is the use of a fixed grid, but they are generally not well suited for the handling of sharp concentration fronts. In contrast, Lagrangian methods use either a deforming grid or a fixed grid in deforming coordinates. For example, O'Neill (123) proposed a moving, deforming coordinate system that deals with steep concentration fronts with no additional nodes \ and time steps. However, for lack of a fixed grid, it is difficult at handling fixed sources, fixed boundary conditions, and nonhomogeneous structures of media. In addition, it requires more computational effort in gen-erating coordinate systems, especially for 3-D field problems. Mixed Eulerian-Lagrangian methods attempt to eliminate such difficulties by combining the simplicity of a fixed Eulerian grid with the computational power of the Lagrangian approach (106). A complete review of these methods was given by Farmer (124).
The method of characteristics combining with FDM (MOC/FDM) was presented and successfully applied to field problems by several researchers (125)(126)(127). Considering the 2-D ADE [Eq. (14)] that can be combined with the GFE [Eq. (13)] to give at = 1 a mDijac -Viaac+F (61)  The unknown concentration is then divided into two parts: the advection concentration and the dispersion concentration. The former is determined by Eq. (63), which can be solved numerically by tracking the movement of a set of moving points or particles (Lagrangian). The latter is determined by Eq. (64), which can be solved by FDM. Note that there is no advection term in Eq. (64) and the problem of numerical dispersion can be avoided. In MOC, an interpolation process is always required to generate mesh concentrations from the concentrations of moving points and vice versa. Generally, the interpolation process decreases the accuracy of numerical solutions significantly. Another disadvantage of MOC is that the sources (or sink) and boundary conditions must be specially treated to generate or eliminate moving points that make the MOC cumbersome to program and execute. The reports on field applications of the MOC are still limited in 2-D problems.
Several modified forms of the MOC have been presented. One of them is the single-step reverse particle tracking technique that is used to replace the forward moving particle tracking procedure of the MOC (105,128,129). A single-step reverse point of a node is the point where a particle leaving it at the time tk will reach the node location exactly at the time tk+, . In the finite element version of the modified MOC, the concentration of a reverse point at tk is determined by the finite element interpolation, and it is used as the advection concentration ofthe node at tk+ 1. The dispersion concentration of the node can be obtained by the solution of Eq. (64) with FEM. For second-order and third-order versions of the method, see Farmer (124).
Neuman (106) presented an adaptive Eulerian-Lagrangian FEM where the forward moving particle tracking technique is used around each sharp front to improve the accuracy ofnumerical solutions. Away from such fronts the single-step reverse particle tracking technique is used to save the computational effort. The number of moving particles is adaptive. When the front sharpens, moving particles are inserted; when the front flattens, they are eliminated. The adaptive scheme maximizes computational efficiency.
A similar method that is named the hybrid movingpoint method (HMPM) was presented by Farmer and Norman (107). Farmer (124) noted that the finite element version of the modified MOC is a good method for smooth problems (fronts spread over three or more mesh lengths), and the HMPM a good method for problems with sharp fronts.
A variation of MPM is the random-walk method (RWM). It is based on the concept that dispersion in porous media is a random process. In this method, advection transport is simulated by a large number of moving particles like the MRM but the particle movement takes place in continuous space. Dispersion transport is simulated by adding a random-walk process to each particle in each time-step, instead of the solution of the dispersion equation. Consequently, it is not necessary to calculate the mesh concentration and update the particle concentration in every time-step. The concentration distribution needs to be calculated only when it is of interest. The main drawback of the RWM is that the computed concentrations are strongly dependent on sample size. If the number of moving particles is not adequate, a coarse solution may be observed. However, the RWM is conceptually simple and relatively easy to program. It has been used to solve 2-D field problems with many thousands of moving particles. A complete report with a computer program was given by Prickett et al. (108).

Other Alternative Methods
There are some alternative methods for the solution of advection-dominated transport problems. The finite analytic method (FAM) was presented by Chen and Li (130) and Hwang et al. (109). The basic idea of the FAM is to invoke the local analytic solution of the governing equation in the numerical solution of the problem. In this method, the flow region is divided into small rectangular elements. For each element, the boundary conditions prescribed on the four sides are approximated by either a second-degree polynomial or by combining a linear and an exponential function. The local analytic solution is then expressed in an algebraic form relating an interior nodal value to its neighboring nodal values. The system of algebraic equations derived from local solutions is then solved to provide the solution of the total problem. Because the discretization equation as-sociated with each node is derived from the analytic solution technique, the numerical difficulties caused by the truncation error can be avoided. Hwang et al. (109) noted that the coefficients in interior nodal expression have upstream shift that make the behavior of numerical solutions of the FAM be like that of UWMs. For the case of a high Peclet number, the FAM also cannot produce sharp fronts as predicted by analytical solutions. An advantage of the FAM is that it does not require artificially designated upstream parameters as in UWMs. However, the FAM has not met much approval as yet.
When the boundary element method (BEM) is used to solve the contaminant transport problem, the governing equation, e.g., 2-D ADE [Eq. (14)] must first be translated into the following form: at ax a The translation can be done by assuming that the coordinate system is defined in the principal axes of the dispersion tensor. Molecular diffusion can be ignored, and dimensionless variables are used in Eqs. (65) and (66). Using Green's second identity, we have XC(rit) = { aG-Gac dF + GFd (67) in which: (F) is the boundary of the flow region (Qi); function G = In (r-ri) where r-ri is the distance from any field point r = (x,y) in (fQ) to a base point ri = (x,,yi); X = 2ni, -Fr, 0, if (xi,yi) in (fQ), on a smooth part of (F), or a kink of (F), respectively, where 0 is the internal angle at the kink. The integrals on the right-hand side of Eq. (67) can be calculated approximately by linear interpolation along (F) and the finite element interpolation in (fQ), respectively, then discretization equations are obtained. Unfortunately, the coefficient matrix of the system of linear equations is fully populated. For details of the method, see Taigbenu and Liggett (110). Numerical examples show that numerical solutions of the BEM display small oscillations and large numerical dispersion when the Peclet number is large. Recently, Sun and Liang (45) presented a new numerical method that is named as advection control MCBM for the solution of 2-D ADEs. The fundamental idea of the method differs from general UWMs and the MOC in such a way that a given concentration term depending upon advection is added to the right-hand side of multiple cell balance equations to control the character of numerical solutions. The flow region is divided into triangular elements. The vertices i,j,k and the center m of each triangular element Aijk are considered as real and invented nodes, respectively (Fig.   7). The concentration Cm of the invented node m is decomposed into two parts: the regular part Cm,r and the advection control part Cm,a. Let Cm = WCm,r + (1-W)Cm,a (68) where 0 < w < 1 is an advection contr9l coefficient that can be determined by an adaptive procedure. Cm,r (Ci + Cj + Ck)13, Ci, Cj, and Ck are the unknown concentrations of real nodes i,j, and k, respectively; Cm,.
is the pure advection concentration of the invented node m that can be obtained by a modified single-step reverse particle tracking technique. By constructing the MCB equations for all nodes and eliminating the unknown concentrations of invented nodes with Eq. (68) one obtains a set of discretization equations that include all unknown concentrations of real nodes on its left-hand side and all pure advection concentrations of invented nodes on its right-hand side. The character of the numerical solution ofthe equation is controlled by w. When w = 1, the method reduces to the general MCB method that is accurate enough for dispersion-dominated cases. When w = 0, the numerical solution is usually determined by pure advection. Several classical problems and a field problem are solved in order to verify and demonstrate the usefulness of the proposed method. The method is very simple both in concept and computation. Also, it is easy to extend to 3-D cases. A finite-volume method with high-resolution upwind schemes is introduced by Putti et al. (131) for the solution of the ADE. Due to an efficient interpolation technique that is employed, the results of this method show good agreement with analytical solutions for a full range of cell Peclet numbers.

Inverse Problems
At the initial stage of constructing a model for a field problem, often the average velocity and dispersion coefficient included in the ADE are unknown. These parameters are difficult to measure directly from the physical point of view and have to be determined from historical observations. The inverse problem of advection-dispersion simulation (ADIP) is determining unknown parameters imbedded in the simulation model by using a set of observations of concentration and head distributions in time and space. Sometimes, the unknown parameters also include the strength of sink (or source) and boundary conditions. Literature published in this area is limited because most concepts and methods used in the inverse problem of groundwater flow (GFIP) can also be used in the ADIP. Another reason is that the ADIP is more difficult to solve than the GFIP and new results are hard to find. In the ADIP, we have to deal with two problems: the nonlinear problem arising from the dispersion coefficient (depending on the velocity of groundwater flow) and the problem of differentiating the effects of advection transport and dispersion transport. These problems increase the uncertainty of the ADIP. Therefore, in-vestigating the ADIP should include two fields: extending methods used in GFIP to the ADIP, and developing new methods for the ADIP.
Note that the ADIP is generally coupled with the GFIP (32,(132)(133)(134)(135). The observed concentrations depend on the porous velocity and the latter depends on the unknown porosity and hydraulic conductivity that have to be identified by solving the GFIP.

Parameter Identification
In order to solve the ADIP numerically, the unknown parameter must be first replaced by a vettor with finite dimension N. This can be done by either using the zonation method or the finite element interpolation method (22,136). As the GFIP, the ADIP also can be translated into an optimization problem, depending upon what kind of performance criterion is used.
When the equation error criterion is used, the discretization equations derived from any numerical method are rewritten as (69) in which: [A]LXN = coefficient matrix and b = L-dimensional column vector with [A] and b dependent on the observed concentration distribution; p = N-dimensional unknown parameter vector; E = residual column vector with components ei, (i = 1,2, ... ,L); and L is the number of observations. Numerical methods based on minimizing the residual E are classified as direct L methods (15). If minE lEil is taken as the minimization i= 1 criterion, then the ADIP is translated into a linear programming problem. If min TE iS taken as the minimization criterion, then the unknown parameter vector can be estimated as (70) Direct methods require an interpolation process to obtain the observed concentration at all nodes; the observation noise and the interpolation error cannot be avoided. Unfortunately, direct methods are generally unstable in the presence of noise.

(p) = ([A]T[A])-1 [A]T{b)
Another approach is the use of output least square error criterion. The objective function of minimization is the weighted sum of squared differences between simulated and observed concentrations: min E = (C* -C(p))T [W (C* -C(p)) (71) in which: C* = observed concentration vector with noise; C(p) = nodal output, when p is used as the nodal parameter; [W] = diagonal matrix of weights, its element wii may be taken as 1/Ci2 (p) (29) (137) used the quadratic programming method, and Wagner and Gorelick (29) used a quasi-Newtonian algorithm. In fact, all methods mentioned in Yeh (25) for the GFIP can be extended to the ADIP directly. A variety of statistical methods can be done in the same way (16,18,138). The identification of the sink or source term is a little easier than that of dispersion parameters. For example, Gorelick (139) used the linear programming method to identify sources of groundwater pollution.
However, the difficulty of inverse problems cannot be moved by the application of numerical methods. It asserts itself either by the ill-condition of the coefficient matrix or by the existence of many local minimizers in the optimization problem. The identifiability of unknown parameters is defined initially as the one-to-one property of mapping from the space of system outputs to the space of parameters (140). However, the uniqueness of such mapping is extremely difficult to establish and often nonexistent. Another definition of identifiability was given by Chavent (141). The parameter is said to be output least square identifiable if and only if a unique solution of the optimization problem Eq. (71) exists and the solution depends continuously on observations. It is also very difficult to achieve in field problems.
Yeh and Sun (21) presented a concept of the extended identifiability. The parameter p is said to be 8-identifiable if the solutions of Eq. (71) must be equivalent when they are used in the prediction model (a small, prescribed error is allowed). For field problems, the 8identifiability generally can be achieved by only point measurements. Chavent (142) reviewed the concepts of identifiability of distributed parameters systems. He summarized and compared five definitions of identifiability.
Recently, Sun and Yeh (135) presented a new concept of identifiability for coupled inverse problems that is named management equivalence identifiability (MEI). An identified parameter is said to be management equivalence indentifiable, if it can satisfy the accuracy requirement of a given management model. Using MEI, the whole procedure of groundwater modeling, from the experimental design, parameter identification to management decision can be considered systemically.

Sensitivity Analysis
In order to solve the optimization problem [Eq. (71)], we have to calculate sensitivity coefficients aCi(p)/apj in Eq. (73) or gradient components aE(p)/hpj, j = 1,2, ..., N. There are three approaches that can be selected: the influence coefficient method (differencing method), the sensitivity equation method, and the variational method (25). If the number of components of the parameter to be identified is greater than the number of observations (N > L) the variational method would be advantageous. The adjoint sensitivity method presented by Carter (143) and Chavent (144) is easy to extend to advectiondispersion equations. If the unknown parameter p is the longitudinal dispersivity aL, we have OA= J fVq(x,y,t- ¶) DLVC(X,Y, ¶) dtdQ (74) where (fi) is the exclusive subdomain of node i; V is the gradient operator; and C(x,y,t) is the solution of the 2-D ADE with additional conditions. q' is the time derivative of q(x,y,t) which is the solution of the following adjoint problem: x[Dap -] + a (Vaq) + i= (75) with final and boundary conditions: In Sun and Yeh (134), the adjoint problem is derived for general coupled problems through introducing a set of adjoint operation rules. An example of parameter identification for coupled groundwater flow and mass transport is given.

Experiment Design
An experiment design X is a decision that includes the following components: a) locations of observation wells (space distribution of observations; b) observation times (time distribution of observations); c) locations of injection and pumping wells (locations of stimulation); and d) pumping or recharge rate and the concentration of injection water (strength of stimulation).
The method of sensitivity analysis can help engineers to design field tracer experiments. Mercer and Faust (149) presented an example where a two-well tracer test was designed to determine the well spacing, the injection flow rate, and the injection period.
The reliability of parameter estimate depends upon the quantity and quality of observations. If the estimated parameter p is located in a neighborhood near the true parameter p* in the parameter space, then we have the theoretical covariance matrix: (80) where a2 is the common variance of the random error in observations. Obviously, the maximization of the determination of ([J]T[J]) will both improve the condition number of Eq. (72) and decrease the uncertainty of the estimated parameter. Therefore, the objective max A = det ([J]T[J]) (81) is often used as a criterion of the optimal experiment design (30). The senstivity coefficients that are elements of [J] in Eq. (81) depend upon the unknown true parameter; Eq. (80) is correct only if p is a good approximation of p*. It is very difficult to make an optimal experiment design by using Eq. (81) for a field problem.
Yeh and Sun (21) presented a new concept named the sufficiency of observations. The observations of an experimental design are said to be sufficient if any two estimations Pi and P2 that are non-equivalent for the objective of model prediction are 5-identifiable when these observations are used in the identification procedure. In other words, these observations are adequate for determining a weak, unique solution of the inverse problem. The sufficiency of design X can be judged by solving a nonlinear programming problem. Obviously, this method is very useful for practically designing a field experiment.
Recently, this concept has been extended to deal with management problems (135). The question of whether the observations are sufficient for the MEI can be answered by solving a constrained optimization problem. Hsu and Yeh (150) presented a design criterion that minimizes experimental cost while meeting the reliability of model prediction. They formulated the optimal experimental design in terms of a nonlinear mixed-in-teger programming problem and solved it through a heuristic approach.

Reliability Estimation
After getting a parameter from an identification process, we should estimate the reliability of the identified parameter. The covariance estimation [Eq. (80)] may not be available because the true parameter is unknown. In general cases, Monte Carlo simulation can be used to evaluate the reliability without any additional assumption. In the Monte Carlo method, each trial is a different realization of observations used to solve the inverse problem. The output set of Monte Carlo runs is then used to estimate the statistical properties of the identified parameter. Several hundreds or thousands of realizations may be necessary to arrive at a stable estimate of statistical results. If the kth realization of parameter estimate is Pk, then the expected value and variance of the parameter are Once the expected value and the variance of the.unknown parameter are obtained, it is necessary to further estimate the uncertainty of the output (concentration) of the constructed model. For this reason, we need to estimate the uncertainty of the solution of an ADE from the uncertainty of its coefficients. Tang and Pinder (151,152), introduced a numerical scheme based on a perturbation expansion technique and its extension to solve the problem. From numerical examples, they concluded that the greatest uncertainty of the solution occurs at the concentration front, and the transport equation appears to dampen the uncertainty associated with the physical parameters: velocity and dispersivity.
In fact, the ADE governing the contaminant transport in groundwater should be seen as a stochastic partial differential equation because the uncertainty is included in its coefficients, initial and boundary conditions, and sink or source terms. If the model is used in a prediction or a decision-making process, quantifying the uncertainty of the model may play an important role. For a complex problem, it is convenient to use the Monte Carlo method to evaluate the uncertainty. Nelson et al. (153) presented an analysis sequence that begins with field data and results in the confidence limits of the model output. A field model is used to explain the technique. When the covariance matrix of parameter estimation error is provided as a part of the model calibration, then the velocity field and fluid flow paths with the associated travel times are generated for each realization. Also, the confidence limits of the model output are obtained by the Monte Carlo simulation.
Note that the reliability of a model prediction is already included in the 8-identifiability (21,135). Provided observations are sufficient for the 8-identifiability; the identified parameter must be equivalent with the true one in the sense of satisfying prediction and management requirements.
To obtain the uncertainty of using an estimated parameter, another approach is the use of stochastic models (20,(154)(155)(156), where the covariance of the estimated parameters can be obtained directly by the use of conditioning or unconditioning geostatistical methods.

Summary
In this paper, we focus our attention on numerical methods for solving 3-D problems, advection-dominated problems, and inverse problems. The solution of these problems is very important, not only in theoretical investigations, but also in practical applications.
In order to describe the shape of contaminant plumes correctly, we need 3-D and multilayer models. The main problem in applying 3-D models is decreasing the computational expenses both on human and machine resources. Thus the Galerkin FEM with special elements and basis functions, the mixed FD-FE method, and the MCBM are introduced. The MCBM keeps the local mass conservative. It is simple in concept and easy to use for field problems. A layer-by-layer iteration process is efficient for solving the discretization equations. A software package with the functions of figure input and figure output was also presented.
The UWM, modified MOC, BEM, FAM, and ACM were introduced and compared for solving advectiondominated transport problems. The advantages and disadvantages of each method were shown. The adaptive Eulerian-Lagrangian method seems to be a good approach for getting a high-quality solution. For field problems, a kind of modified MCBM or ACM is a simple, efficient method to obtain a satisfactory solution without oscillations.
Most methods used for ADIP are those that directly extend the one for GFIP. Several optimization techniques based on different criteria were introduced. In the solution of ADIP, the main difficulty was its illposedness. We do not know how many observations are adequate for getting a unique solution of the inverse problem. However, a weak identifiability, i.e., the 8identifiability, can be determined with limited observations. The calculation of sensitivity coefficients was necessary, both in the numerical methods of the ADIP and in the experimental design. The adjoint sensitivity analysis was discussed for the ADE and general coupled probems. In order to construct a practical model, a systematic procedure from the data collection to the uncertainty estimates of the built model is required, in which the Monte Carlo method plays an important role.