- Email: [email protected]
A method to solve the advection-dispersion equation with a kinetic adsorption isotherm
Advances in Wafer Resources, Vol. 19, No. 4, pp. 193-206, 1996
ELSEVIER
Copyright 0 1996 Elsevier Science Limited Printed in Great Britain. All rights reserved 0309-1708/96/$15.00+0.00
0309-1708(95)00045-3
A method to solve the advection-dispersion equation with a kinetic adsorption isotherm J. J. A. van Kooten* Department of Mathematics, Wageningen Agricultural University, Dreyenlaan 4, 6703 HA Wageningen, The Netherlands (Received 8 May 1995; accepted 23 November 1995)
This study deals with a method to solve the transport equations for a kinetically adsorbing solute in a porous medium with spatially varying velocity field and dispersion coefficients. Making use of the stochastic nature of a first-order kinetic p~:ocess,we show that the advection-dispersion equation and the adsorption isotherm can be decoupled. Once the solution for a non-adsorbing solute is known, the method provides an exact solution for the kinetically adsorbing solute. The method is worked out in four examples. In particular we demonstrate how the method can be applied simultaneously with a numerical transport code: the advective-dispersive transport is computed numerically, whereas kinetic effects are incorporated analytically. The proposed approach may be useful in field scale applications with complex flow patterns. Copyright 0 1996 Elsevier Science Limited Key worak kinetic adsorption, advection, dispersion, multi-dimensional varying flow, residence time distributions, analytical solutions.
1 INTRODUCTION
spatially
mathematics can be understood from the physics of the kinetic process. In this way the practical applicability of the kinetic model may be enlarged. The conventional advection-dispersion model for a solute subject to first-order reversible kinetics is
Solute transport in groundwater may be highly affected by interactions between the solute and the solid. Therefore transport models should contain equilibrium or non-equilibrium isotherms that account for these interactions. The most simple model is obtained by applying the linear equilibrium isotherm. Linear equilibrium adsorption only retards the transport of a solute. The local equilibrium assumption (LEA) may be valid if the rate of mass change due to the adsorption process is much faster than that (due to the flow process. If the LEA is not valid a kinetic model may give a better description of the solute transport. Several authors3T18,36 have presented criteria to decide whether an equilibrium or kinetic model is required. An equilibrium approach is often preferred above a kinetic approach; the retardation factor in the equilibrium model can be estimated more easily from experimental data than the kinetic reaction rates. Moreover, a kinetic model requires more complicated mathematics than the equilibrium model. One of the aims of this study is to show that the
( )
g+g=-$(UiC)+& I
X3 ,,=k,C-k,S,
I
+Y(X),
Dijg
J
(14 (lb)
where x is the vector of space coordinates, C(x, t) (M/L3) is the concentration of the solute in the free phase, S(x, t) (M/L3) is the concentration in the adsorbed phase, Y (L/T) is the fluid velocity vector, Dij (L2/T) is the dispersion tensor, y(x) (M/L3T) is the rate of zero-order production at the point x and kl (T-i) and k2 (T-l) are the forward and backward reaction rate, respectively. With respect to the indices we use the Einstein summation convention. Model (1) is mathematically equivalent to a dual porosity model, i.e. a model for transport in a medium in which two zones can be distinguished: a mobile zone in which flow and transport take place and an immobile zone (e.g. dead end pores) in which the water is nearly stagnant.”
*Currently at: Baan Research, P.O. Box 250,’6710 BG Ede, The Netherlands. 193
194
J. J. A. van Kooten
In the literature analytical solutions of eqn (1) are available for various initial and boundary conditions. Lindstrom and Narasimham26 have presented a onedimensional solution for a previously distributed contaminant. Solutions for migration through a semiinfinite column with different types of solute input at the inlet have been obtained by Ogata,27 Cameron and Klute,’ van Genuchten and Wierenga” and de Smedt solutions for a and Wierenga.32 Three-dimensional point spill of pollution in an infinite domain have been derived by Carnahan and Remer6 and Goltz and Roberts.13 Lassey23 has derived solutions for arbitrary solute input after t = 0. Also, recently solutions have been derived with zero-order production (y # 0); Torido et a1.35 have presented a wide class of onedimensional solutions assuming arbitrary input concentration, initial conditions and production functions. In all papers the velocity field was assumed to be uniform and one-dimensional. The dispersion coefficients were taken as constant. At field scale these conditions will often not be satisfied. All solutions were obtained with the Laplace transform method. From a mathematical point of view this method may be quite powerful. However, it does not always give an insight into the relationship between the structure of a solution and the underlying physical processes. Moreover, the Laplace method fails if model parameters are space dependent. Studies dealing with the transport of a kinetically adsorbing solute in a porous media with heterogeneous hydraulic and/or chemical properties have not primarily focused on the resident concentrations, but on the spatial moments of the solute plume1’7>29’37 or on the arrival time distribution at a plane perpendicular to the mean fluid flow.s>30>31~33~34 The objective of this study is to develop a method to solve the concentration, C, and S from model (1) in a multi-dimensional system with spatially varying w, D and y. The sorption rates kl and k2 are assumed to be constant. All coefficients are assumed to be timeindependent. We interpret the adsorption kinetics as a stochastic process in which a particle is alternately free or sorbed (Section 2). Keller and Giddings” have presented four conditional distributions for the residence time in the free phase during a time interval t. For arbitrary initial and boundary conditions we show that the solution of eqn (1) is a convolution of solutions of the classical ADE (la) without the adsorption term and a combination of these distributions (Section 3). So, in fact, we show that advectiondispersion eqn (la) and the reaction eqn (lb) can be decoupled. Analytical solutions of the classical ADE have been derived for numerous initial and boundary conditions.9124125Now, for all these problems we also have the solution for a kinetically adsorbing contaminant. Some examples are considered in Section 4. Furthermore, we demonstrate how the method may
be applied simultaneously with a numerical transport model. The transport equation is solved numerically, whereas the kinetic adsorption is incorporated analytically.
2 RESIDENCE TIME DISTRIBUTIONS IN THE FREE PHASE The kinetic adsorption process described by eqn (lb) can be considered as a stochastic process in which particles are alternately in the free and sorbed phase. Between transitions the residence times in the free (f) and sorbed (s) phase are exponentially distributed with mean l/k, and l/k2, respectively. During a time interval t a particle may make several transitions from the one state to the other. Assuming that at least one transition has taken place, Keller and Giddings” have presented probability density functions for the time fraction x (0 5 x < 1) that a particle has spent in the free phase (see also Giddings and Eyring”). These distributions depend on the initial and final state of a particle. For the purpose of this study we rewrite these distributions in terms of r (= xt ), the residence time in the fluid phase during a time interval t
$j-(T t) =
c
Co
(klkzT)“(t
-
+-’
e-(k,-k2)T-k2f
n!(n - I)!
n=l
=
A~(‘,
t)
=
kl
2 n=O
=
kl
(k1k2T)n(t
-
T)”
e-(k,-k2)T-k2f
(n!)2
e-_(kl-kzk-kzf
Zo(2&&2+
-
4,
(3)
where the first and second subscript of the functions h denote the initial and final state of a particle, respectively. In eqns (2) and (3) Z. and Zi denote the modified Bessel functions. For the distributions Zzf~and hq the condition that at least one transition has taken place is satisfied automatically. In the two other cases there is also the possibility that a particle remains in the initial state all the time; the probabilities that a particle spends time I in the free or sorbed phase are exp (-k, t) or exp (-k2t), respectively. Some integration and differentiation properties of the distributions are given in the Appendix A. In general, only the expressions in terms of the modified Bessel functions are given (see e.g. Valocchi and Quinodoz38). We have also given the summation representation because it is closer to the underlying
A method to solve the advection-dispersion equation
physical process: e.g. the counter n is the number of times a particle has returned to the initial state. Only in the free-p:hase particles are subject to advection and dispersion. In the sorbed phase particles are immobile. Because the adsorption rates are constant the above-described exchange process between the free and sorbed phase is not affected by the advectiondispersion process, even when the advection and dispersion parameters are space dependent. So, the adsorption process and the advection-dispersion process may be studied separately. Once the temporal evolution of the advective-dispersive displacement is known, the adsorption kinetics can be incorporated naturally by utilizing the fact that the residence time, 7, in the free-space is stochastic; with the aid of the distributions (2), (3), (4) and (5) the expected advective-dispersive displacement can be computed. In Section 3 we use this idea to decouple the transport eqn (la) and the adsorption isotherm (lb). Other authors have applied the same idea in different contexts. Valocchi and Quinodoz38 have developed an efficient random walk algorithm for simulating the spread of a kinetically adsorbing particles. In each time step At they use the above h-distributions to determine the fraction of time that a particle has spent in the free phase. A modification of the algorithm has been described by Andrecivic and Foufoula-Geourgiou.’ They show that for large k2t the distributions approach Gaussian distributions that can be fully characterised by the first and second moment of the residence time. In Quinodoz and Valocchi38 expressions are presented for the first and second spatial moment of an instantaneously injected contaminant in porous formations with randomly heterogeneous hydraulic conductivity. Given the spatial moments of a non-reacting contaminant, they use distribution (2) to derive an expression for the spatial moments of a kinetically adsorbing contaminant. Comparable results have been obtained by Cvetcovic and Shapiro8 (eqn (31)) and Selroos and Cvetcovic3’ (eqn (16)), who analyse the effect of kinetic adsorption on the arrival time distribution of a solute at a control plane in a heterogene0u.s porous formation.
3 DECOUPLING
OF THE TRANSPORT REACTION EQUATION
195
In the sequel we refer to MAD as the advectiondispersion operator. Let 0 (R c RN, N = 1,2 or 3) be the region in which contaminant transport is studied. The boundary dfl may consist of two disjunct parts. At the part dRD the concentration is prescribed. At the part dfi,v the mass flux is prescribed. In mathematical terms dRD is a Dirichlet boundary and d& is a Neumann boundary. A method to solve the kinetic model (1) for arbitrary initial and boundary conditions is developed in three steps. First we consider the case of non-zero initial conditions and zero boundary conditions. Next, we consider the opposite case: zero initial conditions and non-zero boundary conditions. Finally, the solution for arbitrary initial and boundary conditions is obtained by superposition of the two former cases. We also give a physical interpretation of the solutions. 3.1 Transport of a contaminant distributed at or before t=O
In this section we consider the transport of a contaminant which has been distributed at or before t = 0. It is assumed that the distributions of both C and S at t = 0 are known C(X, 0) = Ci”i,(x) and
S(X,
0) = Sinit(
(7)
Furthermore, we assume that no contaminant is added to or extracted from the system so that y = 0 and the boundary conditions are
aciD,
( C(x, t) = 0
at
[ni*(ui-Dij&)C=O
atdoN,
(8)
where n is the outward normal vector on 80. An expression for the solution of model (1) subject to these initial and boundary conditions can be found with physical arguments. We first give this expression. Next, we explain how it is obtained. C(x,
t)
=
C;(x,
t) e+’
S(X,
t)
=
Sinit
e+’
AND
In this section we show that the solutions C and S can be expressed as an integral of a proper combination of the h-distributions and C’, the concentration of a non-adsorbing contaminant satisfying
+ 1;
(c’(X,
t)
+ cs,(X,
~)h,(~>
t))
d7, (9b)
where both C’ and 6;, are solutions of ADE (6) (with y = 0) subject to the initial conditions
with
(6)
4&b,
C;(x7 O) =
Cinit(X)
and
C.f(x, O) =
&nitCX), (10)
and zero boundary conditions.
J. J. A. van Kooten
196
The physical reasoning behind expression (9) is as follows. During the time that particles are in the free phase their transport is governed by the classical ADE (6). The remaining time particles are immobile. The residence time in the free-phase r is stochastic. With the help of the residence time distributions frfl, hf,, h,f and h, the expected advective-dispersive displacement can be computed. The free phase transport of the contaminant with initial state ‘f ’is described by C;, and the free phase transport of the contaminant with initial state ‘s’is described by 6;,. To compute the concentration C in the free-phase at time t the distributions with final state ‘f'should be used and to compute the concentration S the distributions with final state ‘s’ should be used. The term exp (-k,t) in eqn (9a) represents the probability that particles have remained in the free-phase continuously. The term exp (+t ) in (9b) represents the probability that a particle has not moved at all. The mathematical verification is given in Appendix B. It is often possible to simplify expression (9). For example, if all contaminant is injected instantaneously at t = 0 (see Section 4.1) then S = 0, so that C;, = 0. In many applications it is assumed the initial distributions C and S are in equilibrium, which implies that Cs = (k,/k&“, so that only C; has to be determined. A one dimensional example of the latter case is discussed in Lindstrom and Narasimham.26
have the same structure. Generalizing the result to multi-dimensional systems with spatially varying velocity fields and dispersion coefficients and y # 0 we have found that the solution of (1) with zero initial conditions and boundary conditions (12) can be expressed as C(x, t) = C,(x, t) eCkir +~~~o(x,~){h~(T,t)+~h~~(~,t
(134
(13b) where Co is the concentration of a non-adsorbing contaminant satisfying ADE (6) with zero initial condition and boundary conditions given by eqn (12). The mathematical verification is sketched in Appendix B. Because contamination is constantly entering and leaving the porous medium the physical interpretation of expression (13) is less straightforward than that of (9). However, using the differentiation properties of Appendix A one may check that the time derivatives of (13) can be written as
a+, t) -=
aco;,(x, t)
e-k,r
at
3.2 Transport of a continuously injected contaminant If a domain is uncontaminated conditions are C(x, 0) = 0
and
at t = 0, the initial
S(x, 0) = 0.
(11)
We assume that from t = 0 there is a constant input and output of contaminant. This input/output may take place at the boundary of the domain and inside the domain, e.g. at an injection or a withdrawal well. Contaminant input/output at the boundary is described by at dRo,
( CC& t) = C,(x) ~n;‘(,iC-Dii~)
=F(X)
at dRN.
wx, 1) -= at
at
f
+
aco,(x, 4 hff(T,t>dT,
f
lb
f0
t aCo(x,T) h,dT,t) dT. 0 aT
(14b)
The physical interpretation of eqn (14) is straightforward. When particles are in the free phase the change in the concentration is described by eqn (6). Because the residence time in the free phase is stochastic, this change should be weighted with the proper residence time distribution. Because particles enter or leave the medium in the free phase, only the distributions with initial state ‘f' are relevant. The term exp (-k, t ) in eqns (13a) and (14a) again represents the probability that particles have remained in the free-phase continuously.
(12) 3.3 Arbitrary initial and boundary conditions
Production/withdrawal inside the domain is described by the term y in eqn (1). The functions C,, F and y may vary in space. One-dimensional examples of this case with y = 0 are discussed in 0gata,27 Cameron and Klute,’ van Genuchten and Wierenga” and de Smedt and Wierenga.32 We have rewritten their solutions in terms of the residence time distributions hff, hfi, hd and h,. For the solutions that are expressed in terms of the Goldstein J-function’2 we have applied relation A3 and A4 of Appendix A. All solutions turned out to
The solution of model (1) (with y # 0) subject to non-zero initial conditions c(x, O) =
Chit(X)
and
s(x, O) =
GdCX)7 (15)
and non-zero boundary conditions ( C(x, t) =
cb(x)
at f3RD,
(16)
A method to solve the advection-dispersion equation is the superposition of eqns superposition can be written as C(x,
t) =
Cl(x,
(9) and
(13). This
t) e+’
197
neither the kinetic nor the equilibrium model describe the transport of a solute satisfactorily.5 A combined equilibrium and kinetic approach may be required:
r +
J
,@x,+J+,t)
+cz(
x, 4
h&,
t))
d7, (174
S(X,
t) =
Sini,
e+’
+ 1; (c,
(x, +fi(~,
t) + z;(x,
%(7,
t))
d7, (17b)
where C, E Cf + C,, and Cz = Cs + (ki/kz)C,-,, so that Ci and Cz are concentrations of a non-adsorbing contaminant satisfying eqn (6) with y](x) = y(x) and y2(x) = (k,/k&y(x), respectively. The initial and boundary conditions are C‘1(x~O) =
Cinit(X)7
where k3 (T-l) and k4 (T-l) again denote a forward and backward sorption rate, respectively. After introduction of the retarded parameters 2rf = vi/R, DC = DijlR, 7R = y/R, SR = S/R and @ = kJ/R model (23) transforms into a model equivalent to (l), so that the method presented in this study also applies to the combined model. As a point of mathematical interest we mention that the decoupling method is also valid for models in which the advection-dispersion operator MAD (6) is replaced by an arbitrary linear differential operator M=a(x)+bi(x)~+Cij(X)~ I
(i,j,k, at &lo,
(19) A considerable simphfication is achieved if it is assumed that the initial concentration distributions are in equilibrium, i.e. = (kl lk2)Cinit
(20)
which implies that C2 := (ki/k#i to be determined.
so that only Ci has
3.4 Further generalisations If contaminant input only takes place during time z,,, i.e. for t 5 Ti:,, the boundary conditions and/or y are not equal to 0 and for t > T they are 0, then the solution of model (1) is Cr(x, t) = C(x, t)
and
ST(X,
for t 5 Tin,
t) =
S(X,
t)
CT(x, t) = C(X, t) - C(X, t - Tin) ST(X,
t) =
S(X,
t) -
S(X,
t -
Tin)
and
1
J
etc = 1.. .N),
where a, bi, cij, dijk etc. are scalar functions position x.
at a0,v.
&it
(22)
a2
at bRo,
f e’1 (x, t) = cb(x)
as ;;r;-=k3C-k4S
(21)
for t > Tiny
where C and S are given by eqn (17). It is well known that for large kl and k2 a kinetic model approximately predicts the same transport as a linear equilibrium model with retardation factor R = 1 + (k,/k2). On the other hand, it may occur that
(23)
of the
4 SOME EXAMPLES In Section 3 we have shown that once the transport of a non-adsorbing solute is known, then the transport of a kinetically adsorbing solute can be determined. In this section the theory is demonstrated in some examples. In the first example we show how a solution of a kinetic model presented by Carnaham and Remer6 and Goltz and Robertsi can be obtained as a special case of the method described in this study. In the second example the method is applied to predict the breakthrough of a contaminant at a well. We extend the results for the breakthrough of a non-adsorbing contaminant obtained in previous publications.‘4Y21 In examples 3 and 4 we demonstrate how the method can be combined with a numerical transport model. In all examples we assume that the dispersion tensor is given by Dij = aTIu[Sij + (aL - ar) z,
(24)
where aL and aT denote the longitudinal and transversal dispersivities.2 In example 2 we allow aL and aT to vary in space. The integrals with respect to r (see e.g. expressions (9), (13) (17) (27) and (31)) were computed with a Gauss-Kronrod rule of the IMSL-Math Library.” For the modified Bessel functions in the h-distributions we
J. J. A. van Kooten
198 have implemented Blair.4
the Chebyshev
approximations
of
3,5,
-NON-~RBING,S~
. ,
4.1 Example 1: 3-D transport of instantaneously injected solute In three dimensions the kinetic model for transport in an isotropic formation with uniform flow in the xl-direction reads
(254
G-4
At t = 0 an amount Q (M) of solute is injected in the origin (xi, x2, x3) = (0, 0, 0), so that the following initial conditions apply to eqn (25a). C(x, 0) = QW
and
S(x, 0) = 0,
KINETICALLYADSORBINGSOLUTE 3.
(25b)
where S(x) denotes the Dirac delta function. As boundary condition we require that C vanishes at infinity. The concentration of a non-adsorbing solute, i.e. kl = k2 = 0 in eqn (25) is given by
DH
I
3 022
I
x:
s. -8 2P 4 1.5 -
t t) epklt + o zf(x,
S(x, t) = 1 c’(x, $I&, J
J
0.5 o-
(see e.g. Carnahan & Remer6). Now, from eqn (9) we directly obtain the concentrations of the kinetically adsorbing solute in the free and sorbed phase C(x, t) = cf(x,
l-
033
(26)
T)~#(T, t) dr,
J
n- 2.5 -
P 8
(Xl -A2
1
-I
(b)
Fig. 1. Evolution in time of the concentration profiles c;- (eqn (26)), C and S (eqn (27)) along the xl-axis. (Q = 100(M), Dz2 = D33 = 0.05 (L2/T), wI = 1 (M/L), D,, = 0.25 (L’/T), k, = k2 = 0.1 (T-l)). 4.2 Example 2: Contamination of a pumping well
(27)
t) dr.
(Note that in this example CX= 0, because S = 0 at t = 0). Solutions equivalent to eqn (27) were presented by Carnahan and Remer6 and Goltz and Roberts.13 However, they obtained the solutions after lengthy Laplace transforms. Here, the solutions simply follow with a more general theory which, moreover, is closer to the underlying physical process. In Fig. 1 a comparison is made between the evolution of the concentration profiles C’, C and S along the xi-axis (x2 = x3 = 0). For small t the concentration profiles of C and S have a strong non-Gaussian shape. This can be explained from the fact that shortly after the injection the solute concentration close to the origin is high so that a large amount of solute will be adsorbed. The adsorbed amount is released slowly. Furthermore the figures clearly illustrate the enhanced longitudinal spreading and tailing due to kinetic adsorption.
In previous publications we have studied the arrival of conservative or linearly decaying contaminants at a well in a confined aquifer (e.g. van Herwaarden & Grasman,” van der Hoek,16 van Herwaarden,14 van Kooten,21’22). A typical flow pattern for groundwater discharge is shown in Fig. 2. Van Herwaarden and Grasman” have analysed the effect of the transversal dispersion on the arrival rate. If a contaminant is subject to first-order kinetics, then during intervals a particle is adsorbed and not displaced in the transversal direction. Thus, kinetic adsorption only affects the time at which particles pass a certain point and not the ultimate transversal spreading. Therefore, the approximation for the fraction of a contaminant that reaches the well presented by van Herwaarden14 (eqn (4.5)) may also be used for a kinetically adsorbing contaminant. In van Kooten22 we have derived approximations for the arrival time distribution at a well. For example, for particles released in the point x far inside the
A method to solve the advection-dispersion
equation
199 I
Fig. 2. Flow pattern for a well in a uniform background flow. catchment area, where the effect of the transversal dispersion is negligible, the arrival time distribution can be approximated by (28a) where p and rs2 are asymptotic approximations mean and variance of t:he arrival time
/+) =
Tadv(X)
ai,C-4
+ F(x)
for the
&(t) = s,(t) =
JR
conservative cases are Ci,it(x)g(x, t) dx
and (29)
R Git(X)g(X, t) dx,
J and, thus, the arrival time distribution for a kinetically adsorbing contaminant is g(t) = &(t) eFklr
and
+ 1; (&(+,j+,
(28b) 02(x) = 2 In (28b) Tad”is the advlective travelling time to the well and C is a coordinate along the streamline through x, see Fig. 3. Note that for any flow pattern p and c? can easily be determined by particle tracking. For a detailed description of how to take into account transversal dispersion we refer to van Kooten.22 With the residence time distributions of Section 2 kinetic adsorption can be incorporated without making any further approximation. It is assumed that in an aquifer Q at t = 0 the concentrations of the contaminant in the free and sorbed phase are given by Ci,i,(X) and sinit( The arrival time distributions of the
Fig. 3. The coordinate C is taken along the streamline through
the point x.
corresponding
For an instantaneous reduces to
t) +&W&
t)) d7
(30)
spill in the point x eqn (30)
g(x, t) = 2(x, t) eCklt + : g(x, T)hfl(T, t) dz (31) J We have worked out the above for the special case of a well in a uniform background flow (Fig. 2). The corresponding velocity field is given by 211(x,,x2) = 1 - 2 x: + x; ’
-x2 V2(Xl,X2)
=
-. xy + x; (34
We assume that at t = 0 in the point (x1,x2) = (-4,l) an instantaneous spill of contaminant takes place. For various values of kl and k2 the approximation for the arrival time distribution (31) at the well is displayed in Fig. 4. The dispersivities are aL = 0.1 and aT = 0.01. The figures show that the tail due to kinetic effects becomes more pronounced if k, increases and/or k2 decreases, i.e. if the expected time of one stay in the fluid phase (l/k,) decreases and/or the expected time of one stay (1 /k2) in the fluid phase increases. To assess the accuracy of the approximations we have made a comparison with the results of random walk simulations. There exist various random walk algorithms
J. J. A. van Kooten
200
AI~~~~TIMEDISI’RIBUTIONSAT~EI~.
~~~DISIRIBUTIOI’JSAT~~ 0.7b ( 0.6 .
,
0.6 -
wa4dsorbingaln~nant
nm-ads&ingamtaminant
8
30
35
40
tO
(a)
0.7c
Oo
45
5
10
(b)
15
20 25 t0.J
30
35
40
45
Fig. 4. Some arrival time distributions at the well for a contaminant spilled in (-4,l).
In (a) the adsorption rate ki (T-l) is kept fixed. In (b) the desorption rate k2 (T-i) is kept fixed.
that take into account kinetic adsorption.‘>38 Here, it suffices to select the most natural and straightforward method (method 1 of Valocchi and Quinodoz38). In the free-phase the stochastic displacement of a particle is governed by AX=
I3 wl+-&,+-~i2 ax (
d ay
At >
+~~$z,&t+&&$$zr&t v2+&D2i
AY= (
+2022 8Y
At >
+~~$z&t-~z+$z*6t (33)
where At is the time step and Z, and ZT are random deviates from a uniform distribution with mean 0 and variance 1 (see Kinzelbach20 or van Kooten*‘). First the time r- of one stay in the fluid phase is generated from an exponential distribution. During time r- a particle is advected-dispersed according to eqn (34). Next, from an exponential distribution the time rs that a particle stays in the adsorbed phase is generated. The time T$ is added to t, the total time elapsed from the beginning of the simulation. This process is repeated, until the particle enters the well. we carried out 15,000 random From the point (-4,l ) and k2 = 0.2 (T’). The time walks with ki = O-5 (T- 1’ step At was chosen to be 0.001 except if rf was exceeded. In the latter case At was taken equal to the difference between T- and the time that has been spent in the
COMPARISONOF RANDOM WALK AND ASYMFIWI’IC
APPROXIMATION
42 -
! OSB * &
1
O-l0.05 -
OO
Fig. 5. Comparison of the arrival time distributions obtained with the random walk simulation and computed from eqn (31).
A method to solve the advection-dispersion equation
fluid phase yet. From Fig. 5 we may conclude that approximation (31) agrees quite well with the arrivaltime distribution obtained from the simulations. The small differences are mainly caused by the irregular behavior of the simulat:ions, and it is to be expected that the agreement will still become better if the number of random walks is increased.
201
transport is computed numerically whereas the kinetic adsorption is incorporated analytically (in the sequel we refer to the latter one as the mixed numericalanalytical solution). We study the lD-transport of a solute through a semi-infinite column. From t = 0 to t = Tin the concentration at the inlet x = 0 is maintained at some constant value C(0, t) = Cb. Next no more solute enters the domain, i.e. C(0, t) = 0 for t = Ti:,n. Ogata and Banks** and Ogata27 have presented solutions for a non-adsorbing and kinetically adsorbing solute, respectively. The solution for the kinetically adsorbing species can also be obtained by substituting the solution for the non-adsorbing species in expression (13). With the FE-code COLTRAP3’ we have simulated the transport at the interval 0 I x 5 1 with I = 100(L), w = 1 (L/T), D = 1 (L*/L), Cb = 1 (M/L) and Tir, = 10(T)). At the endpoint x = 1 a no-flux condition is imposed. The fact that the interval is finite does not undermine the comparison with the analytical solution for the semi-infinite system; van Genuchten and Alves’ have shown that for Peclet numbers Pe = vi/D greater than 20 the solutions for the finite and infinite system deviate only close to the endpoint. We have run COLTRAP with various values of the
4.3 Example 3: Comparison between an analytic solution,
a numerical solution and a mixed numerical-analytical solution Currently, there exist many finite element (FE) and finite difference (FD) codes to solve the ADE for a conservative contaminant. The aim of the examples, 3 and 4, is to demonstrate blow the theory of Section 3 can be used simultaneous1.y with such a code. Once the transport of a conservative contaminant has been computed numerically, with some post processing the transport of a kinetical1.y adsorbing contaminant can be computed. In this example we make a comparison between an analytical solution of the kinetic model, a full numerical solution and a solution in which the advective-dispersive TIME PROPILE IN x=30
aes.
.
.
TIME PROFILE IN x=30 . . . . .
.
.
= aMlyticai dution ._______ = numerical 8ohtion .... .........
= mixed solution
k,=O.l(ur)
xl
40
60
a0
100
12-O 140 160
160 xl0
tm
tO
(a)
TIME PROFILE IN x=30
.
ai2.
= mixed solution
._______= numerical solution
k2 = 0.02 (l/T)
k, = O.&i (1JQ k2 = 0.02 (l/l)
aaz.
Oo#2040 (4
.
= mixed solution
ao2. I-J.
.
= analytical solution 0.1.
._______ = numericalsolution
A, = 0.M (lm
TIME PROFILE IN x=30 . _ . .
ai2.
analytical solution 0.1 -
-
~==O.Ol(Lm
-
. 60
. eo
-
100
.
l#)
.
.
140 160
.
I
180 #r)
tc0
Fig. 6. Analytical solutions2’, numerical solutions and mixed numerical-analytical solutions for transport through a semi-infinite column. We have displayed the time profiles of the free phase concentration C (M/L) in the point x = 30. (w = 1 (L/T), D = 1 (L/T’), T, = IO(T), gridsize Ax = 2.0(L)). The numerical simulations of (a) and (c) were carried out with time step At = 1.0. In (b) and (d) At = 0.5.
J. J. A. van Kooten
202
flow takes place. Initially the aquifer is unpolluted. From t = 0 at the top of the aquifer for x between 7 and 10 m and y between - 1 and 1 m contaminated water constantly leaches into the aquifer. The concentration of the contaminant in the incoming water is C, = 100 pg/m3. At the remainder of the top surface C, = 0. Because of the symmetry in the y-direction it suffices to solve the transport problem in only one half of the aquifer, with a no-flux condition at the surface y = 0. At the downstream boundary x = 40 we impose a zero concentration condition. The surfaces x = 0, y = 0 and z = 0 are no-flux boundaries. The aquifer was covered by 40 x 20 x 12 = 9600 cubic elements. Because of the zero initial conditions the theory for simulating the transport of a kinetically adsorbing contaminant is described in Section 3.2. With the METROPOL code the concentration profile C0 of a non-adsorbing contaminant was computed from t = 0 to t = 20 days with time step At = 0.2 days and dispersivities aL = 0.6 m and aT = 0.3 m. The computation time at a Sun3 workstation was about 6 h. In Fig. 8(a) the flow pattern and the concentration contours in the y = 0 plane at t = 20 are shown. For each nodal point the time profile of Co was stored in a data file. The time profiles were postprocessed by substitution in eqn (13). In between partition point ti = iAt the concentration was determined by linear interpolation. The concentration contours for a kinetically adsorbing contaminant with kl = 0.1 and k2 = 0.1 (days-‘) are displayed in Fig. 8(b). The post processing could be carried out in about 2min, so that the costs of the post processing are almost negligible compared to the cost of solving the ADE numerically. In Fig. 8(c) we have drawn the concentration contours resulting from an equilibrium model with retardation factor R = 1 + (kl/k2) = 2. The solution of the equilibrium model was obtained from the METROPOL output by replacing the time t by the retarded time tR = t/R. Compared to the kinetic model the downstream spreading is clearly underestimated. The ‘kinetic’ profile is more flat than the ‘equilibrium’ profile. So, the application of a kinetic or equilibrium model may have large consequences on, for example, the delineation of
gridsize Ax and timestep At. We have computed the transport of non-adsorbing and kinetically adsorbing contaminants from t = 0 to t = 200(T). The time profiles of the concentrations in the point x = 30 were stored in data files. The time profile of a non-adsorbing contaminant was post-processed by substitution in eqn (13). In between partition point ti = iAt the concentration was determined by linear interpolation. We applied applied expression (21) to obtain a solution for pulse input. From the simulation results it (perhaps surprisingly) turned out that the accuracy of the full numerical and the mixed numerical-analytical solution is of the same order. Some examples are shown in Fig. 6. One might expect that the mixed solution is slightly better because no error is introduced due to the time discretization of the adsorption isotherm (1 b). However, this error appears to be dominated by the error due to the space-time discretization of the parabolic ADE (la). The main advantage of the mixed (or decoupling) method is that for different adsorption rates only the post-processing has to be repeated, whereas for the full numerical method the simulation code has to be rerun each time. So the mixed method is computationally much more efficient.
4.4 Example 4: Simultaneous application with a numerical model for transport in spatially varying flow In this example we demonstrate how the method presented in this paper can be applied to more complicated flow patterns. The numerical computations were carried out with the FE code METROPOL, which has been developed by the National Institute of Public Health and Environmental Protection in the Netherlands. METROPOL can only deal with non-adsorbing or linearly decaying solutes. With our method the effect of kinetic adsorption is incorporated. We consider contaminant transport in an aquifer which is recharged at the top (Fig. 7). From below the formation is confined by an impervious layer. The effective velocities (m/day) of the incoming and outgoing flow are indicated in the figure. In the y-direction no
V=O.l
V = 0.6 -t
10m I
z
Y k
40m
x Fig. 7.
Sketch of the aquifer considered in Section 4.3.
J. J. A. van Kooten
204
LZ~.~~) or cannot be derived at all. Therefore, the method will probably find its major application when it is used in combination with a numerical transport model. With such a model the transport of a non-adsorbing contaminant is computed first. Next, by substitution of the output in expressions (9), (13) or (17) the effect of the kinetic adsorption can be incorporated analytically. Examples are given in Sections 4.3 and 4.4. The proposed approach may be useful in field-scale applications with complex flow fields. The main advantage above solving the complete model (1) numerically is that no new run with the numerical code has to be performed when the values of the reaction rates k, and k2 are changed. Only the post-processing has to be repeated. This feature may be cost-effective when the model output has to be fitted with experimental data. Another advantage is that the method is computationally more efficient. The costs of solving one differential equation numerically are lower than the cost of solving two coupled equations. Although in the mixed numerical-analytical method no errors are introduced due to time discretization of the reaction isotherm (lb), example 3 (Section 4.3) revealed that a mixed solution has the same order of accuracy as a full numerical solution; the dominating error stems from the space-time discretization of the ADE. A disadvantage of the approach may be that the time profiles of the concentrations in the nodal points have to be stored, so that for large grids much computer memory is required. In practice this disadvantage is not too large, because one is often only interested in the breakthrough of a contaminant at some particular points. In this case, just as for the traditional numerical method, only in the observation points do the time profiles need to be stored. ACKNOWLEDGEMENTS
residence time distributions. Stoch. Hydrol. Hydraul., 5 (1991) 155-71.
2. Bear, J. & Verruyt, A., Modelling Groundwater Flow and Pollution. Reidel, Dordrecht, 1987. 3. Bahr, J. M. & Rubin, J., Direct comparison of kinetic and local equilibrium formulations for solute transport affected by surface reactions. Water Resour. Res., 23(3) (1987) 438-52.
4. Blair, J. M., Rational Chebyshev approximations for the modified Bessel functions IO(x) and 11(x). Maths Comp., 28 (1974) 581-3.
5. Cameron, D. R. and Klute, A., Convective-dispersive solute transport with a combined equilibrium and kinetic adsorption model. Water Resour. Res., 13 (1977) 183-8.
and 6. Carnahan, C. L. & Remer, J. S., Nonequilibrium equilibrium sorption with a linear sorption isotherm during mass transport through an infinite porous medium: some analytical solutions. J. Hydrol., 73 (1984) 227-58.
7. Chrysikopolous, C. V., Kitanidis, P. K. & Roberts, P. V., Generalized Taylor-Aris moment analysis of the transport of sorbing solutes through porous media with spatiallyperiodic retardation factor. Trans. Porous Media, 7 (1992) 163-85.
8. Cvetkovic, V. D. & Shapiro, A. M., Mass arrival of sorptive solute in heterogeneous porous media. Water Resour. Res., 26(9) (1991) 2057-67.
9. van Genuchten, M. Th. & Alves, W. J., Analytical solutions of the one-dimensional convective-dispersive solute transport equations. Tech. Bull. 1661, US Department of Agriculture, Washington, DC, 1982, 151 pp. 0. van Genuchten, M. Th. & Wierenga, P. J., Mass transfer studies in sorbing porous media. 1. Analytical solutions. Soil Sci. Sot. Am. J., 40 (1976) 473-80.
1. Giddings, J. C. & Eyring, H., A molecular dynamic theory of chromatography. J. Phys. Chem., 59 (1995) 416-21. 2. Goldstein, S., On the mathematics of exchange processes in fixed columns. 1. Mathematical solutions and asymptotic expansions. Proc. R. Sot. London, Ser. A, 219 (1953) 151-71.
3. Goltz, M. N. & Roberts, P. V., Three-dimensional solutions for solute transport in an infinite medium with mobile and immobile zones. Water Resour. Res., 22(7) (1986).
The author thanks M. de Gee, J. Grasman, P. A. C. Raats and S. E. A. T. M. van der Zee for valuable discussions and for reading the paper critically. The author is grateful to 0. A. Buse for his help with drawing the figures and also thanks R. van de Weerd and the Institute of Public Health and Environmental Protection (in particular A. Leijnse) for their permission to use the COLTRAP and METROPOL code. I also would like to thank the reviewers for their constructive remarks and their suggestions for improving my paper. The work was funded by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Dutch Organisation for Scientific Research). REFERENCES 1. Andricevic, R. & Foufoula-Georgiou, E., Modeling kinetic
non-equilibrium using the first two moments of the
14. van Herwaarden, 0. A., Spread of pollution by dispersive groundwater flow. SIAM J. Appl. Math., 54 (1) (1994) 26-41. 15. van Herwaarden, 0. A. & Grasman, J., Dispersive groundwater flow and pollution. Math. Model Meth. Appl. Sci., 1 (1991) 61-81. 16. van der Hoek, C. J., Contamination of a well in a uniform background flow. Stoch. Hydrol. Hydraul., 6 (1992) 191-208.
17. IMSL. User’s Manual Math/Library, Fortran Subroutines for Mathematical Applications, IMSL, Houston, Texas, 1987. 18. Jennings, A. A. & Kirkner, D. J. Instantaneous equilibrium approximation analysis. J. Hydraul. Engng, 110 (1984) 1700-1717. 19. Keller, R. A. & Giddings, J. C., Multiple zones and spots in chromatography. J. Chromatog., 3 (1960) 205-20. 20. Kinzelbach, W., The random walk method in pollutant transport simulation. In Groundwater Flow and Quality Modeling, eds E. Custodio et af. Reidel, Dordrecht, 1988, pp. 227-45. 21. van Kooten, J. J. A., Groundwater contaminant transport
A method to solve the advection-dispersion including adsorption amd first-order decay. Stoch. Hydrol. Hydraul., 8 (1994) 185;-205. 22. van Kooten, J. J. A., An asymptotic method to predict the contamination of a pu.mping well. Adv. Water Resour., 18 (5) (1995) 295-313. 23. Lassey, K. R., Unidimensional solute transport incorporating equilibrium an’d rate-limited isotherms with firstorder loss. 1. Model conceptualizations and analytical solutions. Water Resour. Res., 24(3) (1988) 343-50. 24. Leij, F. J. t Dane, J. H., Analytical solutions of the onedimensional advection equation and two- or threedimensional dispersion equation. Water Resour. Res., 26 (1990) 1475-82. 25. Leij, F. J., Skags, ‘T. H., van Genuchten, M. Th., Analytical solutions Ifor transport in three-dimensional semi-infinite porous media. Water Resour. Res., 27 (1991) 2719-34. 26. Lindstrom, F. T. & Narasimham, M. N. L. Mathematical theory of a kinetic model for dispersion of previously distributed chemicals in a sorbing porous medium. SIAM J. Appl. Math., 24 (1973) 496-510. 27. Ogata, A., Mathematics of Dispersion with Linear Adsorption Isotherm. US Geol. Surv. Prof. Pap, 411-H, 1964, 9pp. 28. Ogata, A. & Banks, R. B., A Solution of Differential Equation of Longitudinal Dispersion in Porous Media. US Geological Survey Professional Paper, 41 l-A, Al-A9, 1961. 29. Quinodoz, H. A. M. &: Valocchi, A. J., Stochastic analysis of the transport of kinetically sorbing solutes in aquifers with randomly heterogeneous hydraulic conductivity. Water Resour. Res., 2!)(g) (1993) 3227-40. 30. Sardin, M., Schweich, D., Leij, F. J. & van Genuchten, M. Th., Modeling the nonequilibrium transport of linearly interacting solutes in porous media: a review. Water Resour. Res., 27(9) (1991) 2287-307. 31. Selroos, J. & Cvetkovic, V., Modeling solute advection coupled with sorption kinetics in heterogeneous formations. Water Resour. Res., 28 (5) (1992) 1271-78. 32. Smedt, F. de & Wierenga, P. J., A generalized solution for solute flow in soils with mobile and immobile water. Water Resour. Res., lS(5) (1979) 1137-41. 33. Jury, W. A., Sposito, G. and White, R. E., A transfer function model of solute transport through soil, 1. Fundamental concepts. Water Resour. Res. 22(2) (1986) 243-37. 34. Sposito, G., White, R. E., Darrah, P. R. & Jury, W. A., A transfer function model of solute through soil. 3. The advection-dispersion equation. Water Resour. Res., 22(2) (1986) 255-62. 35. Torido, N., Leij, F. J. & van Genuchten, M. Th., A comprehensive set of .analytical solutions for nonequilibrium solute transport with first-order decay and zero-order production. Water Resour. Res., 29 (1993) 2167-82. 36. Valocchi, A. J., Validity of the local equilibrium assump tion for modeling sorbing solute through homogeneous soils. Water Resour. Res., 6 (1985) 808-20. 37. Valocchi, A. J., Spatial moment analysis of the transport of kinetically adsorbin solutes through stratified aquifers. Water Resour. Res., 2!i (2) (1989) 273-9. 38. Valocchi, A. J. & Quinodoz, H. A. M., Application of the random walk method to simulate the transport of kinetically adsorbing solutes. In: Groundwater Contamination: Proceedings of the Third ScientiJic Assembly of the ZAHS, ed. L. M. Abriola. Baltimore, 1989, pp. 35-42. 39. van de Weerd, R. Ct Leijnse, A., Colloid facilitated contaminant transport in porous media. Model verification, calibration and sensitivity analysis. In preparation (1995).
equation
205
APPENDIX A: SOME MATHEMATICAL PROPERTIES OF THE DISTRIBUTIONS hff, 5,
kf AND ASS
Al Differentiation properties By differentiation
3
dt
of (2), (3), (4) and (5), we obtain
= klhff - kzhfi
and
%
= klha - kzh,. (A21
A2 Integration properties
In the analytical solutions of the kinetic model (1) presented in literature, frequently the Goldstein Jfunction appears. ‘* The J function and the distributions hfl, hf,, hsf and h, are related as follows -
5 J(k,T, kz(t ;
7)) = hr/(T,t> + hfi(T,t),
J(h(t - 4, h4 = h&T, t>+ MT, t>,
(A31 (A41
(see also de Smedt & Wierenga3*). Using J(a, 0) = e-’ and J(O,b) = 1 we find exp(-kit)+
exp (-k2t)
: h,(q t) dr+ ji hd(T, t) dr= 1. (A6) J The determination of the integral of only a single h-function is less straightforward. Let H(t)
+
=
(A7)
Using (Al) and (A5) we find that H satisfies the differential equation i H = -(k,
+ k2)H + k2(l - emkl’).
(AS)
The solution of (AS) is
hn(T,t) dT =
k2
kl
k,
+
kl
e-(kl+k2)t_ e-klra (A9)
From (AS) it follows that
&
(1 _
e-(kl+W)a
(AlO)
In a similar way it can be shown that h,(T, t) dr = & 1
+ --& 2
h
+k2
e -(k,+k,)r _ e-kz’7
(All) h&T, t) dr = & 1
(1 - e-(klfk2)t). 2
6412)
206
J. J. A. van Kooten
APPENDIX B: MATHEMATICAL VERIFICATION OF THE SOLUTIONS OF THE KINETIC MODEL
Applying partial integration we find
Bl Verification of the solution presented in Section 3.1 It is obvious that C (9a) and S (9b) satisfy the initial conditions (7). Because C’ and C, have zero boundary conditions, C and S have zero boundary conditions too. To check whether the differential eqns (la) and (lb) are fulfilled we start with differentiating S to t dS - e- k,j - = -k2Sinit eCkzt + kl Cf at
Note that $_,-(0, t) = 0 and d;,(x, O&(0, kzSinit(x) eeklt. Substituting integration property
W) In (Bl) we used hfi(t, t) = kl exp (-&It) and h,(t, t) = 0. Applying differentiation property (A2) and ordening the terms that belong to C and S, one can see that (Bl) is equivalent to eqn (lb). For C we have on the one hand
E _ %f at -37
e-k,r
- klCfepklt + Cfhff(t,
t) + C&(t,
t)
B2 Verification of the solution presented in Section 3.2 If x E dMN then for any t, ~o(x, t) = Cb(x), so that (13) reduces to C(x, t) = C,(x) exp (-kl t)
kl h,&,t) +-Q&d)
P4 Since C, and C, are solutions of ADE (6) with y # 0 we have on the other hand MADE
MADcf
+
c-e
e-
klf
~{MAD~;(%T)~~((~.
acf
t)
+
MADCsb,+sf(T
-k,l
at
+
2
(x, +Q-(T, t) +
$j (x, T)~&T, t)}dT
f))
d7
t) =
(Al) in (B2) we obtain that (B2) + (Bl) = (B4), so that eqn (la) is satisfied too.
d7.
Since (kl/k2)hs. = hfs it follows from integration property (A5) that C(x, t) = &(x), so that at 8RN C satisfied the right boundary condition. In a similar way it can be shown that at aRo the solute flux equals F(x). To verify that C and S satisfy the differential equations of the kinetic model the steps of (Bl) can be repeated with C’ replaced by Co and Cs by (kl/kz)~o. To include the effect of the production term y, integration property (A5) must be applied in step (B3).
ELSEVIER
Copyright 0 1996 Elsevier Science Limited Printed in Great Britain. All rights reserved 0309-1708/96/$15.00+0.00
0309-1708(95)00045-3
A method to solve the advection-dispersion equation with a kinetic adsorption isotherm J. J. A. van Kooten* Department of Mathematics, Wageningen Agricultural University, Dreyenlaan 4, 6703 HA Wageningen, The Netherlands (Received 8 May 1995; accepted 23 November 1995)
This study deals with a method to solve the transport equations for a kinetically adsorbing solute in a porous medium with spatially varying velocity field and dispersion coefficients. Making use of the stochastic nature of a first-order kinetic p~:ocess,we show that the advection-dispersion equation and the adsorption isotherm can be decoupled. Once the solution for a non-adsorbing solute is known, the method provides an exact solution for the kinetically adsorbing solute. The method is worked out in four examples. In particular we demonstrate how the method can be applied simultaneously with a numerical transport code: the advective-dispersive transport is computed numerically, whereas kinetic effects are incorporated analytically. The proposed approach may be useful in field scale applications with complex flow patterns. Copyright 0 1996 Elsevier Science Limited Key worak kinetic adsorption, advection, dispersion, multi-dimensional varying flow, residence time distributions, analytical solutions.
1 INTRODUCTION
spatially
mathematics can be understood from the physics of the kinetic process. In this way the practical applicability of the kinetic model may be enlarged. The conventional advection-dispersion model for a solute subject to first-order reversible kinetics is
Solute transport in groundwater may be highly affected by interactions between the solute and the solid. Therefore transport models should contain equilibrium or non-equilibrium isotherms that account for these interactions. The most simple model is obtained by applying the linear equilibrium isotherm. Linear equilibrium adsorption only retards the transport of a solute. The local equilibrium assumption (LEA) may be valid if the rate of mass change due to the adsorption process is much faster than that (due to the flow process. If the LEA is not valid a kinetic model may give a better description of the solute transport. Several authors3T18,36 have presented criteria to decide whether an equilibrium or kinetic model is required. An equilibrium approach is often preferred above a kinetic approach; the retardation factor in the equilibrium model can be estimated more easily from experimental data than the kinetic reaction rates. Moreover, a kinetic model requires more complicated mathematics than the equilibrium model. One of the aims of this study is to show that the
( )
g+g=-$(UiC)+& I
X3 ,,=k,C-k,S,
I
+Y(X),
Dijg
J
(14 (lb)
where x is the vector of space coordinates, C(x, t) (M/L3) is the concentration of the solute in the free phase, S(x, t) (M/L3) is the concentration in the adsorbed phase, Y (L/T) is the fluid velocity vector, Dij (L2/T) is the dispersion tensor, y(x) (M/L3T) is the rate of zero-order production at the point x and kl (T-i) and k2 (T-l) are the forward and backward reaction rate, respectively. With respect to the indices we use the Einstein summation convention. Model (1) is mathematically equivalent to a dual porosity model, i.e. a model for transport in a medium in which two zones can be distinguished: a mobile zone in which flow and transport take place and an immobile zone (e.g. dead end pores) in which the water is nearly stagnant.”
*Currently at: Baan Research, P.O. Box 250,’6710 BG Ede, The Netherlands. 193
194
J. J. A. van Kooten
In the literature analytical solutions of eqn (1) are available for various initial and boundary conditions. Lindstrom and Narasimham26 have presented a onedimensional solution for a previously distributed contaminant. Solutions for migration through a semiinfinite column with different types of solute input at the inlet have been obtained by Ogata,27 Cameron and Klute,’ van Genuchten and Wierenga” and de Smedt solutions for a and Wierenga.32 Three-dimensional point spill of pollution in an infinite domain have been derived by Carnahan and Remer6 and Goltz and Roberts.13 Lassey23 has derived solutions for arbitrary solute input after t = 0. Also, recently solutions have been derived with zero-order production (y # 0); Torido et a1.35 have presented a wide class of onedimensional solutions assuming arbitrary input concentration, initial conditions and production functions. In all papers the velocity field was assumed to be uniform and one-dimensional. The dispersion coefficients were taken as constant. At field scale these conditions will often not be satisfied. All solutions were obtained with the Laplace transform method. From a mathematical point of view this method may be quite powerful. However, it does not always give an insight into the relationship between the structure of a solution and the underlying physical processes. Moreover, the Laplace method fails if model parameters are space dependent. Studies dealing with the transport of a kinetically adsorbing solute in a porous media with heterogeneous hydraulic and/or chemical properties have not primarily focused on the resident concentrations, but on the spatial moments of the solute plume1’7>29’37 or on the arrival time distribution at a plane perpendicular to the mean fluid flow.s>30>31~33~34 The objective of this study is to develop a method to solve the concentration, C, and S from model (1) in a multi-dimensional system with spatially varying w, D and y. The sorption rates kl and k2 are assumed to be constant. All coefficients are assumed to be timeindependent. We interpret the adsorption kinetics as a stochastic process in which a particle is alternately free or sorbed (Section 2). Keller and Giddings” have presented four conditional distributions for the residence time in the free phase during a time interval t. For arbitrary initial and boundary conditions we show that the solution of eqn (1) is a convolution of solutions of the classical ADE (la) without the adsorption term and a combination of these distributions (Section 3). So, in fact, we show that advectiondispersion eqn (la) and the reaction eqn (lb) can be decoupled. Analytical solutions of the classical ADE have been derived for numerous initial and boundary conditions.9124125Now, for all these problems we also have the solution for a kinetically adsorbing contaminant. Some examples are considered in Section 4. Furthermore, we demonstrate how the method may
be applied simultaneously with a numerical transport model. The transport equation is solved numerically, whereas the kinetic adsorption is incorporated analytically.
2 RESIDENCE TIME DISTRIBUTIONS IN THE FREE PHASE The kinetic adsorption process described by eqn (lb) can be considered as a stochastic process in which particles are alternately in the free and sorbed phase. Between transitions the residence times in the free (f) and sorbed (s) phase are exponentially distributed with mean l/k, and l/k2, respectively. During a time interval t a particle may make several transitions from the one state to the other. Assuming that at least one transition has taken place, Keller and Giddings” have presented probability density functions for the time fraction x (0 5 x < 1) that a particle has spent in the free phase (see also Giddings and Eyring”). These distributions depend on the initial and final state of a particle. For the purpose of this study we rewrite these distributions in terms of r (= xt ), the residence time in the fluid phase during a time interval t
$j-(T t) =
c
Co
(klkzT)“(t
-
+-’
e-(k,-k2)T-k2f
n!(n - I)!
n=l
=
A~(‘,
t)
=
kl
2 n=O
=
kl
(k1k2T)n(t
-
T)”
e-(k,-k2)T-k2f
(n!)2
e-_(kl-kzk-kzf
Zo(2&&2+
-
4,
(3)
where the first and second subscript of the functions h denote the initial and final state of a particle, respectively. In eqns (2) and (3) Z. and Zi denote the modified Bessel functions. For the distributions Zzf~and hq the condition that at least one transition has taken place is satisfied automatically. In the two other cases there is also the possibility that a particle remains in the initial state all the time; the probabilities that a particle spends time I in the free or sorbed phase are exp (-k, t) or exp (-k2t), respectively. Some integration and differentiation properties of the distributions are given in the Appendix A. In general, only the expressions in terms of the modified Bessel functions are given (see e.g. Valocchi and Quinodoz38). We have also given the summation representation because it is closer to the underlying
A method to solve the advection-dispersion equation
physical process: e.g. the counter n is the number of times a particle has returned to the initial state. Only in the free-p:hase particles are subject to advection and dispersion. In the sorbed phase particles are immobile. Because the adsorption rates are constant the above-described exchange process between the free and sorbed phase is not affected by the advectiondispersion process, even when the advection and dispersion parameters are space dependent. So, the adsorption process and the advection-dispersion process may be studied separately. Once the temporal evolution of the advective-dispersive displacement is known, the adsorption kinetics can be incorporated naturally by utilizing the fact that the residence time, 7, in the free-space is stochastic; with the aid of the distributions (2), (3), (4) and (5) the expected advective-dispersive displacement can be computed. In Section 3 we use this idea to decouple the transport eqn (la) and the adsorption isotherm (lb). Other authors have applied the same idea in different contexts. Valocchi and Quinodoz38 have developed an efficient random walk algorithm for simulating the spread of a kinetically adsorbing particles. In each time step At they use the above h-distributions to determine the fraction of time that a particle has spent in the free phase. A modification of the algorithm has been described by Andrecivic and Foufoula-Geourgiou.’ They show that for large k2t the distributions approach Gaussian distributions that can be fully characterised by the first and second moment of the residence time. In Quinodoz and Valocchi38 expressions are presented for the first and second spatial moment of an instantaneously injected contaminant in porous formations with randomly heterogeneous hydraulic conductivity. Given the spatial moments of a non-reacting contaminant, they use distribution (2) to derive an expression for the spatial moments of a kinetically adsorbing contaminant. Comparable results have been obtained by Cvetcovic and Shapiro8 (eqn (31)) and Selroos and Cvetcovic3’ (eqn (16)), who analyse the effect of kinetic adsorption on the arrival time distribution of a solute at a control plane in a heterogene0u.s porous formation.
3 DECOUPLING
OF THE TRANSPORT REACTION EQUATION
195
In the sequel we refer to MAD as the advectiondispersion operator. Let 0 (R c RN, N = 1,2 or 3) be the region in which contaminant transport is studied. The boundary dfl may consist of two disjunct parts. At the part dRD the concentration is prescribed. At the part dfi,v the mass flux is prescribed. In mathematical terms dRD is a Dirichlet boundary and d& is a Neumann boundary. A method to solve the kinetic model (1) for arbitrary initial and boundary conditions is developed in three steps. First we consider the case of non-zero initial conditions and zero boundary conditions. Next, we consider the opposite case: zero initial conditions and non-zero boundary conditions. Finally, the solution for arbitrary initial and boundary conditions is obtained by superposition of the two former cases. We also give a physical interpretation of the solutions. 3.1 Transport of a contaminant distributed at or before t=O
In this section we consider the transport of a contaminant which has been distributed at or before t = 0. It is assumed that the distributions of both C and S at t = 0 are known C(X, 0) = Ci”i,(x) and
S(X,
0) = Sinit(
(7)
Furthermore, we assume that no contaminant is added to or extracted from the system so that y = 0 and the boundary conditions are
aciD,
( C(x, t) = 0
at
[ni*(ui-Dij&)C=O
atdoN,
(8)
where n is the outward normal vector on 80. An expression for the solution of model (1) subject to these initial and boundary conditions can be found with physical arguments. We first give this expression. Next, we explain how it is obtained. C(x,
t)
=
C;(x,
t) e+’
S(X,
t)
=
Sinit
e+’
AND
In this section we show that the solutions C and S can be expressed as an integral of a proper combination of the h-distributions and C’, the concentration of a non-adsorbing contaminant satisfying
+ 1;
(c’(X,
t)
+ cs,(X,
~)h,(~>
t))
d7, (9b)
where both C’ and 6;, are solutions of ADE (6) (with y = 0) subject to the initial conditions
with
(6)
4&b,
C;(x7 O) =
Cinit(X)
and
C.f(x, O) =
&nitCX), (10)
and zero boundary conditions.
J. J. A. van Kooten
196
The physical reasoning behind expression (9) is as follows. During the time that particles are in the free phase their transport is governed by the classical ADE (6). The remaining time particles are immobile. The residence time in the free-phase r is stochastic. With the help of the residence time distributions frfl, hf,, h,f and h, the expected advective-dispersive displacement can be computed. The free phase transport of the contaminant with initial state ‘f ’is described by C;, and the free phase transport of the contaminant with initial state ‘s’is described by 6;,. To compute the concentration C in the free-phase at time t the distributions with final state ‘f'should be used and to compute the concentration S the distributions with final state ‘s’ should be used. The term exp (-k,t) in eqn (9a) represents the probability that particles have remained in the free-phase continuously. The term exp (+t ) in (9b) represents the probability that a particle has not moved at all. The mathematical verification is given in Appendix B. It is often possible to simplify expression (9). For example, if all contaminant is injected instantaneously at t = 0 (see Section 4.1) then S = 0, so that C;, = 0. In many applications it is assumed the initial distributions C and S are in equilibrium, which implies that Cs = (k,/k&“, so that only C; has to be determined. A one dimensional example of the latter case is discussed in Lindstrom and Narasimham.26
have the same structure. Generalizing the result to multi-dimensional systems with spatially varying velocity fields and dispersion coefficients and y # 0 we have found that the solution of (1) with zero initial conditions and boundary conditions (12) can be expressed as C(x, t) = C,(x, t) eCkir +~~~o(x,~){h~(T,t)+~h~~(~,t
(134
(13b) where Co is the concentration of a non-adsorbing contaminant satisfying ADE (6) with zero initial condition and boundary conditions given by eqn (12). The mathematical verification is sketched in Appendix B. Because contamination is constantly entering and leaving the porous medium the physical interpretation of expression (13) is less straightforward than that of (9). However, using the differentiation properties of Appendix A one may check that the time derivatives of (13) can be written as
a+, t) -=
aco;,(x, t)
e-k,r
at
3.2 Transport of a continuously injected contaminant If a domain is uncontaminated conditions are C(x, 0) = 0
and
at t = 0, the initial
S(x, 0) = 0.
(11)
We assume that from t = 0 there is a constant input and output of contaminant. This input/output may take place at the boundary of the domain and inside the domain, e.g. at an injection or a withdrawal well. Contaminant input/output at the boundary is described by at dRo,
( CC& t) = C,(x) ~n;‘(,iC-Dii~)
=F(X)
at dRN.
wx, 1) -= at
at
f
+
aco,(x, 4 hff(T,t>dT,
f
lb
f0
t aCo(x,T) h,dT,t) dT. 0 aT
(14b)
The physical interpretation of eqn (14) is straightforward. When particles are in the free phase the change in the concentration is described by eqn (6). Because the residence time in the free phase is stochastic, this change should be weighted with the proper residence time distribution. Because particles enter or leave the medium in the free phase, only the distributions with initial state ‘f' are relevant. The term exp (-k, t ) in eqns (13a) and (14a) again represents the probability that particles have remained in the free-phase continuously.
(12) 3.3 Arbitrary initial and boundary conditions
Production/withdrawal inside the domain is described by the term y in eqn (1). The functions C,, F and y may vary in space. One-dimensional examples of this case with y = 0 are discussed in 0gata,27 Cameron and Klute,’ van Genuchten and Wierenga” and de Smedt and Wierenga.32 We have rewritten their solutions in terms of the residence time distributions hff, hfi, hd and h,. For the solutions that are expressed in terms of the Goldstein J-function’2 we have applied relation A3 and A4 of Appendix A. All solutions turned out to
The solution of model (1) (with y # 0) subject to non-zero initial conditions c(x, O) =
Chit(X)
and
s(x, O) =
GdCX)7 (15)
and non-zero boundary conditions ( C(x, t) =
cb(x)
at f3RD,
(16)
A method to solve the advection-dispersion equation is the superposition of eqns superposition can be written as C(x,
t) =
Cl(x,
(9) and
(13). This
t) e+’
197
neither the kinetic nor the equilibrium model describe the transport of a solute satisfactorily.5 A combined equilibrium and kinetic approach may be required:
r +
J
,@x,+J+,t)
+cz(
x, 4
h&,
t))
d7, (174
S(X,
t) =
Sini,
e+’
+ 1; (c,
(x, +fi(~,
t) + z;(x,
%(7,
t))
d7, (17b)
where C, E Cf + C,, and Cz = Cs + (ki/kz)C,-,, so that Ci and Cz are concentrations of a non-adsorbing contaminant satisfying eqn (6) with y](x) = y(x) and y2(x) = (k,/k&y(x), respectively. The initial and boundary conditions are C‘1(x~O) =
Cinit(X)7
where k3 (T-l) and k4 (T-l) again denote a forward and backward sorption rate, respectively. After introduction of the retarded parameters 2rf = vi/R, DC = DijlR, 7R = y/R, SR = S/R and @ = kJ/R model (23) transforms into a model equivalent to (l), so that the method presented in this study also applies to the combined model. As a point of mathematical interest we mention that the decoupling method is also valid for models in which the advection-dispersion operator MAD (6) is replaced by an arbitrary linear differential operator M=a(x)+bi(x)~+Cij(X)~ I
(i,j,k, at &lo,
(19) A considerable simphfication is achieved if it is assumed that the initial concentration distributions are in equilibrium, i.e. = (kl lk2)Cinit
(20)
which implies that C2 := (ki/k#i to be determined.
so that only Ci has
3.4 Further generalisations If contaminant input only takes place during time z,,, i.e. for t 5 Ti:,, the boundary conditions and/or y are not equal to 0 and for t > T they are 0, then the solution of model (1) is Cr(x, t) = C(x, t)
and
ST(X,
for t 5 Tin,
t) =
S(X,
t)
CT(x, t) = C(X, t) - C(X, t - Tin) ST(X,
t) =
S(X,
t) -
S(X,
t -
Tin)
and
1
J
etc = 1.. .N),
where a, bi, cij, dijk etc. are scalar functions position x.
at a0,v.
&it
(22)
a2
at bRo,
f e’1 (x, t) = cb(x)
as ;;r;-=k3C-k4S
(21)
for t > Tiny
where C and S are given by eqn (17). It is well known that for large kl and k2 a kinetic model approximately predicts the same transport as a linear equilibrium model with retardation factor R = 1 + (k,/k2). On the other hand, it may occur that
(23)
of the
4 SOME EXAMPLES In Section 3 we have shown that once the transport of a non-adsorbing solute is known, then the transport of a kinetically adsorbing solute can be determined. In this section the theory is demonstrated in some examples. In the first example we show how a solution of a kinetic model presented by Carnaham and Remer6 and Goltz and Robertsi can be obtained as a special case of the method described in this study. In the second example the method is applied to predict the breakthrough of a contaminant at a well. We extend the results for the breakthrough of a non-adsorbing contaminant obtained in previous publications.‘4Y21 In examples 3 and 4 we demonstrate how the method can be combined with a numerical transport model. In all examples we assume that the dispersion tensor is given by Dij = aTIu[Sij + (aL - ar) z,
(24)
where aL and aT denote the longitudinal and transversal dispersivities.2 In example 2 we allow aL and aT to vary in space. The integrals with respect to r (see e.g. expressions (9), (13) (17) (27) and (31)) were computed with a Gauss-Kronrod rule of the IMSL-Math Library.” For the modified Bessel functions in the h-distributions we
J. J. A. van Kooten
198 have implemented Blair.4
the Chebyshev
approximations
of
3,5,
-NON-~RBING,S~
. ,
4.1 Example 1: 3-D transport of instantaneously injected solute In three dimensions the kinetic model for transport in an isotropic formation with uniform flow in the xl-direction reads
(254
G-4
At t = 0 an amount Q (M) of solute is injected in the origin (xi, x2, x3) = (0, 0, 0), so that the following initial conditions apply to eqn (25a). C(x, 0) = QW
and
S(x, 0) = 0,
KINETICALLYADSORBINGSOLUTE 3.
(25b)
where S(x) denotes the Dirac delta function. As boundary condition we require that C vanishes at infinity. The concentration of a non-adsorbing solute, i.e. kl = k2 = 0 in eqn (25) is given by
DH
I
3 022
I
x:
s. -8 2P 4 1.5 -
t t) epklt + o zf(x,
S(x, t) = 1 c’(x, $I&, J
J
0.5 o-
(see e.g. Carnahan & Remer6). Now, from eqn (9) we directly obtain the concentrations of the kinetically adsorbing solute in the free and sorbed phase C(x, t) = cf(x,
l-
033
(26)
T)~#(T, t) dr,
J
n- 2.5 -
P 8
(Xl -A2
1
-I
(b)
Fig. 1. Evolution in time of the concentration profiles c;- (eqn (26)), C and S (eqn (27)) along the xl-axis. (Q = 100(M), Dz2 = D33 = 0.05 (L2/T), wI = 1 (M/L), D,, = 0.25 (L’/T), k, = k2 = 0.1 (T-l)). 4.2 Example 2: Contamination of a pumping well
(27)
t) dr.
(Note that in this example CX= 0, because S = 0 at t = 0). Solutions equivalent to eqn (27) were presented by Carnahan and Remer6 and Goltz and Roberts.13 However, they obtained the solutions after lengthy Laplace transforms. Here, the solutions simply follow with a more general theory which, moreover, is closer to the underlying physical process. In Fig. 1 a comparison is made between the evolution of the concentration profiles C’, C and S along the xi-axis (x2 = x3 = 0). For small t the concentration profiles of C and S have a strong non-Gaussian shape. This can be explained from the fact that shortly after the injection the solute concentration close to the origin is high so that a large amount of solute will be adsorbed. The adsorbed amount is released slowly. Furthermore the figures clearly illustrate the enhanced longitudinal spreading and tailing due to kinetic adsorption.
In previous publications we have studied the arrival of conservative or linearly decaying contaminants at a well in a confined aquifer (e.g. van Herwaarden & Grasman,” van der Hoek,16 van Herwaarden,14 van Kooten,21’22). A typical flow pattern for groundwater discharge is shown in Fig. 2. Van Herwaarden and Grasman” have analysed the effect of the transversal dispersion on the arrival rate. If a contaminant is subject to first-order kinetics, then during intervals a particle is adsorbed and not displaced in the transversal direction. Thus, kinetic adsorption only affects the time at which particles pass a certain point and not the ultimate transversal spreading. Therefore, the approximation for the fraction of a contaminant that reaches the well presented by van Herwaarden14 (eqn (4.5)) may also be used for a kinetically adsorbing contaminant. In van Kooten22 we have derived approximations for the arrival time distribution at a well. For example, for particles released in the point x far inside the
A method to solve the advection-dispersion
equation
199 I
Fig. 2. Flow pattern for a well in a uniform background flow. catchment area, where the effect of the transversal dispersion is negligible, the arrival time distribution can be approximated by (28a) where p and rs2 are asymptotic approximations mean and variance of t:he arrival time
/+) =
Tadv(X)
ai,C-4
+ F(x)
for the
&(t) = s,(t) =
JR
conservative cases are Ci,it(x)g(x, t) dx
and (29)
R Git(X)g(X, t) dx,
J and, thus, the arrival time distribution for a kinetically adsorbing contaminant is g(t) = &(t) eFklr
and
+ 1; (&(+,j+,
(28b) 02(x) = 2 In (28b) Tad”is the advlective travelling time to the well and C is a coordinate along the streamline through x, see Fig. 3. Note that for any flow pattern p and c? can easily be determined by particle tracking. For a detailed description of how to take into account transversal dispersion we refer to van Kooten.22 With the residence time distributions of Section 2 kinetic adsorption can be incorporated without making any further approximation. It is assumed that in an aquifer Q at t = 0 the concentrations of the contaminant in the free and sorbed phase are given by Ci,i,(X) and sinit( The arrival time distributions of the
Fig. 3. The coordinate C is taken along the streamline through
the point x.
corresponding
For an instantaneous reduces to
t) +&W&
t)) d7
(30)
spill in the point x eqn (30)
g(x, t) = 2(x, t) eCklt + : g(x, T)hfl(T, t) dz (31) J We have worked out the above for the special case of a well in a uniform background flow (Fig. 2). The corresponding velocity field is given by 211(x,,x2) = 1 - 2 x: + x; ’
-x2 V2(Xl,X2)
=
-. xy + x; (34
We assume that at t = 0 in the point (x1,x2) = (-4,l) an instantaneous spill of contaminant takes place. For various values of kl and k2 the approximation for the arrival time distribution (31) at the well is displayed in Fig. 4. The dispersivities are aL = 0.1 and aT = 0.01. The figures show that the tail due to kinetic effects becomes more pronounced if k, increases and/or k2 decreases, i.e. if the expected time of one stay in the fluid phase (l/k,) decreases and/or the expected time of one stay (1 /k2) in the fluid phase increases. To assess the accuracy of the approximations we have made a comparison with the results of random walk simulations. There exist various random walk algorithms
J. J. A. van Kooten
200
AI~~~~TIMEDISI’RIBUTIONSAT~EI~.
~~~DISIRIBUTIOI’JSAT~~ 0.7b ( 0.6 .
,
0.6 -
wa4dsorbingaln~nant
nm-ads&ingamtaminant
8
30
35
40
tO
(a)
0.7c
Oo
45
5
10
(b)
15
20 25 t0.J
30
35
40
45
Fig. 4. Some arrival time distributions at the well for a contaminant spilled in (-4,l).
In (a) the adsorption rate ki (T-l) is kept fixed. In (b) the desorption rate k2 (T-i) is kept fixed.
that take into account kinetic adsorption.‘>38 Here, it suffices to select the most natural and straightforward method (method 1 of Valocchi and Quinodoz38). In the free-phase the stochastic displacement of a particle is governed by AX=
I3 wl+-&,+-~i2 ax (
d ay
At >
+~~$z,&t+&&$$zr&t v2+&D2i
AY= (
+2022 8Y
At >
+~~$z&t-~z+$z*6t (33)
where At is the time step and Z, and ZT are random deviates from a uniform distribution with mean 0 and variance 1 (see Kinzelbach20 or van Kooten*‘). First the time r- of one stay in the fluid phase is generated from an exponential distribution. During time r- a particle is advected-dispersed according to eqn (34). Next, from an exponential distribution the time rs that a particle stays in the adsorbed phase is generated. The time T$ is added to t, the total time elapsed from the beginning of the simulation. This process is repeated, until the particle enters the well. we carried out 15,000 random From the point (-4,l ) and k2 = 0.2 (T’). The time walks with ki = O-5 (T- 1’ step At was chosen to be 0.001 except if rf was exceeded. In the latter case At was taken equal to the difference between T- and the time that has been spent in the
COMPARISONOF RANDOM WALK AND ASYMFIWI’IC
APPROXIMATION
42 -
! OSB * &
1
O-l0.05 -
OO
Fig. 5. Comparison of the arrival time distributions obtained with the random walk simulation and computed from eqn (31).
A method to solve the advection-dispersion equation
fluid phase yet. From Fig. 5 we may conclude that approximation (31) agrees quite well with the arrivaltime distribution obtained from the simulations. The small differences are mainly caused by the irregular behavior of the simulat:ions, and it is to be expected that the agreement will still become better if the number of random walks is increased.
201
transport is computed numerically whereas the kinetic adsorption is incorporated analytically (in the sequel we refer to the latter one as the mixed numericalanalytical solution). We study the lD-transport of a solute through a semi-infinite column. From t = 0 to t = Tin the concentration at the inlet x = 0 is maintained at some constant value C(0, t) = Cb. Next no more solute enters the domain, i.e. C(0, t) = 0 for t = Ti:,n. Ogata and Banks** and Ogata27 have presented solutions for a non-adsorbing and kinetically adsorbing solute, respectively. The solution for the kinetically adsorbing species can also be obtained by substituting the solution for the non-adsorbing species in expression (13). With the FE-code COLTRAP3’ we have simulated the transport at the interval 0 I x 5 1 with I = 100(L), w = 1 (L/T), D = 1 (L*/L), Cb = 1 (M/L) and Tir, = 10(T)). At the endpoint x = 1 a no-flux condition is imposed. The fact that the interval is finite does not undermine the comparison with the analytical solution for the semi-infinite system; van Genuchten and Alves’ have shown that for Peclet numbers Pe = vi/D greater than 20 the solutions for the finite and infinite system deviate only close to the endpoint. We have run COLTRAP with various values of the
4.3 Example 3: Comparison between an analytic solution,
a numerical solution and a mixed numerical-analytical solution Currently, there exist many finite element (FE) and finite difference (FD) codes to solve the ADE for a conservative contaminant. The aim of the examples, 3 and 4, is to demonstrate blow the theory of Section 3 can be used simultaneous1.y with such a code. Once the transport of a conservative contaminant has been computed numerically, with some post processing the transport of a kinetical1.y adsorbing contaminant can be computed. In this example we make a comparison between an analytical solution of the kinetic model, a full numerical solution and a solution in which the advective-dispersive TIME PROPILE IN x=30
aes.
.
.
TIME PROFILE IN x=30 . . . . .
.
.
= aMlyticai dution ._______ = numerical 8ohtion .... .........
= mixed solution
k,=O.l(ur)
xl
40
60
a0
100
12-O 140 160
160 xl0
tm
tO
(a)
TIME PROFILE IN x=30
.
ai2.
= mixed solution
._______= numerical solution
k2 = 0.02 (l/T)
k, = O.&i (1JQ k2 = 0.02 (l/l)
aaz.
Oo#2040 (4
.
= mixed solution
ao2. I-J.
.
= analytical solution 0.1.
._______ = numericalsolution
A, = 0.M (lm
TIME PROFILE IN x=30 . _ . .
ai2.
analytical solution 0.1 -
-
~==O.Ol(Lm
-
. 60
. eo
-
100
.
l#)
.
.
140 160
.
I
180 #r)
tc0
Fig. 6. Analytical solutions2’, numerical solutions and mixed numerical-analytical solutions for transport through a semi-infinite column. We have displayed the time profiles of the free phase concentration C (M/L) in the point x = 30. (w = 1 (L/T), D = 1 (L/T’), T, = IO(T), gridsize Ax = 2.0(L)). The numerical simulations of (a) and (c) were carried out with time step At = 1.0. In (b) and (d) At = 0.5.
J. J. A. van Kooten
202
flow takes place. Initially the aquifer is unpolluted. From t = 0 at the top of the aquifer for x between 7 and 10 m and y between - 1 and 1 m contaminated water constantly leaches into the aquifer. The concentration of the contaminant in the incoming water is C, = 100 pg/m3. At the remainder of the top surface C, = 0. Because of the symmetry in the y-direction it suffices to solve the transport problem in only one half of the aquifer, with a no-flux condition at the surface y = 0. At the downstream boundary x = 40 we impose a zero concentration condition. The surfaces x = 0, y = 0 and z = 0 are no-flux boundaries. The aquifer was covered by 40 x 20 x 12 = 9600 cubic elements. Because of the zero initial conditions the theory for simulating the transport of a kinetically adsorbing contaminant is described in Section 3.2. With the METROPOL code the concentration profile C0 of a non-adsorbing contaminant was computed from t = 0 to t = 20 days with time step At = 0.2 days and dispersivities aL = 0.6 m and aT = 0.3 m. The computation time at a Sun3 workstation was about 6 h. In Fig. 8(a) the flow pattern and the concentration contours in the y = 0 plane at t = 20 are shown. For each nodal point the time profile of Co was stored in a data file. The time profiles were postprocessed by substitution in eqn (13). In between partition point ti = iAt the concentration was determined by linear interpolation. The concentration contours for a kinetically adsorbing contaminant with kl = 0.1 and k2 = 0.1 (days-‘) are displayed in Fig. 8(b). The post processing could be carried out in about 2min, so that the costs of the post processing are almost negligible compared to the cost of solving the ADE numerically. In Fig. 8(c) we have drawn the concentration contours resulting from an equilibrium model with retardation factor R = 1 + (kl/k2) = 2. The solution of the equilibrium model was obtained from the METROPOL output by replacing the time t by the retarded time tR = t/R. Compared to the kinetic model the downstream spreading is clearly underestimated. The ‘kinetic’ profile is more flat than the ‘equilibrium’ profile. So, the application of a kinetic or equilibrium model may have large consequences on, for example, the delineation of
gridsize Ax and timestep At. We have computed the transport of non-adsorbing and kinetically adsorbing contaminants from t = 0 to t = 200(T). The time profiles of the concentrations in the point x = 30 were stored in data files. The time profile of a non-adsorbing contaminant was post-processed by substitution in eqn (13). In between partition point ti = iAt the concentration was determined by linear interpolation. We applied applied expression (21) to obtain a solution for pulse input. From the simulation results it (perhaps surprisingly) turned out that the accuracy of the full numerical and the mixed numerical-analytical solution is of the same order. Some examples are shown in Fig. 6. One might expect that the mixed solution is slightly better because no error is introduced due to the time discretization of the adsorption isotherm (1 b). However, this error appears to be dominated by the error due to the space-time discretization of the parabolic ADE (la). The main advantage of the mixed (or decoupling) method is that for different adsorption rates only the post-processing has to be repeated, whereas for the full numerical method the simulation code has to be rerun each time. So the mixed method is computationally much more efficient.
4.4 Example 4: Simultaneous application with a numerical model for transport in spatially varying flow In this example we demonstrate how the method presented in this paper can be applied to more complicated flow patterns. The numerical computations were carried out with the FE code METROPOL, which has been developed by the National Institute of Public Health and Environmental Protection in the Netherlands. METROPOL can only deal with non-adsorbing or linearly decaying solutes. With our method the effect of kinetic adsorption is incorporated. We consider contaminant transport in an aquifer which is recharged at the top (Fig. 7). From below the formation is confined by an impervious layer. The effective velocities (m/day) of the incoming and outgoing flow are indicated in the figure. In the y-direction no
V=O.l
V = 0.6 -t
10m I
z
Y k
40m
x Fig. 7.
Sketch of the aquifer considered in Section 4.3.
J. J. A. van Kooten
204
LZ~.~~) or cannot be derived at all. Therefore, the method will probably find its major application when it is used in combination with a numerical transport model. With such a model the transport of a non-adsorbing contaminant is computed first. Next, by substitution of the output in expressions (9), (13) or (17) the effect of the kinetic adsorption can be incorporated analytically. Examples are given in Sections 4.3 and 4.4. The proposed approach may be useful in field-scale applications with complex flow fields. The main advantage above solving the complete model (1) numerically is that no new run with the numerical code has to be performed when the values of the reaction rates k, and k2 are changed. Only the post-processing has to be repeated. This feature may be cost-effective when the model output has to be fitted with experimental data. Another advantage is that the method is computationally more efficient. The costs of solving one differential equation numerically are lower than the cost of solving two coupled equations. Although in the mixed numerical-analytical method no errors are introduced due to time discretization of the reaction isotherm (lb), example 3 (Section 4.3) revealed that a mixed solution has the same order of accuracy as a full numerical solution; the dominating error stems from the space-time discretization of the ADE. A disadvantage of the approach may be that the time profiles of the concentrations in the nodal points have to be stored, so that for large grids much computer memory is required. In practice this disadvantage is not too large, because one is often only interested in the breakthrough of a contaminant at some particular points. In this case, just as for the traditional numerical method, only in the observation points do the time profiles need to be stored. ACKNOWLEDGEMENTS
residence time distributions. Stoch. Hydrol. Hydraul., 5 (1991) 155-71.
2. Bear, J. & Verruyt, A., Modelling Groundwater Flow and Pollution. Reidel, Dordrecht, 1987. 3. Bahr, J. M. & Rubin, J., Direct comparison of kinetic and local equilibrium formulations for solute transport affected by surface reactions. Water Resour. Res., 23(3) (1987) 438-52.
4. Blair, J. M., Rational Chebyshev approximations for the modified Bessel functions IO(x) and 11(x). Maths Comp., 28 (1974) 581-3.
5. Cameron, D. R. and Klute, A., Convective-dispersive solute transport with a combined equilibrium and kinetic adsorption model. Water Resour. Res., 13 (1977) 183-8.
and 6. Carnahan, C. L. & Remer, J. S., Nonequilibrium equilibrium sorption with a linear sorption isotherm during mass transport through an infinite porous medium: some analytical solutions. J. Hydrol., 73 (1984) 227-58.
7. Chrysikopolous, C. V., Kitanidis, P. K. & Roberts, P. V., Generalized Taylor-Aris moment analysis of the transport of sorbing solutes through porous media with spatiallyperiodic retardation factor. Trans. Porous Media, 7 (1992) 163-85.
8. Cvetkovic, V. D. & Shapiro, A. M., Mass arrival of sorptive solute in heterogeneous porous media. Water Resour. Res., 26(9) (1991) 2057-67.
9. van Genuchten, M. Th. & Alves, W. J., Analytical solutions of the one-dimensional convective-dispersive solute transport equations. Tech. Bull. 1661, US Department of Agriculture, Washington, DC, 1982, 151 pp. 0. van Genuchten, M. Th. & Wierenga, P. J., Mass transfer studies in sorbing porous media. 1. Analytical solutions. Soil Sci. Sot. Am. J., 40 (1976) 473-80.
1. Giddings, J. C. & Eyring, H., A molecular dynamic theory of chromatography. J. Phys. Chem., 59 (1995) 416-21. 2. Goldstein, S., On the mathematics of exchange processes in fixed columns. 1. Mathematical solutions and asymptotic expansions. Proc. R. Sot. London, Ser. A, 219 (1953) 151-71.
3. Goltz, M. N. & Roberts, P. V., Three-dimensional solutions for solute transport in an infinite medium with mobile and immobile zones. Water Resour. Res., 22(7) (1986).
The author thanks M. de Gee, J. Grasman, P. A. C. Raats and S. E. A. T. M. van der Zee for valuable discussions and for reading the paper critically. The author is grateful to 0. A. Buse for his help with drawing the figures and also thanks R. van de Weerd and the Institute of Public Health and Environmental Protection (in particular A. Leijnse) for their permission to use the COLTRAP and METROPOL code. I also would like to thank the reviewers for their constructive remarks and their suggestions for improving my paper. The work was funded by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Dutch Organisation for Scientific Research). REFERENCES 1. Andricevic, R. & Foufoula-Georgiou, E., Modeling kinetic
non-equilibrium using the first two moments of the
14. van Herwaarden, 0. A., Spread of pollution by dispersive groundwater flow. SIAM J. Appl. Math., 54 (1) (1994) 26-41. 15. van Herwaarden, 0. A. & Grasman, J., Dispersive groundwater flow and pollution. Math. Model Meth. Appl. Sci., 1 (1991) 61-81. 16. van der Hoek, C. J., Contamination of a well in a uniform background flow. Stoch. Hydrol. Hydraul., 6 (1992) 191-208.
17. IMSL. User’s Manual Math/Library, Fortran Subroutines for Mathematical Applications, IMSL, Houston, Texas, 1987. 18. Jennings, A. A. & Kirkner, D. J. Instantaneous equilibrium approximation analysis. J. Hydraul. Engng, 110 (1984) 1700-1717. 19. Keller, R. A. & Giddings, J. C., Multiple zones and spots in chromatography. J. Chromatog., 3 (1960) 205-20. 20. Kinzelbach, W., The random walk method in pollutant transport simulation. In Groundwater Flow and Quality Modeling, eds E. Custodio et af. Reidel, Dordrecht, 1988, pp. 227-45. 21. van Kooten, J. J. A., Groundwater contaminant transport
A method to solve the advection-dispersion including adsorption amd first-order decay. Stoch. Hydrol. Hydraul., 8 (1994) 185;-205. 22. van Kooten, J. J. A., An asymptotic method to predict the contamination of a pu.mping well. Adv. Water Resour., 18 (5) (1995) 295-313. 23. Lassey, K. R., Unidimensional solute transport incorporating equilibrium an’d rate-limited isotherms with firstorder loss. 1. Model conceptualizations and analytical solutions. Water Resour. Res., 24(3) (1988) 343-50. 24. Leij, F. J. t Dane, J. H., Analytical solutions of the onedimensional advection equation and two- or threedimensional dispersion equation. Water Resour. Res., 26 (1990) 1475-82. 25. Leij, F. J., Skags, ‘T. H., van Genuchten, M. Th., Analytical solutions Ifor transport in three-dimensional semi-infinite porous media. Water Resour. Res., 27 (1991) 2719-34. 26. Lindstrom, F. T. & Narasimham, M. N. L. Mathematical theory of a kinetic model for dispersion of previously distributed chemicals in a sorbing porous medium. SIAM J. Appl. Math., 24 (1973) 496-510. 27. Ogata, A., Mathematics of Dispersion with Linear Adsorption Isotherm. US Geol. Surv. Prof. Pap, 411-H, 1964, 9pp. 28. Ogata, A. & Banks, R. B., A Solution of Differential Equation of Longitudinal Dispersion in Porous Media. US Geological Survey Professional Paper, 41 l-A, Al-A9, 1961. 29. Quinodoz, H. A. M. &: Valocchi, A. J., Stochastic analysis of the transport of kinetically sorbing solutes in aquifers with randomly heterogeneous hydraulic conductivity. Water Resour. Res., 2!)(g) (1993) 3227-40. 30. Sardin, M., Schweich, D., Leij, F. J. & van Genuchten, M. Th., Modeling the nonequilibrium transport of linearly interacting solutes in porous media: a review. Water Resour. Res., 27(9) (1991) 2287-307. 31. Selroos, J. & Cvetkovic, V., Modeling solute advection coupled with sorption kinetics in heterogeneous formations. Water Resour. Res., 28 (5) (1992) 1271-78. 32. Smedt, F. de & Wierenga, P. J., A generalized solution for solute flow in soils with mobile and immobile water. Water Resour. Res., lS(5) (1979) 1137-41. 33. Jury, W. A., Sposito, G. and White, R. E., A transfer function model of solute transport through soil, 1. Fundamental concepts. Water Resour. Res. 22(2) (1986) 243-37. 34. Sposito, G., White, R. E., Darrah, P. R. & Jury, W. A., A transfer function model of solute through soil. 3. The advection-dispersion equation. Water Resour. Res., 22(2) (1986) 255-62. 35. Torido, N., Leij, F. J. & van Genuchten, M. Th., A comprehensive set of .analytical solutions for nonequilibrium solute transport with first-order decay and zero-order production. Water Resour. Res., 29 (1993) 2167-82. 36. Valocchi, A. J., Validity of the local equilibrium assump tion for modeling sorbing solute through homogeneous soils. Water Resour. Res., 6 (1985) 808-20. 37. Valocchi, A. J., Spatial moment analysis of the transport of kinetically adsorbin solutes through stratified aquifers. Water Resour. Res., 2!i (2) (1989) 273-9. 38. Valocchi, A. J. & Quinodoz, H. A. M., Application of the random walk method to simulate the transport of kinetically adsorbing solutes. In: Groundwater Contamination: Proceedings of the Third ScientiJic Assembly of the ZAHS, ed. L. M. Abriola. Baltimore, 1989, pp. 35-42. 39. van de Weerd, R. Ct Leijnse, A., Colloid facilitated contaminant transport in porous media. Model verification, calibration and sensitivity analysis. In preparation (1995).
equation
205
APPENDIX A: SOME MATHEMATICAL PROPERTIES OF THE DISTRIBUTIONS hff, 5,
kf AND ASS
Al Differentiation properties By differentiation
3
dt
of (2), (3), (4) and (5), we obtain
= klhff - kzhfi
and
%
= klha - kzh,. (A21
A2 Integration properties
In the analytical solutions of the kinetic model (1) presented in literature, frequently the Goldstein Jfunction appears. ‘* The J function and the distributions hfl, hf,, hsf and h, are related as follows -
5 J(k,T, kz(t ;
7)) = hr/(T,t> + hfi(T,t),
J(h(t - 4, h4 = h&T, t>+ MT, t>,
(A31 (A41
(see also de Smedt & Wierenga3*). Using J(a, 0) = e-’ and J(O,b) = 1 we find exp(-kit)+
exp (-k2t)
: h,(q t) dr+ ji hd(T, t) dr= 1. (A6) J The determination of the integral of only a single h-function is less straightforward. Let H(t)
+
=
(A7)
Using (Al) and (A5) we find that H satisfies the differential equation i H = -(k,
+ k2)H + k2(l - emkl’).
(AS)
The solution of (AS) is
hn(T,t) dT =
k2
kl
k,
+
kl
e-(kl+k2)t_ e-klra (A9)
From (AS) it follows that
&
(1 _
e-(kl+W)a
(AlO)
In a similar way it can be shown that h,(T, t) dr = & 1
+ --& 2
h
+k2
e -(k,+k,)r _ e-kz’7
(All) h&T, t) dr = & 1
(1 - e-(klfk2)t). 2
6412)
206
J. J. A. van Kooten
APPENDIX B: MATHEMATICAL VERIFICATION OF THE SOLUTIONS OF THE KINETIC MODEL
Applying partial integration we find
Bl Verification of the solution presented in Section 3.1 It is obvious that C (9a) and S (9b) satisfy the initial conditions (7). Because C’ and C, have zero boundary conditions, C and S have zero boundary conditions too. To check whether the differential eqns (la) and (lb) are fulfilled we start with differentiating S to t dS - e- k,j - = -k2Sinit eCkzt + kl Cf at
Note that $_,-(0, t) = 0 and d;,(x, O&(0, kzSinit(x) eeklt. Substituting integration property
W) In (Bl) we used hfi(t, t) = kl exp (-&It) and h,(t, t) = 0. Applying differentiation property (A2) and ordening the terms that belong to C and S, one can see that (Bl) is equivalent to eqn (lb). For C we have on the one hand
E _ %f at -37
e-k,r
- klCfepklt + Cfhff(t,
t) + C&(t,
t)
B2 Verification of the solution presented in Section 3.2 If x E dMN then for any t, ~o(x, t) = Cb(x), so that (13) reduces to C(x, t) = C,(x) exp (-kl t)
kl h,&,t) +-Q&d)
P4 Since C, and C, are solutions of ADE (6) with y # 0 we have on the other hand MADE
MADcf
+
c-e
e-
klf
~{MAD~;(%T)~~((~.
acf
t)
+
MADCsb,+sf(T
-k,l
at
+
2
(x, +Q-(T, t) +
$j (x, T)~&T, t)}dT
f))
d7
t) =
(Al) in (B2) we obtain that (B2) + (Bl) = (B4), so that eqn (la) is satisfied too.
d7.
Since (kl/k2)hs. = hfs it follows from integration property (A5) that C(x, t) = &(x), so that at 8RN C satisfied the right boundary condition. In a similar way it can be shown that at aRo the solute flux equals F(x). To verify that C and S satisfy the differential equations of the kinetic model the steps of (Bl) can be repeated with C’ replaced by Co and Cs by (kl/kz)~o. To include the effect of the production term y, integration property (A5) must be applied in step (B3).