- Email: [email protected]
Chapter 6 Regional Pollutant Transport Models
188
Chapter 6
REGIONAL POLLUTANT TRANSPORT MODELS SCOPE OF APPLICATION Groundwater pollution due t o human activities occurs in a number of ways. Several situations are shown i n figure 6.1. Polluted water may infiltrate into the aquifer from a polluted surface water body, from leaking waste-water pipes, ponds or cesspools (bacterial pollution). The pollutant may also be leached from the top soil by rain and be carried down into the saturat _1 zone by seepage
6.1
(nitrates, pesticides).
Similarly, seepage through landfills and waste deposits
is a source of pollution. Pollutants may enter the soil in a form immiscible with water (e.g. oil, chlorinated hydrocarbons). Their gradual dissolution by seepage water or the passing groundwater flow causes groundwater pollution. While the transport of pollutants in the unsaturated zone is essentially limited to vertical transport between groundlevel and the top of the saturated zone, long-range spreading of pollutants is only possible in the saturated zone where dissolved pollutants are carried along by the prevailing flow, mainly in horizontal directions. The distance that can be covered within a certain time depends on the velocity of flow and the persistance of the pollutant. We consider the often complicated processes in the unsaturated zone on the whole as a source term for the transport in the saturated zone and limit the further discussion to the saturated zone. Depending on the density of solute, the pollutant may influence the flow field (hydrodynamically active solute). We only consider hydrodynamically inactive solutes, this means pollutant concentrations that are so small that density induced flow can be ignored. We further restrict our discussion to regional transport which means the horizontal scale of transport is much larger than the thickness of the aquifer. I n these cases the flow field can essentially be considered as horizontally two-dimensional. In the following chapters we assume that regional pollutant transport can also be described by horizontally two-dimensional models as long as vertical averages of concentrations are considered. In the vicinity of infiltrations or in strongly stratified aquifers, three-dimensional aspects transport may, however, be important. Such situations are shown i n figure 6.2. Any transport model requires as an essential input the velocity field of flow. We assume that the flow field is either known a priori or is modelled parallel of
to transport by a flow model. The flow model yields heads, from which specific flowrates (filter velocities) are calculated by means of the Darcy-law. The transport model requires pore velocities. These are obtained from specific flowrates by division by the effective porosity ne' The derivations given below follow the work of Bear.
preci pi tat ion
P
Fig.
6.1: Sources
of
groundwater pollution
00
(0
190
1 1 1 1 1 1 1 waste deposit
surface table
streamline
concentration
a) Strong vertical flow component
b) Both vertical velocity distribution and concentration distribution are non-homogeneous
Fig.
6.2: Examples o f situations where depth-averaging is not reasonable
191
PHENOMENA TO BE CONSIDERED AND BASIC EQUATIONS
6.2
We first discuss the behaviour of an ideal tracer which experiences neither adsorption nor chemical transformation. On the microscopic level of a single pore space, there are then only two basic processes of transport. These are the convective transport by the water flow and molecular diffusion. The convective pollutant flow through an infinitesimal area element dA at a point T3 inside a pore is given by
- ( U ~ . ~ ~ , U ~n3) =. (nx,ny,nz),
G3.E3 c dA
u c dA
=
with
G3
=
- Jconv' "3 dA
(6.1) =
(x.y.2).
-
and jconv= u3c
z3
un is the velocity component normal to the area element, the 3-dimensional velocity vector at location y3, and G3 the unit normal vector of the area element dA. (Figure 6 . 3 ) The diffusive pollutant flow through the same area element is given by Fick's law. It is proportional to the negative gradient of concentration normal to the area element dA.
The total flux at any point in space is given by
For
all practical purposes we have to look at macroscopic volumes when taking
pollutant balances. Therefore we must calculate averages of the flux vector ( 6 . 3 ) over volumes comprising many pore spaces. Averages can be taken over a geometrical volume or over the portion of pore space contained in it, which is effectively taking part in flow. Averages over total volume are obtained from averages over the effective pore space by multiplication by the effective porosity ne'
I n the averaging process it is useful to decompose local concentration and velocity into a macroscopic average over the effective pore space and the local variations from that average.
c = c +6c
192
F i g . 6.3: Nomenclature for the definition o f infinitesimal convective
pollutant flux
variable poresize
velocity profile within pore
bending of streamlines around grains
Fig. 6.4: Causes of microscopic variations i n velocity
193 I n s e r t i n g (6.4)
i n t o (6.3).
m u l t i p l y i n g b y t h e e f f e c t i v e p o r o s i t y and a v e r a g i n g
-V
we o b t a i n t h e v o l u m e t r i c a v e r a g e f l u x j 3
(6.5) Terms
-
r 3 6 c and
s i d e o f (6.5)
c6G3
vanish d u r i n g averaging.
The f i r s t t e r m o n t h e r i g h t - h a n d
i s c a l l e d t h e c o n v e c t i o n ( o r a d v e c t i o n ) term.
the molecular diffusion
term.
The s e c o n d t e r m i s
caused by t h e v e l o c i t y p r o f i l e w i t h i n a pore,
These a r e
t h e b e n d i n g of s t r e a m l i n e s a r o u n d
t h e g r a i n s and t h e v a r y i n g t h i c k n e s s o f p o r e s ( F i g u r e 6.4). volume grows l a r g e r ,
It
The t h i r d t e r m i s c a l l e d t h e d i s p e r s i o n t e r m .
t a k e s i n t o a c c o u n t t h e f l u x due t o m i c r o s c o p i c v a r i a t i o n s i n v e l o c i t y .
I f t h e averaging
inhomogeneities i n t h e p e r m e a b i l i t y c o n t r i b u t e t o t h e varia-
b i l i t y o f v e l o c i t y as w e l l .
On a r e g i o n a l s c a l e ,
we u s u a l l y a v e r a g e o v e r t h e d e p t h
I n t h a t case c l a y lenses and l a y e r s o f d i f f e r e n t p e r m e a b i l i t y
o f the aquifer.
c r e a t e t h e m a j o r d e v i a t i o n s o f l o c a l v e l o c i t i e s f r o m t h e average.
The d i s p e r s i o n
due to a v e r a g i n g o v e r l a r g e r v o l u m e s c o n t a i n i n g i n h o m o g e n e i t i e s o f t h e m a t r i x i s c a l l e d macro-dispersion. s t r a t e d i n f i g u r e s 6.5
B o t h d i s p e r s i o n o n t h e m i c r o - and m a c r o l e v e l a r e i l l u and 6.6.
The a c t i o n o f t h e l o c a l v a r i a b i l i t y o f v e l o c i t y
i s a m i x i n g and s p r e a d i n g o f an i n i t i a l c o n c e n t r a t i o n d i s t r i b u t i o n .
As t h i s p r o -
c e s s u n d e r c e r t a i n c i r c u m s t a n c e s l o o k s v e r y much l i k e a d i f f u s i v e p r o c e s s t h e d i s p e r s i o n t e r m i s w r i t t e n i n f o r m a l a n a l o g y t o t h e d i f f u s i o n t e r m as
The a n a l o g o f t h e d i f f u s i o n c o e f f i c i e n t i s t h e d i s p e r s i o n t e n s o r I D 3 .
It has
t h e f o r m o f a s e c o n d r a n k t e n s o r t o b e a b l e t o accommodate a n i s o t r o p y o f d i s p e r sion.
N o t e t h a t d i s p e r s i o n i s a l w a y s a n i s o t r o p i c - e v e n i f f l o w i s i s o t r o p i c - as
t h e d i s p e r s i o n i n d i r e c t i o n o f f l o w i s u s u a l l y an o r d e r o f m a g n i t u d e l a r g e r t h a n t h a t i n transverse direction. E q u a t i o n (6.6)
i s a working hypothesis which i s o n l y j u s t i f i e d i f t h e v e l o c i t y
v a r i a t i o n s i n t h e a v e r a g i n g v o l u m e e x h i b i t e n o u g h randomness i n t h e s e n s e t h a t any t r a c e r m o l e c u l e t r a n s v e r s i n g t h e a v e r a g i n g volume w i l l e x p e r i e n c e t h e whole spectrum o f v e l o c i t i e s . The a v e r a g i n g p r o c e d u r e c a n b e p e r f o r m e d o v e r a m a c r o s c o p i c v o l u m e a r o u n d a n y p o i n t i n space. f l u x variable. A
j3
= neu3c
-
U s i n g (6.7)
A s s i g n i n g t h e a v e r a g e v a l u e t o t h a t p o i n t , we o b t a i n a c o n t i n u o u s D r o p p i n g a l l a v e r a g i n g b a r s i t i s w r i t t e n as
ne(Dm
+
(6.7)
a s l o c a l p o l l u t a n t f l u x i m p l i e s t h a t t h e p o r o u s medium h a s been
r e p l a c e d b y an e q u i v a l e n t c o n t i n u u m .
194
Average translation of initial concentration distribution (convective transport 1
Actual pathlines
F i g . 6.5:
Schematic r e p r e s e n t a t i o n o f d i s p e r s i o n
Actual transport (convection+ dispersion)
195
POLLUTANT DISTRIBUTION AT TIME t = O permeability k,
I----
-
top
7 -I
direction of flow aquifer
bottom depthaveraged concentration C
distance x
POLLUTANT DISTRIBUTION AT TIME
ti > O
pollutant
aquifer
Y/////////////// ////////////// bottom depthaveraged concentration C
distonce x Fig. 6.6: Schematlc representation o f the macrodisperslon process i n a layered aauifer
196 Our g o a l i s t h e t r a n s p o r t e q u a t i o n .
F i r s t t h e 3-dimensional
i s d e r i v e d f r o m a mass b a l a n c e a r o u n d a n a r b i t r a r y volume.
transport equation
Then t h e two-dimen-
s i o n d l e q u a t i o n i s o b t a i n e d b y a v e r a g i n g f l u x e s o v e r t h e a q u i f e r d e p t h and t a k i n g a mass b a l a n c e o v e r a c o n t r o l v o l u m e e x t e n d i n g f r o m t h e b o t t o m t o t h e t o p o f t h e aquifer. The p o l l u t a n t mass b a l a n c e o v e r an a r b i t r a r y v o l u m e V c o n t a i n e d i n a b o u n d a r y s u r f a c e S demands t h a t t h e t o t a l f l u x o v e r t h e s u r f a c e S and t h e p r o l i f e r a t i o n r a t e o f i n t e r n a l p o l l u t i o n s o u r c e s and s i n k s i n s i d e V b a l a n c e t h e s t o r a g e o f p o l l u t a n t i n s i d e V p e r u n i t time.
( F i g u r e 6.7)
?i3 i s t h e u n i t n o r m a l v e c t o r on t h e s u r f a c e S d i r e c t e d t o w a r d s t h e o u t s i d e o f volume V.
0 i s t h e volume-specific
are counted p o s i t i v e .
source term.
I n p u t s i n t o t h e c o n t r o l volume
By means o f t h e Gauss i n t e g r a l t h e o r e m t h e s u r f a c e i n t e g r a l
i s t r a n s f o r m e d i n t o a volume i n t e g r a l .
Inserting
n t o e q u a t i o n (6.8)
t h e volume i n t e g r a l s i s z e r o , dimensiona
and d e m a n d i n g t h a t t h e sum o f t h e i n t e g r a n d s u n d e r we o b t a i n t h e p a r t a 1 d i f f e r e n t i a l e q u a t i o n o f 3-
transport
(6.10)
I n t h e f o l l o w i n g c h a p t e r s we a r e o n l y i n t e r e s t e d i n t h e t r a n s p o r t e q u a t i o n i n two h o r i z o n t a l dimensions.
F o r m a l l y i t c o u l d be o b t a i n e d b y a s s u m i n g u z = 0 and
a c / a z = 0 and r e p l a c i n g t h e v e r t i c a l b o u n d a r y c o n d i t i o n s b y e x t e r n a l s o u r c e s a n d sinks.
T h i s p r o c e d u r e w o u l d e s s e n t i a l l y assume c o m p l e t e v e r t i c a l h o m o g e n e i t y a n d
instantaneous v e r t i c a l mixing.
The d e p t h - a v e r a g e d e q u a t i o n i s ,
however,
t o a c e r t a i n d e g r e e t o v e r t i c a l l y inhomogeneous a q u i f e r s as w e l l .
applicable
To s e e t h i s we
average t h e f l u x v e c t o r over t h e depth o f t h e a q u i f e r . We c o n s i d e r f i r s t a c o n f i n e d a q u i f e r w i t h b o t t o m and t o p i m p e r v i o u s t o t h e p o l lutant.
I n t h i s case t h e d e p t h average of j z must be zero.
The a v e r a g e s o v e r j x
y i e l d a two dimensional f l u x v e c t o r I n a n a l o g y t o e q u a t i o n s (6.4). the Y l o c a l v e l o c i t y and t h e l o c a l c o n c e n t r a t i o n a t d e p t h z c a n b e decomposed i n t o t h e i r
and j
r e s p e c t i v e d e p t h averages and d e v i a t i o n s f r o m those. c o n v e c t i v e f l u x c a n t h e n be w r i t t e n as
The d e p t h a v e r a g e o f t h e
197
w
Fig. 6.7: Nomenclature for the pollutant mass balance over an arbitrary control volume
Fig.
top of aquifer
bottom of aquifer
6.8: Nomenclature for the pollutant mass balance over a control volume extending from the bottom to the top of an aquifer
198
~
n Uc
_
=
n GC
+
_
6(n G)6c
e
(6.11)
with I;= (ux,uy). The bars indicate a depth average.The second term on the righthand side of (6.11) represents the macrodispersion due to vertical inhomogeneities. It is again modelled by the Fickian approach. Note that the Fickian model is only reasonable if mixing over the depth is fast compared to the time scale of the plume to be modelled. This problem will be discussed further later on. We write (6.12)
We now assume that the depth average of the dispersive part of jx and j can be Y written as (6.13)
While ID' expresses the dispersive effect due to velocity variations from layer to layer, ID" describes the average dispersion within one layer. With an effective porosity which is reasonably constant over depth, the depth averaged dispersive flux i s finally obtained as
All averaging bars are omitted in the following. If the assumption of an impervious top of the aquifer is relaxed, the depth integral of jz may be non-zero. It must in that case balance an external pollutant flux, e.g.
pollutant inflow by infiltration or pollutant outflow by abstraction
of polluted water. Neglecting the diffusive flux over the boundary, we can write
199
(6.15)
q i s t h e f l o w r a t e per u n i t area i n v e r t i c a l d i r e c t i o n from o r t o t h e aquifer.
c. in
i s t h e c o n c e n t r a t i o n o f t h e i n f l o w i n g w a t e r i n t h e case o f i n f i l t r a t i o n o r t h e average c o n c e n t r a t i o n i n t h e a q u i f e r i n t h e case o f a b s t r a c t i o n o f water.
This
i m p l i e s t h e assumption t h a t e i t h e r t h e c o n c e n t r a t i o n a t t h e t o p o f t h e a q u i f e r i s a p p r o x i m a t e l y equal t o t h e depth-averaged perfect wells.
c o n c e n t r a t i o n o r water i s abstracted by
As i n t h e f l o w m o d e l s , q i s t a k e n p o s i t i v e f o r i n f l o w t o t h e a q u i f e r
The t r a n s p o r t e q u a t i o n i n t w o h o r i z o n t a l d i m e n s i o n s i s o b t a i n e d b y t a k i n g a mass b a l a n c e a r o u n d a c o n t r o l v o l u m e e x t e n d i n g f r o m t h e t o p t o t h e b o t t o m o f t h e a q u i f e r as shown i n f i g u r e 6.8.
The h o r i z o n t a l and v e r t i c a l f l u x e s a c r o s s t h e
p e r i p h e r y o f t h e v o l u m e m u s t b a l a n c e i n t e r n a l s i n k s a n d s o u r c e s as w e l l as t h e p o l l u t a n t s t o r a g e i n s i d e t h e volume. (6.16)
= ( n x . n ) w h i c h h a s n o component i n z - d i r e c t i o n . w i t h a normal v e c t o r Y The s u r f a c e i n t e g r a l i s c o n v e r t e d i n t o a v o l u m e i n t e g r a l b y t h e Gauss-theorem.
P-*-
n JdS =
S
I-
V-TdV
V
From ( 6 . 1 6 )
(6.17)
=
A and ( 6 . 1 7 )
t h e d i f f e r e n t i a l form o f t h e t r a n s p o r t equation i s obtained
b y equating t h e integrands.
(6.18)
I n s e r t i n g e q u a t i o n (6.14),
we f i n a l l y o b t a i n t h e B e a r e q u a t i o n .
- qcin
= 0
With m = h - b the equation i s also v a l i d f o r the phreatic aquifer. d i n g i n a d s o r p t i o n and c h e m i c a l r e a c t i o n ,
(6.19)
Before b u i l -
the size of the dispersion coefficients
i s discussed. The d i s p e r s i o n t e n s o r i s i n d i a g o n a l f o r m i f one o f t h e c o o r d i n a t e a x e s i s
200
aligned with the velocity vector. Assuming alignment with the x-axis, it has the form
(6.20)
The dispersion coefficients can be written as the product of the absolute value of
velocity and a length scale, called dispersivity (Scheidegger, 1957).
longitudinal dispersion coefficient: D L transverse dispersion coefficient : D,
=
(6.21)
c1 u
L
= CY u
T
A
with u = / u / . The tensor elements in an arbitrarily oriented coordinate system are obtained by rotation. (6.22)
with u
=
IT/.
I n the application of models numerical values of dispersivities are needed: unfortunately they are rarely known a priori. In laboratory experiments longitudinal dispersivities of 0.01 to 10 cm were found for different types of granular materials. They reflect the microdispersivity due to the interaction of flow and grain. Their size varies with porosity, grain diameter, shape of grains, and grading of grains (Klotz, 1973). Longitudinal dispersivities found in tracer experiments in the field are usually much larger (e.g. Lenda. Zuber. 1970) (Table 6.1).
This is only partly due to
the larger heterogeneity o f natural soil material by itself. The main reason is the influence of small-scale inhomogeneities in permeability of the aquifer: this means the onset of macrodispersion. Even larger dispersivities are found in regional pollution transport. Due to this scale behaviour of dispersion no generally valid values for longitudinal dispersivities can be given. In the literature values between 0.1 m and 500 rn can be found for porous aquifers (Figure 6.9).
201
Wikon.1971
F i g . 6.9:
Wilson3978 Fried, 1975 'Pinder, 1973 Fried.19Ey
ibnikow. 1971
Scale dependence o f l o n g i t u d i n a l d i s p e r s i v i t y ( a f t e r Beims,
There a r e two reasons f o r t h e observed s c a l e dependence. s i z e o f a t r a c e r cloud,
1983)
F i r s t , w i t h growing
l a r g e r and l a r g e r i n h o m o g e n e i t i e s can c o n t r i b u t e t o i t s
As
l o n g as t h e i r e x t e n s i o n i s t o o s m a l l t o be d e s c r i b e d i n
d e t a i l by t h e f l o w f i e l d ,
t h e i r average e f f e c t w i l l show up i n t h e a p p a r e n t d i s -
d i s p e r s i v e spreading. persivity.
An a s y m p t o t i c s t a t e w i t h c o n s t a n t d i s p e r s i v i t i e s i s reached when t h e
s i z e o f t h e plume i s l a r g e a g a i n s t t h e l a r g e s t unknown randomly d i s t r i b u t e d i n homogeneities.
A second reason f o r t h e a p p a r e n t g r o w t h o f d i s p e r s i v i t i e s i s t h e i n a d e q u a t e use o f t h e F i c k i a n d i s p e r s i o n model on s i t u a t i o n s which show t o o l i t t l e randomness t o resemble a d i f f u s i o n process. by t h e l a y e r e d a q u i f e r .
A
s i m p l e c o n c e p t u a l example i s p r e s e n t e d
The m a j o r m a c r o d i s p e r s i o n e f f e c t i n h o r i z o n t a l l y P-dimen-
s i o n a l t r a n s p o r t models o f a l a y e r e d a q u i f e r stems f r o m a v e r a g i n g o v e r t h e depth,
I f an a q u i f e r was composed o f p a r a l l e l h o r i z o n t a l l a y e r s w i t h z e r o mass t r a n s f e r and d i f f e r i n g h o r i z o n t a l v e l o c i t i e s ,
t h e d i f f e r e n t i a l c o n v e c t i o n would l e a d t o a
s p r e a d i n g w h i c h corresponds t o a d i s p e r s i v i t y g r o w i n g t o i n f i n i t y w i t h d i s t a n c e (Mercado. 1967).
A F i c k i a n model which uses a c o n s t a n t d i s p e r s i v i t y does n o t de-
s c r i b e t h i s s i t u a t i o n adequately. exists,
I f transverse mixing i n v e r t i c a l d i r e c t i o n
t h e observed d i s p e r s i v i t y w i l l n o t grow i n d e f i n i t e l y b u t w i l l approach a
f i n i t e a s y m p t o t i c v a l u e a f t e r some f l o w d i s t a n c e when t h e r e i s an e q u i l i b r i u m b e t ween l o n g i t u d i n a l s p r e a d i n g due t o d i f f e r e n t i a l c o n v e c t i o n and t r a n s v e r s e d i s p e r s i v e mixing (Taylor,
1953). The l a r g e r t h e v e r t i c a l m i x i n g , t h e f a s t e r t h e asymp-
202
totic Fickian behaviour of longitudinal macrodispersion will be reached (Gelhar et al., 1979). For very small vertical mixing it may occur that asymptotic behaviour is never reached within a whole groundwater basin. A possibility to account for the changing scale of dispersivity, as long as we
consider a momentary point-like discharge, is the introduction of a time- or distance-dependent aL that grows from a very small value at the origin of the tracer cloud to an asymptotic value (Matheron. De Marsily, 1982). The hor3zontally transverse dispersivity is generally an order of magnitude smaller than the longitudinal one. In laboratory experiments a ratio a,/ a L = 0.1 was found (Klotz, Seiler, 1980). Ratios between 0.01 and 0.3 are reported from field studies (Pickens. Grisak, 1980). The vertically transverse dispersivity may be still smaller. Yet, it effectively governs the longitudinal macrodispersion process as was discussed above. In the following it is assumed that macrodispersion can be described by constant dispersivities a L and aT throughout the aquifer. With a molecular diffusion coefficient of Dm = lo-’ m2/s at a temperature of 10°C. we can neglect the diffusion term against the dispersion term for typical groundwater velocities u in the range of 0.10 m/d to 10 m/d. In stagnant groundwater, however, diffusion may be the only effective transport process. Up to now we considered a conservative pollutant; this means a pollutant which does not decay (e.g. NaC1).
We now introduce a first-order reaction into equation (6.19). This is done via the volume-source/sink-term as reaction constitutes a sink of pollutant mass. In a first order reaction the rate of decay is proportional to the concentration present djcn,)
~-
dt
-
Acne
~
=
a
(6.23)
where X is the decay constant. The same type o f law holds for radioactive decay. I n the case where a pollutant is adsorbed by (or desorbed from) the matrix o f grains, the mass balance must include not only the dissolved pollutant mass but also the adsorbed pollutant mass. While the dissolved concentrations c are usually measured as mass of pollutant per water volume, the adsorbed concentrations c are measured as mass of pollutant per mass of dry matrix material. To compare the two on a geometrical volume basis, a normalization factor is required. The total pollutant mass in a unit volume cut out of the aquifer is given by
AM
=
c n
+
ca(l-ne) p
(6.24)
where p i s the density of the dry matrix material.Including adsorption and first order reaction in equation (6.19) and neglecting molecular diffusion against dis-
203
persion, we obtain
a (mn
c+m( l-ne) pea) at
(6.25)
where CI now denotes volume sources and sinks other than first order reactions. The adsorbed mass appears only in the storage term and in the decay term, the latter being due to the fact that the adsorbed pollutant mass will also decay. To complete equation (6.25) a further equation for the development of the adsorbed concentration is needed. If the adsorption process is fast compared to the typical time scale of flow, we can assume that the adsorbed concentration ca i s always in equilibrium with the dissolved concentration c; this means that ca
=
f(c)
(6.26)
The function f(c) i s called isothermal as it describes the equilibrium at constant temperature. In the simplest case the isothermal i s linear c
=
KC
(6.27)
This is usually the case for very small concentrations (Henry's law). If dissolved and adsorbed concentrations are not in equilibrium, a separate differential equation for ca must be given. The simplest model states that the exchange between adsorbed and dissolved phases is proportional to the difference in respective concentrations.
Equations (6.25) and (6.28) form a system of partial differential equations for c and ca. Up to now, we assumed that the total porosity n i s identical with the effective porosity as we neglected any pollutant mass stored in the part of pore space not accessible to convective transport. In reality, n and ne may differ considerably. The immobile water accounting for the difference n-ne in porosities can also receive pollutant mass by diffusion. The retention effect due to dead-end pore volume can be taken into account in the same way as non-equilibrium adsorption. Writing balance equations for the mobile and the immobile pore water sepa-
204
r a t e l y w i t h an e x c h a n g e t e r m o f t h e f o r m ( 6 . 2 8 ) ,
we o b t a i n ( C o a t s ,
Smith,
1961)
(6.29)
where c1 and c 2 a r e t h e c o n c e n t r a t i o n s i n t h e m o b i l e a n d i m m o b i l e p o r e w a t e r r e spectively.
nl
and n2 a r e t h e f r a c t i o n s o f p o r e s p a c e o c c u p i e d b y m o b i l e and i m -
m o b i l e water. The n o n - e q u i l i b r i u m p o l l u t a n t e x c h a n g e b e t w e e n m a t r i x and w a t e r o r b e t w e e n i m m o b i l e and m o b i l e w a t e r p r o d u c e s a b a c k w a r d t a i l i n g i n p l u m e s w h i c h i s o b s e r v e d
I n t h e f o l l o w i n g we n e g l e c t t h i s e f f e c t as i t i s u s u a l l y n o t a domina-
i n nature.
ting effect.
Further,
it i s n o t e a s i l y separable from dispersion.
m o b i l e w a t e r i s n o t o n l y t h e w a t e r i n a c t u a l dead-end
pores,
Note t h a t im-
but also the water
i n l i t t l e pervious lenses through which groundwater h a r d l y flows. I n t h e f o l l o w i n g we c o n s i d e r o n l y e q u i l i b r i u m a d s o r p t i o n . i s o t h e r m a l i n t o e q u a t i o n (6.25)
with
R
=
1
+
Inserting the linear
we o b t a i n
pi<(l-ne)/ne
Assuming t h a t m,
ne and
R
show o n l y v e r y s m a l l s p a t i a l g r a d i e n t s a n d t h a t
t h e saturated thickness m i s constant over time,
we c a n d i v i d e e q u a t i o n (6.30)
b y t h e f a c t o r m neR.
+ c-(gc) R
- 1 V-(-DVc)
R
+
hc
~
(6.31)
' n 1 mR , in-
We n o t i c e t h a t t h e i n t r o d u c t i o n o f l i n e a r a d s o r p t i o n i s e s s e n t i a l l y e q u i v a l e n t t o a r e t a r d a t i o n o f t h e t r a n s p o r t p r o c e s s as t h e p o r e v e l o c i t y a d i m i n i s h e d v e l o c i t y ;/R
i s replaced by
b o t h i n t h e c o n v e c t i v e t e r m and t h e d i s p e r s i o n t e n s o r .
The d i v i s i o n o f s o u r c e t e r m s b y R r e f l e c t s t h e f a c t t h a t p a r t o f t h e p o l l u t a n t mass i n j e c t e d w i l l be a d s o r b e d o n t o t h e m a t r i x and n o t c o n t r i b u t e t o d i s s o l v e d c o n c e n t r a t i o n a t t h e t i m e o f i n j e c t i o n . The c o e f f i c i e n t R i s c a l l e d r e t a r d a t i o n factor. I n a p h r e a t i c a q u i f e r t h e s a t u r a t e d t h i c k n e s s m=h-b may v a r y w i t h t i m e . f o r e t h e s t o r a g e t e r m i n e q u a t i o n (6.30)
m u s t b e w r i t t e n as
There-
205
a ( nemRc )
___-
at
-
ac n mRe at
+
(6.32)
am
n Rce a
t
The s e c o n d t e r m o n t h e r i g h t - h a n d s i d e d e s c r i b e s s t o r a g e o f p o l l u t a n t mass assoc i a t e d w i t h movement o f t h e w a t e r t a b l e . The v o l u m e s o u r c e s 0 a r e c o n v e r t e d i n t o a r e a l s o u r c e s b y t h e d e f i n i t i o n
u
(6.33)
Sint/m
=
T h e y c o v e r a l l p o l l u t a n t i n p u t s t h a t a r e n o t c o u p l e d t o an i n p u t o f w a t e r . d i s s o l u t i o n o f substances i m m i s c i b l e w i t h water,
The
such as c h l o r o h y d r o c a r b o n s w h i c h
h a v e p e n e t r a t e d t h e s a t u r a t e d z o n e i n u n d i s s o l v e d form,
c a n be d e s c r i b e d b y t h i s
term. The c o n v e c t i o n t e r m c a n be t r a n s f o r m e d b y means o f t h e c o n t i n u i t y e q u a t i o n o f water f l o w i n t o
We f i n a l l y end u p w i t h t h e common e q u a t i o n (e.9. ac
+
ri-
R'VC
-
-
ID-
V * ( - VRc )
+
hc
-
* e
'lnt
(C&
- __ n mR
-
Konikow,
__ nemR
Grove,
(S%-n
at
R -a )m= e
at
1977) 0 (6.35)
f o r b o t h t h e p h r e a t i c and c o n f i n e d a q u i f e r . The f o r m o f t h e t r a n s p o r t e q u a t i o n u s e d i n t h e f o l l o w i n g c h a p t e r s assumes steady-state
flow.
Therefore,
a l l t e r m s i n (6.35)
containing time derivatives o f
h and m a r e o m i t t e d . I n some f i e l d a p p l i c a t i o n s as w e l l as i n s o i l c o l u m n e x p e r i m e n t s t h e one-dimens i o n a l form o f t h e t r a n s p o r t equation i s applicable. b y assuming t h a t u tion.
Also,
=
Y l a t e r a l i n f l o w s and i n t e r n a l s o u r c e s a r e o m i t t e d .
ac + ~u a c - -(--) a DLac a t R ~ X ax R a x
+ hc
=
(6.36)
0
The t r a n s p o r t e q u a t i o n i s a s e c o n d - o r d e r
p a r t i a l d i f f e r e n t i a l equation which
r e q u i r e s i n i t i a l c o n d i t i o n s and b o u n d a r y c o n d i t i o n s . by t h e c o n c e n t r a t i o n d i s t r i b u t i o n c(x.y.tO) As i n t h e c a s e o f t h e f l o w e q u a t i o n , conditions.
We o b t a i n i t f r o m ( 6 . 3 5 )
0 and t h a t t h e c o n c e n t r a t i o n c h a s n o g r a d i e n t s i n y - d i r e c -
I n i t i a l conditions are given
a t t h e s t a r t i n g t i m e to o f s i m u l a t i o n .
t h e r e a r e t h r e e p o s s i b l e t y p e s o f boundary
F i r s t k i n d ( D i r i c h l e t ) b o u n d a r y c o n d i t i o n s s p e c i f y p r e s c r i b e d con-
c e n t r a t i o n s on t h e b o u n d a r y .
I n n e r b o u n d a r i e s o f t h e f i r s t k i n d can d e s c r i b e p o l -
206
lution sources. Second kind (Neumann) boundary conditions specify the concentration gradient normal to a boundary. This means they prescribe the dispersive flux. The convective flux at boundaries cannot be prescribed by second-kind boundary conditions. It is a result of the interplay of concentrations on the boundary, as given, for example, by first-kind boundary conditions, and velocities. At impervious boundaries the flow model has to yield zero water flow across the boundary. Then we need only prescribe a vanishing dispersive flux, i.e. &/an = 0. Prescribing the total flux, i. e. the sum of convective and dispersive flux, on the boundary corresponds to a boundary condition of the third kind (Cauchy boundary condition), as the total flux is a linear combination of the concentration on the boundary, c, and the normal derivative of the concentration on the boundary, adan. (6.37) If the dispersive flux is small compared to the convective flux, the third-kind boundary condition approaches a first-kind boundary condition. I n the following, we consider first-kind and second-kind boundary conditions only. The combined action of convection, diffusion-dispersion, adsorption, and decay is shown for one-dimensional transport i n figure 6.10. Initially, a homogeneous distribution of dissolved concentration between x=O and x=a is given. While convection alone causes a spatial translation of the concentration distribution, diffusion-dispersion leads to a spreading. When switching on linear adsorption, the dissolved concentration is diminished if the pollutant mass in the distribution remains unchanged, as part of the dissolved pollutant mass will go to the matrix, For convenience of comparison we leave the dissolved pollutant concentrations unchanged, adding the adsorbed pollutant mass to the system. Adsorption causes a delay of the transport process, the spreading included. While convection and diffusion-dispersion leave the pollutant mass unchanged, decay reduces the pollutant mass and consequently the area under the distribution.
207
POLLUTANT DISTRIBUTION AT TIME
W
X
Concentration c
t =O
-
direction of flow
____e
xlo
distance x
POLLUTANT DISTRIBUTION AT TIME
t, > 0
X=Ufr
X
X
CA
*
c 4 Effect of convection, dispersion, adsorption. degradation
F i g . 6.10:
Schematic d e s c r i p t i o n o f t h e e f f e c t s o f c o n v e c t i o n , d i s p e r s i o n , a d s o r p t i o n , and chemical d e g r a d a t i o n on p o l l u t a n t t r a n s p o r t
Chapter 6
REGIONAL POLLUTANT TRANSPORT MODELS SCOPE OF APPLICATION Groundwater pollution due t o human activities occurs in a number of ways. Several situations are shown i n figure 6.1. Polluted water may infiltrate into the aquifer from a polluted surface water body, from leaking waste-water pipes, ponds or cesspools (bacterial pollution). The pollutant may also be leached from the top soil by rain and be carried down into the saturat _1 zone by seepage
6.1
(nitrates, pesticides).
Similarly, seepage through landfills and waste deposits
is a source of pollution. Pollutants may enter the soil in a form immiscible with water (e.g. oil, chlorinated hydrocarbons). Their gradual dissolution by seepage water or the passing groundwater flow causes groundwater pollution. While the transport of pollutants in the unsaturated zone is essentially limited to vertical transport between groundlevel and the top of the saturated zone, long-range spreading of pollutants is only possible in the saturated zone where dissolved pollutants are carried along by the prevailing flow, mainly in horizontal directions. The distance that can be covered within a certain time depends on the velocity of flow and the persistance of the pollutant. We consider the often complicated processes in the unsaturated zone on the whole as a source term for the transport in the saturated zone and limit the further discussion to the saturated zone. Depending on the density of solute, the pollutant may influence the flow field (hydrodynamically active solute). We only consider hydrodynamically inactive solutes, this means pollutant concentrations that are so small that density induced flow can be ignored. We further restrict our discussion to regional transport which means the horizontal scale of transport is much larger than the thickness of the aquifer. I n these cases the flow field can essentially be considered as horizontally two-dimensional. In the following chapters we assume that regional pollutant transport can also be described by horizontally two-dimensional models as long as vertical averages of concentrations are considered. In the vicinity of infiltrations or in strongly stratified aquifers, three-dimensional aspects transport may, however, be important. Such situations are shown i n figure 6.2. Any transport model requires as an essential input the velocity field of flow. We assume that the flow field is either known a priori or is modelled parallel of
to transport by a flow model. The flow model yields heads, from which specific flowrates (filter velocities) are calculated by means of the Darcy-law. The transport model requires pore velocities. These are obtained from specific flowrates by division by the effective porosity ne' The derivations given below follow the work of Bear.
preci pi tat ion
P
Fig.
6.1: Sources
of
groundwater pollution
00
(0
190
1 1 1 1 1 1 1 waste deposit
surface table
streamline
concentration
a) Strong vertical flow component
b) Both vertical velocity distribution and concentration distribution are non-homogeneous
Fig.
6.2: Examples o f situations where depth-averaging is not reasonable
191
PHENOMENA TO BE CONSIDERED AND BASIC EQUATIONS
6.2
We first discuss the behaviour of an ideal tracer which experiences neither adsorption nor chemical transformation. On the microscopic level of a single pore space, there are then only two basic processes of transport. These are the convective transport by the water flow and molecular diffusion. The convective pollutant flow through an infinitesimal area element dA at a point T3 inside a pore is given by
- ( U ~ . ~ ~ , U ~n3) =. (nx,ny,nz),
G3.E3 c dA
u c dA
=
with
G3
=
- Jconv' "3 dA
(6.1) =
(x.y.2).
-
and jconv= u3c
z3
un is the velocity component normal to the area element, the 3-dimensional velocity vector at location y3, and G3 the unit normal vector of the area element dA. (Figure 6 . 3 ) The diffusive pollutant flow through the same area element is given by Fick's law. It is proportional to the negative gradient of concentration normal to the area element dA.
The total flux at any point in space is given by
For
all practical purposes we have to look at macroscopic volumes when taking
pollutant balances. Therefore we must calculate averages of the flux vector ( 6 . 3 ) over volumes comprising many pore spaces. Averages can be taken over a geometrical volume or over the portion of pore space contained in it, which is effectively taking part in flow. Averages over total volume are obtained from averages over the effective pore space by multiplication by the effective porosity ne'
I n the averaging process it is useful to decompose local concentration and velocity into a macroscopic average over the effective pore space and the local variations from that average.
c = c +6c
192
F i g . 6.3: Nomenclature for the definition o f infinitesimal convective
pollutant flux
variable poresize
velocity profile within pore
bending of streamlines around grains
Fig. 6.4: Causes of microscopic variations i n velocity
193 I n s e r t i n g (6.4)
i n t o (6.3).
m u l t i p l y i n g b y t h e e f f e c t i v e p o r o s i t y and a v e r a g i n g
-V
we o b t a i n t h e v o l u m e t r i c a v e r a g e f l u x j 3
(6.5) Terms
-
r 3 6 c and
s i d e o f (6.5)
c6G3
vanish d u r i n g averaging.
The f i r s t t e r m o n t h e r i g h t - h a n d
i s c a l l e d t h e c o n v e c t i o n ( o r a d v e c t i o n ) term.
the molecular diffusion
term.
The s e c o n d t e r m i s
caused by t h e v e l o c i t y p r o f i l e w i t h i n a pore,
These a r e
t h e b e n d i n g of s t r e a m l i n e s a r o u n d
t h e g r a i n s and t h e v a r y i n g t h i c k n e s s o f p o r e s ( F i g u r e 6.4). volume grows l a r g e r ,
It
The t h i r d t e r m i s c a l l e d t h e d i s p e r s i o n t e r m .
t a k e s i n t o a c c o u n t t h e f l u x due t o m i c r o s c o p i c v a r i a t i o n s i n v e l o c i t y .
I f t h e averaging
inhomogeneities i n t h e p e r m e a b i l i t y c o n t r i b u t e t o t h e varia-
b i l i t y o f v e l o c i t y as w e l l .
On a r e g i o n a l s c a l e ,
we u s u a l l y a v e r a g e o v e r t h e d e p t h
I n t h a t case c l a y lenses and l a y e r s o f d i f f e r e n t p e r m e a b i l i t y
o f the aquifer.
c r e a t e t h e m a j o r d e v i a t i o n s o f l o c a l v e l o c i t i e s f r o m t h e average.
The d i s p e r s i o n
due to a v e r a g i n g o v e r l a r g e r v o l u m e s c o n t a i n i n g i n h o m o g e n e i t i e s o f t h e m a t r i x i s c a l l e d macro-dispersion. s t r a t e d i n f i g u r e s 6.5
B o t h d i s p e r s i o n o n t h e m i c r o - and m a c r o l e v e l a r e i l l u and 6.6.
The a c t i o n o f t h e l o c a l v a r i a b i l i t y o f v e l o c i t y
i s a m i x i n g and s p r e a d i n g o f an i n i t i a l c o n c e n t r a t i o n d i s t r i b u t i o n .
As t h i s p r o -
c e s s u n d e r c e r t a i n c i r c u m s t a n c e s l o o k s v e r y much l i k e a d i f f u s i v e p r o c e s s t h e d i s p e r s i o n t e r m i s w r i t t e n i n f o r m a l a n a l o g y t o t h e d i f f u s i o n t e r m as
The a n a l o g o f t h e d i f f u s i o n c o e f f i c i e n t i s t h e d i s p e r s i o n t e n s o r I D 3 .
It has
t h e f o r m o f a s e c o n d r a n k t e n s o r t o b e a b l e t o accommodate a n i s o t r o p y o f d i s p e r sion.
N o t e t h a t d i s p e r s i o n i s a l w a y s a n i s o t r o p i c - e v e n i f f l o w i s i s o t r o p i c - as
t h e d i s p e r s i o n i n d i r e c t i o n o f f l o w i s u s u a l l y an o r d e r o f m a g n i t u d e l a r g e r t h a n t h a t i n transverse direction. E q u a t i o n (6.6)
i s a working hypothesis which i s o n l y j u s t i f i e d i f t h e v e l o c i t y
v a r i a t i o n s i n t h e a v e r a g i n g v o l u m e e x h i b i t e n o u g h randomness i n t h e s e n s e t h a t any t r a c e r m o l e c u l e t r a n s v e r s i n g t h e a v e r a g i n g volume w i l l e x p e r i e n c e t h e whole spectrum o f v e l o c i t i e s . The a v e r a g i n g p r o c e d u r e c a n b e p e r f o r m e d o v e r a m a c r o s c o p i c v o l u m e a r o u n d a n y p o i n t i n space. f l u x variable. A
j3
= neu3c
-
U s i n g (6.7)
A s s i g n i n g t h e a v e r a g e v a l u e t o t h a t p o i n t , we o b t a i n a c o n t i n u o u s D r o p p i n g a l l a v e r a g i n g b a r s i t i s w r i t t e n as
ne(Dm
+
(6.7)
a s l o c a l p o l l u t a n t f l u x i m p l i e s t h a t t h e p o r o u s medium h a s been
r e p l a c e d b y an e q u i v a l e n t c o n t i n u u m .
194
Average translation of initial concentration distribution (convective transport 1
Actual pathlines
F i g . 6.5:
Schematic r e p r e s e n t a t i o n o f d i s p e r s i o n
Actual transport (convection+ dispersion)
195
POLLUTANT DISTRIBUTION AT TIME t = O permeability k,
I----
-
top
7 -I
direction of flow aquifer
bottom depthaveraged concentration C
distance x
POLLUTANT DISTRIBUTION AT TIME
ti > O
pollutant
aquifer
Y/////////////// ////////////// bottom depthaveraged concentration C
distonce x Fig. 6.6: Schematlc representation o f the macrodisperslon process i n a layered aauifer
196 Our g o a l i s t h e t r a n s p o r t e q u a t i o n .
F i r s t t h e 3-dimensional
i s d e r i v e d f r o m a mass b a l a n c e a r o u n d a n a r b i t r a r y volume.
transport equation
Then t h e two-dimen-
s i o n d l e q u a t i o n i s o b t a i n e d b y a v e r a g i n g f l u x e s o v e r t h e a q u i f e r d e p t h and t a k i n g a mass b a l a n c e o v e r a c o n t r o l v o l u m e e x t e n d i n g f r o m t h e b o t t o m t o t h e t o p o f t h e aquifer. The p o l l u t a n t mass b a l a n c e o v e r an a r b i t r a r y v o l u m e V c o n t a i n e d i n a b o u n d a r y s u r f a c e S demands t h a t t h e t o t a l f l u x o v e r t h e s u r f a c e S and t h e p r o l i f e r a t i o n r a t e o f i n t e r n a l p o l l u t i o n s o u r c e s and s i n k s i n s i d e V b a l a n c e t h e s t o r a g e o f p o l l u t a n t i n s i d e V p e r u n i t time.
( F i g u r e 6.7)
?i3 i s t h e u n i t n o r m a l v e c t o r on t h e s u r f a c e S d i r e c t e d t o w a r d s t h e o u t s i d e o f volume V.
0 i s t h e volume-specific
are counted p o s i t i v e .
source term.
I n p u t s i n t o t h e c o n t r o l volume
By means o f t h e Gauss i n t e g r a l t h e o r e m t h e s u r f a c e i n t e g r a l
i s t r a n s f o r m e d i n t o a volume i n t e g r a l .
Inserting
n t o e q u a t i o n (6.8)
t h e volume i n t e g r a l s i s z e r o , dimensiona
and d e m a n d i n g t h a t t h e sum o f t h e i n t e g r a n d s u n d e r we o b t a i n t h e p a r t a 1 d i f f e r e n t i a l e q u a t i o n o f 3-
transport
(6.10)
I n t h e f o l l o w i n g c h a p t e r s we a r e o n l y i n t e r e s t e d i n t h e t r a n s p o r t e q u a t i o n i n two h o r i z o n t a l dimensions.
F o r m a l l y i t c o u l d be o b t a i n e d b y a s s u m i n g u z = 0 and
a c / a z = 0 and r e p l a c i n g t h e v e r t i c a l b o u n d a r y c o n d i t i o n s b y e x t e r n a l s o u r c e s a n d sinks.
T h i s p r o c e d u r e w o u l d e s s e n t i a l l y assume c o m p l e t e v e r t i c a l h o m o g e n e i t y a n d
instantaneous v e r t i c a l mixing.
The d e p t h - a v e r a g e d e q u a t i o n i s ,
however,
t o a c e r t a i n d e g r e e t o v e r t i c a l l y inhomogeneous a q u i f e r s as w e l l .
applicable
To s e e t h i s we
average t h e f l u x v e c t o r over t h e depth o f t h e a q u i f e r . We c o n s i d e r f i r s t a c o n f i n e d a q u i f e r w i t h b o t t o m and t o p i m p e r v i o u s t o t h e p o l lutant.
I n t h i s case t h e d e p t h average of j z must be zero.
The a v e r a g e s o v e r j x
y i e l d a two dimensional f l u x v e c t o r I n a n a l o g y t o e q u a t i o n s (6.4). the Y l o c a l v e l o c i t y and t h e l o c a l c o n c e n t r a t i o n a t d e p t h z c a n b e decomposed i n t o t h e i r
and j
r e s p e c t i v e d e p t h averages and d e v i a t i o n s f r o m those. c o n v e c t i v e f l u x c a n t h e n be w r i t t e n as
The d e p t h a v e r a g e o f t h e
197
w
Fig. 6.7: Nomenclature for the pollutant mass balance over an arbitrary control volume
Fig.
top of aquifer
bottom of aquifer
6.8: Nomenclature for the pollutant mass balance over a control volume extending from the bottom to the top of an aquifer
198
~
n Uc
_
=
n GC
+
_
6(n G)6c
e
(6.11)
with I;= (ux,uy). The bars indicate a depth average.The second term on the righthand side of (6.11) represents the macrodispersion due to vertical inhomogeneities. It is again modelled by the Fickian approach. Note that the Fickian model is only reasonable if mixing over the depth is fast compared to the time scale of the plume to be modelled. This problem will be discussed further later on. We write (6.12)
We now assume that the depth average of the dispersive part of jx and j can be Y written as (6.13)
While ID' expresses the dispersive effect due to velocity variations from layer to layer, ID" describes the average dispersion within one layer. With an effective porosity which is reasonably constant over depth, the depth averaged dispersive flux i s finally obtained as
All averaging bars are omitted in the following. If the assumption of an impervious top of the aquifer is relaxed, the depth integral of jz may be non-zero. It must in that case balance an external pollutant flux, e.g.
pollutant inflow by infiltration or pollutant outflow by abstraction
of polluted water. Neglecting the diffusive flux over the boundary, we can write
199
(6.15)
q i s t h e f l o w r a t e per u n i t area i n v e r t i c a l d i r e c t i o n from o r t o t h e aquifer.
c. in
i s t h e c o n c e n t r a t i o n o f t h e i n f l o w i n g w a t e r i n t h e case o f i n f i l t r a t i o n o r t h e average c o n c e n t r a t i o n i n t h e a q u i f e r i n t h e case o f a b s t r a c t i o n o f water.
This
i m p l i e s t h e assumption t h a t e i t h e r t h e c o n c e n t r a t i o n a t t h e t o p o f t h e a q u i f e r i s a p p r o x i m a t e l y equal t o t h e depth-averaged perfect wells.
c o n c e n t r a t i o n o r water i s abstracted by
As i n t h e f l o w m o d e l s , q i s t a k e n p o s i t i v e f o r i n f l o w t o t h e a q u i f e r
The t r a n s p o r t e q u a t i o n i n t w o h o r i z o n t a l d i m e n s i o n s i s o b t a i n e d b y t a k i n g a mass b a l a n c e a r o u n d a c o n t r o l v o l u m e e x t e n d i n g f r o m t h e t o p t o t h e b o t t o m o f t h e a q u i f e r as shown i n f i g u r e 6.8.
The h o r i z o n t a l and v e r t i c a l f l u x e s a c r o s s t h e
p e r i p h e r y o f t h e v o l u m e m u s t b a l a n c e i n t e r n a l s i n k s a n d s o u r c e s as w e l l as t h e p o l l u t a n t s t o r a g e i n s i d e t h e volume. (6.16)
= ( n x . n ) w h i c h h a s n o component i n z - d i r e c t i o n . w i t h a normal v e c t o r Y The s u r f a c e i n t e g r a l i s c o n v e r t e d i n t o a v o l u m e i n t e g r a l b y t h e Gauss-theorem.
P-*-
n JdS =
S
I-
V-TdV
V
From ( 6 . 1 6 )
(6.17)
=
A and ( 6 . 1 7 )
t h e d i f f e r e n t i a l form o f t h e t r a n s p o r t equation i s obtained
b y equating t h e integrands.
(6.18)
I n s e r t i n g e q u a t i o n (6.14),
we f i n a l l y o b t a i n t h e B e a r e q u a t i o n .
- qcin
= 0
With m = h - b the equation i s also v a l i d f o r the phreatic aquifer. d i n g i n a d s o r p t i o n and c h e m i c a l r e a c t i o n ,
(6.19)
Before b u i l -
the size of the dispersion coefficients
i s discussed. The d i s p e r s i o n t e n s o r i s i n d i a g o n a l f o r m i f one o f t h e c o o r d i n a t e a x e s i s
200
aligned with the velocity vector. Assuming alignment with the x-axis, it has the form
(6.20)
The dispersion coefficients can be written as the product of the absolute value of
velocity and a length scale, called dispersivity (Scheidegger, 1957).
longitudinal dispersion coefficient: D L transverse dispersion coefficient : D,
=
(6.21)
c1 u
L
= CY u
T
A
with u = / u / . The tensor elements in an arbitrarily oriented coordinate system are obtained by rotation. (6.22)
with u
=
IT/.
I n the application of models numerical values of dispersivities are needed: unfortunately they are rarely known a priori. In laboratory experiments longitudinal dispersivities of 0.01 to 10 cm were found for different types of granular materials. They reflect the microdispersivity due to the interaction of flow and grain. Their size varies with porosity, grain diameter, shape of grains, and grading of grains (Klotz, 1973). Longitudinal dispersivities found in tracer experiments in the field are usually much larger (e.g. Lenda. Zuber. 1970) (Table 6.1).
This is only partly due to
the larger heterogeneity o f natural soil material by itself. The main reason is the influence of small-scale inhomogeneities in permeability of the aquifer: this means the onset of macrodispersion. Even larger dispersivities are found in regional pollution transport. Due to this scale behaviour of dispersion no generally valid values for longitudinal dispersivities can be given. In the literature values between 0.1 m and 500 rn can be found for porous aquifers (Figure 6.9).
201
Wikon.1971
F i g . 6.9:
Wilson3978 Fried, 1975 'Pinder, 1973 Fried.19Ey
ibnikow. 1971
Scale dependence o f l o n g i t u d i n a l d i s p e r s i v i t y ( a f t e r Beims,
There a r e two reasons f o r t h e observed s c a l e dependence. s i z e o f a t r a c e r cloud,
1983)
F i r s t , w i t h growing
l a r g e r and l a r g e r i n h o m o g e n e i t i e s can c o n t r i b u t e t o i t s
As
l o n g as t h e i r e x t e n s i o n i s t o o s m a l l t o be d e s c r i b e d i n
d e t a i l by t h e f l o w f i e l d ,
t h e i r average e f f e c t w i l l show up i n t h e a p p a r e n t d i s -
d i s p e r s i v e spreading. persivity.
An a s y m p t o t i c s t a t e w i t h c o n s t a n t d i s p e r s i v i t i e s i s reached when t h e
s i z e o f t h e plume i s l a r g e a g a i n s t t h e l a r g e s t unknown randomly d i s t r i b u t e d i n homogeneities.
A second reason f o r t h e a p p a r e n t g r o w t h o f d i s p e r s i v i t i e s i s t h e i n a d e q u a t e use o f t h e F i c k i a n d i s p e r s i o n model on s i t u a t i o n s which show t o o l i t t l e randomness t o resemble a d i f f u s i o n process. by t h e l a y e r e d a q u i f e r .
A
s i m p l e c o n c e p t u a l example i s p r e s e n t e d
The m a j o r m a c r o d i s p e r s i o n e f f e c t i n h o r i z o n t a l l y P-dimen-
s i o n a l t r a n s p o r t models o f a l a y e r e d a q u i f e r stems f r o m a v e r a g i n g o v e r t h e depth,
I f an a q u i f e r was composed o f p a r a l l e l h o r i z o n t a l l a y e r s w i t h z e r o mass t r a n s f e r and d i f f e r i n g h o r i z o n t a l v e l o c i t i e s ,
t h e d i f f e r e n t i a l c o n v e c t i o n would l e a d t o a
s p r e a d i n g w h i c h corresponds t o a d i s p e r s i v i t y g r o w i n g t o i n f i n i t y w i t h d i s t a n c e (Mercado. 1967).
A F i c k i a n model which uses a c o n s t a n t d i s p e r s i v i t y does n o t de-
s c r i b e t h i s s i t u a t i o n adequately. exists,
I f transverse mixing i n v e r t i c a l d i r e c t i o n
t h e observed d i s p e r s i v i t y w i l l n o t grow i n d e f i n i t e l y b u t w i l l approach a
f i n i t e a s y m p t o t i c v a l u e a f t e r some f l o w d i s t a n c e when t h e r e i s an e q u i l i b r i u m b e t ween l o n g i t u d i n a l s p r e a d i n g due t o d i f f e r e n t i a l c o n v e c t i o n and t r a n s v e r s e d i s p e r s i v e mixing (Taylor,
1953). The l a r g e r t h e v e r t i c a l m i x i n g , t h e f a s t e r t h e asymp-
202
totic Fickian behaviour of longitudinal macrodispersion will be reached (Gelhar et al., 1979). For very small vertical mixing it may occur that asymptotic behaviour is never reached within a whole groundwater basin. A possibility to account for the changing scale of dispersivity, as long as we
consider a momentary point-like discharge, is the introduction of a time- or distance-dependent aL that grows from a very small value at the origin of the tracer cloud to an asymptotic value (Matheron. De Marsily, 1982). The hor3zontally transverse dispersivity is generally an order of magnitude smaller than the longitudinal one. In laboratory experiments a ratio a,/ a L = 0.1 was found (Klotz, Seiler, 1980). Ratios between 0.01 and 0.3 are reported from field studies (Pickens. Grisak, 1980). The vertically transverse dispersivity may be still smaller. Yet, it effectively governs the longitudinal macrodispersion process as was discussed above. In the following it is assumed that macrodispersion can be described by constant dispersivities a L and aT throughout the aquifer. With a molecular diffusion coefficient of Dm = lo-’ m2/s at a temperature of 10°C. we can neglect the diffusion term against the dispersion term for typical groundwater velocities u in the range of 0.10 m/d to 10 m/d. In stagnant groundwater, however, diffusion may be the only effective transport process. Up to now we considered a conservative pollutant; this means a pollutant which does not decay (e.g. NaC1).
We now introduce a first-order reaction into equation (6.19). This is done via the volume-source/sink-term as reaction constitutes a sink of pollutant mass. In a first order reaction the rate of decay is proportional to the concentration present djcn,)
~-
dt
-
Acne
~
=
a
(6.23)
where X is the decay constant. The same type o f law holds for radioactive decay. I n the case where a pollutant is adsorbed by (or desorbed from) the matrix o f grains, the mass balance must include not only the dissolved pollutant mass but also the adsorbed pollutant mass. While the dissolved concentrations c are usually measured as mass of pollutant per water volume, the adsorbed concentrations c are measured as mass of pollutant per mass of dry matrix material. To compare the two on a geometrical volume basis, a normalization factor is required. The total pollutant mass in a unit volume cut out of the aquifer is given by
AM
=
c n
+
ca(l-ne) p
(6.24)
where p i s the density of the dry matrix material.Including adsorption and first order reaction in equation (6.19) and neglecting molecular diffusion against dis-
203
persion, we obtain
a (mn
c+m( l-ne) pea) at
(6.25)
where CI now denotes volume sources and sinks other than first order reactions. The adsorbed mass appears only in the storage term and in the decay term, the latter being due to the fact that the adsorbed pollutant mass will also decay. To complete equation (6.25) a further equation for the development of the adsorbed concentration is needed. If the adsorption process is fast compared to the typical time scale of flow, we can assume that the adsorbed concentration ca i s always in equilibrium with the dissolved concentration c; this means that ca
=
f(c)
(6.26)
The function f(c) i s called isothermal as it describes the equilibrium at constant temperature. In the simplest case the isothermal i s linear c
=
KC
(6.27)
This is usually the case for very small concentrations (Henry's law). If dissolved and adsorbed concentrations are not in equilibrium, a separate differential equation for ca must be given. The simplest model states that the exchange between adsorbed and dissolved phases is proportional to the difference in respective concentrations.
Equations (6.25) and (6.28) form a system of partial differential equations for c and ca. Up to now, we assumed that the total porosity n i s identical with the effective porosity as we neglected any pollutant mass stored in the part of pore space not accessible to convective transport. In reality, n and ne may differ considerably. The immobile water accounting for the difference n-ne in porosities can also receive pollutant mass by diffusion. The retention effect due to dead-end pore volume can be taken into account in the same way as non-equilibrium adsorption. Writing balance equations for the mobile and the immobile pore water sepa-
204
r a t e l y w i t h an e x c h a n g e t e r m o f t h e f o r m ( 6 . 2 8 ) ,
we o b t a i n ( C o a t s ,
Smith,
1961)
(6.29)
where c1 and c 2 a r e t h e c o n c e n t r a t i o n s i n t h e m o b i l e a n d i m m o b i l e p o r e w a t e r r e spectively.
nl
and n2 a r e t h e f r a c t i o n s o f p o r e s p a c e o c c u p i e d b y m o b i l e and i m -
m o b i l e water. The n o n - e q u i l i b r i u m p o l l u t a n t e x c h a n g e b e t w e e n m a t r i x and w a t e r o r b e t w e e n i m m o b i l e and m o b i l e w a t e r p r o d u c e s a b a c k w a r d t a i l i n g i n p l u m e s w h i c h i s o b s e r v e d
I n t h e f o l l o w i n g we n e g l e c t t h i s e f f e c t as i t i s u s u a l l y n o t a domina-
i n nature.
ting effect.
Further,
it i s n o t e a s i l y separable from dispersion.
m o b i l e w a t e r i s n o t o n l y t h e w a t e r i n a c t u a l dead-end
pores,
Note t h a t im-
but also the water
i n l i t t l e pervious lenses through which groundwater h a r d l y flows. I n t h e f o l l o w i n g we c o n s i d e r o n l y e q u i l i b r i u m a d s o r p t i o n . i s o t h e r m a l i n t o e q u a t i o n (6.25)
with
R
=
1
+
Inserting the linear
we o b t a i n
pi<(l-ne)/ne
Assuming t h a t m,
ne and
R
show o n l y v e r y s m a l l s p a t i a l g r a d i e n t s a n d t h a t
t h e saturated thickness m i s constant over time,
we c a n d i v i d e e q u a t i o n (6.30)
b y t h e f a c t o r m neR.
+ c-(gc) R
- 1 V-(-DVc)
R
+
hc
~
(6.31)
' n 1 mR , in-
We n o t i c e t h a t t h e i n t r o d u c t i o n o f l i n e a r a d s o r p t i o n i s e s s e n t i a l l y e q u i v a l e n t t o a r e t a r d a t i o n o f t h e t r a n s p o r t p r o c e s s as t h e p o r e v e l o c i t y a d i m i n i s h e d v e l o c i t y ;/R
i s replaced by
b o t h i n t h e c o n v e c t i v e t e r m and t h e d i s p e r s i o n t e n s o r .
The d i v i s i o n o f s o u r c e t e r m s b y R r e f l e c t s t h e f a c t t h a t p a r t o f t h e p o l l u t a n t mass i n j e c t e d w i l l be a d s o r b e d o n t o t h e m a t r i x and n o t c o n t r i b u t e t o d i s s o l v e d c o n c e n t r a t i o n a t t h e t i m e o f i n j e c t i o n . The c o e f f i c i e n t R i s c a l l e d r e t a r d a t i o n factor. I n a p h r e a t i c a q u i f e r t h e s a t u r a t e d t h i c k n e s s m=h-b may v a r y w i t h t i m e . f o r e t h e s t o r a g e t e r m i n e q u a t i o n (6.30)
m u s t b e w r i t t e n as
There-
205
a ( nemRc )
___-
at
-
ac n mRe at
+
(6.32)
am
n Rce a
t
The s e c o n d t e r m o n t h e r i g h t - h a n d s i d e d e s c r i b e s s t o r a g e o f p o l l u t a n t mass assoc i a t e d w i t h movement o f t h e w a t e r t a b l e . The v o l u m e s o u r c e s 0 a r e c o n v e r t e d i n t o a r e a l s o u r c e s b y t h e d e f i n i t i o n
u
(6.33)
Sint/m
=
T h e y c o v e r a l l p o l l u t a n t i n p u t s t h a t a r e n o t c o u p l e d t o an i n p u t o f w a t e r . d i s s o l u t i o n o f substances i m m i s c i b l e w i t h water,
The
such as c h l o r o h y d r o c a r b o n s w h i c h
h a v e p e n e t r a t e d t h e s a t u r a t e d z o n e i n u n d i s s o l v e d form,
c a n be d e s c r i b e d b y t h i s
term. The c o n v e c t i o n t e r m c a n be t r a n s f o r m e d b y means o f t h e c o n t i n u i t y e q u a t i o n o f water f l o w i n t o
We f i n a l l y end u p w i t h t h e common e q u a t i o n (e.9. ac
+
ri-
R'VC
-
-
ID-
V * ( - VRc )
+
hc
-
* e
'lnt
(C&
- __ n mR
-
Konikow,
__ nemR
Grove,
(S%-n
at
R -a )m= e
at
1977) 0 (6.35)
f o r b o t h t h e p h r e a t i c and c o n f i n e d a q u i f e r . The f o r m o f t h e t r a n s p o r t e q u a t i o n u s e d i n t h e f o l l o w i n g c h a p t e r s assumes steady-state
flow.
Therefore,
a l l t e r m s i n (6.35)
containing time derivatives o f
h and m a r e o m i t t e d . I n some f i e l d a p p l i c a t i o n s as w e l l as i n s o i l c o l u m n e x p e r i m e n t s t h e one-dimens i o n a l form o f t h e t r a n s p o r t equation i s applicable. b y assuming t h a t u tion.
Also,
=
Y l a t e r a l i n f l o w s and i n t e r n a l s o u r c e s a r e o m i t t e d .
ac + ~u a c - -(--) a DLac a t R ~ X ax R a x
+ hc
=
(6.36)
0
The t r a n s p o r t e q u a t i o n i s a s e c o n d - o r d e r
p a r t i a l d i f f e r e n t i a l equation which
r e q u i r e s i n i t i a l c o n d i t i o n s and b o u n d a r y c o n d i t i o n s . by t h e c o n c e n t r a t i o n d i s t r i b u t i o n c(x.y.tO) As i n t h e c a s e o f t h e f l o w e q u a t i o n , conditions.
We o b t a i n i t f r o m ( 6 . 3 5 )
0 and t h a t t h e c o n c e n t r a t i o n c h a s n o g r a d i e n t s i n y - d i r e c -
I n i t i a l conditions are given
a t t h e s t a r t i n g t i m e to o f s i m u l a t i o n .
t h e r e a r e t h r e e p o s s i b l e t y p e s o f boundary
F i r s t k i n d ( D i r i c h l e t ) b o u n d a r y c o n d i t i o n s s p e c i f y p r e s c r i b e d con-
c e n t r a t i o n s on t h e b o u n d a r y .
I n n e r b o u n d a r i e s o f t h e f i r s t k i n d can d e s c r i b e p o l -
206
lution sources. Second kind (Neumann) boundary conditions specify the concentration gradient normal to a boundary. This means they prescribe the dispersive flux. The convective flux at boundaries cannot be prescribed by second-kind boundary conditions. It is a result of the interplay of concentrations on the boundary, as given, for example, by first-kind boundary conditions, and velocities. At impervious boundaries the flow model has to yield zero water flow across the boundary. Then we need only prescribe a vanishing dispersive flux, i.e. &/an = 0. Prescribing the total flux, i. e. the sum of convective and dispersive flux, on the boundary corresponds to a boundary condition of the third kind (Cauchy boundary condition), as the total flux is a linear combination of the concentration on the boundary, c, and the normal derivative of the concentration on the boundary, adan. (6.37) If the dispersive flux is small compared to the convective flux, the third-kind boundary condition approaches a first-kind boundary condition. I n the following, we consider first-kind and second-kind boundary conditions only. The combined action of convection, diffusion-dispersion, adsorption, and decay is shown for one-dimensional transport i n figure 6.10. Initially, a homogeneous distribution of dissolved concentration between x=O and x=a is given. While convection alone causes a spatial translation of the concentration distribution, diffusion-dispersion leads to a spreading. When switching on linear adsorption, the dissolved concentration is diminished if the pollutant mass in the distribution remains unchanged, as part of the dissolved pollutant mass will go to the matrix, For convenience of comparison we leave the dissolved pollutant concentrations unchanged, adding the adsorbed pollutant mass to the system. Adsorption causes a delay of the transport process, the spreading included. While convection and diffusion-dispersion leave the pollutant mass unchanged, decay reduces the pollutant mass and consequently the area under the distribution.
207
POLLUTANT DISTRIBUTION AT TIME
W
X
Concentration c
t =O
-
direction of flow
____e
xlo
distance x
POLLUTANT DISTRIBUTION AT TIME
t, > 0
X=Ufr
X
X
CA
*
c 4 Effect of convection, dispersion, adsorption. degradation
F i g . 6.10:
Schematic d e s c r i p t i o n o f t h e e f f e c t s o f c o n v e c t i o n , d i s p e r s i o n , a d s o r p t i o n , and chemical d e g r a d a t i o n on p o l l u t a n t t r a n s p o r t