- Email: [email protected]
A new numerical method for groundwater fiow and solute transport using velocity field*
356
2008,20(3):356-364
A NEW NUMERICAL METHOD FOR GROUNDWATER FIOW AND SOLUTE TRANSPORT USING VELOCITY FIELD* ZHANG Qian-fei School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200030, China Shanghai Municipal Engineering Design Institute, Shanghai 200092, China,E-mail: [email protected] LAN Shou-qi Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China WANG Yan-ming Shanghai Municipal Engineering Design Institute, Shanghai 200092, China XU Yong-fu School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200030, China
(Received July 18, 2007, Revised March 9, 2008)
Abstract:A new numerical method for groundwater flow analysis was introduced to estimate simultaneously velocity vectors and water pressure head. The method could be employed to handle the vertical flow under variably saturated conditions and for horizontal flow as well. The method allows for better estimation of velocities at the element nodes which can be used as direct input to transport models. The advection-dispersion process was treated by the Eulerian-Lagrangian approach with particle tracking technique using the velocities at FEM nodes. The method was verified with the classical one dimensional model and applied to simulate contaminant transport process through a slurry wall as a barrier to prevent leachate pollution from a sanitary landfill. Key words:groundwater flow, solute transport, velocity vectors, finite element method, Eulerian-Lagrangian approach, particle tracking technique
1. Introduction Groundwater flow and transport analysis have been an important research topic in the last three decades, due to the fact that many geo-environmental engineering problems have a direct or indirect impact by groundwater flow and solute transport. The transport of dissolved or suspended materials by flowing water is of broad practical significance in relation to environmental protection and resource utilization [1-3]. The transport and fate of chemical or biological species is a major concern in relation to groundwater contamination [4,5]. The process of * Project supported by the National Natural Science Foundation of China (Grant No. 40201024). Biography: ZHANG Qian-fei (1975- ), Male, Ph. D., Senior Engineer
displacement of salts and nutrients in soils has an important bearing on agricultural production [6,7]. Environmental engineering problems in the recent two decades forced the research interest to more sophisticated numerical modeling, in order to predict the contaminant migration in the geological formation more accurately. Many techniques for solving numerically the advection-dispersion equation have been reported, using the most popular techniques, such as the Finite Element Method (FEM) [8, 9], Finite Difference Method (FDM) [10,11] , Boundary Element Method (BEM) [12], Fuzzy Sets Approach (FSA) [13], Artificial Neural Networks method (ANN) [14] , etc. Conventional analysis of groundwater flow was generally made by assuming that the total or pressure head is the only unknown variable. The velocity vectors were calculated after obtaining the pressure head distribution by applying Darcy’s law. However,
357
velocity vectors have to be calculated directly with the numerical method due to the following reasons: (1) to be used as direct input for transport model analysis, (2) to be used for inverse analysis when the measured velocity data from the field is available, and (3) to estimate nodal velocities more accurately especially in case of heterogeneous media. The advection-dispersion equation has been widely used to describe the Fickian transport of pollutants in the atmosphere, oceans, lakes, rivers and subsurface bodies of water. In recent years, there has been much interest in the possibility of using this equation to predict the effect of groundwater flow on the migration of radionuclides from geological repositories of nuclear waste to the biosphere [15]. This latter problem poses a special challenge to the numerical analyst because subsurface transport is often controlled by complex three-dimensional flow patterns which may vary with time [16], in an anisotropic dispersion process whose tonsorial description depends on the velocity decay, sorption, and chemical reactions. presented the Neuman and Sorek[17] Eulerian-Lagrangian method for solving the advection-dispersion equation, where the advection was solved independently at each time step, which they named it conventional method of continuous particle tracking. The residual dispersion was treated by the Eulerian finite element method. The aim of the Eulerian-Lagrangian methods is to combine the simplicity of the fixed Eulerian grid with the computational power of the Lagrangian approach. On the other hand, the Eulerian-Lagrangian method proposed by many researchers is based on interpolating the flow velocities at the element nodes which are estimated from the water flow model [18]. The technique based on the continuous forward particle tracking method which has to obtain the values of concentration, changed by the advection during a time interval, and projected onto the fixed finite element grid [19]. This projection process may yield numerical errors in the solutions. Therefore, the numerical method for estimating the velocity field directly at the element nodes should be used as the direct input to the transport model in order to make a better estimation of concentration caused by advection-dominated transport especially in highly heterogeneous porous media. For conducting transport analysis more accurately and inverting the hydraulic parameters with the water head and velocity vectors, a groundwater flow model is developed herein, which directly estimate the water pressure head as well as the velocity vector distribution in saturated and unsaturated porous media.
2. Groundwater flow analysis considering both pressure head and velocity vectors 2.1 Governing equations The continuity equation can be written in the form
Vi ,i = S s
wh wt
(1)
where Vi (i = 1, 2,3) is Darcy’s velocity components along the Cartesian coordinate axes, h the total head, S s the storage coefficient, and t the time. On the other hand, according to Darcy’s law, the equation of motion can be written as
Vi = K ij (T ) h, j
(2)
where K ij (i, j = 1, 2,3) is the tensor of the variably saturated permeability coefficients and it is a function of volumetric moisture content ( T ). Equations (1) and (2) are the governing equations for flows in saturated-unsaturated porous media with the pressure head and velocity vectors as unknown variables. Solving these two equations simultaneously require the boundary and initial conditions to be specified, given as Initial conditions:
h( x , t )
= h( x, 0)
t =0
Vi ( x, t )
t =0
for total head
= Vi ( x, 0) for velocity vectors
(3a) (3b)
Boundary conditions:
h( x , t )
*1
= hˆ( x , t )
for total head
Vi ( x , t )
*2
= Vˆi ( x , t )
for velocity vectors
K ij (T ) h, j ni
*3
= Qˆ ( x , t )
for flow rate
(3c)
(3d) (3e)
where x is the position vector, ni the unit normal vector along the Cartesian coordinate axes, hˆ the head, Vˆi the known velocity components, and Qˆ the prescribed flow rate. known
total
358
2.2 FEM technique The Galerkin-type finite element technique was implemented to formulate the finite element discretization. Linear hexahedral isoparametric element was used to represent the behavior of the total head h , and quadratic hexahedral isoparametric element was used to express the velocity vectors Vi . According to the Galerkin procedure, to approximately solve the governing equations, let the solution domain be subdivided into a number of finite elements and let a trial solution (h,Vi ) be expressed in the form:
(4b)
wc = <( DC vc) + q wt
20
J =1
where N J is the linear shape function for total head, and N J is the quadratic shape function for the velocity vector components. The formulation with the Galerkin-type FEM produces the following matrix form:
0
ª A « « 0 « 0 « x ¬ 'tC
0 0
A 0 'tC
A y
'tC z
B x º Vx ½ »° ° B y » ° Vy ° ® ¾= B z » ° Vz ° » D ¼ °¯ht 't °¿
0 ½ ° ° 0 ° ° ¾ ® 0 ° ° °¯ht D + 'tF °¿ where
AIJ = ³ ( N I N J )dV , V
BIJi = ³ ( N I Kij (T ) N J , j )dV , V
where V is the elemental volume, S the surface area of the element, ni the normal vector to the surface, and the others are the same as defined above.
(4a)
J =1
Vi ( x , t ) = ¦ N J ( x )VJi (t )
S
3. Advection-dispersion transport analysis 3.1 Governing equations In the theory presented here, we will consider only the main processes of manly controling the solute transport within the geological formation, namely, the advection and dispersion processes. Consider the advection-dispersion equation
8
h( x , t ) = ¦ N J ( x ) hJ (t )
FIJ = ³ ( N I N JVi ni )dV
(5)
where c is the concentration, t the time, the gradient operator, D the diffusion tensor, v the seepage velocity vector, and q the source term. Here, D is a symmetric positive-semidefinite tensor. The equation is subject to the initial and boundary conditions
c( x, 0) = C0 ( x )
(7)
( DC vc)
(8)
where * is the boundary, n the unit vector normal to * , pointing outwards, C0 , C and Q are the initial concentration, varying concentration and mass flux respectively, and D a controling parameter, which determines the type of prevailing boundary condition , that is, if D o f , Eq.(8) is the boundary condition for the concentration, if D = 0 , it is the boundary condition for the mass flux, in other cases, it is a mixed condition. Equation (8) applies to inflow and no-flow boundaries. Along outflow boundaries, it is common to assume that
DC
D w = + v Dt wt
DIJ = ³ ( N I Ss N J )dV ,
we can rewrite Eq.(6) in the Lagrangian form as
V
(9)
Using the material derivative operator
CIJi = ³ ( N I ,i N J )dV , V
(6)
(10)
359
Dc = <( Dc ) fc + q Dt
(11)
where f = v . Here c no longer represents the concentration at a point in space-time, but rather the concentration of a fluid particle moving at the velocity v . The path of this particle is described by the material derivative of x , which leads to the characteristic equation
Dx =v Dt
(13)
(14)
(15)
the Cauchy condition
vc
(16)
along inflow boundaries, and
vc
(17)
along the no-flow boundaries. The conditions along the outflow boundaries have no effect on c and are thus irrelevant to it. Clearly, the value of c for a given fluid particle remains constant as the latter is advected through the flow field. Thus, the “advection problem” defined by Eqs.(14)-(17) can be solved independently of cˆ . The residual concentration, cˆ , must satisfy
Dcˆ = <( Dcˆ) fcˆ + q + g Dt where
cˆ( x , 0) = 0
(20)
and the conditions along the inflow boundaries read (21)
where h = Dc
subjected to the initial condition
c ( x, 0) = C0 ( x )
subjected to the homogeneouse initial condition
(12)
One way to perform such a decomposition is to let c satisfy the homogeneous first-order partial differential equation
Dc wc = + vc = 0 Dt wt
(19)
( Dcˆ + vcˆ)
Equations (11)-(12) subjected to Eqs.(7)-(9) can be replaced by two sets of equations, one in terms of c , and the other in terms of cˆ , where
c = c + cˆ
g = <( Dc ) fc
(18)
Dcˆ + D (cˆ C ) * = u
(22)
where u = Q D c + h , and along the outflow boundaries
Dcˆ
(23)
Note that g , h , u and c play the role of prescribed source functions as c is known. We will refer to Eqs. (18)-(23) as the residual “dispersion problem”. It can be seen that the hyperbolic-parabolic advection-dispersion problem defined by Eqs.(6)-(9) can be formally decoupled into a purely hyperbolic “advection problem” defined in terms of c , and another predominantly parabolic residual “dispersion problem” defined in terms of cˆ . The approach is to first solve the advection problem for c , as the latter is dependent of cˆ and then solve the residual dispersion problem for cˆ . Another alternative is to solve Eq.(11) in its Lagrangian form without decomposing c into c and cˆ , and this approach is suitable for the regions with steep concentration gradients. The numerical scheme described below utilizes both of these alternatives jointly in an adaptive manner. 3.2 Eulerian-Lagrangian FEM approach Suppose that c is known at time tk , and we wish to compute it at tk +1 = tk + 't . First we set
c ( x , t k ) { c( x , tk ) ,
cˆ( x, tk ) { 0
(24)
We then solve the advection problem for c ( x , tk +1 ) by using the continuous forward particle tracking in the vicinity of steep fronts, and single-step reverse particle tracking away from such fronts. The residual dispersion problem is solved for c or cˆ by the finite element method.
360
The finite element equations in terms of c will be given. Let c ( x , t ) is approximated by a finite
the node, n , of the finite element grid. The projected values of c are denoted by k Cn such that
element function, c N ( x, t ) , defined as Mn
N
N
c( x , t ) | c ( x , t ) = ¦ cn (t )[ n ( x )
(25)
k
n =1
where N is the number of grid points or nodes, cn the concentration at node n , and [ n a basis function satisfying
[ n ( xm ) = G nm
(26)
where xm is x at node m , and G nm is the Kronceker delta. 3.2.1 Continuous forward particle tracking If there is a steep concentration front, the approach is to introduce moving particles at strategic points around the front and continuously track their positions along the paths. For any existing particle, p , located at point x p at time tk , the value of c is assigned to be c pk = c ( x p , tk ) for the duration of the time interval (tk , tk +1 ) . By virtue of Eqs.(24)-(25), we have N
c pk = c( x p , tk ) | ¦ cmk [ m ( x p )
(27)
m =1
For any new particle, r , at ( xr , tk ) the value of
c is assigned to be crk = c ( xr , tk ) . For a new particle along an inflow boundary, the value of crk is assigned as
crk =
DC + Q v
(28) xr ,tk
according to Eq.(16) for the duration of the time interval. At the end of the time step, each particle, p , reaches a new position
x
k p
=x +
³ vDt
k +1 p
, tk +1 )
p =1
Mn
(30)
Clearly, k Cn can be viewed as c at tk of a fictitious moving particle which reaches the node n at time tk +1 . 3.2.2 Single-step reverse particle tracking In the areas where the concentration gradients are mild, k Cn is computed by sending a fictitious particle from each node, n , backward to the point tk +1 k
xn { xn
³ vDt
(31)
tk
in each time step, which means that a particle leaving from k xn at tk will reach the grid point location, xn , exactly at tk +1 . If the velocity field does not change with time and 't is fixed, then k xn remains constant for each n and needs to be computed only once. By the virtue of Eqs.(24)-(25), k Cn can be expressed as k
N
Cn { c ( k xn , tk ) | ¦ cmk [ m ( k xn )
(32)
m =1
This requires much less computer time and storage than the procedure with the continuous forward particle tracking. However, the method cannot be used near sharp fronts because there exhibits numerical dispersion near such fronts. 3.2.3 Dispersion computation by finite element method Applying the Galerkin orthogonalization procedure to Eq.(11) leads to
ª Dc º <( Dc N ) + fc N q » [ n dR = 0, « R Dt ¬ ¼
³
tk +1 k +1 p
Cn { c ( xn , tk +1 ) =
¦c(x
(29)
tk
At this stage, there is a need to project c onto
n = 1, 2," , N
(33)
where R is the flow region bounded by * (for the
361
moment, we will treat the time derivatives as known functions). Application of Green”s first identity gives
ª Dc º N N «¬ Dt + fc q »¼ [ n dR + ³R Dc <[ n dR
³
R
³* Dc
N
n = 1, 2," , N
Bnm = ³ (v
at time t k 1 . This unconventional approach renders the problem purely parabolic. The first term in Eq.(34) can be approximated as
Dc Dc [ n dR | n R Dt Dt
³ [ dR
if n and m are on the boundaries, and
inflow or noflow
Bnm = 0
(42)
if n and m are on the outflow boundary. As D o f at x n , cn is known and Bnm is not needed. F is a symmetric matrix of order N defined as
(35)
n
R
(41)
*
In what follows, we treat each node of the fixed finite element grid as a moving particle having reached x n
³
B is a symmetric “boundary matrix” of order N , with the elements
Fnm = ³ f [ n[ m dR
(43)
R
Note that this is analogous to what one does in most conventional finite difference schemes. To approximate the time derivatives by finite difference, we write
Dcn cnk +1 k Cn | , Dt 't
dcn cnk +1 cnk | dt 't
W is a diagonal “capacity matrix” of order N defined as
Wnm = G mn ³ f [ n dR
(44)
R
(36)
Since n is viewed as a particle reaching xn at tk +1 , we must use a backward difference scheme. From Eqs.(8)-(9),
Q is an N-dimensional “source vector” defined as
Qn = < n + ³ q[ n dR
(45)
R
where
³* Dc
N
N
*
(46)
*
(c N C ) Q ][ n d*
If n is on the inflow or no-flow boundaries, and
(37)
along the inflow and no-flow boundaries, and
³* Dc
N
(38)
along the outflow boundaries. Substituting these together with Eqs.(25) and (35)-(36) into Eq.(34) leads to the matrix equation
1 · k +1 1 k § ¨ A+ B + F + W ¸c = Q + W C 't ¹ 't ©
(39)
if n is on the outflow or noflow boundaries. As before, < n is not needed if D o f at xn . c k is the N-dimensional vector of c k and corresponding vector of k Cn .
C is the
The final step is to project c k +1 onto moving particles if such particles exist in the flow field. Let p be such a particle, then c kp +1 is computed
c kp +1 = c pk + cˆ N ( x p , tk +1 ) = c pk + N
R
k
according to
where A is a symmetric positive-semidefinite “dispersion matrix” of order N , defined as
Anm = ³ D[n <[ m dR
(47)
(40)
k +1 m
¦ (c m =1
k Cm )[ m ( x p )
(48)
362
where cmk +1 k Cm is equivalent to cˆmk +1 by the virtue of Eq.(32). Equation (48) has been derived with a finite element interpolation scheme for cˆ which is similar to that used for c in Eq.(25). For the next time step, c p is set to have the above value of c p in accordance with Eq.(24), this supersedes the use of Eq.(27). 4. Verification and application of model The verification and application of the proposed method are described as follows. Firstly, the classical one dimensional model was used to compare the numerical results using the proposed method with the analytical solutions. Then, the method was applied to simulate contaminant transport process through a slurry wall as a barrier to prevent leachate pollution from a sanitary landfill. 4.1 Verification of the proposed method The finite element discretization of the one dimensional model domain is shown in Fig.1, which consists of 100 elements and 202 nodes. The velocities were calculated first using the developed method, and the nodal velocities were directly input to the transport model. The analytical solutions for the one dimensional model of advection-dispersion equation in Ref.[20].
Fig.1 Verification of one dimensional model
Different calculation examples were conducted in order to verify the method under different Peclet numbers. Figures 2(a) to 2(c) show the primary results of the one dimensional model compared to the analytical solutions with the variation of the Peclet number ranging from 0.2 to 20. It is clear that the method is able to achieve good estimate of the concentration profile for the wide range of the Peclet number. The method does not cause any over-shooting or under-shooting as well as negative concentration, and only very small numerical dispersion existed at very high Peclet numbers. 4.2 Application of model The slurry walls are in common use as vertical cut-off barriers to allow the containment of liquids in the storage lagoons. Trenches are excavated under cement-bentonite slurry which initially supports the trench during construction and later self-hardens in
situ, hence turns into single-phase. In 2007, in order to contain leachates arising from “piggy-backing” of an existing landfill called the Qizishan landfill located in Suzhou, China, a slurry wall was installed in the downstream of the landfill site. The new wall was installed as a remedial measure for the early wall installed in 1993 due to the implementation of a new phase of landfill to be operated under different working conditions. A sketch for the section across the installed wall can be seen in Fig.3.
Fig.2(a)
Breakthrough cover curve for low Peclet number (Pe = 0.2)
Fig.2(b)
Breakthrough cover curve for Peclet number (Pe = 2)
Fig.2(c) Breakthrough cover curve for high Peclet number (Pe = 20)
Assume that a steady leachate concentration has been generated and keeps conservative in the upstream side of the wall. The models of dvective
363
Fig. 3 Cross section through cement-bentonite slurry wall
and dispersive transport are conducted to calculate the breakthrough curves and predicated the breakthrough times of the worst case. During the simulation process, it is assumed that the barrier retains its chemical integrity with time and no retardation is taken into account for the pollutant transport. It is shown in Fig.4 that it will be 50 years before breakthrough at 0.1% of the source concentration and 80 years before breakthrough at 5% of that for a 600 mm slurry wall. For the demand of contamination c control = 0.001 , it means that the breakthrough c0 time of the slurry wall will be 50 years which exceeds the 40 years’ life cycle of the landfill. Therefore, it is clear that the 600 mm slurry wall is effective to prevent the leachate pollution from landfill inside.
Fig.4
Predicated plot for leachate transport
From the models of advective-dispersive transport, it is also shown that the advective flux acts as the dominant process for contaminant transport across the slurry wall. All the conclusions are similar to the in situ investigation of the contaminant transport through slurry wall in a landfill situated in
Yorkshire, UK by Philip[21] as well as the simulation of solute transport across low-permeability barrier walls by Philip[22]. 5. Conclusions The main conclusions of the results of this study can be reached as follow: (1) A new numerical groundwater flow method has been introduced for estimating simultaneously the velocity vectors and pressure head . It can be used for estimating the velocities directly at nodal elements which can be used as direct input to transport models without need to interpolation process and therefore avoiding the error which may occur especially in heterogeneous media. (2) The advection-dispersion phenomena have been simulated using the velocities at the nodes directly, and the Eulerian-Lagrangian FEM for the governing equations have been established. (3) The proposed method has been verified with the classical one dimnsional model, and the comparison with analytical solutions shows that the method is able to obtain good estimate of the concentration profile for the wide range of the Peclet number. (4) The simulation of contaminant transport process through a slurry wall as a barrier has conducted to evaluate the effect of contaminant control of the slurry wall. The conclusions show that the advective flux acts as the dominant process for contaminant transport, and it at least needs 50 years to result in the breakthrough at 0.1% of the source concentration for a 600 mm slurry wall.
364
[12]
References [1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
GRIFT B. V., GRIFFIOEN J. Modelling assessment of regional groundwater contamination due to historic smelter emissions of heavy metals [J]. Journal of Contaminant Hydrology, 2008, 96: 48-68. ZHANG Qian-fei, WANG Yan-ming and XU Yong-fu. Seepage analysis of landfill foundations in Shanghai Laogang landfill phase IV [J]. Journal of Hydrodynamics, Ser. B, 2006, 18(5): 613-619. ADIL Elkrail, SHU Long-cang and HAO Zhen-chun. Numerical simulation of groundwater dynamics for Songhuajiang River valley in China [J]. Journal of Hydrodynamics, Ser. B, 2004, 16(3): 332-335. LI Yong, WANG Chao and YANG Lin-zhang et al. Influence of seepage face obliquity on discharge of groundwater and its pollutant into lake from a typical unconfined aquifer[J]. Journal of Hydrodynamics, Ser. B, 2007, 19(6): 756-761. OLAF A. C., ALBERT J. V. Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state [J]. Advances in Water Resources, 2007, 30: 1668-1679. CARL C. H., PETER B. and METTE D. et al. Groundwater flow and transport of nutrients through a riparian meadow – Field data and modeling [J]. Journal of Hydrology, 2006, 331: 315-335. YIN Hai-long, XU Zu-xin and LI Huai-zheng et al. Numerical modeling of wastewater transport and degradation in soil aquifer [J], Journal of Hydrodynamics, Ser. B, 2006, 18(5): 597-605. MAZZIA A., PUTTI M. Three-dimensional mixed finite element-finite volume approach for the solution of density-dependent flow in porous media [J]. Journal of Computational and Applied Mathematics, 2006, 185: 347-359. JAVADI A. A., NAJJAR M. M. Finite element modeling of contaminant transport in soils including the effect of chemical reactions [J]. Journal of Hazardous Materials, 2007, 143: 690-701. MARK M. M., CHARLES T. Finite difference approximations for fractional advection ̢ dispersion flow equations[J]. Journal of Computational and Applied Mathematics, 2004, 172: 65-77. JAMES R. C., ALAN J. R. Finite difference modeling of contaminant transport using analytic element flow
[13]
[14]
[15]
[16]
[17] [18]
[19]
[20]
[21]
[22]
solutions [J]. Advances in Water Resources, 2006, 29: 1075-1087. TODD C. R., YU G. Q. Determination of groundwater flownets, fluxes, velocities, and travel times using the complex variable boundary element method [J]. Engineering Analysis with Boundary Elements, 2006, 30: 1030-1044. DOU C. H., WOLDT W. and BOGARDI I. et al. Numerical solute transport simulation using fuzzy sets approach[J]. Journal of Contaminant Hydrology, 1997, 27: 107-126. YOON H. S., HYUN Y. J. and LEE K. K. Forecasting solute breakthrough curves through the unsaturated zone using artificial neural networks [J]. Journal of Hydrology, 2007, 335: 68-77. JEFFREY M. M., CLIFFORD I. V. and DONALD I. S. Groundwater flow with energy transport and water̢ice phase change: Numerical simulations, benchmarks, and application to freezing in peat bogs [J]. Advances in Water Resources, 2007, 30: 966-983. THOMAS L., LI S. G. A finite analytic method for solving the 2-D time-dependent advection-diffusion equation with time-invariant coefficients [J]. Advances in Water Resources, 2005, 28: 117-133. NEUMAN S. P., SOREK S. Eulerian-Lagrangisan methods for advection-dispersion [J]. Water Resources Research, 1982, (18): 1441-1468. BINNING P., CELIA M. A. A forward particle tracking Eulerian-Lagrangisan localized adjoint method for solution of the contaminant transport equation in three dimensions [J]. Advances in Water Resources, 2002, 25: 147-157. THOMAS F. R., MICHAEL A. C. An overview of research on Eulerian ̢ Lagrangian localized adjoint methods[J]. Advances in Water Resources, 2002, 25: 1215-1231. CHEN Chong-xi, LI Guo-min. Theory and model of groundwater flow and solute transport [M]. Wuhan: Press of China University of Geosciences, 1996(in Chinese ). PHILIP L. K. An investigation into contaminant transport processes through single–phase cement-bentonite slurry walls [J]. Engineering Geology, 2001, (60): 209-221. PHILIP T. H., LEONARD F. K. and GEORGE Z. H. Simulation of solute transport across low-permeability barrier walls [J].Journal of Contaminant Hydrology, 2006, 85: 247-270.
2008,20(3):356-364
A NEW NUMERICAL METHOD FOR GROUNDWATER FIOW AND SOLUTE TRANSPORT USING VELOCITY FIELD* ZHANG Qian-fei School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200030, China Shanghai Municipal Engineering Design Institute, Shanghai 200092, China,E-mail: [email protected] LAN Shou-qi Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China WANG Yan-ming Shanghai Municipal Engineering Design Institute, Shanghai 200092, China XU Yong-fu School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200030, China
(Received July 18, 2007, Revised March 9, 2008)
Abstract:A new numerical method for groundwater flow analysis was introduced to estimate simultaneously velocity vectors and water pressure head. The method could be employed to handle the vertical flow under variably saturated conditions and for horizontal flow as well. The method allows for better estimation of velocities at the element nodes which can be used as direct input to transport models. The advection-dispersion process was treated by the Eulerian-Lagrangian approach with particle tracking technique using the velocities at FEM nodes. The method was verified with the classical one dimensional model and applied to simulate contaminant transport process through a slurry wall as a barrier to prevent leachate pollution from a sanitary landfill. Key words:groundwater flow, solute transport, velocity vectors, finite element method, Eulerian-Lagrangian approach, particle tracking technique
1. Introduction Groundwater flow and transport analysis have been an important research topic in the last three decades, due to the fact that many geo-environmental engineering problems have a direct or indirect impact by groundwater flow and solute transport. The transport of dissolved or suspended materials by flowing water is of broad practical significance in relation to environmental protection and resource utilization [1-3]. The transport and fate of chemical or biological species is a major concern in relation to groundwater contamination [4,5]. The process of * Project supported by the National Natural Science Foundation of China (Grant No. 40201024). Biography: ZHANG Qian-fei (1975- ), Male, Ph. D., Senior Engineer
displacement of salts and nutrients in soils has an important bearing on agricultural production [6,7]. Environmental engineering problems in the recent two decades forced the research interest to more sophisticated numerical modeling, in order to predict the contaminant migration in the geological formation more accurately. Many techniques for solving numerically the advection-dispersion equation have been reported, using the most popular techniques, such as the Finite Element Method (FEM) [8, 9], Finite Difference Method (FDM) [10,11] , Boundary Element Method (BEM) [12], Fuzzy Sets Approach (FSA) [13], Artificial Neural Networks method (ANN) [14] , etc. Conventional analysis of groundwater flow was generally made by assuming that the total or pressure head is the only unknown variable. The velocity vectors were calculated after obtaining the pressure head distribution by applying Darcy’s law. However,
357
velocity vectors have to be calculated directly with the numerical method due to the following reasons: (1) to be used as direct input for transport model analysis, (2) to be used for inverse analysis when the measured velocity data from the field is available, and (3) to estimate nodal velocities more accurately especially in case of heterogeneous media. The advection-dispersion equation has been widely used to describe the Fickian transport of pollutants in the atmosphere, oceans, lakes, rivers and subsurface bodies of water. In recent years, there has been much interest in the possibility of using this equation to predict the effect of groundwater flow on the migration of radionuclides from geological repositories of nuclear waste to the biosphere [15]. This latter problem poses a special challenge to the numerical analyst because subsurface transport is often controlled by complex three-dimensional flow patterns which may vary with time [16], in an anisotropic dispersion process whose tonsorial description depends on the velocity decay, sorption, and chemical reactions. presented the Neuman and Sorek[17] Eulerian-Lagrangian method for solving the advection-dispersion equation, where the advection was solved independently at each time step, which they named it conventional method of continuous particle tracking. The residual dispersion was treated by the Eulerian finite element method. The aim of the Eulerian-Lagrangian methods is to combine the simplicity of the fixed Eulerian grid with the computational power of the Lagrangian approach. On the other hand, the Eulerian-Lagrangian method proposed by many researchers is based on interpolating the flow velocities at the element nodes which are estimated from the water flow model [18]. The technique based on the continuous forward particle tracking method which has to obtain the values of concentration, changed by the advection during a time interval, and projected onto the fixed finite element grid [19]. This projection process may yield numerical errors in the solutions. Therefore, the numerical method for estimating the velocity field directly at the element nodes should be used as the direct input to the transport model in order to make a better estimation of concentration caused by advection-dominated transport especially in highly heterogeneous porous media. For conducting transport analysis more accurately and inverting the hydraulic parameters with the water head and velocity vectors, a groundwater flow model is developed herein, which directly estimate the water pressure head as well as the velocity vector distribution in saturated and unsaturated porous media.
2. Groundwater flow analysis considering both pressure head and velocity vectors 2.1 Governing equations The continuity equation can be written in the form
Vi ,i = S s
wh wt
(1)
where Vi (i = 1, 2,3) is Darcy’s velocity components along the Cartesian coordinate axes, h the total head, S s the storage coefficient, and t the time. On the other hand, according to Darcy’s law, the equation of motion can be written as
Vi = K ij (T ) h, j
(2)
where K ij (i, j = 1, 2,3) is the tensor of the variably saturated permeability coefficients and it is a function of volumetric moisture content ( T ). Equations (1) and (2) are the governing equations for flows in saturated-unsaturated porous media with the pressure head and velocity vectors as unknown variables. Solving these two equations simultaneously require the boundary and initial conditions to be specified, given as Initial conditions:
h( x , t )
= h( x, 0)
t =0
Vi ( x, t )
t =0
for total head
= Vi ( x, 0) for velocity vectors
(3a) (3b)
Boundary conditions:
h( x , t )
*1
= hˆ( x , t )
for total head
Vi ( x , t )
*2
= Vˆi ( x , t )
for velocity vectors
K ij (T ) h, j ni
*3
= Qˆ ( x , t )
for flow rate
(3c)
(3d) (3e)
where x is the position vector, ni the unit normal vector along the Cartesian coordinate axes, hˆ the head, Vˆi the known velocity components, and Qˆ the prescribed flow rate. known
total
358
2.2 FEM technique The Galerkin-type finite element technique was implemented to formulate the finite element discretization. Linear hexahedral isoparametric element was used to represent the behavior of the total head h , and quadratic hexahedral isoparametric element was used to express the velocity vectors Vi . According to the Galerkin procedure, to approximately solve the governing equations, let the solution domain be subdivided into a number of finite elements and let a trial solution (h,Vi ) be expressed in the form:
(4b)
wc = <( DC vc) + q wt
20
J =1
where N J is the linear shape function for total head, and N J is the quadratic shape function for the velocity vector components. The formulation with the Galerkin-type FEM produces the following matrix form:
0
ª A « « 0 « 0 « x ¬ 'tC
0 0
A 0 'tC
A y
'tC z
B x º Vx ½ »° ° B y » ° Vy ° ® ¾= B z » ° Vz ° » D ¼ °¯ht 't °¿
0 ½ ° ° 0 ° ° ¾ ® 0 ° ° °¯ht D + 'tF °¿ where
AIJ = ³ ( N I N J )dV , V
BIJi = ³ ( N I Kij (T ) N J , j )dV , V
where V is the elemental volume, S the surface area of the element, ni the normal vector to the surface, and the others are the same as defined above.
(4a)
J =1
Vi ( x , t ) = ¦ N J ( x )VJi (t )
S
3. Advection-dispersion transport analysis 3.1 Governing equations In the theory presented here, we will consider only the main processes of manly controling the solute transport within the geological formation, namely, the advection and dispersion processes. Consider the advection-dispersion equation
8
h( x , t ) = ¦ N J ( x ) hJ (t )
FIJ = ³ ( N I N JVi ni )dV
(5)
where c is the concentration, t the time, the gradient operator, D the diffusion tensor, v the seepage velocity vector, and q the source term. Here, D is a symmetric positive-semidefinite tensor. The equation is subject to the initial and boundary conditions
c( x, 0) = C0 ( x )
(7)
( DC vc)
(8)
where * is the boundary, n the unit vector normal to * , pointing outwards, C0 , C and Q are the initial concentration, varying concentration and mass flux respectively, and D a controling parameter, which determines the type of prevailing boundary condition , that is, if D o f , Eq.(8) is the boundary condition for the concentration, if D = 0 , it is the boundary condition for the mass flux, in other cases, it is a mixed condition. Equation (8) applies to inflow and no-flow boundaries. Along outflow boundaries, it is common to assume that
DC
D w = + v Dt wt
DIJ = ³ ( N I Ss N J )dV ,
we can rewrite Eq.(6) in the Lagrangian form as
V
(9)
Using the material derivative operator
CIJi = ³ ( N I ,i N J )dV , V
(6)
(10)
359
Dc = <( Dc ) fc + q Dt
(11)
where f = v . Here c no longer represents the concentration at a point in space-time, but rather the concentration of a fluid particle moving at the velocity v . The path of this particle is described by the material derivative of x , which leads to the characteristic equation
Dx =v Dt
(13)
(14)
(15)
the Cauchy condition
vc
(16)
along inflow boundaries, and
vc
(17)
along the no-flow boundaries. The conditions along the outflow boundaries have no effect on c and are thus irrelevant to it. Clearly, the value of c for a given fluid particle remains constant as the latter is advected through the flow field. Thus, the “advection problem” defined by Eqs.(14)-(17) can be solved independently of cˆ . The residual concentration, cˆ , must satisfy
Dcˆ = <( Dcˆ) fcˆ + q + g Dt where
cˆ( x , 0) = 0
(20)
and the conditions along the inflow boundaries read (21)
where h = Dc
subjected to the initial condition
c ( x, 0) = C0 ( x )
subjected to the homogeneouse initial condition
(12)
One way to perform such a decomposition is to let c satisfy the homogeneous first-order partial differential equation
Dc wc = + vc = 0 Dt wt
(19)
( Dcˆ + vcˆ)
Equations (11)-(12) subjected to Eqs.(7)-(9) can be replaced by two sets of equations, one in terms of c , and the other in terms of cˆ , where
c = c + cˆ
g = <( Dc ) fc
(18)
Dcˆ + D (cˆ C ) * = u
(22)
where u = Q D c + h , and along the outflow boundaries
Dcˆ
(23)
Note that g , h , u and c play the role of prescribed source functions as c is known. We will refer to Eqs. (18)-(23) as the residual “dispersion problem”. It can be seen that the hyperbolic-parabolic advection-dispersion problem defined by Eqs.(6)-(9) can be formally decoupled into a purely hyperbolic “advection problem” defined in terms of c , and another predominantly parabolic residual “dispersion problem” defined in terms of cˆ . The approach is to first solve the advection problem for c , as the latter is dependent of cˆ and then solve the residual dispersion problem for cˆ . Another alternative is to solve Eq.(11) in its Lagrangian form without decomposing c into c and cˆ , and this approach is suitable for the regions with steep concentration gradients. The numerical scheme described below utilizes both of these alternatives jointly in an adaptive manner. 3.2 Eulerian-Lagrangian FEM approach Suppose that c is known at time tk , and we wish to compute it at tk +1 = tk + 't . First we set
c ( x , t k ) { c( x , tk ) ,
cˆ( x, tk ) { 0
(24)
We then solve the advection problem for c ( x , tk +1 ) by using the continuous forward particle tracking in the vicinity of steep fronts, and single-step reverse particle tracking away from such fronts. The residual dispersion problem is solved for c or cˆ by the finite element method.
360
The finite element equations in terms of c will be given. Let c ( x , t ) is approximated by a finite
the node, n , of the finite element grid. The projected values of c are denoted by k Cn such that
element function, c N ( x, t ) , defined as Mn
N
N
c( x , t ) | c ( x , t ) = ¦ cn (t )[ n ( x )
(25)
k
n =1
where N is the number of grid points or nodes, cn the concentration at node n , and [ n a basis function satisfying
[ n ( xm ) = G nm
(26)
where xm is x at node m , and G nm is the Kronceker delta. 3.2.1 Continuous forward particle tracking If there is a steep concentration front, the approach is to introduce moving particles at strategic points around the front and continuously track their positions along the paths. For any existing particle, p , located at point x p at time tk , the value of c is assigned to be c pk = c ( x p , tk ) for the duration of the time interval (tk , tk +1 ) . By virtue of Eqs.(24)-(25), we have N
c pk = c( x p , tk ) | ¦ cmk [ m ( x p )
(27)
m =1
For any new particle, r , at ( xr , tk ) the value of
c is assigned to be crk = c ( xr , tk ) . For a new particle along an inflow boundary, the value of crk is assigned as
crk =
DC + Q v
(28) xr ,tk
according to Eq.(16) for the duration of the time interval. At the end of the time step, each particle, p , reaches a new position
x
k p
=x +
³ vDt
k +1 p
, tk +1 )
p =1
Mn
(30)
Clearly, k Cn can be viewed as c at tk of a fictitious moving particle which reaches the node n at time tk +1 . 3.2.2 Single-step reverse particle tracking In the areas where the concentration gradients are mild, k Cn is computed by sending a fictitious particle from each node, n , backward to the point tk +1 k
xn { xn
³ vDt
(31)
tk
in each time step, which means that a particle leaving from k xn at tk will reach the grid point location, xn , exactly at tk +1 . If the velocity field does not change with time and 't is fixed, then k xn remains constant for each n and needs to be computed only once. By the virtue of Eqs.(24)-(25), k Cn can be expressed as k
N
Cn { c ( k xn , tk ) | ¦ cmk [ m ( k xn )
(32)
m =1
This requires much less computer time and storage than the procedure with the continuous forward particle tracking. However, the method cannot be used near sharp fronts because there exhibits numerical dispersion near such fronts. 3.2.3 Dispersion computation by finite element method Applying the Galerkin orthogonalization procedure to Eq.(11) leads to
ª Dc º <( Dc N ) + fc N q » [ n dR = 0, « R Dt ¬ ¼
³
tk +1 k +1 p
Cn { c ( xn , tk +1 ) =
¦c(x
(29)
tk
At this stage, there is a need to project c onto
n = 1, 2," , N
(33)
where R is the flow region bounded by * (for the
361
moment, we will treat the time derivatives as known functions). Application of Green”s first identity gives
ª Dc º N N «¬ Dt + fc q »¼ [ n dR + ³R Dc <[ n dR
³
R
³* Dc
N
n = 1, 2," , N
Bnm = ³ (v
at time t k 1 . This unconventional approach renders the problem purely parabolic. The first term in Eq.(34) can be approximated as
Dc Dc [ n dR | n R Dt Dt
³ [ dR
if n and m are on the boundaries, and
inflow or noflow
Bnm = 0
(42)
if n and m are on the outflow boundary. As D o f at x n , cn is known and Bnm is not needed. F is a symmetric matrix of order N defined as
(35)
n
R
(41)
*
In what follows, we treat each node of the fixed finite element grid as a moving particle having reached x n
³
B is a symmetric “boundary matrix” of order N , with the elements
Fnm = ³ f [ n[ m dR
(43)
R
Note that this is analogous to what one does in most conventional finite difference schemes. To approximate the time derivatives by finite difference, we write
Dcn cnk +1 k Cn | , Dt 't
dcn cnk +1 cnk | dt 't
W is a diagonal “capacity matrix” of order N defined as
Wnm = G mn ³ f [ n dR
(44)
R
(36)
Since n is viewed as a particle reaching xn at tk +1 , we must use a backward difference scheme. From Eqs.(8)-(9),
Q is an N-dimensional “source vector” defined as
Qn = < n + ³ q[ n dR
(45)
R
where
³* Dc
N
N
*
(46)
*
(c N C ) Q ][ n d*
If n is on the inflow or no-flow boundaries, and
(37)
along the inflow and no-flow boundaries, and
³* Dc
N
(38)
along the outflow boundaries. Substituting these together with Eqs.(25) and (35)-(36) into Eq.(34) leads to the matrix equation
1 · k +1 1 k § ¨ A+ B + F + W ¸c = Q + W C 't ¹ 't ©
(39)
if n is on the outflow or noflow boundaries. As before, < n is not needed if D o f at xn . c k is the N-dimensional vector of c k and corresponding vector of k Cn .
C is the
The final step is to project c k +1 onto moving particles if such particles exist in the flow field. Let p be such a particle, then c kp +1 is computed
c kp +1 = c pk + cˆ N ( x p , tk +1 ) = c pk + N
R
k
according to
where A is a symmetric positive-semidefinite “dispersion matrix” of order N , defined as
Anm = ³ D[n <[ m dR
(47)
(40)
k +1 m
¦ (c m =1
k Cm )[ m ( x p )
(48)
362
where cmk +1 k Cm is equivalent to cˆmk +1 by the virtue of Eq.(32). Equation (48) has been derived with a finite element interpolation scheme for cˆ which is similar to that used for c in Eq.(25). For the next time step, c p is set to have the above value of c p in accordance with Eq.(24), this supersedes the use of Eq.(27). 4. Verification and application of model The verification and application of the proposed method are described as follows. Firstly, the classical one dimensional model was used to compare the numerical results using the proposed method with the analytical solutions. Then, the method was applied to simulate contaminant transport process through a slurry wall as a barrier to prevent leachate pollution from a sanitary landfill. 4.1 Verification of the proposed method The finite element discretization of the one dimensional model domain is shown in Fig.1, which consists of 100 elements and 202 nodes. The velocities were calculated first using the developed method, and the nodal velocities were directly input to the transport model. The analytical solutions for the one dimensional model of advection-dispersion equation in Ref.[20].
Fig.1 Verification of one dimensional model
Different calculation examples were conducted in order to verify the method under different Peclet numbers. Figures 2(a) to 2(c) show the primary results of the one dimensional model compared to the analytical solutions with the variation of the Peclet number ranging from 0.2 to 20. It is clear that the method is able to achieve good estimate of the concentration profile for the wide range of the Peclet number. The method does not cause any over-shooting or under-shooting as well as negative concentration, and only very small numerical dispersion existed at very high Peclet numbers. 4.2 Application of model The slurry walls are in common use as vertical cut-off barriers to allow the containment of liquids in the storage lagoons. Trenches are excavated under cement-bentonite slurry which initially supports the trench during construction and later self-hardens in
situ, hence turns into single-phase. In 2007, in order to contain leachates arising from “piggy-backing” of an existing landfill called the Qizishan landfill located in Suzhou, China, a slurry wall was installed in the downstream of the landfill site. The new wall was installed as a remedial measure for the early wall installed in 1993 due to the implementation of a new phase of landfill to be operated under different working conditions. A sketch for the section across the installed wall can be seen in Fig.3.
Fig.2(a)
Breakthrough cover curve for low Peclet number (Pe = 0.2)
Fig.2(b)
Breakthrough cover curve for Peclet number (Pe = 2)
Fig.2(c) Breakthrough cover curve for high Peclet number (Pe = 20)
Assume that a steady leachate concentration has been generated and keeps conservative in the upstream side of the wall. The models of dvective
363
Fig. 3 Cross section through cement-bentonite slurry wall
and dispersive transport are conducted to calculate the breakthrough curves and predicated the breakthrough times of the worst case. During the simulation process, it is assumed that the barrier retains its chemical integrity with time and no retardation is taken into account for the pollutant transport. It is shown in Fig.4 that it will be 50 years before breakthrough at 0.1% of the source concentration and 80 years before breakthrough at 5% of that for a 600 mm slurry wall. For the demand of contamination c control = 0.001 , it means that the breakthrough c0 time of the slurry wall will be 50 years which exceeds the 40 years’ life cycle of the landfill. Therefore, it is clear that the 600 mm slurry wall is effective to prevent the leachate pollution from landfill inside.
Fig.4
Predicated plot for leachate transport
From the models of advective-dispersive transport, it is also shown that the advective flux acts as the dominant process for contaminant transport across the slurry wall. All the conclusions are similar to the in situ investigation of the contaminant transport through slurry wall in a landfill situated in
Yorkshire, UK by Philip[21] as well as the simulation of solute transport across low-permeability barrier walls by Philip[22]. 5. Conclusions The main conclusions of the results of this study can be reached as follow: (1) A new numerical groundwater flow method has been introduced for estimating simultaneously the velocity vectors and pressure head . It can be used for estimating the velocities directly at nodal elements which can be used as direct input to transport models without need to interpolation process and therefore avoiding the error which may occur especially in heterogeneous media. (2) The advection-dispersion phenomena have been simulated using the velocities at the nodes directly, and the Eulerian-Lagrangian FEM for the governing equations have been established. (3) The proposed method has been verified with the classical one dimnsional model, and the comparison with analytical solutions shows that the method is able to obtain good estimate of the concentration profile for the wide range of the Peclet number. (4) The simulation of contaminant transport process through a slurry wall as a barrier has conducted to evaluate the effect of contaminant control of the slurry wall. The conclusions show that the advective flux acts as the dominant process for contaminant transport, and it at least needs 50 years to result in the breakthrough at 0.1% of the source concentration for a 600 mm slurry wall.
364
[12]
References [1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
GRIFT B. V., GRIFFIOEN J. Modelling assessment of regional groundwater contamination due to historic smelter emissions of heavy metals [J]. Journal of Contaminant Hydrology, 2008, 96: 48-68. ZHANG Qian-fei, WANG Yan-ming and XU Yong-fu. Seepage analysis of landfill foundations in Shanghai Laogang landfill phase IV [J]. Journal of Hydrodynamics, Ser. B, 2006, 18(5): 613-619. ADIL Elkrail, SHU Long-cang and HAO Zhen-chun. Numerical simulation of groundwater dynamics for Songhuajiang River valley in China [J]. Journal of Hydrodynamics, Ser. B, 2004, 16(3): 332-335. LI Yong, WANG Chao and YANG Lin-zhang et al. Influence of seepage face obliquity on discharge of groundwater and its pollutant into lake from a typical unconfined aquifer[J]. Journal of Hydrodynamics, Ser. B, 2007, 19(6): 756-761. OLAF A. C., ALBERT J. V. Two-dimensional concentration distribution for mixing-controlled bioreactive transport in steady state [J]. Advances in Water Resources, 2007, 30: 1668-1679. CARL C. H., PETER B. and METTE D. et al. Groundwater flow and transport of nutrients through a riparian meadow – Field data and modeling [J]. Journal of Hydrology, 2006, 331: 315-335. YIN Hai-long, XU Zu-xin and LI Huai-zheng et al. Numerical modeling of wastewater transport and degradation in soil aquifer [J], Journal of Hydrodynamics, Ser. B, 2006, 18(5): 597-605. MAZZIA A., PUTTI M. Three-dimensional mixed finite element-finite volume approach for the solution of density-dependent flow in porous media [J]. Journal of Computational and Applied Mathematics, 2006, 185: 347-359. JAVADI A. A., NAJJAR M. M. Finite element modeling of contaminant transport in soils including the effect of chemical reactions [J]. Journal of Hazardous Materials, 2007, 143: 690-701. MARK M. M., CHARLES T. Finite difference approximations for fractional advection ̢ dispersion flow equations[J]. Journal of Computational and Applied Mathematics, 2004, 172: 65-77. JAMES R. C., ALAN J. R. Finite difference modeling of contaminant transport using analytic element flow
[13]
[14]
[15]
[16]
[17] [18]
[19]
[20]
[21]
[22]
solutions [J]. Advances in Water Resources, 2006, 29: 1075-1087. TODD C. R., YU G. Q. Determination of groundwater flownets, fluxes, velocities, and travel times using the complex variable boundary element method [J]. Engineering Analysis with Boundary Elements, 2006, 30: 1030-1044. DOU C. H., WOLDT W. and BOGARDI I. et al. Numerical solute transport simulation using fuzzy sets approach[J]. Journal of Contaminant Hydrology, 1997, 27: 107-126. YOON H. S., HYUN Y. J. and LEE K. K. Forecasting solute breakthrough curves through the unsaturated zone using artificial neural networks [J]. Journal of Hydrology, 2007, 335: 68-77. JEFFREY M. M., CLIFFORD I. V. and DONALD I. S. Groundwater flow with energy transport and water̢ice phase change: Numerical simulations, benchmarks, and application to freezing in peat bogs [J]. Advances in Water Resources, 2007, 30: 966-983. THOMAS L., LI S. G. A finite analytic method for solving the 2-D time-dependent advection-diffusion equation with time-invariant coefficients [J]. Advances in Water Resources, 2005, 28: 117-133. NEUMAN S. P., SOREK S. Eulerian-Lagrangisan methods for advection-dispersion [J]. Water Resources Research, 1982, (18): 1441-1468. BINNING P., CELIA M. A. A forward particle tracking Eulerian-Lagrangisan localized adjoint method for solution of the contaminant transport equation in three dimensions [J]. Advances in Water Resources, 2002, 25: 147-157. THOMAS F. R., MICHAEL A. C. An overview of research on Eulerian ̢ Lagrangian localized adjoint methods[J]. Advances in Water Resources, 2002, 25: 1215-1231. CHEN Chong-xi, LI Guo-min. Theory and model of groundwater flow and solute transport [M]. Wuhan: Press of China University of Geosciences, 1996(in Chinese ). PHILIP L. K. An investigation into contaminant transport processes through single–phase cement-bentonite slurry walls [J]. Engineering Geology, 2001, (60): 209-221. PHILIP T. H., LEONARD F. K. and GEORGE Z. H. Simulation of solute transport across low-permeability barrier walls [J].Journal of Contaminant Hydrology, 2006, 85: 247-270.