M I G R A T I O N OF P O L L U T A N T S

IN G R O U N D W A T E R .

II. A D S O R B A B L E P O L L U T A N T S A N D N U M E R I C A L DISPERSION REDUCTION KENNETH

N. C A R T E R ,

Jr,, MELIO

SAENZ,

and D A V I D 1. W I L S O N *

Departments of Chemistry and of Civil and Environmental Engineering, Vanderbilt University, Nashville, TN37235, U.S.A. and P H I L I P W, R O S T E N , Jr.

ERM-Southeast lnc,. Brentwood, TN 37027, U~S.A.

(Received August, 1983) Abstract. We discuss here the partial differential equations governing the migration of a decomposing pollutant adsorbing according to a Langmuir isotherm and undergoing 2-dimensional flow in a saturated aquifer. The equation governing the mass transfer of the pollutant to the surfaces within the aquifer are solved in.closed form, permitting the use of larger values of the time increment At in the numerical integration of the dispersion-advection equation governing the behavior of the dissolved pollutant. In this numerical integration transverse numerical dispersion is eliminated by using conformal coordinates (velocity potential and stream function), and longitudinal numerical dispersion is very substantially reduced by use of an asymmetrical 4-point formula to represent the advection term. Some representative results are given as contour maps. The mass transfer rate coefficient is estimated as the least positive eigenvalue of a diffusion problem.

I. Introduction

The discovery of a number of sites in the U.S. at which groundwater has been contaminated by toxic pollutants has spurred interest in the modeling of the dynamics of pollutants in aquifers. The optimal siting of monitoring wells around hazardous waste landfills requires that we be able to model the evolution of plumes of various solutes in the aquifers of interest. Generally these systems cannot be modeled analytically if we take dispersion into account. If we have a constant, uniform flow field, a non-adsorbing pollutant, and a diagonal dispersion tensor, the partial differential equation describing the movement of the solute can be solved analytically (Ogata, 1970). This provides a useful first approximation from which to get insight about the movement of a non-adsorbing pollutant in an aquifer. Usually, however, one resorts to numerical methods. Anderson's (1979) detailed and lucid review discusses the methods which have been employed and notes the liabilities and benefits of each. The direct attack on the partial differential equation (1) by finite difference methods was found tO lead to rather severe numerical dispersion, an artifact * To whom correspondence should be addressed.

EnvironmentalMonitoring andAssessmem 4 (1984) 171-202. 9 1984 by D. Reidel Publishing Company.

0167-6369/84/0042-0171504.80,

172

K. N. CARTER ET AL.

of the computational method which can be of the same magnitude as the physical dispersion (Oster et al., 1970; Bredehoeft, 1971). In a search for a solution ~r -

-

= v. (BVc) - V. (vc) -

~t

kc

(1)

c(r, t) = pollutant concentration, moles cm- 3 at r = ix + jy at time t, D=

dispersion tensor, cm 2 sec-1,

v=

velocity of liquid, cm sec- 1,

k=

decomposition rate constant, sec- ~,

to this difficulty, Gardner et al. (1964) used the method of characteristics to get approximate solutions to Equation (1); Anderson (1979) notes that, while this approach does eliminate numerical dispersion, it is cumbersome to program and execute. Finite element techniques have also been employed for solving this problem. Pinder (1974) and Grove (1977) have discussed this approach, finite difference methods, and the method of characteristics. In this paper we first present the equations describing the movement of a pollutant in an aquifer, the pollutant being capable of decomposition and of adsorption on the pore surfaces within the aquifer. We then rewrite the dispersion-advection equation in a suitably chosen pair of conformal coordinates, the velocity potential and the stream function which eliminates transverse numerical dispersion (Shamir and Harleman, 1967). We then note that the equation governing mass transfer of solute between the bulk solution and the surface phase adsorbed on the porous aquifer medium can be solved in closed form by making a very reasonable approximation. This replaces half of the differential equations of the system by algebraic equations and also permits the use of larger time increments in the numerical integration of the dispersion-advection equation than would otherwise be possible. We then examine a number of finite difference formulas for representing advection (only) in one dimension. Some of these have the disadvantage of leading to excessive longitudinal numerical dispersion. Some yield rather sharp plume fronts, but also yield objectionable oscillatory behavior on either side of the plume front. Two 4-point asymmetrical formulas are found to yield rather sharp plume fronts and to be remarkably free of oscillatory behavior on either side of the plume front. These formulas are then used to represent the advection term in the conformal representation of the dispersionadvection equation. The system is then solved for some representative situations and the results are displayed as line printer contour maps or conventional contour maps. We then turn to the problem of estimating the mass transfer rate coefficient for the movement of solute from the liquid flowing through the interstices of the aquifer medium across the boundary layer to the phase adsorbed at the surfaces of the medium. This is formulated as an eigenvalue problem, and the mass transfer rate coefficient is identified with the smallest non-zero eigenvalue of the system.

MIGRATION OF POLLUTANTS IN GROUNDWATER

173

2. The Model Equations We take as our starting point a slight generalization of Equation (1) which describes the dynamics of a solute undergoing first-order decomposition and mass-transfer-limited adsorption on the surfaces available within the porous aquifer, af = v. (DVc) - v. Vc - kc -

~ (~

3t

r/

- ~).

(2)

Here 0 since we are assuming that the velocity v is derived from a potential K=

mass transfer rate constant for solute adsorption, sec- l

r/~-

porosity of aquifer concentration of adsorbed solute, mole c m - 3 ~max/(1 + C1/2/Ce), adsorbed solute concentration in equilibrium with solution concentration G (3)

~m~ = Langmuir isotherm capacity parameter, mole c m - 3

cl/2 = Langrnuir isotherm concentration parameter, mole cm - 3, the solute concentration in water at which ~e -----!~ 2 max" We also require an equation describing the change of r

t) with time; we take

a~ - - = X ( ~ e -- ~) -- k ' s ~ '

(4)

at

Here k' = decomposition rate constant for adsorbed solute, sec- ~. In choosing the form of the dispersion tensor we follow Scheidegger (1961), 2

Do: Z m=l

2

Z %m. vmv" n=l

(5)

