- Email: [email protected]
A boundary element and particle tracking model for advective transport in zoned aquifers
Journal of Hydrology, 124 (1991) 159~176
159
Elsevier Science Publishers B.V., Amsterdam
[2]
A B O U N D A R Y E L E M E N T A N D PARTICLE TRACKING MODEL FOR ADVECTIVE T R A N S P O R T IN ZONED AQUIFERS
P, LATINOPOULOS and K. KATSIFARAKIS
Department of Civil Engineering, School of Technology, Aristotle University of Thessaloniki, Thessaloniki (Greece) (Received J a n u a r y 4, 1990; accepted after revision July 29, 1990)
ABSTRACT Latinopoulos, P. and Katsifarakis, K., 1991. A boundary element and particle tracking model for advective transport in zoned aquifers. J. Hydrol., 124:159 176. One of the most crucial factors t h a t affects both the accuracy and efficiency of a groundwater transport model is the flow simulation method which produces the velocity field. This paper presents a boundary element technique t h a t incorporates the zoning approach to allow aquifer heterogeneity. First, the accuracy of the flow model is verified by comparing it with analytic solutions for well-pumping and injection problems. The boundary element flow model is then coupled with a particle tracking procedure for advection-based transport simulation. The combined model also calculates streamlines, travel times and b r e a k t h r o u g h curves and is therefore a powerful tool for simulation and management of advection-dominated contamination problems.
INTRODUCTION
Current experience in the field of water resources management shows that the contamination of groundwater reserves is a widespread problem. This fact accounts for the remarkable increase in the development and use of groundwater simulation models in recent years. The main concern of a groundwater modeller is to produce a reliable numerical code for the investigation of a usually complex system in order to understand the physical processes that govern its response to past, present and future stresses and to evaluate the efficiency of remedial actions involved in proposed restoration schemes. The case of solute transport modelling is much more difficult than modelling groundwater flow alone. Moreover, the flow problem is an element of the transport problem, which means that possible inconsistencies arising from the former would inevitably be retained, if not magnified, within the simulation of the latter. Apart from the scarcity of solute concentration data and the uncertainties with regard to both the pollution sources and the transport parameters, it is well-known that the numerical models of solute transport suffer from inaccuracies, due mainly to numerical dispersion. The last is a particularly acute problem in advection-dominated transport and a common
0022-1694/91/$03.50
cC'. 1991
Elsevier Science Publishers B.V.
160
P. LATINOPOULOS AND K. KATSIFARAKIS
practice to tackle it is to separate the advective and dispersive parts. Nevertheless, a crucial factor in these types of simulation is the accurate estimation of the water velocity, as it directly influences not only the dispersive flux, but also the dominating advective transport. In addition, there is a class of problems where advection is assumed to be the sole driving transport mechanism, and therefore dispersion is ignored in the simulation procedure. In most of the existing numerical codes the calculation of the velocity field is based on solution of the flow equation using either the finite difference method or the finite element method. Experience shows that direct use of these domain-type methods often results in inaccurate and unstable solutions of the transport problem (Galeati and Gambolati, 1989). To overcome the accuracy problem various techniques have been proposed, principally for the finite element method, such as solving the flow problem directly for velocities (Yeh, 1981) or using interpolation functions that yield continuous velocities. Within the same context Frind and Matanga (1985) proposed a different approach by using stream functions as well as potentials in a finite element 'dual formulation' model for steady flow in the vertical plane of the homogeneous, anisotropic aquifer. Another critical point is the accuracy of the numerical solutions near singularities, such as pumping or injection wells. This is very important in groundwater quality management problems since at least one water supply or injection well is always under operation. The strong dependence of the solution accuracy near the wells upon the domain mesh discretization has led developers of domain-type models to incorporate local adjustments and corrections by employing complex numerical algorithms either for the finite difference method (Lerner, 1989) or for the finite element method (Charbeneau and Street, 1979a,b). Finally, when new locations of wells are required the mesh has to be changed in order to retain the same degree of accuracy for all operating wells. A method that can overcome the above problems is the boundary element method. Firstly it produces a continuous velocity field in the interior of the flow domain because it is grid-independent. That holds even for piecewise homogeneous aquifers, provided that physical discontinuities in transmissivity are discretized through interzonal boundaries. Secondly, it is very accurate within the domain and especially so near wells, where the natural representation of the solution reproduces the corresponding analytic ones (Latinopoulos et al., 1982; Latinopoulos, 1986). Thirdly, it requires much less input data and computational time than the domain-type solutions and, particularly for solute transport simulation, a minimal amount of bookkeeping (Kemblowski, 1986). These advantages led us to develop a boundary element model that accurately and efficiently calculates the water-pore velocities; it is therefore very useful when advection is either the sole or the major driving mechanism of solute transport. The model is an extension of a previous one (Latinopoulos, 1986) mainly in that it handles aquifer heterogeneity by employing the zoning technique. One essential characteristic of the model is that all integrals
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECT1VE TRANSPORT
]61
included in head and velocity calculations are evaluated by analytic expressions, leading to very accurate results, as will be shown later. Another feature of the model is that it can be efficiently combined with a particle-tracking numerical code that considers advection-based transport. The ability of the combined algorithm to simulate solute transport in heterogeneous aquifers makes it a useful tool in many applications concerned with groundwater quality management. In principle there are two limitations to the presented boundary elementbased approach. The first concerns the steady-state assumption for the flow problem. A possible extension to handle unsteady flow would heavily tax the overall efficiency of the specific boundary element technique. Moreover, in the class of transport problems for which the method is proposed the travel times of contaminant transport are large (in the order of years), so that the steadystate assumption is generally acceptable. The second drawback of the model refers to the overall reliability of the zoning technique. With this approach, which is very popular in boundary element solutions for both porous (Brebbia and Chang, 1979; Chugh and Falvey, 1984; Lafe et al., 1981) and fractured aquifers (Liggett and Medina, 1987; Shapiro and Andersson, 1983) the heterogeneous flow domain is divided into a number of zones which are each assumed to be homogeneous. The reported limitation arises from the fact that, as the number of zones increases, the efficiency of the boundary element technique decreases. This was the motivation for other boundary element-based approaches to heterogeneity, such as the introduction of specific functional permeability distributions (Cheng, 1984) and the expansion of the solution into perturbation serie.s (Lafe and Cheng, 1987). In conclusion, taking into consideration the facts that that in real-world applications the permeability information is in most cases insufficient and that, in general, uncertainties with regard to pollution histories and parameters are the rule rather than the exception, a restriction to a small number of zones with constant parameters may sometimes be dictated even by the model's calibration procedure (Sophocleous, 1984). THEORETICAL BACKGROUND
We consider the case of a plane, horizontal steady-state groundwater flow in a confined heterogeneous aquifer. The governing equation is obtained by combining the equation of conservation of mass and Darcy's law (Bear, 1979) and is
~ f c~¢\
O(T~
+ q = 0
(1)
where 4) is the piezometric head, T is the aquifer's transmissivity (generally a function of space) and q is a source function. To formulate our model, we assume that the non-homogeneous aquifer can be divided into a number of zones of constant transmissivity. For the sake of simplicity we also assume a
162
P. LATINOPOULOS AND K. KATSIFARAKIS
generally isotropic material, although anisotropy can be easily treated by a geometric scaling procedure, provided t hat zonal boundaries can be equally scaled. Under these assumptions and considering that the source term represents only point singularities, such as wells, eqn. (1) becomes
The above equation is valid for each single zone characterized by its own constant transmissivity value T. In this equation Qw is the pumping rate of the wth well (which is negative if it is an injection well) and 5w, which equals 5(x - x~:, y - Yw), is the Dirac delta function; x~, and Yw are the location coordinates of the operating wells (w = 1, W). As stated above the problem of a zoned aquifer can be solved by satisfying eqn. (2) in each subregion (zone), thus decomposing it into as many subproblems as the total number of zones with different transmissivities (Bear, 1979). To obtain simultaneous solutions for ¢ in all zones both external and interzonal boundary conditions must be specified. The former need no further comment while the latter are obtained by considering: (a) A compatability equation which expresses the fact t hat the piezometric heads are the same when a point j at an interzonal boundary, dividing two adjacent zones k and k + 1, is approached from both sides, i.e. Cj.k =
Oj,k~l
(3)
(b) The continuity equation, expressing the equality of the seepage rates across the boundary, i.e.
T~ (c~¢/~n)j,~ -
T~_l(80/On)2,k÷l
(4)
where Th and Tk+ 1 are the constant transmissivities of the adjacent homogeneous zones k and k + 1, respectively, and n is the distance measured along the normal to the boundary. The pore-water velocities of the flow within each homogeneous and isotropic zone are given by expressions
E
~. c~x
(5)
V~, -
K (?0
(6)
where Vx and V~ are the components of the velocity vector in the directions of the cartesian coordinates x and y, respectively, K is the hydraulic conductivity of the specific porous material of each zone, and e is its porosity. By definition, the transmissivity of each zone is T which is equal to K d where d is the aquifer's thickness. In what follows all K, d and ~ are assumed to be constant within each zone (taken usually equal to their average values), but varying between different zones.
BOUNDARY ELEMENT AND PART1CLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
163
If the variation of the piezometric head is known, normally after the solution of eqn. (2), then the velocity field can be calculated by employing eqns. (5) and (6). The velocity field can next be used to trace streamlines and to calculate travel times between specific points along them. This is very useful in simulating pure advective transport of a non-dispersive and non-reactive solute, but is also required in the general case where the two mechanisms of solute transport, i.e. advective transport and dispersion, are both important. In such problems one must solve the flow equation to produce the water velocities which are needed to calculate both the advective terms and the dispersion coefficients (Bear, 1979). INTEGRAL FORMULATIONOF THE FLOWMODEL In order to solve the governing partial differential equation (eqn. (2)) the boundary element method (BEM) is used. After incorporating the boundary conditions for the external boundaries in any combination of the three possible types, namely Dirichlet, Neumann and Cauchy, the integral formulation of the problem can be obtained by applying the standard BEM approach (Brebbia, 1978; Liggett and Liu, 1983). Although applications of the BEM with linear elements in zoned problems are reported in the literature (Liggett and Medina, 1987), in the present model constant-boundary elements, i.e. node-centred elements, are used, mainly for three reasons: (a) it is easier to develop and simpler to employ the analytical integration formulae for the heads, but especially for the flow velocities, using constant elements rather than higherapproximation ones, (b) specific terms of these formulae are incorporated in the particle-tracking pkocedure, presented later in this paper, and (c) the accuracy obtained even with this simple discretization has proven to be very satisfactory in previous applications (Latinopoulos, 1986) as well as in the present case. In brief, the procedure followed in solving the groundwater flow problem in zoned aquifers is based on the simple fact that each zone is treated separately as a single homogeneous unit in which the BEM is applied. Thus, the boundary of each zone is discretized in linear segments along which the values of the piezometric head ¢ and its derivative 3~/~n, out and normal to the boundary, are assumed to be constant and equal to the respective values Cj and O~bj/?nat the midpoint j of each segment. The boundary conditions for the external boundary of each zone are either prescribed flux or prescribed head or a combination of the two (for Cauchy-type boundary conditions), whereas for the interzonal boundary elements neither the head nor its gradient is prescribed. The final solution for the whole domain can be obtained by a proper coupling of the different zones. Performing the standard procedure of the boundary element method the following set of equations is constructed for each individual zone (Latinopoulos, 1986).
[h] {b}
= [D] {3~/?n} + {P].
(7)
164
P. LATINOPOULOS AND K. KATSIFARAK[S
\ j+l
/'
/~
5
(a)
Cb)
Fig. 1. Notation for the boundary element formulation.
where
A,j
-
! ~iJ'fi ~dFj
z: Dij
=
f ln(~)dYj
ifi
¢ j
if i
j
(8)
(9)
In the above three equations (eqns. (8)-(10)) i is the node of the ith element, j is a point moving along the j t h element, and w is a well location. The corresponding distances rij and riw are shown in Fig. l(a). The system of eqns. (7) is next rearranged in such a way that all the unknowns are on the left-hand side. The final set of equations for the whole domain is constructed by assembling all similar sets of BEM equations of the existing zones into a single global set. In addition, the individual zones are coupled within this global set of equations by incorporating the compatibility equations (3) and (4) for the interzonal boundary elements. The set of equations so assembled is finally solved by Gaussian elimination. It should be noted here that, although more efficient structures of the linear system are reported in the literature, we choose the above simple method because, as will be shown later, our problem requires the solution of only one system, while the time-consuming part in the simulations is the particle-tracking module.Therefore the relevant reduction in efficiency of the full model is a minor issue. When, through the solution of the system of equations, the values of ~bj and c~¢j/~n are known for all the external and interface boundary elements, the piezometric head distribution inside any zone can be calculated. Assuming that
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
165
point i is now within a specific zone which is bounded by N elements (where N stands for the total number of external and interface boundary elements of that zone only), its piezometric head is calculated in exactly the same way as in the case of a single homogeneous aquifer, i.e. through solution of the relation 2uq~i
= ~ __~ in ( 1 )
+ 2N ~b2!rij'n dVj
w: !
j=l
drj
3 In r~
(11)
In the above eqn. (11) the index k characterizing the specific zone for which it applies is omitted for the sake of simplicity (e.g. W and N should read Wk, Nh and so on). The line integrals appearing in eqns. (8), (9) and (11) have been analytically evaluated in a previous work (Ingham et al., 1981) and they have also been applied to various groundwater modelling studies (Latinopoulos et al., 1982; Latinopoulos, 1986). Using the relations given in eqns. (8) and (9) and according to the notation shown in Fig. l(b) these integrals are A°
=
)'
l if i ¢ j
Di~ = 13cos~(1 - ln/a) + /2cosfl(1 - lnl2) - /3~'sin~
)
(12)
and
Ao = 7c Dij
= Ii 1 - In
~))I if i
= j
(13)
The next step in the development of the model is the calculation of the velocity field within the different zones of the aquifer. This can be achieved by evaluating the gradient of the piezometric head from eqn. (11) and substituting it into eqns. (5) and (6). Using the analytical expressions of eqn. (12) (holding for internal points only) the following analytical formulae for the velocity field within each particular zone are obtained (Latinopoulos, 1986).
g S ~ Qwcw Vx = 2~n[ .... L1 Tr~
N ~ 1 Cj
2 1 ~-n I_/asin~ + ~
ll(12c2cosot + laclcosfl ) (/2/a) 2 sin7 ';c°t~ + In
(14)
166
P. LATINOPOULOS AND K. KATSIFARAKIS
K
Vv "
f ~
Qwd w
N
2nn [ w:l z_, Tr~
I1(I,fl2COS ~ + 1 3 d l c o s f l )
~ 1 ~j j=
' (12l:J ~ sin7
- j l~[_13sin~ + ~ -
:,cot~ * In
(15)
where cl
=
xj
c2 cu,
=
xi
dl
-
Yj-
(16a)
Yi
xj÷l -
xi
d2
-
Y~I
Yl
x~: -
xi
d~.
-
Yu - Yi
(16b) (16c)
As mentioned in that work (Latinopoulos, 1986), the main advantages of the straightforward procedure presented in calculating heads and velocities are high accuracy and computational efficiency. Especially in the presence of singularities, such as pumping and injection wells, their natural representation by the logarithmic term in the head relation and the derived gradient expressions in the velocity ones, is an exact reproduction of the corresponding analytical solutions for infinite confined aquifers. This can be easily seen by inspecting the well terms in eqns. (11), (14) and (15) which are the exact analtyical solutions for homogeneous aquifers of infinite extent, whereas the remaining terms approximate the boundary effects. This specific property of the boundary element solutions is a major advantage over the domain-type ones, which, as mentioned in the introduction, require a lot of extra development and computational effort to simulate near-well behaviour accurately. VERIFICATION OF THE NUMERICALFLOWMODEL In this part the boundary element model is tested for accuracy against an analytic solution. Having already verified the model for the resulting head (Latinopoulos et al., 1982) and for velocity values (Latinopoulos, 1986) in homogeneous aquifers, we are now interested in its performance in a zoned aquifer. Looking in the available literature for analytic solutions, the most suitable that we found is a generalized solution for drawdown due to pumping in a rectangular aquifer having two zones with distinct properties (Chan et al., 1977). This solution requires some mathematical manipulations to be tailored to our case and we now briefly describe it. We consider two-dimensional steady flow in a rectangular aquifer consisting of two rectangular zones of constant but different transmissivities T1 and T2, as shown in Fig. 2 Three sides of this aquifer are assumed to be impermeable while the fourth is kept at a constant head q~. A single well, located in Zone 1 at a point with coordinates ~ and ~ in a cartesian system, also shown in Fig. 2, is pumping continuously at a rate Q0. Thus, the steady-state analytic solution of this flow problem can be obtained by considering eqn. (2) as the governing equation. From the general, transient solution of Chan et al. (1977) we obtain the following analytic expressions for the steady-state head distribution.
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
167
,y Ii/////11/I/////I//////////
ZONE
2
ZONE1
(~,~) b
= const
c),p/8x-o1
i
~L
X 0",1/////I
/ / / t z ///////////////////
O ~--c
~'
o~O/Oy=o 4
Fig. 2. The idealized two-zone aquifer with a single pumping well.
For 0
~< x
~< c
¢(x,s)
= ¢¢-~
Oo[
a-C
+ 2T1 ~ cosflmy COSflm~coshflmx sinhflm(a - ~)1 For c
~< x
¢(x.y)
=
~<
4o--~1
For ~
~< x
c~(x,s)
=
+
Oo[
a -
cosflmy c o s f l m , sinhfl.~(a - ~) A ( x - c.c;flm)]
+2
(18)
~< a
~o -
2
(17)
Oo[
~
a -
x
~-. eosflmy Cosfl.,~ sinhflm(a - x) A(~ - c, c;flm),~=~L
fl,.A(a - c. c; tim)
(19)
168
P. LATINOPOULOS AND K. KATSIFARAKIS
where tim = m~/b and A ( X , Y ; f l m ) = TlcoshflmX c o s h f l ~ Y + T2 s i n h f l m X sinh tim Y B e c a u s e we are primarily i n t e r e s t e d in c h e c k i n g the a c c u r a c y of the velocity field r e s u l t i n g from our b o u n d a r y element model, in addition to the above, we need to o b t a i n the velocity c o m p o n e n t s in the x and y directions. T a k i n g the g r a d i e n t of the head from eqns. (17) (19) and s u b s t i t u t i n g into eqns. (5) and (6) we o b t a i n the following expressions. For 0 ~< x ~< c
2QoT2 ~
Vv(x,y)
cosflmy COSflm~ sinhflmX sinhflm(a - ~)
(20)
2QoT e ~ sinflmy cosflm~ coshflmx sinhflm(a - ~) b~2d2 ,, ~ A(a C ~,, ~', -[~7,)
-
(21)
For c ~< x ~<
Vx (x,y)
-
~(x,y)
-
A'(x - C,C;flm) 2Q0 .... ~ 1 COSflmy COSflm~l sinhfl~(a - ~) A(a - C,C;flm) beld~ A(x c,c;flm) b~ld12Q°~:1 ~ sinflmy c o s f l ~ sinhfl~(a - ~) A(a - c,c;flm)
(22) (23)
F o r ~ <~ x ~< a
Q0
V~(x,y)
~ cosfl~y c o s f l ~ coshflm(a - x) b~ldl I1 + 2 m=l
-
(24)
A(( - c, c; tim)7 × A(a c, c; ~ - ) J" a(~ -
(x,y) -
c,c;~)
b~.~ 2Q°dl ~=1 sinflmy COSflm~ sinhflm(a - x) A(a - c,c;flm)
(25)
w h e r e ~k and d~ (k = 1,2) are the porosity and t h i c k n e s s of each zone and A'(x - c, c; tim) - T1 sinhfim(X - c) coshflmc + T2coshflm(x - c) sinhflmC In order to c h e c k the b o u n d a r y element model against the above a n a l y t i c solution, an aquifer of the type shown in Fig. 2 is assumed with a = 600m, b - 900m, and c = 300m. These small h o r i z o n t a l dimensions are chosen in order to m a k e the influence of the b o u n d a r i e s u p o n the solution substantial, t a k i n g into a c c o u n t the fact t h a t the m a g n i t u d e of the well's pumping r a t e is 432m a day 1. A wide r a n g e of transmissivity ratios, i.e. T~/Te, is used as the v a r y i n g i n p u t to both the a n a l y t i c and n u m e r i c a l models. F o r the B E M model two b o u n d a r y discretizations are examined: (a) one with 30 e x t e r n a l and 9 i n t e r z o n a l elements of equal length, n a m e l y 100m, and (b) a n o t h e r with the same 30 e x t e r n a l elements but a finer i n t e r z o n a l discretization with 18 elements, each 50 m long. Also, two possible single-well locations are assumed, one in the middle of Zone 1 (~ - 450 m, ~ - 450m) and a n o t h e r closer to the
169
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR AI)VECTIVE TRANSPORT
interface (~ = 350 m, ~ = 450 m). In all simulations both the piezometric head and the velocity field are calculated and checked. In general, the results of this verification procedure are very satisfactory, and are summarized in the following. (a) The solutions for the primary variable, the head, are more accurate than those of the velocity, the secondary variable. Analytical and numerical solutions for heads are almost identical, with a maximum relative error, in most cases, of less than 1% (the exception being the near-boundary areas as explained below). (b) The numerical solutions are sensitive to mesh discretization close to the domain and interzonal boundaries. This is a well-known characteristic and is intrinsic to the method: high accuracy is not retained within a boundary zone of a width approximately equal to half the length of the element (Brebbia, 1978). In fact, the errors calculated at points near the interface are almost halved in value when the fine r a t h e r than the coarser discretization is employed. As an example, the relative error for the velocities at points along a line located 25 m from the densely discretized interface is, on average, in the order of 4% while the same average is around 9% for points equally distant from the external (coarser) boundary elements. (c) For points located off those boundary zones the increase in accuracy, relative to the near-boundary ones, is significant. Especially in the vicinity of the pumping well, almost identical values for both head and velocities are obtained, owing to the natural representation of the well's behaviour in the numerical solution, as explained earlier. It should be noted that the same degree of accuracy is found for the well's locations both distant and close to the interface boundary. This is because the head gradients due to pumping only are very large in the vicinity of the well and therefore the influence of the boundary conditions is comparatively small.
,y'
y
!
1.o
i.o
.~
09
0.9
0.8
0.8
~O'}Ok~,?.6
"' ~"
0.7
,,
-,
06 , #
0.50.4
04
03
O.3 --
-
0.2
0.2 --
0.1
0.1 --
,
v~ I i I ' ] I I -1.0-0.8-0.6-0.4-q2
'
i
[
~ 0.0
Q2
0.4
0.6
0.8
1.0
v~ 0.0
i 0.0
I I I t I ~ I I i 0.2 0.4 0.6 0 . 8 1.0
Fig. 3. Dimensionless velocity components distribution along x analytic; (z~)numerical BEM solution.
275m for
T1/T2
: 104. (--)
170
P. LATINOPOULOSAND K. KATSIFARAKIS ,y"
y,
i
7.0
1.0
:~
0.9
0.9
O.8
0 8 •,
0.7
07 -
"~6
0.6--
~(~
0.5 0.4
04
02
0.2 --
o.1 I
i
I
q
l
'
I
i
-10-08-0.6-0.4-0.2
[
v~ ~
ik
0.0 0.2
0.4
06
0.8
7.0
o.~
v;
0.0
l
0.0
i
I
0.2 0.4
i
I
i
r
0.6 0.8
i
i B
7.0
Fig. 4. D i m e n s i o n l e s s velocity c o m p o n e n t s d i s t r i b u t i o n a l o n g x - 3 2 5 m for T 1/ T 2 = 10 -4. (--) a n a l y t i c ; (a) n u m e r i c a l B E M solution.
(d) The effect of transmissivity contrasts between the two zones upon the numerical model's accuracy was also investigated. The fact t hat no significant differences were found in the error distribution for extreme cases, e.g. those with T1/T2 = 10 4 and T1/T2 = 104, shows t hat the model is very accurate, at least for realistic situations. As an example of the whole verification procedure the plots in Figs. 3 and 4 are presented for the distant well location and fine interzonal discretization case. In both figures only the upper part of the symmetric velocity distribution is presented along the lines 25m from the interface. The dimensionless variables used in these figures read y' = (y - b/2)/b, V'~ = Vx/max Vx and V~ = VSmax Vx, where max Vx is the analytical value of Vx along the central line, y = b/2. The close agreement of the two solutions verifies the above comments regarding the model's accuracy. (e) The well-known phenomenon of refraction was also observed in the numerical results. In fact, the ratio of the velocity components in the ydirection for points located on either side of and very close to the interface, e.g. 20 cm, was found to be equal to the transmissivity ratio in all cases examined, as expected. (Of course, these values suffer from systematic errors due to the BEM itself but, nevertheless, they indicate the satisfaction of the compatibility criteria.) FORMULATION AND VERIFICATIONOF THE ADVECTIVE TRANSPORTMODEL In the present work the boundary element methodology of Latinopoulos (1986) is extended in that the simulation code is now the zoned aquifer boundary element model. In a similar way to the homogeneous case, the simulation model is coupled with a particle tracking program which is based on the numerical code RESSQ (Javandel et al., 1984). The basic assumption of RESSQ is that steady-state two-dimensional flow is considered in a horizontal homogeneous and isotropic aquifer of infinite (or semi-infinite) extent. Despite
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
171
this lack of high sophistication, RESSQ has been successfully used not only in pure simulation problems (e.g. Latinopoulos, 1986) but also in conjunction with non-linear optimization in groundwater quality management (Greenwald and Gorelick, 1989). By modifying the code of RESSQ in such a way that the calculation of the velocity field is performed through our BEM code, we obtain a handy, yet powerful tool for aquifer management because now the flow domain is bounded as well as heterogeneous. In brief, the application procedure of the combined model includes the following stages. (a) Solution of the boundary element system and calculation of the steadystate heads and head gradients for all boundary elements. (b) Generation of the streamlines by particle tracking. In this stage fluid particles are tracked along their streamlines by using local values of the velocity vector, as calculated using eqns. (14) and (15). For every new particle location a check is made in order to find out whether it is still in the same zone as the previous one. This can be achieved in the following way (Katsifarakis and Latinopoulos, 1990). The sum of all angles ~, (see Fig. l(b)) formed around the new location of the particle, in relation to the specific zone's boundary is compared with 27r. If it is less than 2:r then the particle has entered another zone, otherwise it is still in the same zone as the previous location. A similar procedure is followed to identify the new zone if the former condition is true. This calculation adds very little extra work, since the values of 7 are already known from the velocities evaluation at the same location (see eqns. (14) and (15)). The method's sole drawback is that, as presented, it applies to convex zones only. Alternative solutions for non-convex domains and a further discussion of those techniques can be found in the work of Katisfarakis and Latinopoulos (1990). (c) Particle travel times between any two points on a streamline are easily calculated in the process of streamline generation. (d) Contaminant fronts around injection wells and the time variations of concentrations at production wells or at outflow boundaries are also calculated simultaneously with the streamlines (Javandel et al., 1984). The accuracy of the combined model for zoned aquifers is checked against analytical solutions for infinite aquifers. This is because at this stage we wish to evaluate the performance of the model's output, which includes streamline generation and particle travel times, and it is known that it is very difficult to derive analytic expressions for the related variables if the aquifer is finite in extent. Thus, we consider a rather large aquifer, symmetric in space and with no natural flow (Fig. 5). This can be achieved by assuming equal head values (e.g. qS~, = 0) at all boundaries. By considering a pair of wells, an injection well in Zone 1 and a pumping well in Zone 2, located close to each other in the central part of the aquifer, as the sole stresses to the system, we have an approximate representation of flow within an infinite aquifer (at least in a large area around the wells). The numerical data for a series of example problems based on the above aquifer configuration are: T~ = 10 2m2s 1 ~ _ varying, d l = d~ (= d, the
e. LATINOPOULOSAND K. KATSIFARAKIS
172 l"
1500m
"1
ZONE 2
ZONE 1
/
3000m
150rn 150m
]-
3000 rn -
-
"t
Fig. 5. The idealized two-zone aquifer with a well doublet.
aquifer thickness) = 14m, ~1 = ~2 (= ~, its porosity) = 0.20, Q~n = 10m3h-1 and Qp = varying, where Qi, and Qp are the injection and production well flow rates, respectively. As far as the boundary element discretization is concerned, the specific needs of this study and especially the requirement for accurate description of the crossing of streamlines with the interzonal boundary lead us to apply a fine discretization. Thus we employed 40 external and 20 internal boundary elements with a proportionally increasing spacing between the interface nodes, the smallest two being 50 m long on either side of the intersection of the interface with the line connecting the two wells, and the largest two equal to 200 m near the constant head external boundaries. The first case examined is that of a homogeneous aquifer for which the two zones, still divided by the interface, have the same value of transmissivity (i.e. T1 = T2). We also assume that Qp equals Qin, the so-called 'well doublet' condition. With this data we run a semi-analytic model for infinite aquifers (Javandel et ah, 1984), the numerical model for bounded and homogeneous aquifers (Latinopoulos, 1986) and the present model for heterogeneous aquifers. The results from the three models are identical in the central part of the aquifer, and deviate near the boundaries, where displacement of curved streamlines is expected owing firstly to the existence of the boundaries themselves and secondly to the explicit scheme of streamline calculations. As in all three models above the time-marching scheme is essentially the same (Javandel et al., 1984); the numerical errors are equally affected by the value of the length step A1. Nevertheless, for a step length between I and 5 m all the
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
173
central streamlines showed a deviation of less than 5% from the corresponding analytic paths. As an example, the initial breakthrough to the production well, i.e. the shortest travel time between the two wells, is found to equal 3.03 years using all three models for Al = 5m. Application of the corresponding analytical relation (Bear, 1979)
4 nedx2o
tb
3
(26)
Qo
in which Q0 = Qin = Qp and x0 is half the distance between the two wells, leads to a value of tb = 3.01 years, representing a very small relative error of 0.6%. As this result refers to the shortest straight streamline, several arrival times for longer and apparently curved streamlines have also been calculated. The results from the two BEM models are again the same and differ from the semi-analytic result by up to 5% for travel times up to 20 years and up to 25% for travel times in the order of 100 years. Considering the fact that the streamlines correponding to the last case deviate from the semi-analytic ones mainly because of the boundary effects and should therefore not be taken into account, the final conclusion from this part of the study is that the combined BEMparticle tracking model is very accurate in generating streamlines. In a second numerical example using the same aquifer (Fig. 5), we assume that only the injection well in Zone 1 is operating, with Qin = 10m3h 1. (In reality, a negligible discharge rate of Qp = 3.6 × 10-4m3h 1 is retained in the production well in order to serve as a monitoring well in the numerical model). Leaving the rest of the data unchanged the particle-tracking model is applied for different transmissivity ratios, by changing the value of T2, and the travel time from the injection well to the monitoring well (equally distant from the interface) is calculated. The output of this application example is summarized in Fig. 6, where the numerical results are successfuly compared with the corresponding analytic expression
1000 5OO
IO0 5O
10 5
'/v 1
f 2
I 5
l 10
I 2O
I 5O
I 100
~
Fig. 6. V a r i a t i o n of d i m e n s i o n l e s s i n i t i a l b r e a k t h r o u g h w i t h t r a n s m i s s i v i t y r a t i o for t w o - z o n e 'infinite' aquifer. (--) analytic; (4) numerical BEM solution.
174
P. LATINOPOULOS AND K. KATSIFARAKIS
tbQ~,
-
~dx~
_
~(1 + ~) [ ~' - 1 + 3 + (72 - 1) In 2
?~7]
(27)
w h e r e r is the dimensionless initial b r e a k t h r o u g h and ~ = T 2 / T 1, also dimensionless. E q u a t i o n (27) has been derived by using the steady-state head expressions for the case of a single well n e a r an i n t e r f a c e line which separates two semi-infinite zones with different transmissivities (Bear, 1979). A l t h o u g h not clearly seen in Fig. 6 (mainly on a c c o u n t of the logarithmic scaling of the variables), t h e r e are some small differences in the results o b t a i n e d by the two approaches. Thus, the r e l a t i v e e r r o r of the n u m e r i c a l l y c a l c u l a t e d t r a v e l time for the h o m o g e n e o u s case ( T 1 / T 2 = 1) is a b o u t 1.5% and increases to a v a l u e of 6% for T I / T 2 = 100. These larger errors t h a t are o b t a i n e d w h e n the transmissivity r a t i o is i n c r e a s i n g are probably a t t r i b u t a b l e to the fact t h a t w h e n T2 ~ T1 the velocities in Zone 2 become v e r y small and thus more sensitive to r o u n d off errors. APPLICATION OF THE ADVECTIVE TRANSPORT MODEL As an example of application to a real case, the n u m e r i c a l model is applied to the complex aquifer of Fig. 7, which consists of four zones of different ~=0
o
50 too
2oom
~,=o
Fig. 7. Pathlines and pollution fronts in an aquifer with four zones of different transmissivity.
BOUNDARYELEMENTAND PARTICLETRACKINGMODELFOR ADVECTIVETRANSPORT
175
transmissivity. An i n j e c t i o n well operates at a point w i t h i n Zone 3. P o l l u t a n t particles, which e n t e r the aquifer t h r o u g h this well are t r a c k e d to obtain pathlines and p o l l u t i o n fronts. The n u m e r i c a l d a t a for this application are: T 1 = 0.001m 2s 1, T.~ = 0.01m 2s ', T~ = 0.005m 2s ', T4 = 0.0005m 2s ', d l = d2 = d3 - d4 20m, ~1 = ~2 = ~3 = ~4 = 0.20, Q = 0.01m3s 1 P a t h l i n e s a p p e a r in Fig. 7 as solid lines, while pollution fronts after 1, 2 and 3 years of well o p e r a t i o n a p p e a r as b r o k e n lines. It can be seen t h a t p o l l u t a n t particles move more slowly towards areas of lower transmissivity. It can also be seen t h a t pathlines are r e f r a c t e d as they cross interfaces between zones of different transmissivity and t h a t t h e y i n t e r s e c t p e r p e n d i c u l a r l y with the c o n s t a n t head boundaries. CONCLUSIONS
A b o u n d a r y element model for the c a l c u l a t i o n of heads and velocities in zoned aquifers has been developed and verified. A p a r t from its overall a c c u r a c y and simplicity, the model produces a c o n t i n u o u s velocity field, which is part i c u l a r l y i m p o r t a n t in g r o u n d w a t e r flows driven m a i n l y by well pumping and injection. The m a j o r f e a t u r e of this work is the successful coupling of the B E M model with a p a r t i c l e - t r a c k i n g technique. The combined model has p r o v e n to be an efficient and a c c u r a t e tool for advective t r a n s p o r t s i m u l a t i o n in zoned aquifers: a series of tests against a n a l y t i c a l solutions, especially n e a r i n t e r z o n a l boundaries, shows t h a t the n u m e r i c a l errors produced are not significant, are explicable, and are easily c o n t r o l l a b l e by simple means such as local grid r e f i n e m e n t and use of small time steps. Finally, a realistic application of the coupled model emphasizes its utility in solving real-life e n g i n e e r i n g problems. REFERENCES Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York, 569 pp. Brebbia, C.A., 1978. The Boundary Element Method for Engineers. Pentech, London, 255 pp. Brebbia, C.A. and Chang, O.V., 1979. Boundary elements applied to seepage problems in zoned anisotropic soils. Adv. Eng. Software, 1(3): 95-105. Chan, Y.K., Mullineux, N. and Reed, J.R., 1977. Analytic solution for drawdown in an unconfinedconfined rectangular aquifer. J. Hydrol., 34:287 296. Cheng, A. H-D., 1984. Darcy's flow with variable permeability: a boundary integral solution. Water Resour. Res., 20(7): 98(~984. Charbeneau, R.J. and Street, R.L., 1979a. Modeling groundwater flow fields containing point singularities: a technique for singularity removal. Water Resour. Res., 15(3): 583 594. Gharbeneau, R.J. and Street, R.L., 1979b. Modeling groundwater flow fields containing point singularities: streamlines, travel times, and breakthrough curves. Water Resour. Res., 15(6): 1445-1450. Chugh, A.K. and Falvey, H.T., 1984. Seepage analysis using the boundary element method. Report No. REC-ERC-83-11, Bureau of Reclamation, U.S. Department of the Interior, Engineering and Research Center, Denver, CO, 55 pp. Frind, E.O. and Matanga, G.B., 1985. The dual formulation of flow for contaminant transport modeling, 1, Review of theory and accuracy aspects. Water Resour. Res., 21(2): 159 169.
176
P. LATINOPOULOSAND K. KATSIFARAKIS
Galeati, G. and Gambolati, G., 1989. On the boundary conditions and point sources in the finite element integration of the transport equation. Water Resour. Res., 25(5): 847-856. Greenwald, R.M. and Gorelick, S.M., 1989. Particle travel times of contaminants incorporated into a planning model for groundwater plume capture. J. Hydrol., 107: 73-98. Ingham, D.B., Heggs, P.J. and Manzoor, M., 1981. The numerical solution of plane potential problems by improved boundary integral equation methods. J. Comput. Phys., 42: 77-98. Javandel, I., Doughty, C. and Tsang, C.F., 1984. Groundwater Transport: Handbook of Mathematical Models. Am. Geophys. Union, Water Resour. Monogr. 10, Washington, DC, 228 pp. Katsifarakis, K. and Latinopoulos, P., 1990. Two simple numerical techniques for use in particle tracking models of zoned aquifers. Hydrosoft, 3(4): 161 163. Kemblowski, M., 1986. A boundary element random walk model of mass transport in groundwater. d. Hydrol., 85:305 318. Lafe, O.E. and Cheng, A. H-D., 1987. A perturbation boundary element code for steady state groundwater flow in heterogeneous aquifers. Water Resour. Res., 23(6): 1079-1084. Lafe, O.E., Liggett, J.A. and Liu, P. L-F., 1981. BIEM solutions to combinations of leaky, layered, confined, unconfined, nonisotropic aquifers. Water Resour. Res., 17(5): 1431 1444. Latinopoulos, P., 1986. A boundary element approach for modeling groundwater movement. Adv. Water Resour., 9:171 177. Latinopoulos, P., Ganoulis, J. and Tolikas, D., 1982. Efficiency of integral equation method for modelling regional groundwater flows. Proc. Int. Conf. on Modern Approaches to Groundwater Resources Management, Capri, Italy, October 1982, Instituto di Idraulica e Costruzioni Idrauliche, Politecnico di Milano, pp. 207 215. Lerner, D.N., 1989. Predicting pumping water levels in single and multiple wells using regional groundwater models. J. Hydrol., 105: 39-55. Liggett, J.A. and Liu, P. L-F., 1983. The Boundary Integral Equation Method for Porous Media Flow. Allen and Unwin, London, 255 pp. Liggett, J.A. and Medina, D.E., 1987. Flow in three-dimensional fracture networks using a discrete approach. NATO Advanced Research Workshop on Advances in Analytical and Numerical Groundwater Flow and Quality Modeling, Lisbon, Portugal. Shapiro, A.M. and Andersson, J., 1983. Steady state fluid response in fractured rock: a boundary element solution for a coupled discrete fracture continuum model. Water Resour. Res., 19(4): 959-969. Sophocleous, M.A., 1984. Groundwater-flow parameter estimation and quality modeling of the Equus Beds aquifers in Kansas, U.S.A.J. Hydrol., 69:197 222. Yeh, G.T., 1981. On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow. Water Resour. Res., 17(5): 1529-1534.
159
Elsevier Science Publishers B.V., Amsterdam
[2]
A B O U N D A R Y E L E M E N T A N D PARTICLE TRACKING MODEL FOR ADVECTIVE T R A N S P O R T IN ZONED AQUIFERS
P, LATINOPOULOS and K. KATSIFARAKIS
Department of Civil Engineering, School of Technology, Aristotle University of Thessaloniki, Thessaloniki (Greece) (Received J a n u a r y 4, 1990; accepted after revision July 29, 1990)
ABSTRACT Latinopoulos, P. and Katsifarakis, K., 1991. A boundary element and particle tracking model for advective transport in zoned aquifers. J. Hydrol., 124:159 176. One of the most crucial factors t h a t affects both the accuracy and efficiency of a groundwater transport model is the flow simulation method which produces the velocity field. This paper presents a boundary element technique t h a t incorporates the zoning approach to allow aquifer heterogeneity. First, the accuracy of the flow model is verified by comparing it with analytic solutions for well-pumping and injection problems. The boundary element flow model is then coupled with a particle tracking procedure for advection-based transport simulation. The combined model also calculates streamlines, travel times and b r e a k t h r o u g h curves and is therefore a powerful tool for simulation and management of advection-dominated contamination problems.
INTRODUCTION
Current experience in the field of water resources management shows that the contamination of groundwater reserves is a widespread problem. This fact accounts for the remarkable increase in the development and use of groundwater simulation models in recent years. The main concern of a groundwater modeller is to produce a reliable numerical code for the investigation of a usually complex system in order to understand the physical processes that govern its response to past, present and future stresses and to evaluate the efficiency of remedial actions involved in proposed restoration schemes. The case of solute transport modelling is much more difficult than modelling groundwater flow alone. Moreover, the flow problem is an element of the transport problem, which means that possible inconsistencies arising from the former would inevitably be retained, if not magnified, within the simulation of the latter. Apart from the scarcity of solute concentration data and the uncertainties with regard to both the pollution sources and the transport parameters, it is well-known that the numerical models of solute transport suffer from inaccuracies, due mainly to numerical dispersion. The last is a particularly acute problem in advection-dominated transport and a common
0022-1694/91/$03.50
cC'. 1991
Elsevier Science Publishers B.V.
160
P. LATINOPOULOS AND K. KATSIFARAKIS
practice to tackle it is to separate the advective and dispersive parts. Nevertheless, a crucial factor in these types of simulation is the accurate estimation of the water velocity, as it directly influences not only the dispersive flux, but also the dominating advective transport. In addition, there is a class of problems where advection is assumed to be the sole driving transport mechanism, and therefore dispersion is ignored in the simulation procedure. In most of the existing numerical codes the calculation of the velocity field is based on solution of the flow equation using either the finite difference method or the finite element method. Experience shows that direct use of these domain-type methods often results in inaccurate and unstable solutions of the transport problem (Galeati and Gambolati, 1989). To overcome the accuracy problem various techniques have been proposed, principally for the finite element method, such as solving the flow problem directly for velocities (Yeh, 1981) or using interpolation functions that yield continuous velocities. Within the same context Frind and Matanga (1985) proposed a different approach by using stream functions as well as potentials in a finite element 'dual formulation' model for steady flow in the vertical plane of the homogeneous, anisotropic aquifer. Another critical point is the accuracy of the numerical solutions near singularities, such as pumping or injection wells. This is very important in groundwater quality management problems since at least one water supply or injection well is always under operation. The strong dependence of the solution accuracy near the wells upon the domain mesh discretization has led developers of domain-type models to incorporate local adjustments and corrections by employing complex numerical algorithms either for the finite difference method (Lerner, 1989) or for the finite element method (Charbeneau and Street, 1979a,b). Finally, when new locations of wells are required the mesh has to be changed in order to retain the same degree of accuracy for all operating wells. A method that can overcome the above problems is the boundary element method. Firstly it produces a continuous velocity field in the interior of the flow domain because it is grid-independent. That holds even for piecewise homogeneous aquifers, provided that physical discontinuities in transmissivity are discretized through interzonal boundaries. Secondly, it is very accurate within the domain and especially so near wells, where the natural representation of the solution reproduces the corresponding analytic ones (Latinopoulos et al., 1982; Latinopoulos, 1986). Thirdly, it requires much less input data and computational time than the domain-type solutions and, particularly for solute transport simulation, a minimal amount of bookkeeping (Kemblowski, 1986). These advantages led us to develop a boundary element model that accurately and efficiently calculates the water-pore velocities; it is therefore very useful when advection is either the sole or the major driving mechanism of solute transport. The model is an extension of a previous one (Latinopoulos, 1986) mainly in that it handles aquifer heterogeneity by employing the zoning technique. One essential characteristic of the model is that all integrals
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECT1VE TRANSPORT
]61
included in head and velocity calculations are evaluated by analytic expressions, leading to very accurate results, as will be shown later. Another feature of the model is that it can be efficiently combined with a particle-tracking numerical code that considers advection-based transport. The ability of the combined algorithm to simulate solute transport in heterogeneous aquifers makes it a useful tool in many applications concerned with groundwater quality management. In principle there are two limitations to the presented boundary elementbased approach. The first concerns the steady-state assumption for the flow problem. A possible extension to handle unsteady flow would heavily tax the overall efficiency of the specific boundary element technique. Moreover, in the class of transport problems for which the method is proposed the travel times of contaminant transport are large (in the order of years), so that the steadystate assumption is generally acceptable. The second drawback of the model refers to the overall reliability of the zoning technique. With this approach, which is very popular in boundary element solutions for both porous (Brebbia and Chang, 1979; Chugh and Falvey, 1984; Lafe et al., 1981) and fractured aquifers (Liggett and Medina, 1987; Shapiro and Andersson, 1983) the heterogeneous flow domain is divided into a number of zones which are each assumed to be homogeneous. The reported limitation arises from the fact that, as the number of zones increases, the efficiency of the boundary element technique decreases. This was the motivation for other boundary element-based approaches to heterogeneity, such as the introduction of specific functional permeability distributions (Cheng, 1984) and the expansion of the solution into perturbation serie.s (Lafe and Cheng, 1987). In conclusion, taking into consideration the facts that that in real-world applications the permeability information is in most cases insufficient and that, in general, uncertainties with regard to pollution histories and parameters are the rule rather than the exception, a restriction to a small number of zones with constant parameters may sometimes be dictated even by the model's calibration procedure (Sophocleous, 1984). THEORETICAL BACKGROUND
We consider the case of a plane, horizontal steady-state groundwater flow in a confined heterogeneous aquifer. The governing equation is obtained by combining the equation of conservation of mass and Darcy's law (Bear, 1979) and is
~ f c~¢\
O(T~
+ q = 0
(1)
where 4) is the piezometric head, T is the aquifer's transmissivity (generally a function of space) and q is a source function. To formulate our model, we assume that the non-homogeneous aquifer can be divided into a number of zones of constant transmissivity. For the sake of simplicity we also assume a
162
P. LATINOPOULOS AND K. KATSIFARAKIS
generally isotropic material, although anisotropy can be easily treated by a geometric scaling procedure, provided t hat zonal boundaries can be equally scaled. Under these assumptions and considering that the source term represents only point singularities, such as wells, eqn. (1) becomes
The above equation is valid for each single zone characterized by its own constant transmissivity value T. In this equation Qw is the pumping rate of the wth well (which is negative if it is an injection well) and 5w, which equals 5(x - x~:, y - Yw), is the Dirac delta function; x~, and Yw are the location coordinates of the operating wells (w = 1, W). As stated above the problem of a zoned aquifer can be solved by satisfying eqn. (2) in each subregion (zone), thus decomposing it into as many subproblems as the total number of zones with different transmissivities (Bear, 1979). To obtain simultaneous solutions for ¢ in all zones both external and interzonal boundary conditions must be specified. The former need no further comment while the latter are obtained by considering: (a) A compatability equation which expresses the fact t hat the piezometric heads are the same when a point j at an interzonal boundary, dividing two adjacent zones k and k + 1, is approached from both sides, i.e. Cj.k =
Oj,k~l
(3)
(b) The continuity equation, expressing the equality of the seepage rates across the boundary, i.e.
T~ (c~¢/~n)j,~ -
T~_l(80/On)2,k÷l
(4)
where Th and Tk+ 1 are the constant transmissivities of the adjacent homogeneous zones k and k + 1, respectively, and n is the distance measured along the normal to the boundary. The pore-water velocities of the flow within each homogeneous and isotropic zone are given by expressions
E
~. c~x
(5)
V~, -
K (?0
(6)
where Vx and V~ are the components of the velocity vector in the directions of the cartesian coordinates x and y, respectively, K is the hydraulic conductivity of the specific porous material of each zone, and e is its porosity. By definition, the transmissivity of each zone is T which is equal to K d where d is the aquifer's thickness. In what follows all K, d and ~ are assumed to be constant within each zone (taken usually equal to their average values), but varying between different zones.
BOUNDARY ELEMENT AND PART1CLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
163
If the variation of the piezometric head is known, normally after the solution of eqn. (2), then the velocity field can be calculated by employing eqns. (5) and (6). The velocity field can next be used to trace streamlines and to calculate travel times between specific points along them. This is very useful in simulating pure advective transport of a non-dispersive and non-reactive solute, but is also required in the general case where the two mechanisms of solute transport, i.e. advective transport and dispersion, are both important. In such problems one must solve the flow equation to produce the water velocities which are needed to calculate both the advective terms and the dispersion coefficients (Bear, 1979). INTEGRAL FORMULATIONOF THE FLOWMODEL In order to solve the governing partial differential equation (eqn. (2)) the boundary element method (BEM) is used. After incorporating the boundary conditions for the external boundaries in any combination of the three possible types, namely Dirichlet, Neumann and Cauchy, the integral formulation of the problem can be obtained by applying the standard BEM approach (Brebbia, 1978; Liggett and Liu, 1983). Although applications of the BEM with linear elements in zoned problems are reported in the literature (Liggett and Medina, 1987), in the present model constant-boundary elements, i.e. node-centred elements, are used, mainly for three reasons: (a) it is easier to develop and simpler to employ the analytical integration formulae for the heads, but especially for the flow velocities, using constant elements rather than higherapproximation ones, (b) specific terms of these formulae are incorporated in the particle-tracking pkocedure, presented later in this paper, and (c) the accuracy obtained even with this simple discretization has proven to be very satisfactory in previous applications (Latinopoulos, 1986) as well as in the present case. In brief, the procedure followed in solving the groundwater flow problem in zoned aquifers is based on the simple fact that each zone is treated separately as a single homogeneous unit in which the BEM is applied. Thus, the boundary of each zone is discretized in linear segments along which the values of the piezometric head ¢ and its derivative 3~/~n, out and normal to the boundary, are assumed to be constant and equal to the respective values Cj and O~bj/?nat the midpoint j of each segment. The boundary conditions for the external boundary of each zone are either prescribed flux or prescribed head or a combination of the two (for Cauchy-type boundary conditions), whereas for the interzonal boundary elements neither the head nor its gradient is prescribed. The final solution for the whole domain can be obtained by a proper coupling of the different zones. Performing the standard procedure of the boundary element method the following set of equations is constructed for each individual zone (Latinopoulos, 1986).
[h] {b}
= [D] {3~/?n} + {P].
(7)
164
P. LATINOPOULOS AND K. KATSIFARAK[S
\ j+l
/'
/~
5
(a)
Cb)
Fig. 1. Notation for the boundary element formulation.
where
A,j
-
! ~iJ'fi ~dFj
z: Dij
=
f ln(~)dYj
ifi
¢ j
if i
j
(8)
(9)
In the above three equations (eqns. (8)-(10)) i is the node of the ith element, j is a point moving along the j t h element, and w is a well location. The corresponding distances rij and riw are shown in Fig. l(a). The system of eqns. (7) is next rearranged in such a way that all the unknowns are on the left-hand side. The final set of equations for the whole domain is constructed by assembling all similar sets of BEM equations of the existing zones into a single global set. In addition, the individual zones are coupled within this global set of equations by incorporating the compatibility equations (3) and (4) for the interzonal boundary elements. The set of equations so assembled is finally solved by Gaussian elimination. It should be noted here that, although more efficient structures of the linear system are reported in the literature, we choose the above simple method because, as will be shown later, our problem requires the solution of only one system, while the time-consuming part in the simulations is the particle-tracking module.Therefore the relevant reduction in efficiency of the full model is a minor issue. When, through the solution of the system of equations, the values of ~bj and c~¢j/~n are known for all the external and interface boundary elements, the piezometric head distribution inside any zone can be calculated. Assuming that
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
165
point i is now within a specific zone which is bounded by N elements (where N stands for the total number of external and interface boundary elements of that zone only), its piezometric head is calculated in exactly the same way as in the case of a single homogeneous aquifer, i.e. through solution of the relation 2uq~i
= ~ __~ in ( 1 )
+ 2N ~b2!rij'n dVj
w: !
j=l
drj
3 In r~
(11)
In the above eqn. (11) the index k characterizing the specific zone for which it applies is omitted for the sake of simplicity (e.g. W and N should read Wk, Nh and so on). The line integrals appearing in eqns. (8), (9) and (11) have been analytically evaluated in a previous work (Ingham et al., 1981) and they have also been applied to various groundwater modelling studies (Latinopoulos et al., 1982; Latinopoulos, 1986). Using the relations given in eqns. (8) and (9) and according to the notation shown in Fig. l(b) these integrals are A°
=
)'
l if i ¢ j
Di~ = 13cos~(1 - ln/a) + /2cosfl(1 - lnl2) - /3~'sin~
)
(12)
and
Ao = 7c Dij
= Ii 1 - In
~))I if i
= j
(13)
The next step in the development of the model is the calculation of the velocity field within the different zones of the aquifer. This can be achieved by evaluating the gradient of the piezometric head from eqn. (11) and substituting it into eqns. (5) and (6). Using the analytical expressions of eqn. (12) (holding for internal points only) the following analytical formulae for the velocity field within each particular zone are obtained (Latinopoulos, 1986).
g S ~ Qwcw Vx = 2~n[ .... L1 Tr~
N ~ 1 Cj
2 1 ~-n I_/asin~ + ~
ll(12c2cosot + laclcosfl ) (/2/a) 2 sin7 ';c°t~ + In
(14)
166
P. LATINOPOULOS AND K. KATSIFARAKIS
K
Vv "
f ~
Qwd w
N
2nn [ w:l z_, Tr~
I1(I,fl2COS ~ + 1 3 d l c o s f l )
~ 1 ~j j=
' (12l:J ~ sin7
- j l~[_13sin~ + ~ -
:,cot~ * In
(15)
where cl
=
xj
c2 cu,
=
xi
dl
-
Yj-
(16a)
Yi
xj÷l -
xi
d2
-
Y~I
Yl
x~: -
xi
d~.
-
Yu - Yi
(16b) (16c)
As mentioned in that work (Latinopoulos, 1986), the main advantages of the straightforward procedure presented in calculating heads and velocities are high accuracy and computational efficiency. Especially in the presence of singularities, such as pumping and injection wells, their natural representation by the logarithmic term in the head relation and the derived gradient expressions in the velocity ones, is an exact reproduction of the corresponding analytical solutions for infinite confined aquifers. This can be easily seen by inspecting the well terms in eqns. (11), (14) and (15) which are the exact analtyical solutions for homogeneous aquifers of infinite extent, whereas the remaining terms approximate the boundary effects. This specific property of the boundary element solutions is a major advantage over the domain-type ones, which, as mentioned in the introduction, require a lot of extra development and computational effort to simulate near-well behaviour accurately. VERIFICATION OF THE NUMERICALFLOWMODEL In this part the boundary element model is tested for accuracy against an analytic solution. Having already verified the model for the resulting head (Latinopoulos et al., 1982) and for velocity values (Latinopoulos, 1986) in homogeneous aquifers, we are now interested in its performance in a zoned aquifer. Looking in the available literature for analytic solutions, the most suitable that we found is a generalized solution for drawdown due to pumping in a rectangular aquifer having two zones with distinct properties (Chan et al., 1977). This solution requires some mathematical manipulations to be tailored to our case and we now briefly describe it. We consider two-dimensional steady flow in a rectangular aquifer consisting of two rectangular zones of constant but different transmissivities T1 and T2, as shown in Fig. 2 Three sides of this aquifer are assumed to be impermeable while the fourth is kept at a constant head q~. A single well, located in Zone 1 at a point with coordinates ~ and ~ in a cartesian system, also shown in Fig. 2, is pumping continuously at a rate Q0. Thus, the steady-state analytic solution of this flow problem can be obtained by considering eqn. (2) as the governing equation. From the general, transient solution of Chan et al. (1977) we obtain the following analytic expressions for the steady-state head distribution.
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
167
,y Ii/////11/I/////I//////////
ZONE
2
ZONE1
(~,~) b
= const
c),p/8x-o1
i
~L
X 0",1/////I
/ / / t z ///////////////////
O ~--c
~'
o~O/Oy=o 4
Fig. 2. The idealized two-zone aquifer with a single pumping well.
For 0
~< x
~< c
¢(x,s)
= ¢¢-~
Oo[
a-C
+ 2T1 ~ cosflmy COSflm~coshflmx sinhflm(a - ~)1 For c
~< x
¢(x.y)
=
~<
4o--~1
For ~
~< x
c~(x,s)
=
+
Oo[
a -
cosflmy c o s f l m , sinhfl.~(a - ~) A ( x - c.c;flm)]
+2
(18)
~< a
~o -
2
(17)
Oo[
~
a -
x
~-. eosflmy Cosfl.,~ sinhflm(a - x) A(~ - c, c;flm),~=~L
fl,.A(a - c. c; tim)
(19)
168
P. LATINOPOULOS AND K. KATSIFARAKIS
where tim = m~/b and A ( X , Y ; f l m ) = TlcoshflmX c o s h f l ~ Y + T2 s i n h f l m X sinh tim Y B e c a u s e we are primarily i n t e r e s t e d in c h e c k i n g the a c c u r a c y of the velocity field r e s u l t i n g from our b o u n d a r y element model, in addition to the above, we need to o b t a i n the velocity c o m p o n e n t s in the x and y directions. T a k i n g the g r a d i e n t of the head from eqns. (17) (19) and s u b s t i t u t i n g into eqns. (5) and (6) we o b t a i n the following expressions. For 0 ~< x ~< c
2QoT2 ~
Vv(x,y)
cosflmy COSflm~ sinhflmX sinhflm(a - ~)
(20)
2QoT e ~ sinflmy cosflm~ coshflmx sinhflm(a - ~) b~2d2 ,, ~ A(a C ~,, ~', -[~7,)
-
(21)
For c ~< x ~<
Vx (x,y)
-
~(x,y)
-
A'(x - C,C;flm) 2Q0 .... ~ 1 COSflmy COSflm~l sinhfl~(a - ~) A(a - C,C;flm) beld~ A(x c,c;flm) b~ld12Q°~:1 ~ sinflmy c o s f l ~ sinhfl~(a - ~) A(a - c,c;flm)
(22) (23)
F o r ~ <~ x ~< a
Q0
V~(x,y)
~ cosfl~y c o s f l ~ coshflm(a - x) b~ldl I1 + 2 m=l
-
(24)
A(( - c, c; tim)7 × A(a c, c; ~ - ) J" a(~ -
(x,y) -
c,c;~)
b~.~ 2Q°dl ~=1 sinflmy COSflm~ sinhflm(a - x) A(a - c,c;flm)
(25)
w h e r e ~k and d~ (k = 1,2) are the porosity and t h i c k n e s s of each zone and A'(x - c, c; tim) - T1 sinhfim(X - c) coshflmc + T2coshflm(x - c) sinhflmC In order to c h e c k the b o u n d a r y element model against the above a n a l y t i c solution, an aquifer of the type shown in Fig. 2 is assumed with a = 600m, b - 900m, and c = 300m. These small h o r i z o n t a l dimensions are chosen in order to m a k e the influence of the b o u n d a r i e s u p o n the solution substantial, t a k i n g into a c c o u n t the fact t h a t the m a g n i t u d e of the well's pumping r a t e is 432m a day 1. A wide r a n g e of transmissivity ratios, i.e. T~/Te, is used as the v a r y i n g i n p u t to both the a n a l y t i c and n u m e r i c a l models. F o r the B E M model two b o u n d a r y discretizations are examined: (a) one with 30 e x t e r n a l and 9 i n t e r z o n a l elements of equal length, n a m e l y 100m, and (b) a n o t h e r with the same 30 e x t e r n a l elements but a finer i n t e r z o n a l discretization with 18 elements, each 50 m long. Also, two possible single-well locations are assumed, one in the middle of Zone 1 (~ - 450 m, ~ - 450m) and a n o t h e r closer to the
169
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR AI)VECTIVE TRANSPORT
interface (~ = 350 m, ~ = 450 m). In all simulations both the piezometric head and the velocity field are calculated and checked. In general, the results of this verification procedure are very satisfactory, and are summarized in the following. (a) The solutions for the primary variable, the head, are more accurate than those of the velocity, the secondary variable. Analytical and numerical solutions for heads are almost identical, with a maximum relative error, in most cases, of less than 1% (the exception being the near-boundary areas as explained below). (b) The numerical solutions are sensitive to mesh discretization close to the domain and interzonal boundaries. This is a well-known characteristic and is intrinsic to the method: high accuracy is not retained within a boundary zone of a width approximately equal to half the length of the element (Brebbia, 1978). In fact, the errors calculated at points near the interface are almost halved in value when the fine r a t h e r than the coarser discretization is employed. As an example, the relative error for the velocities at points along a line located 25 m from the densely discretized interface is, on average, in the order of 4% while the same average is around 9% for points equally distant from the external (coarser) boundary elements. (c) For points located off those boundary zones the increase in accuracy, relative to the near-boundary ones, is significant. Especially in the vicinity of the pumping well, almost identical values for both head and velocities are obtained, owing to the natural representation of the well's behaviour in the numerical solution, as explained earlier. It should be noted that the same degree of accuracy is found for the well's locations both distant and close to the interface boundary. This is because the head gradients due to pumping only are very large in the vicinity of the well and therefore the influence of the boundary conditions is comparatively small.
,y'
y
!
1.o
i.o
.~
09
0.9
0.8
0.8
~O'}Ok~,?.6
"' ~"
0.7
,,
-,
06 , #
0.50.4
04
03
O.3 --
-
0.2
0.2 --
0.1
0.1 --
,
v~ I i I ' ] I I -1.0-0.8-0.6-0.4-q2
'
i
[
~ 0.0
Q2
0.4
0.6
0.8
1.0
v~ 0.0
i 0.0
I I I t I ~ I I i 0.2 0.4 0.6 0 . 8 1.0
Fig. 3. Dimensionless velocity components distribution along x analytic; (z~)numerical BEM solution.
275m for
T1/T2
: 104. (--)
170
P. LATINOPOULOSAND K. KATSIFARAKIS ,y"
y,
i
7.0
1.0
:~
0.9
0.9
O.8
0 8 •,
0.7
07 -
"~6
0.6--
~(~
0.5 0.4
04
02
0.2 --
o.1 I
i
I
q
l
'
I
i
-10-08-0.6-0.4-0.2
[
v~ ~
ik
0.0 0.2
0.4
06
0.8
7.0
o.~
v;
0.0
l
0.0
i
I
0.2 0.4
i
I
i
r
0.6 0.8
i
i B
7.0
Fig. 4. D i m e n s i o n l e s s velocity c o m p o n e n t s d i s t r i b u t i o n a l o n g x - 3 2 5 m for T 1/ T 2 = 10 -4. (--) a n a l y t i c ; (a) n u m e r i c a l B E M solution.
(d) The effect of transmissivity contrasts between the two zones upon the numerical model's accuracy was also investigated. The fact t hat no significant differences were found in the error distribution for extreme cases, e.g. those with T1/T2 = 10 4 and T1/T2 = 104, shows t hat the model is very accurate, at least for realistic situations. As an example of the whole verification procedure the plots in Figs. 3 and 4 are presented for the distant well location and fine interzonal discretization case. In both figures only the upper part of the symmetric velocity distribution is presented along the lines 25m from the interface. The dimensionless variables used in these figures read y' = (y - b/2)/b, V'~ = Vx/max Vx and V~ = VSmax Vx, where max Vx is the analytical value of Vx along the central line, y = b/2. The close agreement of the two solutions verifies the above comments regarding the model's accuracy. (e) The well-known phenomenon of refraction was also observed in the numerical results. In fact, the ratio of the velocity components in the ydirection for points located on either side of and very close to the interface, e.g. 20 cm, was found to be equal to the transmissivity ratio in all cases examined, as expected. (Of course, these values suffer from systematic errors due to the BEM itself but, nevertheless, they indicate the satisfaction of the compatibility criteria.) FORMULATION AND VERIFICATIONOF THE ADVECTIVE TRANSPORTMODEL In the present work the boundary element methodology of Latinopoulos (1986) is extended in that the simulation code is now the zoned aquifer boundary element model. In a similar way to the homogeneous case, the simulation model is coupled with a particle tracking program which is based on the numerical code RESSQ (Javandel et al., 1984). The basic assumption of RESSQ is that steady-state two-dimensional flow is considered in a horizontal homogeneous and isotropic aquifer of infinite (or semi-infinite) extent. Despite
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
171
this lack of high sophistication, RESSQ has been successfully used not only in pure simulation problems (e.g. Latinopoulos, 1986) but also in conjunction with non-linear optimization in groundwater quality management (Greenwald and Gorelick, 1989). By modifying the code of RESSQ in such a way that the calculation of the velocity field is performed through our BEM code, we obtain a handy, yet powerful tool for aquifer management because now the flow domain is bounded as well as heterogeneous. In brief, the application procedure of the combined model includes the following stages. (a) Solution of the boundary element system and calculation of the steadystate heads and head gradients for all boundary elements. (b) Generation of the streamlines by particle tracking. In this stage fluid particles are tracked along their streamlines by using local values of the velocity vector, as calculated using eqns. (14) and (15). For every new particle location a check is made in order to find out whether it is still in the same zone as the previous one. This can be achieved in the following way (Katsifarakis and Latinopoulos, 1990). The sum of all angles ~, (see Fig. l(b)) formed around the new location of the particle, in relation to the specific zone's boundary is compared with 27r. If it is less than 2:r then the particle has entered another zone, otherwise it is still in the same zone as the previous location. A similar procedure is followed to identify the new zone if the former condition is true. This calculation adds very little extra work, since the values of 7 are already known from the velocities evaluation at the same location (see eqns. (14) and (15)). The method's sole drawback is that, as presented, it applies to convex zones only. Alternative solutions for non-convex domains and a further discussion of those techniques can be found in the work of Katisfarakis and Latinopoulos (1990). (c) Particle travel times between any two points on a streamline are easily calculated in the process of streamline generation. (d) Contaminant fronts around injection wells and the time variations of concentrations at production wells or at outflow boundaries are also calculated simultaneously with the streamlines (Javandel et al., 1984). The accuracy of the combined model for zoned aquifers is checked against analytical solutions for infinite aquifers. This is because at this stage we wish to evaluate the performance of the model's output, which includes streamline generation and particle travel times, and it is known that it is very difficult to derive analytic expressions for the related variables if the aquifer is finite in extent. Thus, we consider a rather large aquifer, symmetric in space and with no natural flow (Fig. 5). This can be achieved by assuming equal head values (e.g. qS~, = 0) at all boundaries. By considering a pair of wells, an injection well in Zone 1 and a pumping well in Zone 2, located close to each other in the central part of the aquifer, as the sole stresses to the system, we have an approximate representation of flow within an infinite aquifer (at least in a large area around the wells). The numerical data for a series of example problems based on the above aquifer configuration are: T~ = 10 2m2s 1 ~ _ varying, d l = d~ (= d, the
e. LATINOPOULOSAND K. KATSIFARAKIS
172 l"
1500m
"1
ZONE 2
ZONE 1
/
3000m
150rn 150m
]-
3000 rn -
-
"t
Fig. 5. The idealized two-zone aquifer with a well doublet.
aquifer thickness) = 14m, ~1 = ~2 (= ~, its porosity) = 0.20, Q~n = 10m3h-1 and Qp = varying, where Qi, and Qp are the injection and production well flow rates, respectively. As far as the boundary element discretization is concerned, the specific needs of this study and especially the requirement for accurate description of the crossing of streamlines with the interzonal boundary lead us to apply a fine discretization. Thus we employed 40 external and 20 internal boundary elements with a proportionally increasing spacing between the interface nodes, the smallest two being 50 m long on either side of the intersection of the interface with the line connecting the two wells, and the largest two equal to 200 m near the constant head external boundaries. The first case examined is that of a homogeneous aquifer for which the two zones, still divided by the interface, have the same value of transmissivity (i.e. T1 = T2). We also assume that Qp equals Qin, the so-called 'well doublet' condition. With this data we run a semi-analytic model for infinite aquifers (Javandel et ah, 1984), the numerical model for bounded and homogeneous aquifers (Latinopoulos, 1986) and the present model for heterogeneous aquifers. The results from the three models are identical in the central part of the aquifer, and deviate near the boundaries, where displacement of curved streamlines is expected owing firstly to the existence of the boundaries themselves and secondly to the explicit scheme of streamline calculations. As in all three models above the time-marching scheme is essentially the same (Javandel et al., 1984); the numerical errors are equally affected by the value of the length step A1. Nevertheless, for a step length between I and 5 m all the
BOUNDARY ELEMENT AND PARTICLE TRACKING MODEL FOR ADVECTIVE TRANSPORT
173
central streamlines showed a deviation of less than 5% from the corresponding analytic paths. As an example, the initial breakthrough to the production well, i.e. the shortest travel time between the two wells, is found to equal 3.03 years using all three models for Al = 5m. Application of the corresponding analytical relation (Bear, 1979)
4 nedx2o
tb
3
(26)
Qo
in which Q0 = Qin = Qp and x0 is half the distance between the two wells, leads to a value of tb = 3.01 years, representing a very small relative error of 0.6%. As this result refers to the shortest straight streamline, several arrival times for longer and apparently curved streamlines have also been calculated. The results from the two BEM models are again the same and differ from the semi-analytic result by up to 5% for travel times up to 20 years and up to 25% for travel times in the order of 100 years. Considering the fact that the streamlines correponding to the last case deviate from the semi-analytic ones mainly because of the boundary effects and should therefore not be taken into account, the final conclusion from this part of the study is that the combined BEMparticle tracking model is very accurate in generating streamlines. In a second numerical example using the same aquifer (Fig. 5), we assume that only the injection well in Zone 1 is operating, with Qin = 10m3h 1. (In reality, a negligible discharge rate of Qp = 3.6 × 10-4m3h 1 is retained in the production well in order to serve as a monitoring well in the numerical model). Leaving the rest of the data unchanged the particle-tracking model is applied for different transmissivity ratios, by changing the value of T2, and the travel time from the injection well to the monitoring well (equally distant from the interface) is calculated. The output of this application example is summarized in Fig. 6, where the numerical results are successfuly compared with the corresponding analytic expression
1000 5OO
IO0 5O
10 5
'/v 1
f 2
I 5
l 10
I 2O
I 5O
I 100
~
Fig. 6. V a r i a t i o n of d i m e n s i o n l e s s i n i t i a l b r e a k t h r o u g h w i t h t r a n s m i s s i v i t y r a t i o for t w o - z o n e 'infinite' aquifer. (--) analytic; (4) numerical BEM solution.
174
P. LATINOPOULOS AND K. KATSIFARAKIS
tbQ~,
-
~dx~
_
~(1 + ~) [ ~' - 1 + 3 + (72 - 1) In 2
?~7]
(27)
w h e r e r is the dimensionless initial b r e a k t h r o u g h and ~ = T 2 / T 1, also dimensionless. E q u a t i o n (27) has been derived by using the steady-state head expressions for the case of a single well n e a r an i n t e r f a c e line which separates two semi-infinite zones with different transmissivities (Bear, 1979). A l t h o u g h not clearly seen in Fig. 6 (mainly on a c c o u n t of the logarithmic scaling of the variables), t h e r e are some small differences in the results o b t a i n e d by the two approaches. Thus, the r e l a t i v e e r r o r of the n u m e r i c a l l y c a l c u l a t e d t r a v e l time for the h o m o g e n e o u s case ( T 1 / T 2 = 1) is a b o u t 1.5% and increases to a v a l u e of 6% for T I / T 2 = 100. These larger errors t h a t are o b t a i n e d w h e n the transmissivity r a t i o is i n c r e a s i n g are probably a t t r i b u t a b l e to the fact t h a t w h e n T2 ~ T1 the velocities in Zone 2 become v e r y small and thus more sensitive to r o u n d off errors. APPLICATION OF THE ADVECTIVE TRANSPORT MODEL As an example of application to a real case, the n u m e r i c a l model is applied to the complex aquifer of Fig. 7, which consists of four zones of different ~=0
o
50 too
2oom
~,=o
Fig. 7. Pathlines and pollution fronts in an aquifer with four zones of different transmissivity.
BOUNDARYELEMENTAND PARTICLETRACKINGMODELFOR ADVECTIVETRANSPORT
175
transmissivity. An i n j e c t i o n well operates at a point w i t h i n Zone 3. P o l l u t a n t particles, which e n t e r the aquifer t h r o u g h this well are t r a c k e d to obtain pathlines and p o l l u t i o n fronts. The n u m e r i c a l d a t a for this application are: T 1 = 0.001m 2s 1, T.~ = 0.01m 2s ', T~ = 0.005m 2s ', T4 = 0.0005m 2s ', d l = d2 = d3 - d4 20m, ~1 = ~2 = ~3 = ~4 = 0.20, Q = 0.01m3s 1 P a t h l i n e s a p p e a r in Fig. 7 as solid lines, while pollution fronts after 1, 2 and 3 years of well o p e r a t i o n a p p e a r as b r o k e n lines. It can be seen t h a t p o l l u t a n t particles move more slowly towards areas of lower transmissivity. It can also be seen t h a t pathlines are r e f r a c t e d as they cross interfaces between zones of different transmissivity and t h a t t h e y i n t e r s e c t p e r p e n d i c u l a r l y with the c o n s t a n t head boundaries. CONCLUSIONS
A b o u n d a r y element model for the c a l c u l a t i o n of heads and velocities in zoned aquifers has been developed and verified. A p a r t from its overall a c c u r a c y and simplicity, the model produces a c o n t i n u o u s velocity field, which is part i c u l a r l y i m p o r t a n t in g r o u n d w a t e r flows driven m a i n l y by well pumping and injection. The m a j o r f e a t u r e of this work is the successful coupling of the B E M model with a p a r t i c l e - t r a c k i n g technique. The combined model has p r o v e n to be an efficient and a c c u r a t e tool for advective t r a n s p o r t s i m u l a t i o n in zoned aquifers: a series of tests against a n a l y t i c a l solutions, especially n e a r i n t e r z o n a l boundaries, shows t h a t the n u m e r i c a l errors produced are not significant, are explicable, and are easily c o n t r o l l a b l e by simple means such as local grid r e f i n e m e n t and use of small time steps. Finally, a realistic application of the coupled model emphasizes its utility in solving real-life e n g i n e e r i n g problems. REFERENCES Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York, 569 pp. Brebbia, C.A., 1978. The Boundary Element Method for Engineers. Pentech, London, 255 pp. Brebbia, C.A. and Chang, O.V., 1979. Boundary elements applied to seepage problems in zoned anisotropic soils. Adv. Eng. Software, 1(3): 95-105. Chan, Y.K., Mullineux, N. and Reed, J.R., 1977. Analytic solution for drawdown in an unconfinedconfined rectangular aquifer. J. Hydrol., 34:287 296. Cheng, A. H-D., 1984. Darcy's flow with variable permeability: a boundary integral solution. Water Resour. Res., 20(7): 98(~984. Charbeneau, R.J. and Street, R.L., 1979a. Modeling groundwater flow fields containing point singularities: a technique for singularity removal. Water Resour. Res., 15(3): 583 594. Gharbeneau, R.J. and Street, R.L., 1979b. Modeling groundwater flow fields containing point singularities: streamlines, travel times, and breakthrough curves. Water Resour. Res., 15(6): 1445-1450. Chugh, A.K. and Falvey, H.T., 1984. Seepage analysis using the boundary element method. Report No. REC-ERC-83-11, Bureau of Reclamation, U.S. Department of the Interior, Engineering and Research Center, Denver, CO, 55 pp. Frind, E.O. and Matanga, G.B., 1985. The dual formulation of flow for contaminant transport modeling, 1, Review of theory and accuracy aspects. Water Resour. Res., 21(2): 159 169.
176
P. LATINOPOULOSAND K. KATSIFARAKIS
Galeati, G. and Gambolati, G., 1989. On the boundary conditions and point sources in the finite element integration of the transport equation. Water Resour. Res., 25(5): 847-856. Greenwald, R.M. and Gorelick, S.M., 1989. Particle travel times of contaminants incorporated into a planning model for groundwater plume capture. J. Hydrol., 107: 73-98. Ingham, D.B., Heggs, P.J. and Manzoor, M., 1981. The numerical solution of plane potential problems by improved boundary integral equation methods. J. Comput. Phys., 42: 77-98. Javandel, I., Doughty, C. and Tsang, C.F., 1984. Groundwater Transport: Handbook of Mathematical Models. Am. Geophys. Union, Water Resour. Monogr. 10, Washington, DC, 228 pp. Katsifarakis, K. and Latinopoulos, P., 1990. Two simple numerical techniques for use in particle tracking models of zoned aquifers. Hydrosoft, 3(4): 161 163. Kemblowski, M., 1986. A boundary element random walk model of mass transport in groundwater. d. Hydrol., 85:305 318. Lafe, O.E. and Cheng, A. H-D., 1987. A perturbation boundary element code for steady state groundwater flow in heterogeneous aquifers. Water Resour. Res., 23(6): 1079-1084. Lafe, O.E., Liggett, J.A. and Liu, P. L-F., 1981. BIEM solutions to combinations of leaky, layered, confined, unconfined, nonisotropic aquifers. Water Resour. Res., 17(5): 1431 1444. Latinopoulos, P., 1986. A boundary element approach for modeling groundwater movement. Adv. Water Resour., 9:171 177. Latinopoulos, P., Ganoulis, J. and Tolikas, D., 1982. Efficiency of integral equation method for modelling regional groundwater flows. Proc. Int. Conf. on Modern Approaches to Groundwater Resources Management, Capri, Italy, October 1982, Instituto di Idraulica e Costruzioni Idrauliche, Politecnico di Milano, pp. 207 215. Lerner, D.N., 1989. Predicting pumping water levels in single and multiple wells using regional groundwater models. J. Hydrol., 105: 39-55. Liggett, J.A. and Liu, P. L-F., 1983. The Boundary Integral Equation Method for Porous Media Flow. Allen and Unwin, London, 255 pp. Liggett, J.A. and Medina, D.E., 1987. Flow in three-dimensional fracture networks using a discrete approach. NATO Advanced Research Workshop on Advances in Analytical and Numerical Groundwater Flow and Quality Modeling, Lisbon, Portugal. Shapiro, A.M. and Andersson, J., 1983. Steady state fluid response in fractured rock: a boundary element solution for a coupled discrete fracture continuum model. Water Resour. Res., 19(4): 959-969. Sophocleous, M.A., 1984. Groundwater-flow parameter estimation and quality modeling of the Equus Beds aquifers in Kansas, U.S.A.J. Hydrol., 69:197 222. Yeh, G.T., 1981. On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow. Water Resour. Res., 17(5): 1529-1534.