IY[

For an isotropic medium we have O~iiii = O~L

(6)

~i~ = aT 1

and all other %m, vanish. Here ~L = longitudinal dispersivity, cm c~r = transverse dispersivity, cm.

174

K. N. CARTER ET AL.

It is commonly assumed that eL is about three times larger than e r. The dispersion tensor calculated in cartesian coordinates by this prescription is

1 ( eLv2x + e r e O = ~

(eL-er)VxVY~

(7)

\(O~L -- eT) Vx Vy e T ,~2 nt" eL ]): ] '

In an earlier paper (Wilson, 1982) we noted the advantages of working in a coordinate system obtained by a conformal mapping appropriate to the velocity field. If our mapping is given by

w = u(x, y) + iv(x, y) = f(z),

(8)

where u(x, y) is the velocity potential (such that 7u = v) and v(x, y) = const, gives the streamlines, Equation (2) becomes

0t

L0u _

Ah_2

~u + ~vv .ar . kc . . K. (~e

Ou

(9)

~)"

n

Here

h= [(Oxl= + (aYl=l'/= = [(axl= + (o'l=l Lkau/

\a-u/_1

Lkav/

\~-v/J

(i0)

andA is a scale factor for the velocity. We note that Equation (9) is free of mixed second derivatives in the dispersion term, and that the advection term involves only Oc/Ou. The first substantially simplifies the representation of the dispersion by finite differences. The second eliminates transverse numerical dispersion from the advection term and also simplifies its representation by finite differences.

3. The Adsorption Equation We recall Equation (4), which, together with Equation (3), we take as controlling the adsorption and desorption of solute onto and from the surface of the aquifer medium. If we solve Equation (4) simultaneously with Equation (2) by numerical integration, we require rather small values of At in the numerical integration, which requires large amounts of computer time. An alternative approach is as follows. Let us assume that during the time interval At the change in total solute mass in the (i,j)th region (centered at (u;, vj), and of area h2AuAv) is small. We can then neglect the effects of advection, dispersion, and decomposition in examining the adsorption or desorption of solute. This gives us a conservation of solute mass equation (for the (i, j)th region), ~e + nee = ~ + ~ c=-

m

(11)

MIGRATION OF POLLUTANTS IN GROUNDWATER

175

where ce and r are the solution concentration and adsorbed concentration of solute which would pertain at equilibrium, and c and ~ are the initial concentrations. The quantity m is assumed to remain constant during the time interval At. Solving Equation (3) for ce yields C e

C1/2 ~e

--

(12)

This is substituted into Equation (11), to give, after rearrangement, ~e2 -- ~(m + ~max + r/Cl/2) + m~ma~ = 0.

(13)

The desired root of this equation is ~e = 89

+ ~max "b qCl/2) -- [(m + ~max + /~Cl/2)2 -- 4mCmax]l/2) 9

(14)

We are now "in position to solve Equation (4): 0~ + (x + k')~ = ~c~e Ot

(15)

integrates to give r

= r

(x +

k')t] +

~e K

~:+k'

{1 - e x p [ - ( x + k ' ) t ] ) .

(16)

The change in adsorbed solute concentration during the time interval At is given by ~(At) - ~(0), which results in a change of solute concentration in the liquid phase in this region of Ac = r/-'[~(0) - ~(At)].

(17)

So the contribution of adsorption to Equation (2) for the concentration of solute in the moving liquid in the (i,j)th region is

(Oco~

~o[pAt]

~t/=d~

-

~/j[(P + 1)At] r/At

=rl-l[1-exp(-KAt)][~o(pAt)-~e,o(pAt)

x--~].

(18)

Here we calculate Ce,,j from Equation (14), noting that in that equation we must set

m = ~o(pAt) + qco(pAt). 4. Representation of Advection It has been frequently noted (see Anderson, 1979) that the advection term in Equation (2) causes severe problems when it is represented by finite difference ex-

176

K . N . CARTER ET AL.

pressions; a particularly clear and quite entertaining exposition of the folklore associated with this has been presented by Leonard (1979b). Because of our transformation of the problem into the conformal coordinates u and v, we shall need to deal .only with advection in one dimension. The simplest and most physically self-evident representation for the one-dimensional advection term 9 ( v c ) / ~ x in a finite difference scheme is certainly -- [(~(VC)/t~X]n ~-- ( V n - 1 C n - I

(19)

- VnCn)/Ax

which conserves solute and does not exhibit oscillatory instabilities if v A t / A x < 1. If v is constant, as it was in all of our one-dimensional tests, it may be factored out. Equation (19) suffers an acute disadvantage, however, in that it leads to severe numerical dispersion, as illustrated in Figure la, in which we display an approximate solution to Oc v Oc/Ox. (v constant) (20) Ot -

In the absence of numerical dispersion the plume shown would have a perfectly sharp front, since the system generating the plume includes no physical dispersion terms. We find, instead, that this algorithm results in a progressively spreading front as the plume evolves with time, which makes it unsatisfactory for the representation of advection. It is also possible to represent the advection term by v ( c . _ i/2 - c~ + I / 2 ) / A x , and determine the c's by linear interpolation. This yields the same formula one gets as the simplest symmetric central difference formula for - v ~c/Sx,

-

v

---2A-xx ( c . _ , - c ~

(21)

n

If we use this formula for obtaining an approximate solution to Equation (20) we find that we do reduce the extent of numerical dispersion as compared to that resulting from Equation (19), but this algorithm is plagued with oscillatory behavior to an extent which makes it useless for our purposes, as seen in Figure lb. Such oscillatory behavior occurred with all of the central difference algorithms we explored. By a standard Taylor's series expansion method, as discussed by Collatz (1960), one can show that -

v

----n

(-~c._=+Sc._,

~c.+l+~c.+2).

(22)

AX

Approximate solution of Equation (20) by use of this algorithm yielded results of the type shown in Figure lc. We can also write - ( v Oc/dx). as - v ( c . _ 1/2 - c.+ 1/2)/kx and then calculate the c. +_1/e by Taylor's series expansion in terms of c._ 2, c. _ 1, c., c~ + 1, and c. + 2; the result is --

I~

~' AX

n

( -

8

n - 2 +

4

n -

1 -

4On

+1 +

8On +

2)"

(23)

177

MIGRATION OF POLLUTANTS IN GROUNDWATER

Cl

C

e

g h

0

o

03: "6

I

0

500

~

,h

I000 m

x~

Fig. 1.

(a)

Numerical solutions to the equation ((9c/8t) = - v(Oc/Ox) + 3(x); v = 1.0 m d a y - 1, Ax = l0 m, L = I000 m, At = 1.0 day, all solutions plotted at t = 800 days. Oc -

--

~

(c._,

-

c.)lax

;

dx

(b)

(c)

6~C

1

- - - ~ ~ ( c . _ ~ - c . + ~f l a x ; 8x

_ _ _~C ~(

I

(d)

2

2

- - ' 2 C n - 2 + 3 C n - 1 -- 3Cn+ I + ~ 2 C n + 2 ) / A X ;

~x

-__~(--1

c

~x

6~C

(e)

- - ~ (-~c._,

(f)

--8x

(g)

- - - - ~ ( - - 8 C n -'2

t~C

( - i I c, _ 2 + 2 c ~ _ 1

~C

8x

(h)

+ ~c._~, - ~ c . + ,

8x

-

&

--

Ox~

l

( - ~._

-3c,)/Ax;

+ Zc 8 n-I

~ +

r

+ ~c.+,)/ax;

- ~c. - ]c. + ,)lAx ;

-

~c.

-

~c. +,)1~.

178

K. N. CARTER ET A L

Solution of Equation (20) with Equation (23) used to approximate the advection term leads to oscillatory behavior of the type shown in Figure ld. Another variation is to write - (v Oc/Ox),, as - v(c,, _ ,/2 -- Cn+ 1/2)/Ax, where %_ ,/2 and c, + 1/2 are calculated by symmetrical four-point interpolation formulas. We find - •16

Cn-1/2"~

n-2

-t- i ~ C n _ l

(24)

"+" ~ r C n -- l ~ C n + l

and C n + ,/2 ~ ' - - l ~ r

- , + --~Cn + ~ C n + ,

(25)

-- l ~ C n + 2

yielding

--

"V

t3xc)

v

__lcn _ 2 -b 58 cn _

'~ ~X

( -- , 6

5

, -- ~ r

+ l "1" ~ 6 C n + 2 ) .

(26)

n

This yields highly oscillatory results, as shown in Figure le. We conclude on the basis of these results, and in agreement with Leonard (1979a, b) that symmetrical difference formulas for the advection term are not profitable, in that they lead to oscillatory instabilities. These results and Leonard's (1979a) remarks suggest that we explore asymmetrical difference representations for the advection term. The three-point formulas for - ( O c / O x ) . and for (c n _ 1/2 Cn+ 1/2)/AX in terms of Cn_ 2, Cn-1, a n d c~ are identical; the equation is -

-v

-

,,---Axx ( - 2 v-2 + 2c,,-1

9

(27)

The use of this algorithm to obtain an approximate solution to Equation (20) yields results like that shown in Figure If. This shows a front which spreads much more slowly than that obtained by using Equation (19) to represent the advective term, as seen by comparing Figures la and If. It is also seen to be nearly free of the oscillatory behavior associated with the central difference representations which do not contain %. This finding encourages us to further investigate asymmetrical difference representations for the advective term. We look at two, both of the form a _ 2 c , _ 2 + a _ 1 % - 1 + aoC,, + a l c n + ,. For the first we use Taylor's series expansions to show that

_v(

q

V,a x l .

v

~- A x

--

( c . _ ,/2 - c . + 1/2)

v (

Ax

!c

7

3

-- 8 n -- 2 "t- ~ C n _ 1 -- 8 C n --

3c,,+,).

(28)

We used this algorithm to obtain approximate solutions to Equation (20), including in our program an instruction to set equal to zero any small negative concentrations ahead of the advancing front. We find virtually no oscillations, but we do find a hump just

MIGRATION OF POLLUTANTS IN GROUNDWATER

179

behind the advancing front which should not be there and amounts to an error of about 8 ~o, as seen in Figure lg. We obtain our second asymmetrical four-point formula by using Taylor's series expansions to find -

v

\Ox/.

~

--

-

Ax

( - gc._

2 + c._

1 -

~c.

-

~c.

+ 1).

(29)

This algorithm was used to solve Equation (20) approximately, with the results shown in Figure lh. We find that the advancing front remains rather sharply defined; that oscillatory behavior behind the front is essentially absent; that the hump behind the advancing front amounts to about a 3 ~o error; and that the oscillations leading the front are of amplitude less than or equal to 3 ~o of the height of the front. (This is about what we find using Equation (28), also.) Since negative concentrations are physically impossible, we again included an instruction in the program to set any such equal to zero. The results are shown in Figure 2.

1.0

~

0.5

x"

"6

0 Fig. 2.

500

Numerical solution to the equation

- ( O c / O x ) ~ ( - gc.l

_ 2 + c . _ 1 - 89

IO00m

-v(Oc/Ox)+ 3(x); parameters as in Figure 1. and negative values of c ahead of the advancing front are replaced by zero. (Oc/~3t)=

- 89 + i ) / A x ,

We note that the width of the advancing front increases negligibly as the front progresses when one uses either Equation (28) or Equation (29) to represent advection, which is a highly desirable feature. We also note that the midpoint of the front (C = 2dmax) moves from left to right at exactly the speed it should for all of the algorithms examined.

180

K. N. CARTER ET AL.

On the basis of our investigation we conclude that Equation (29) is the most satisfactory of the algorithms examined for representing one-dimensional advection by means of finite differences. All of the plots which we have examined have had the source at the boundary point at the left end of the interval at which point of necessity a lower-order formula was used. If we move the source to a normal interior point with most these algorithms, however, we get an unpleasant surprise. At the source algorithms such as Equation (29) are unable to accomodate to the abrupt change in concentration with position, and erroneously generate on the order of 30 to 4 0 ~ more solute than should be present. Also, the onset of the plume, which should be an abrupt step, is rounded by the algorithm. Some typical results are shown in Figure 3a, in which a 37~o excess of solute is spuriously generated by Equation (29).

1.5

1.0

"ID

o o co 0.5 x" O

0

500

" I O'0 0 m

X~

Fig. 3. Numerical solutions to the equation (c~c/&) = - v ( c ~ c / a x ) + 5(x - 40 m); parameters as in Figure 1. (a) In this plot - ( ~ c , , / a x ) ~ ( - ~c,, _ 2 + c,, _ 1 - 89 - ~c,, + l ) except at the end points of the interval. (b) In this plot the above approximation for (Oc,,/c~x) is used except at the points at the source and immediately on either side of it, at which - ( d c , , / a x ) ~ ( c , , _ ~ - c , , ) / A x .

To eliminate this spurious generation of solute in the vicinity of the source we modify our program to calculate - (v Oc/Ox),, from Equation (19) for n - ns (the source point) and ns + 1. This algorithm has a major advantage over most of the others in that it conserves solute exactly under all conditions if v is positive. The results of this modification are shown in Figure 3b. We find that this modification gives us the desired abrupt

MIGRATION OF POLLUTANTS IN GROUNDWATER

181

step increase in solute concentration, and that no solute is spuriously generated in the vicinity of the source. Nor does this modification of our method interfere with the performance of Equation (29) in carrying out the rest of the numerical integration. Oscillations remain small behind the front, and the relatively narrow width of the front does not increase to any extent as the front advances. We note that the introduction of dispersion terms should improve the performance of our algorithm by reducing the cell Peclet numbers ( v A x / D ) from oo to finite values. We note that the utilization of our recipe for representing one-dimensional advection in two-dimensional calculations is very easy if one works in the conformal coordinate systems mentioned above, in which the velocity of the liquid is always tangent to curves of constant v. One can also test the stabilities of our algorithms by examining the propagation forward in time of an error in one of the values c ( x n , t~), as described, for instance, by Collatz (1960)or Roache (1976). For some of our algorithms actual construction of the propagation table quickly becomes extremely laborious, but one can obtain the same result numerically by letting c ( x n , O) = O, n ~ m , c ( x , , , 0) = 1, and then integrating Equation (20) forward in time using the various finite difference representations for the advection term - vOc/Ox. When this was done with Equation (19) we found that the solution remained stable (did not develop extremely oscillatory behavior), but that the peak broadened proportional to t 1/2. All of the central difference formulas tested by this perturbation method - Equations (21), (22), (23), and (26) - quickly yielded highly oscillatory solutions, and the oscillations rapidly propagated in both directions along the x-axis, confirming our previous conclusion that central difference formulas are useless for representing pure advection. On the other hand, the asymmetrical difference representations for the advective term - Equations (27), (28), and (29) - yielded solutions exhibiting an initial localized distortion, followed by negligible further broadening and negligible oscillatory behavior for the full range of the test, 1000 At. In all cases the Courant number ( v A t / A x ) was 0.1. The perturbation test indicates, in agreement with Leonard, that these asymmetrical difference formulas are uniquely well-suited for representing advection in explicit schemes.

5. Example: A Right Angle Boundary with Adsorption and Decomposition We use the potential flow field occurring within an impermeable right angle boundary as our first illustration of the approach described above. The conformal mapping needed is (Wilson, 1982) w = u + iv = z z ,

z = x + iy

(30)

which yields u = x 2 - y2

(31)

v = 2xy.

(32)

182

K.N. CARTER ET AL.

The metric h of the transformation can be written as ,, =

+

L\ax/

(33)

kay] J

which, on substitution and change of variables gives l

h-

(34)

2(u z + 132)TM

The velocity is given by v = A (dw/dz)*

(35)

where A is a real positive scale factor and * denotes the complex conjugate. This yields l vl = v,, = A / h .

(36)

Substitution of these results into Equation (2), and inclusion of a unit source of pollutant at the point (us,/3,) gives a~C at = 4A(u2..{- /32)1[2 {~-u l ~r" 2( u2 +/)2)

~[

+ - - orr" 2(u 2 + Or

v:)'. ~]}

- k c - 4 ( u z + vZ) " 2 5 ( u

TM

a~U] "l-

- 4A(u: + v:) ' / : a-c- au

- u25(v

- v,) - ~ (~, - 4).

rl

(37)

The dispersion terms (in the curly brackets) were represented by (OC~ disp = ! A h 2 dt/i,j 2 i,j ~AU--~

h -1

[( i+ 1,j + hi,-j 1)

(c,+..j

- el,j)

+

+ (hi,) 1 + h i -1 - 1,j) (Ci- l,j -- el, j)] +

+ O~T [(hi,j--l+1 + h,7j ') (c,,j + x - c,,.,.) + Av 2 + (h,-j ~ + h.-) ~) ( % _ ~ "J -

where

c,,j)]}

(38)

h~.j = h(u i, v i). The advection term was represented by {0C~ adv

/\ ~ t / e/ . j

A

1

Au h~'-)2[ -gce_ 2.j + c;_ l.j

1

2c;.j

1

~c;+ 1.j]

(39)

M I G R A T I O N O F P O L L U T A N T S IN G R O U N D W A T E R

183

which is our Equation (29), the best of the algorithms for one-dimensional advection which we have tested. Adsorption was represented by Equation (18), together with Equations (11) and (14). The source term was written as ......

(Oc~ ~t/e.j

= h_ Z. F. bi . 5j,js/(Au. Av ) t.j .s

(40)

where ~pq is 0 ifp :/: q and is 1 ifp = q. The magnitude of the source (kg d a y - 1) is given by F. The program to solve this system required about 389min to carry out a typical run. These were 50 • 50 arrays (1000 • 1000 m) with the time integration carried out for 500 days. A DEC 1099 computer was used. The vast quantities of results generated were most conveniently presented as contour maps; for reasons of economy and speed these were made on the line printer. The uneven spacing of the points in the contour maps is due to the nonlinear nature of the mapping, Equation (30). We used three types of plots - linear, log2, and logm. In all, M represents the point at which the concentration is a maximum. In the linear plots, 9 represents concentrations between 90 and 100~o of the maximum, 8 represents concentrations between 80 and 90 ~o of the maximum, etc. In the log2 plots, 9 represents concentrations between 50 and 100~o of the maximum, 8 represents concentrations between 25 and 50~o of the maximum, and so on. In the loglo plots, 9 represents concentrations between 10 and 100~/o of the maximum, 8 represents concentrations between 1 and 10~ of the maximum, etc. We also remark that the relative sizes of the line printer vertical and horizontal spaces make the x- and y-axes of our plots have scales differing by about 20 ~ ; the plots should be square. Obviously it is not possible to demonstrate the dependence of the model's results on all of the parameters of the model. We illustrate here the effects of dispersion, adsorption, and decomposition. Figure 4 is a linear plot of a run with the dispersion parameters ~L and c~v set equal to zero. The decomposition rate constants and the capacity factor ~max for the adsorption of solute were also set equal to zero. This run should exhibit a string of 9's clear out to the end of the plume, from which point on we should see only blanks. The use of conformal coordinates eliminates transverse numerical dispersion, but we see that longitudinal dispersion has resulted in a diffuse front to the plume. The thickness of this front is roughly 19Y/o of the length of the plume, which represents about a 40-50~o reduction in longitudinal dispersion compared to results obtained using Equation (19) to represent the advection term. The log2 plot of the same plume is shown in Figure 5. We next introduce physical longitudinal and transverse dispersion by setting ~L = C~T= 0.3. Adsorption and decomposition are still excluded. The linear plot is shown in Figure 6. The plume is still rather narrow, but we see that the effects of these dispersive terms are such as to swamp the substantially smaller effect of numerical longitudinal dispersion. The elapsed time since initiation of the source in this run is, as before, 500 days. The log2 plot of this run, shown in Figure 7, permits us to see the lower levels of pollutant which border the plume.

~LkTtVE CON~£kT~tXldX~ H | G H E S r CO~CN = . 1 J l b [ - O l L O i E S T COdCN; . ~ 0 0 ~ * 0 0 .50U0£*03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

L I ~ E X k PLOT u v TIH[= *5000~OJ

TOCAL POLLUTkNT :

1 43 7 b Y 9

2 2

5

9

9

9 9 9

9

9 Y 9 9 9

9 9 W 9 ~s

V

Fig. 4. Movement of pollutant within a right angle baKzier. The field is 1000 x 1000 m, the aquifer thickness is 1.0m, A =0.002; eL = ~ r = 0 ; k = k ' = 0 ; r/=0.3; ~ = 0 ; xs = 100; y , = 9 0 0 m ; A t = 2 d a y s ; nx = ny = n u = n o = 50, flux = 1.0 kg d a y - 1. This plot is on a linear scale. 500 days have elapsed since the source was initiated. The total quantity of pollutant present is 500 kg. LOC B A $ ~ 2 P&Of Or" ~ C L A T I V ~ C O N C ~ T R A T X O N ~ TIME= .5000~*03 ~I~HEST CO~C~ = . 1 ~ 1 6 ~ - 0 1 LOdEST CO~Cg= .O00UL÷O0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7e55

9

431

9

g

9

9

9 9 9 9 9

9 9 9 9 9 9 9 9

9

9

9 9

9

Fig. 5.

Movement of pollutant within right angle barrier. This is the same plume as shown in Figure 4; here it is plotted on a log2 scale.

LIN~A~ PLUT 0¥ k~LA'H~E CUNC~&T~*T|ON~ T I H L : ,5000[+0d dI~HLST CONC~ = . I ~ l P B - 0 1 L~EST CD8¢~ .O00U~+00 YOT~L P~LLU~AN~ = ,49~2~+03 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

b

1

b

6

1

I

I i

I i

i I

1

h 6 6

1 7

1

I

I I

* I

1 1

7 L;

I

7

1

1 t

'1

I * 1

7 1 T

1 71

U 9

Fig. 6.

M o v e m e n t of pollutant within a right angle barrier; effect of dispersion. Parameters for this run as in Figure 4 except that c~r = c~r = 0.3. The concentration scale is linear.

LOG ~*SE 2 PL~T U~" kLLATIYE CONG[~THATIOkS TIHE= o5000E÷03 hIGHEST CO~CN = .1216£-01 LD~ES~ CO,Ca= ,0000~+00 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

3 76

3 4

9 ?

9 7

3

7

7

9 9

7

9

4

4

4

4 4

3

7

7 7

7

9H 44

5

7

~ 1

666

4443331-

IL2 1

1 1 l

1

4 7

Y

7 7

7

4 4

9 7

3

4

9 7 9

3 7

3

3

7 7

3

7

9

3

3

73 9

7

36

3

3

96 a

3

b

26 2 b

6 63 96 2 962

269 2 696 1 b 1

62 2

9oI

6

~61

159 5~5

51

5~5

1 3

Fig. 7.

4 3

M o v e m e n t of pollutant within a right angle barrier; effect of dispersion. This is the same plume as s h o w n in Figure 6, plotted on a log 2 scale.

186

K. N. CARTER

ET AL.

In both these runs (Figures 4-5 and Figures 6-7) solute is not being destroyed or removed by adsorption, so it is possible to check the algorithm by doing a mass balance. Early in the runs (after 100 days) we find a 2~1 oYo discrepancy; this decreases as time increases, suggesting that the algorithm has some difficulty with the initial very sharply peaked plume, but conserves matter adequately after the plume has become smoothed out somewhat. Let us next permit the aquifer to adsorb pollutant according to a Langmuir isotherm with ~max = 1.0 k g m - 3 and e l / 2 ----- 10 - 5 kgm -3. The pollutant mass transfer rate coefficient, ~, we have set equal to 0.1 days-]. We set e L and er = 0.3, as before. In Figure 8 the linear plot shows a very markedly curtailed plume at 500 days as compared to those for the cases without adsorption. The lOglo plot of Figure 9 shows a similarly reduced plume, as expected. The loglo plot of this plume, shown in Figure 10, shows curtailment of the extent of the plume, too, but indicates that in those regions where the water velocity is relatively high (away from the upper left comer of the figure), the small mass transfer rate coefficient ~ and the high flow velocity of the water permit the 'leakage' of low concentrations of pollutant for considerable distances in front of the main body of the plume. The effect of chemical decomposition of the solute is seen by comparing Figures 11 and 12 with Figures 6 and 7. The run illustrated in Figures 11 and 12 describes the migration of a solute with a first-order decomposition rate constant of 0.01 day - L. The plume at 500 days is drastically reduced compared to that in the absence of reaction, klkEkk PLOT d~" K~L~IV~ C0NC~II(ATIdlI~

TI~,[= . 5 0 0 0 ~ * 0 3 nlGabS! CONCK : . 1 2 1 1 s LOWEST CONCH= .0000~,*00 TOT &L PI3.LL_UT./,NT = 9 L3UOs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 4 U e

il

u o 9 9

"~ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig. 8. Movement of pollutant within a right angle barrier; effect of adsorption. Parameters for this run are as in Figure 6 except that ~c= 0.1 day-1, ~ma~ = 1.0 kg m - 3 , c~/2 = 10-5 kg m - 3 The concentration scale is linear.

LOG BASE 2 PkUT OF E E L ~ T I V E C U U C ~ T ~ a T I O N $ TI~ ,5000E*03 ~XGH~ST CONCa 9 . 1 2 1 1 s LOW~$t CO~C~= .OUOOs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

o

2 9

2

2 2

9

3 4

94

4 494 39

4

Fig. 9. Movement of pollutant within a right angle barrier; effect of adsorption. This is plume shown in Figure 8. LOG B A $ ~ 10 PLdT OF N L ~ A T I Y E C O m C E b T R I ] ] O N S TIME= .500~+OJ ~ I G N L S I ~ONCN = . 1 2 1 1 [ - 0 1 .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

~0~EST .

.

.

CONCk= . 0 o U 0 ~ * 0 0 . . . . . . . . .

.

1

1

.

.

.

.

.

.

.

I 1 1 1 1

11

1

1, 11 I I

I

I

1

4

3 6 7 6 "1

3

5

3

5

4

4

46

~ B b

8 ~ uO 4689 06 46 B 9 e

-

.

.

I 1 1 11111 I I 1

Z 11111 1 1 1 1 1 1 1 1 1 -

-

111111~I~1111, ~1 ,111 111

1 z I i 1 1 1 1 1 I I i

11

1

I I I

I

I 1

I

1 1

1

2 4 I

6 4 4 4

.

1

1

b4

~ a 6 ~b

.

1

1

I

.

I I

J2

51

68

o e9

L 1

3

S

~7 7 V

4

3 3

I 1 I

.

I

1

1

I I

I

1

I 2

5 ? 5 7

5 "/

5

-

4

9 9

3

1

1

4

9 7

7 ?

1 2

5

4

-

I

1 1

1 1

I I 1

1 4 2

11

.

I

I

l I

.

plot of the

a log 2

1 1

1 1

31

.~ 4 2 i

I

461J ~9 6421 40 U o 4~ 357 u 64 ~ ~ 4 S7 b J I b7o

,J2

I~Uk COMPLs I'L b

Fig. 10. Movement of pollutant within a right angle barrier; effect of adsorption. This is a ]ogto plot of the

plume shown in Figure 8.

L I [ ~ A ~ PMUT d f MG~ATIbE ~ U N U ; k t ~ A t I d ~ $ TANEs . 5 0 0 0 L + 0 ~ HIGH~ST C0kCk 9 . 1 1 4 e s ~U~LST C U ~ C ~ . 0 0 0 0 [ + 0 0 ~O~hL P U L L U T A ~ 9 , v ? b ~ g * ~ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

I I I

2 3

t

V

Fig. 11. Movement ofpoUutant within a right angle barrier; effect of chemical decomposition. Parameters for this run are as in Figure 4, except that ~L = ~r = 0.3, k = k' = 0.01 day- 1 There is no adsorption, and the concentration scale is linear. LO~ ~ASL ~ P L d l u r ~ L a T A d L C~AC~nTKATIONS TAeE= * b b 0 ~ L - 0 ~ h I G , ~ S ~ COdex - * 1 1 4 0 ~ - 0 1 LDWL$1 ~O~CX= , ~ V ~ L * 0 ~ . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . .

L[ I

11 22

1 t 1 1

1 dl

2

111 22 1

22 111

222 11

~2222 1 1 1 1 1 1 ~

11

3 3

1 ~

2

1

J

2 11

2

4

9t 5 5

2 -"

3 3

3

b b

4 b

,t

1

'. 4 ~, 5

1

1

"/

1 4 1

"t I

1

?

b

b

1

~

"t 5 1

1 . 15o 1 ~'J5 1 5 5

51 51 1

~*bl

b y~ 3

4 e 3

.

Fig. 12.

Movement of pollutant within a right angle barrier; effect of chemical decomposition. This is a log 2 plot of the plume shown in Figure 11.

LINEA~ r

P ~ U T dP ~ L ~ l l v g

Cr

pO~koTa~ = .622~-b2 . . . . . . . . . . . . . . . . . . . . . . . . . . .

~

. . . . . . . . . . . . . . . . . . . . . . .

1

Fig. 13. Movement of pollutant within a right angle barrier; effect of adsorption and decomposition. Parameters for this run are as in Figure 4, except that ~tL -- ctr = 0.3, k = k' = 0.01 d a y - i, x = 0.I d a y - 1, ~,~a~ = 1.0kg m - 3 , and cl/2 = 1 0 - s k g m -3.

LO~ r 2 PLOT u~ ~L~TIV~ Ca&Cs163 71~= .5000~oJ HI~N~sr (OXCN = ,IAr163 LO*LST . . . . . . . . . . . . . . . . . . . . . . . . . . .

CU~C&: .OUOOL*O0 - . . . . . . . . . . . . . . . . . . . . . . . .

2 4

l

2 6 2

4 4

0 ~ 4

~ 5 4 9 4 4 ~ 3 .i

Fig. 14.

Movement of pollutant within a right angle barrier; effect of adsorption and desorption. This is a log 2 plot of the plume shown in Figure 13.

190

K . N . CARTER ET AL.

shown in Figures 6 and 7. The linear plot of the plume (Figure 11) also shows the expected exponential attenuation with distance from the source. Other data (not shown) indicate that this plume has essentially reached a steady state after 500 days. Both chemical decomposition (k = 0.01 d a y - t) and solute adsorption were included in the run producing the plume shown in Figures 13 and 14. As expected, the extent of the plume after 500 days is much less than that of the non-reacting, non-adsorbing solute examined in Figures 6 and 7. It is also less than that of the plume for the run with a solute which reacts but is not adsorbed (Figures 11 and 12), again as one would expect. 6. E x a m p l e : F l o w around a Cylindrical O b s t a c l e

For a second illustration we examine the movement of a decomposing non-adsorbing pollutant in a two-dimensional potential flow field in the vicinity of a circular obstacle. The function a 2

w = z +--

(41)

= u + iv Z

maps that portion of the z-plane outside a circle of radius a into the w-plane. We start with & = 7 . (DVc) - v. 7c - k c + ~ s~b(r - ri) (42) Ot

t= 1

in which we include the possibility of several point sources of pollutant fluxes s t located at r;. In the u, v coordinate system Equation (42) is essentially the same as Equation (9); the source terms are added and the adsorption term is deleted from Equation (9). The metric h for this transformation is more complex than that for our earlier case (the right-angle comer). It is possible to obtain h in terms of u and v directly, but the resulting expression is objectionably cumbersome. We therefore express h in terms of x and y, noting that these are functions of u and v. We find 2(y__22_xZ)]}(x-2 + y2)Z

aZ[a 2 +

h=

1+

1/2

u 1 {[(u 2 - v 2 - 4a2) 2 + x = - +_ 2 23/2

(43)

4 u 2 v 2 ] 1/2 + ( u 2 - v 2 + 4 a Z ) } 1/2

(44)

4 u Z v 2 ] 1/2 - ( u 2 - v 2 - 4 a 2 ) } 1/z .

(45)

(+ i f u > 0 - if u < 0)

1 {[(u 2 - v z - 4aZ)2 + y = ~v + 23/2 (+ i f v > O - if v < O)

191

M I G R A T I O N O F P O L L U T A N T S IN G R O U N D W A T E R

We locate our grid so that the center of the circular obstacle is equidistant from the four nearest grid points, thus avoiding problems which arise if u or v = 0 in the grid in the w-plane. The discrete approximation to the advection-dispersion equation which we used is given by

'~i~ h-={~ & 2 ,

_

i,j

O~L

-- 1

[(hi_ t,j + hi./l ) (c, _ t,j - c,.j) +

+ (hi.31 + hi+- 1 ,.j) (c,+ ,., - cij)] + 0{T

+ (Av) 2 [(h,5'- i + h,5 l) (c;.j_ i - ci.j) +

+ (h,5' + hi.i+ i ) ( c i . j + l

+

Ah~

2

[

Au - kc~ + ~

1 - - g C i _ 2, j "1- C i _ l , j

Sljh~ZSitbjj

) J

- c,.j)D +

1 9 -- 2ci, J -

1 3ci+

1,j]

-

(46)

l,J

where Slj = Fij/(Au. Av ). As before, the discrete equations must be modified at the boundaries. In addition, the equations must be modified along the line segment in the w-plane into which the circular obstacle is mapped to prohibit diffusion across this line segment. The physical parameters of the model are like those examined m the flow field within a right angle. Therefore, we take this opportunity to examine some of the computational aspects. We will discuss restrictions on the Courant number, variations of Courant number, and effect of grid spacing in the longitudinal and transverse directions. The geometry of this problem requires some care in controlling the value of the Courant number, Nc = I vl A t / A x . A uniform spacing of points along the u-axis dictates nonuniform spacing of points in the x-direction, and the spacing is smallest precisely in the regions of highest v~ for x = 0, y = + a. From these considerations comes the requirement that Nc(] rl --, oo) _< 0.25 if the maximum value of N c (as defined for the corrector step of the integration) is to be _< 1.0. This requirement is necessary for stability in the time integration. The predictor step of the predictor-corrector algorithm conceals a further constraint, for here a value of 2At is used. Nc (predictor) greater than 1.0 may lead to practical instability, but does not always. In Table I are summarized the results for advection-dispersion in one dimension with constant velocity. The amplitudes of the oscillations are small.

192

K. N. CARTER ET AL.

J

Fig. 15, Movement of pollutant around an impermeable obstacle of circular cross-section. Arrows indicate magnitude direction of flow field. Field is 49 by 29. Vx[~ = 1.0. NU = 50, NV = 59, eL = ctr = 0.015, mass flow rate (w) = 2.0 (two sources, of 1.0 each, placed symmetrically on grid lines closest to the obstacle, to simulate a single, on center, source), At = 0.248, time = 39.68, total pollutant = 78,55, concentration at source = 1,925. Solid lines are isograms of concentration. Plot is on log 2 scale with highest displayed value = 0.9875. Fifth level is missing (see text). Sharp bends in contours a short way from the source, and any asymmetries with respect to axis, are artifacts of the contour interpolation algorithm, not of the numerical data. In Figures 15-19 and 21-22 the solute is assumed to be non-adsorbing and not to undergo decomposition. TABLE I Dependence of solution stability on Courant number N C(corr.)

1.00 1,02 0.50 0.42 0.50

N c (pred.)

2.00 2,04 1.00 0.84 1.00

Maximum amplitude of oscillation

Pa

~o of average local value

33.3 33.3 33.3 5.0 10.0

3.4 3.4 1,2 no oscillation no oscillation

Pa, the cell Peclet number, is defined as v(Ax)/D. Our t w o - d i m e n s i o n a l f l o w field a r o u n d an i m p e r m e a b l e obstacle is less forgiving. For the run with results s h o w n in Figure 15, the m a x i m u m value o f the C o u r a n t n u m b e r in the predictor a p p r o a c h e d 2. Strong oscillations in calculated concentration, a m o u n t i n g to a p p r o x i m a t e l y 5 0 ~

a b o v e and b e l o w the averaged local values, were e n c o u n t e r e d in

the regime o f rapid acceleration occuring near the m a x i m a o f velocity at the top and b o t t o m o f the obstacle. This w a s severe e n o u g h to fool the c o n t o u r plotting routine;

MIGRATION OF POLLUTANTS IN GROUNDWATER

193

Fig. 16. Plume of pollutant. Similar to run in Figure 15, but with Courant number halved by reducing At to 0.124, total pollutant = 78.38, concentration at source = 1.926. Log 2 scale with highest displayed value = 0.9875.

Fig. 17. Plume of pollutant. Similar to run in Figure 15, but with tighter longitudinal grid spacing. N U = 99, N V = 59, At = 0.124, total pollutant = 78.58, concentration at source = 1.975. Log 2 scale with highest displayed value = 0.9875.

194

K . N . CARTER ET AL.

contortions of the contour lines are seen near the top and bottom of the obstacle, the second contour level incorrectly intercepts the obstacle, and the fifth contour level does not even appear except as a small dot! We therefore reduced At by a factor of two, thereby likewise halving the Courant numbers. The result is shown in Figure 16, where the spurious oscillations have been eliminated, and all specified contour levels appear. Downstream from the obstacle, the outermost four contour lines match those of Figure 15 to within the width of the pen fine. Two gratifying conclusions may be drawn the localized unphysical oscillations of Figure 15 did not corrupt the solution downstream from their occurrence, and the results of the calculation appear insensitive to the At chosen, provided stability criteria are adequately satisfied. The effect of exceeding N~ < 1.0 in the predictor is not always so severe as in Figure 15. Figure 17 shows a plot with a denser grid in the longitudinal direction, and the maximum values for Nc in the predictor and corrector approaching 2.0 and 1.0 respectively. The wiggles in the contour lines are small. Nevertheless, keeping all Nc < 1.0 might be preferable. Comparison of Figures 17 and 16 reveals the effect of changing the grid spacing in the longitudinal direction. With coarser grid spacing, one expects more longitudinal numerical dispersion on account of truncation error in the advection term, and the plume of Figure 16 is 4~o longer than that in Figure 17 in accord with this expectation. Uncertainties in input parameters render such a small discrepancy irrelevant for practical groundwater modelling concerns. The effect of a change in transverse spacing is dramatically illustrated by contrasting the linear scale line printer plots presented in Figures 18 and 19. We note the relatively longer persistence of 9's in Figure 18 along the streamlines containing the source, in contrast to the more rapid attenuation in Figure 19. Such persistence along these streamlines is a direct consequence of the coarser grid spacing, as made clear by the following argument. If we regard the transverse cross section of the plume to be roughly approximated by a spreading Gaussian, we note that the central difference operator underestimates the absolute value of the second derivative at the peak of a Gaussian, and this underestimation becomes increasingly severe with larger transverse node spacing. Such an underestimation inhibits transverse diffusion from the streamline of highest concentration, just as we have observed. The operator is conservative, so mass is not lost, but is less accurately distributed. Yet another consideration may be invoked. For this it may be helpful to regard the control volume formulation of the finite difference scheme as giving values at the nodepoints that compare, not to the values of the exact solution at that point, but rather to the averaged value of the exact solution within a superimposed control volume. This latter will approach the former in the limit of infinitesimal spacing (except for problems involving the source, as discussed below!). The ratio of the concentration along the source streamline to that at any gridpoint along the line will approach the ratio of the peak heights in the limit of small domains of integration, while in the limit of long intervals over the transverse coordinate, the ratio approaches 1.0, for a conserved substance. It can be contended that this effect, too, may be related to the persistence of 9's (high concentrations) in Figure 18. In Figure 20a and b, this effect of spacing is illustrated for integration over spreading Gaussians. -

MIGRATION

OF

POLLUTANTS

IN GROUNDWATER

195

LInEA~ PLbT Of ;fLATIYE CO~CE~r~ATIO~S

TIMb= .3968L-U2

11777622

,99~9999,9~ Z19 9 9999999

lib8 1188 118 8 ~8

bb22 bb22 622

o

222

'11

b~

B~

b~

i

6 B~

~ 32

2 2 ~ I I

bb22

llTTTb22

Fig. 18. Plumeof pollutant. Linear scale plot related to concentration at the maximum, diminishing in increments of I/10 the maximum concentration. Run similar to that portrayed in Figure 17, but with sparser transverse grid spacing. N U = 99, N V = 30, N X = 50, N Y = 30, At = 0.124, total pollutant = 78.53, concentration at source = 1.027.

blNEAR PLDT OP ~s TZ~s .3968~*02

C~NCENTRATIO~$

1344442~ 1555

265

S5

21 4121

2~2221 2112~T

12

111

1 !

~S 16s 265 )555

t

332 4323 55)21

I I

I

13444421

Fig. 19. Plume of pollutant. Linear scale plot related to concentration at the maximum, diminishing in increments of 1/10 the maximum concentration. Same run as portrayed in Figure t7.

Another feature in this connection concerns the concentration at the sourcepoint(s). Because we are modelling the plume due to a given mass flow rate (mass per unit time) o f contaminant, the resulting concentration at the source is not physically significant, but rather a consequence o f the size and shape o f the control volume around that point and the velocity. For the simple 2 point upwind advection formula, the following

196

K.

'

,

.

,

.

.

,

.

C A R T E R

I

,

.

N.

.

,

.

,

.

,

ET

AL.

4

,

~

t-t-t

t

,

,

.

i

9

.

9

.

.

I

.

.

i

.

i

I

i

.

Fig. 20. Effectof spacing on integrationof normalizedGaussians with standard deviationsof 1.0,2.0, and 3.0. The points show averages over domains of integration of 1.0 and 4.0 respectivelyin a and b. rigorous expression was derived to give the steady state concentration at the source for flow in the u, v coordinate system with no dispersion, c, = w/(v. h" Av), where w is the mass flow rate. The analogous expression for one dimension is c, = w/v. Tests in one dimension revealed that the presence of dispersion lowers the actual values slightly below this. While the total mass is virtually identical for Figures 18 and 19, the concentration at the source in the latter is nearly twice that in the former, in qualitative agreement with the prediction of our equation. Nevertheless, while values near the source are not physically meaningful, values at some distance from the source must be nearly independent of spacing (for sufficiently dense spacing) if the model is to possess predictive meaning in a practical sense. Comparison must be at the same absolute concentrations in each case, not to concentration at the source. For variation in longitudinal spacing such comparison has already been presented in Figures 16 and 17, and the criterion of precision was satisfied. Figure 21 is based on the same data as Figure 18, but with log2 contours at absolute levels identical to those of Figure 17. Comparison of Figures 21 and 17 reveals severe discrepancies near the source, but the contours at lower concentrations once again agree to within practical limits. The length of the two plumes at the lowest contour level plotted differ by only 2.5~/o. Thus we have demonstrated the adequate precision of the model for the longitudinal and transverse spacings employed. The physical character of the solutions is mostly self-explanatory, but a few aspects deserve mention. The plumes narrow slightly near the top and bottom of the obstacle, even though the dispersion coefficients are largest there. This is a consequence of the compression of streamlines in this area. The increased values for dispersion coefficients

MIGRATION OF POLLUTANTS 1N GROUNDWATER

Fig. 21.

197

Plume of pollutant. Log2 plot of same data as in Figure 18, with highest displayed value = 0.9875, to facilitate comparison with Figures 16 and 17.

-...

.2,

---"

Fig. 22. Plume of pollutant. Single source, w = 1.0, placed above axis of symmetry. N U = 99, N V = 59, Al = 0.062; other parameters as in Figure 15. Time = 39.68, total pollutant = 39.21, concentration at source = 1.918. Plot is on log 2 scale with highest displayed value = 0.959. Lowest level has reached edge of computation grid.

198

K. N. CARTER ET AL.

in this region counteract the narrowing somewhat. The notched appearance of the leading edge of the plume is a consequence of the reduced velocities on the streamlines nearest the obstacle at points close to the front and rear of the circle. It was not obvious until the model had been run that this effect would outweight the increase in velocity near the top and bottom of the circular obstacle. One could check the supposition by a comparison of average velocities along streamlines over the distance of interest. We note that our conformal mapping approach treats only inviscid flow; the notch would be still more pronounced for no-slip boundary conditions at the obstacle. Figure 22 presents a plume from a source placed off of the axis of symmetry of the flow field. Note the convergence and subsequent divergence of the plume as it swings around the obstacle. Aside from this, the results are really not too different from those obtained with a uniform flow field. 7. E s t i m a t e

of Mass

Transfer Rate Constant

In the previous sections of this paper we have taken ~, the constant controlling the rate of solute transfer from the moving bulk liquid to the stationary adsorbed phase, as a constant to be determined empirically. In this section we examine a means of estimating this mass transfer rate coefficient. The solute must diffuse through a boundary layer of comparatively motionless liquid in order to be absorbed; we examine this diffusion process in planar and spherical geometries. The planar case is analyzed as follows. The diffusion equation is ~C 02s -- = D -Ot Ox 2

(47)

where D is the diffusion constant, c the solute concentration, and x is the distance measured from the boundary between the boundary layer and the bulk liquid toward the liquid-solid interface, at which x = b. The boundary condition at x = 0 is c(0, t) = co~.

(48)

The boundary condition at b is obtained as follows. From Fick's Law we have D --Oc (b, t ) - dF , c~x dt

(49)

where F is the solute surface concentration. If we assume a Langmuir isotherm we have F =

Fmax 1 + Cl/2/c

(50)

which must be linearized to make further progress, yielding I~ = Fmax C = KC.

Cl/2

(51)

199

MIGRATION OF POLLUTANTS IN GROUNDWATER

Differentiating Equation (51) and substituting in Equation (49) gives x; ac (b, t) = D --~c (b, t)

~t

(52)

Ox

as our second boundary condition. The diffusion equation is readily solved by separation of variables; the result is

+ B,~cos ( ~ / ~ x ) l exp(-2/) .

(53,

The boundary condition at x = 0 and the requirement that

c(x, ~) = coo field

c(x,t)= c~ + ~ Aa sin(~/~ x)exp(-

(54)

The boundary condition at x = b, Equation (52), gives

_ K ~ A asin D

b = -Aa

cos

b .

(55)

This must yield non-zero values of Aa, which requires that tan ( ~ / f b) = , J 2 K - 1 T A B L E II Least root o f t a n u - o~/u= 0 for various

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0,9 1.0 1.5 2.0 3,0

0.31105 0.43284 0.52179 0.59324 0.65327 0.70507 0.75056 0,79103 0.82740 0.86033 0.98824 1.07687 1.19246 1.57080

(56)

200

K . N . CARTER ET AL.

for the non-zero values of 2. (We note from Equation (55) that 2 --- 0 is an eigenvalue of the system, too.) We then take as our value of ~ the smallest positive value of ~ satisfying Equation (56). These can readily be calculated from the first root of Equation (57). tan u - ~/u = 0

(57)

where ~ = b/K and u = (2/1)) 1/2 b; this root is given as a function of ~ in Table II. We next turn to the case with spherical geometry. The diffusion equation is -

dt

r2

.

(58)

r 2 dr

Let us assume that (59)

c(a, t) = c o

gives us the condition at the junction of the bulk solution (r < a) and the boundary layer (a < r < b). At r = b, the liquid-solid interface, we have (as before) K --acat(b, t) = - O _aCor(r, t)r=b

(60)

"

The diffusion equation, Equation (58), is readily solved by the substitution f(r, t) = rc(r, t), followed by separation of variables. The result is

=--+Bo+r

E

~sin

r +

r

The requirement that c(r, oo) = coo gives Ao = 0, Bo = coo. The inner boundary condition, Equation (59), then gives (62)

a-11A,~ sin ( ~ f ~ a ) + B,~cos(~/~ a ) l = 0.

The outer boundary condition, Equation (60), yields, on substitution of Equation (61) and simplification, 0 = Aa

I( 1

K 2 ~ s i n ~ / r - l N / ~-

~ + Db]

+B'~[(12+K2)c~

b

cos ~/~ b1 +

b+lb

sin

b].

(63)

MIGRATION OF POLLUTANTS IN GROUNDWATER

201

In order to have non-zero values of the Az and B~, the determinant of the coefficients of Equation (62) and (63) must vanish; this gives us the requirement that 0 = (~2 + K'~sin(-\N /~b~c~ D / \~t~D a)

c~

-

- ~ ~ / ~ [ c o s (X/r~ b)cos (N/~ a ) + sin (Xf~ b)sin ( X / ~ a ) l . (64) This in turn leads to the eigenvalue equation 2x/~

b- ~+ (K~/D)

= tan [ x ~ ~- ( b - a ) l .

(65)

The smallest positive eigenvalue of Equation (65) we take as our mass transfer rate coefficient ~c. We see that, if we let b ~ ~, b-a being held fixed, this result approaches Equation (56), as expected. It is relatively easy to show that the solution to the planar case, Equation (57), gives an upper bound to the solution of Equation (65). We let fl = b - a, u = x//~D and solve Equation (65) for b- l, obtaining b-1

= u c o t (~u) - K u a - f ( u ) .

(66)

Then, after differentiating and rearranging

df-cscOSu)fc~

sf~-~ul -2Ku"

(67)

In the range 0 < u < rc/fldf/du is negative, so f decreases with increasing u in this range. Therefore b = 1If increases with increasing u and u increases with increasing b in this range. Since f(u) is certainly negative for u > n/2fland positive for u = 0 § we conclude that the desired least positive root of Equation (66) does indeed lie in the range 0 < u < n/t, and that therefore the desired least positive value of ;t increases with increasing b. This makes the least positive root of Equation (56) an upper bound to the least positive root of Equation (65). Acknowledgment

This work was supported by the National Science Foundation. References Anderson, M. P.: 1979, 'Using Models to Simulate the Movement of Contaminants through Groundwater Flow Systems', CRC Crit. Rev. Environ. Control 9, 97-156. Bredehoeft, J. C.: 1971, 'Comment on "Numerical Solution to the Convective Diffusion Equation"', by C. A. Oster, J. C. Sonnichsen, and R. T. Jaske, Water Resour. Res. 7, 755-756.

202

K. N. CARTERET AL.

Collatz, L.: 1960, 'The Numerical Treatment of Differential Equations', 3rd ed., Springer-Verlag, Berlin. Gardner, A. O., Peaceman, D. W., and Pozzi, A. L., Jr.: 1964, 'Numerical Calculation of Multidimensional Miscible Displacement by the Method of Characteristics', J. Soc. Petrol. Eng. 4, 26-36. Grove, D. B.: 1977, 'The Use of Galerkin Finite-Element Methods to Solve Mass-Transport Equations', U.S. Geological Survey, Water Resources Investigations 77-49, 61p., NTIS, PB-277532. Leonard, B. P.: 1979a, 'A Stable and Accurate Convective Modelling Procedure Based on Quadratic Upstream Interpolation', Computer Meth. in AppL Mech. and Eng. 19, 59-98. Leonard, B. P.: 1979b, 'A Survey of Finite Differences of Opinion on Numerical Muddling of the Incomprehensible Defective Confusion Equation ~, in Finite Element Methods for Convection Dominated Flows, T. J. R. Hughes (ed.), American Society of Mechanical Engineers, New York, N.Y. Ogata, A.: 1970, 'Theory of Dispersion in a Granular Medium', U.S. Geol. Surv. Profi Paper 411-I, 34 pp. Oster, C. A., Sonnichsen, J. C., and Jaske, R. T.: 1970, 'Numerical Solution to the Convective Diffusion Equation', Water Resour. Res. 6, 1746-1752. Pinder, G. F.: 1974, 'Progress in Simulation of Contaminant Transport in Porous Media', in A Decade of Progress in Water Resources, S. C. Csollany, Z. A. Saleem, and W. J. Roberts (eds.), American Water Research Association, Urbana, I11., p. 223. Roache, P. J.: 1976, 'Computational Fluid Mechanics', Hermosa Publishers, Albuquerque, N.M. Shamir, U. Y. and Harleman, D. R. F.: 1967, 'Numerical Solutions for Dispersion in Porous Media', Water Resour. Res. 3, 557-581. Scheidegger, A. E.: 1961, 'General Theory of Dispersion in Porous Media',J. Geophys. Res. 66, 3273-3278. Wilson, D. J.: 1982, 'Migration of Pollutants in Groundwater. I. Reduction of Numerical Dispersion by Conformal Mapping', Environ. Monit. Assess. 2, 309-330.