Transport in Porous Media 18: 65-85,1995.
@ 1995 Kluwer Academic Publishers. Printed in the Netherlands.

65

Approximate Analytical Solutions for Solute
Transport in Two-layer Porous Media*
F. J. LEIJ and M. TH. VAN GENUCHTEN
U.S. Salinity Laboratory, USDA-ARS, 4500 Glenwood Drive, Riverside, CA 92501, U.S.A.
(Received: 18 January 1994; in final form: 10 May 1994)

Abstract. Mathematical models for transport in layered media are important for investigating how
restricting layers affect rates of solute migration in soil profiles; they may also improve the analysis
of solute displacement experiments. This study reports an (approximate) analytical solution forsolute transport during steady-state flow in a two-layer medium requiring continuity of solute fluxes
and resident concentrations at the interface. The solutions were derived with Laplace transformations
making use of the binomial theorem. Results based on this solution were found to be in relatively good
agreement with those obtained using numerical inversion of the Laplace transform. An expression
for the flux-averaged concentration in the second layer was also obtained. Zero- and first-orderapproximations for the solute distribution in the second layer were derived for a thin first layer
representing a water film or crust on top of the medium. These thin-layer approximations did not
perform as well as the ‘binomial’ solution, except for the first-order approximation when the Peclet
number, P, of the first layer, was low (P < 5). Results of this study indicate that the ordering of
two layers will affect the predicted breakthrough curves at the outlet of the medium. The two-layer
solution was used to illustrate the effects of dispersion in the inlet or outlet reservoirs using previously
published data on apparatus-induced dispersion.
Key words: Solute transport, composite media, boundary conditions.

1. Introduction
Knowledge of solute transport through composite or layered porous media is ofimportance to better manage and describe the movement of chemicals in naturaland artificial media. Interest in porous media transport is increasingly motivatedby concerns over the presence of a wide variety of chemical substances and wastes
in the subsurface environment. Mathematical models are necessary to assess thefate and movement of such chemicals.
Porous media are seldom homogeneous and the transport properties of these
media will vary spatially and sometimes also temporally. Accurate mathematicalanalyses of transport in heterogeneous media are not easily carried out. However,formulation and mathematical solution of the transport problem becomes possible
if the medium is assumed, somewhat simplistically, to be composed of a seriesof homogeneous layers. In soil science, composite media have been used for
representing stratified soil profiles in which horizons parallel to the soil surface* The U.S. Government right to retain a non-exclusive, royalty free licence in and to any copyrightis acknowledged.

66

F. J. LEIJ AND M. TH. VAN GENUCHTEN

form layers with different transport properties (Selim et al., 1977; Jacobsen etal., 1991). Nakayama et al. (1984) investigated the movement of radionuclides
through a medium consisting of soil and granite. Layered media are often also
created artificially to slow down or prevent chemical movement (clay liners).
Solute transport in porous media has traditionally been described with the deterministic advection-dispersion equation (ADE). Taylor (1953) demonstrated that theconcept of dispersion, although physically quite different, is mathematically equivalent to Fickian diffusion for sufficiently large travel times. This does not always
hold and dispersion will not follow traditional Fickian behavior near interfaces orboundaries, the use of a macroscopically constant dispersion coefficient is formaliyincorrect in such cases (Dagan and Bresler, 1985). Given the lack of alternative
concepts of dispersion that can be conveniently used in relatively simple transportmodels, the ADE with distinct dispersion coefficients for each layer still seems
attractive for modeling transport in composite media (Parker and van Genuchten,1985). This is particularly true if the ADE is used in a consistent manner for data
analysis or prediction purposes.
Mathematical descriptions of transport in layered media have also been employed to evaluate the appropriateness of inlet and outlet boundary conditions forhomogeneous systems. For example, the boundary conditions suggested by Danckwerts (1953) follow from the more general conditions formulated by Wehner andWilhelm (1956) by ignoring dispersion in the influent and effluent reservoirs. Theseassumptions can most easily be evaluated if appropriate mathematical solutions
and parameter values are available for multi-layer transport. Porous media may beviewed as assemblies of ‘independent’ homogeneous layers. Such a concept has
been helpful for formulating interface conditions for transport during steady flow_
perpendicular to the layers. For each layer, a first- or third-type condition is used forthe upper boundary while the lower boundary condition is formulated by means of azero gradient at the outlet or at infinity (e.g., Shamir and Harleman, 1967; Al-Niami
and Rushton, 1979). This approach usually assumes that the influent concentration
for the second layer follows directly from the concentration predicted at the outletof the first layer. The concentration in any layer is now independent of the transport
properties of all downstream layers; the mathematical solution procedure is thengreatly simplified while the problem of formulating more complicated alternativeinterface conditions is circumvented.
An alternative condition at the interface, first used by Wehner and Wilhelm(1956), requires that both the solute flux and the concentration be continuous.
Conditions of this type have been routinely used for heat flow problems involving
composite media(Carslaw and Jaeger, 1959; Gzigik, 1980). Given the similar natureof the ADE and the diffusion or heat flow equation, it is intuitively appealing toformulate the interface condition for solute transport in this manner although werealize that no ‘perfect’ boundary conditionsexist. Despite additional mathematicalcomplications, several authors seem to prefer the combined interface conditions(e.g., van der Laan, 1958; Kreft, 1981b; Barry and Parker, 1987). Such conditions

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

67

may also be appropriate for the inlet and outlet conditions of homogeneous porous
media. Unfortunately, the concentration in influent and effluent reservoirs cantypically not be accurately characterized as a function of time and/or position
except for some special cases (Novakowski, 1992b).
Analytical solutions for one-dimensional transport in composite media areoften derived with Laplace transforms (Carslaw and Jaeger, 1959) and sometimes
with Green’s functions, adjoint solution techniques, and finite integral transforms(Mikhailov and Gzigik, 1984). Solutions can be readily obtained for an arbitrary
number of layers if each layer is viewed as being an effectively semi-infinite medium. The use of Laplace transforms becomes more complicated if the concentrationof a certain upstream layer depends on properties of its downstream layers. Thissituation arises when both concentration and solute flux are required to be continuous at the interfaces. The resulting mathematical problem may be simplified byconsidering only two layers, or by limiting the solution to only the steady-statecase (Kreft, 1981 b). Frequently, the Laplace transform has been inverted numerically. For instance, Barry and Parker (1987) obtained flux-averaged concentrationsin this manner, while Leij et al. (1991) used this approach to predict volumeaveraged concentrations. Novakowski (1992a) similarly used numerical inversion
to predict concentrations in the porous medium as well as in the upstream anddownstream reservoirs. The two- or multi-layered transport problem can also be
conveniently analyzed by means of transfer functions or temporal moments (cf.
Kreft, 198 1 a; Barry and Parker, 1987). Transfer function in the regular time domain(i.e., travel time distributions) could be obtained experimentally, or mathematically by using a Dirac input boundary conditions and Duhamel’s theorem. Jury andUtermann (1992) formulated a joint probability density function for the travel timethrough a two-layer soil by assuming that the travel times in the individual
layers were either perfectly correlated or independent. These authors exploredthe influence of transverse variations in transport due to different paths for waterflow.
To the best of our knowledge, no explicit solutions are available for the ADE
requiring simultaneous continuity in flux and concentration at the interfaces.
The objective of this study is to derive an (approximate) solution for transport in a two-layer medium during steady flow using the method of Laplace transforms. The solution will be evaluated through comparisons with results obtainedby numerical inversion of the Laplace transform. Alternative approximate solutions may be appropriate in some cases. An approximation for the case of a thinfirst layer will be presented. Several examples of resident concentration profiles
versus position or time, and flux-averaged concentration distributions versus
time, will be presented for cases where the two-layer solution may clarify
certain theoretical and experimental aspects of solute transport in layered porous
media.

F. J. LEIJ AND M. TH. VAN GENUCHTEN

68
2. Problem Formulation

The porous medium is assumed to consist of two homogeneous layers subject tosteady water flow perpendicular to the layer interface. The transport and flow properties of both layers are macroscopically uniform in time and space, while the ADE
is assumed to describe transport in each layer. These simplifying assumptions are
made to facilitate the derivation of an analytical solution. The governing transport
equations are
O<X<L,

ac,a%5
ac2- ~237
at = D2 dx2

L<x<co,

t>o,

t>o,

(2)

where C is the volume-averaged (resident) solute concentration, t is time, x is
distance in the direction of flow, L is the position of the interface, D is the dispersion
coefficient, u is the mean pore-water velocity, while the subscripts 1 and 2 referto the first and second layer, respectively. The second layer is chosen to be semiinfinite for a convenient formulation of the outlet condition, this formulation canalso be applied to finite systems.
The partial differential equations are augmented by a zero initial condition,a third-type inlet condition involving a step input, first- and third-type interface
conditions, and a zero gradient at infinity as follows(3)

Cl(X, 0 ) = C2(x, 0 ) = 0,

ac1

u,C1 - DIda:

z=o+

7

(4)

Cl Id- = C21z=Lt,

(5)

-

&DI

(6)

=

(7)

e,u~Cl
dC2
dX

0,

Z-00

where B denotes the volumetric water content and Co is the concentration of
the influent solution at x = 0. As previously stated, conditions (4) and (7) areapproximate because of difficulties to characterize the concentration in the influent
and effluent regions.
The third- or flux-type condition is selected on physical grounds to ensure
mass conservation whereas the first- or concentration-type condition is invoked on
intuitive grounds (Kreft, 1981b). A combined first- and third-type condition at the

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

69

interface has been widely employed, on physical grounds, for heat conduction incomposite media with perfect thermal contact between layers. The analogy between
the governing equations for solute transport and heat conduction, and the fact thatsolute diffusion is one of the processes contributing to solute dispersion justify theuse of a first-type condition at the interface. Note that Equation (6) can also be castin terms of a second-type condition by utilizing continuity in concentration andwater flow (i.e., Ct = C;! and 61~1 = 02~2, respectively).
The multi-layer heat conduction problem has already been solved using Laplacetransforms with respect to time (Carslaw and Jaeger, 1959; Luikov, 1968). Unfortunately, the results for heat flow cannot be applied directly to solute transport
although a similar solution procedure can be adapted, as will be done in the next
section.
3. Solution Procedure
3.1. GENERAL APPROXIMATION
The previously defined problem was solved by applying the Laplace transformwith respect to time. The transformed concentration is represented as

Ci =

J0

O3 Ciexp(-St) dt,(i = 1, 2),

where s is the transformation variable and the overbar denotes a transformed variable. Application of Equation (8) leads to a pair of ordinary differential equations
with accompanying boundary conditions in the Laplace domain. These can besolved using standard methods; the solution for Cr and 6’2 in the Laplace domain
are given by (cf. Leij et al., 1991):
41 cosh[Xt(z - L)] - 42 sinh[Xr(z - L)]41(q2 + i) cosh(XlL) + (qf + iq2) sinh(XrL)’(9)co
C2 = -exp[atlir+ a2(kz - 112)]S~1 exd-X2(z - L)lX

ql(q2 + i) cosh(XtL) + (4: + ia) Sinh(XlL)’

with auxiliary variables

(10)

70

F. J. LEIJ AND M. TH. VAN GENUCHTEN

We were unable to obtain the solution in the real time domain through application of
the inversion theorem. The solutions are first rewritten as (cf. Carslaw and Jaeger,
1959; Brenner, 1962)c, = coexp(e1 k,)s (QI + $)(41 + q2)(9 - q2)exp[Xi(z - 2L)] + (41 + q2) exp(-Xtz) 9x

1 _ (41-$)(91--42)(12)(41+!)(q,+q2)exp(-2X1~)2

c2 = 2C0 ev[adGS+

a2(k2 - 1c2)](9 + 5)(41 + 42)

x 9 exp[-Ad- A2(5 - L)](13)

1 _ (41-$)(91-d(ql+i)(q,+q*)exp(-2X1L)’
2

According to the binomial theorem:
1-=Fwn,
j?fIl<l.
l - w

(14)

n=O

The denominators of the last terms on the right-hand sides of Equations (12) and(13) can be expanded to obtain infinite series:ClM?exp(atkt)g(9 - q2)n+1(q, - t,”( n=o ( 4 1 t 42)n+*(Q1xft i)n+’

exp{--h[2(n-t l)L-- x]}

5 (Ql - q2)n(ql - $)”n=O(41 + q2)n(g, + $)n+l exp[-‘1(2nl t z)1 ) 7

(15)

Ql co
c2 w -exp[utKtt az(lc2 - Ii,)]x

5 (9 - Q2)yq, - $)”n=O(Ql + 42)"+'(ql + $),+l

x exp[-(2n + l)A,L - A2(2 - L)].

(16)

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

71

Such series solutions converge rapidly for small times (Carslaw and Jaeger, 1959;
p. 309). Only the first term of each series is considered here:

x exp[-(2Kr - kr)Js+aT1 ,
)
_

c2 M

4ata&exp[atKt - Kt

(17)

s + a, + (IQ - IIiz)(u2 - &%I
F

These expressions are amenable to analytical inversion. Upon expansion and application of the shifting and convolution theorems, the approximate concentrations in
the real time domain can be written as

exp[-(2Kt - Ict)fi
fi exp[-(2Kt - kr)Js7

c2 M‘tn,a&~ t
exp[atlit - a:r + az(k2 - K2) - Uz(t - r)]a; - UTJ0

72

F. J. LEIJ AND M. TH. VAN GENUCHTEN

x L-l

{

2

(

exp[-(b27~~2Wls - ai(fi(+ ata&’-al,&’-)

exp( --Ii1 &)(3 - a:)(& + al)
exp(-ICI&(ex$Fij@)

[a&l

)I

(S - a:>(&+ aI> ix LF-l, fi exP[-(~;! - W&l
dr
7
s
ai
(
)I

(20)

where ,Cl_?, and C;’ denote the inverse Laplace transform with s ---f t - T ands --f 7, respectively. The following final result can be obtained from entries in thetable of Laplace transforms by van Genuchten and Alves (1982):

Clx 2alCo exp(alkl - &> tA(k~, t) + -&@I, t)

11-$--(1 + 2alkl + 4a$)G(k1, t) + 1
ai - a:tA(2K1 - h, t) + -&{B1(2K1 - h, t)-[l + 2a1(2K1 - k1) + 4a:t]E1(21<1 - kl, t ) } )+a;f[l + al(2K1 - kl) + 2aft]A(2#1 - kl, t)+${[4ait - 2a1(2K1 - k1) - l]B1(21r’l - kl, t)-[4alft - 1 + 2af(2Kl - Icl + 2alt)2]E1(2Kl - 1;1, t)})

21C* - ICI-2ala2 O’ expK4 - 44
I
CC 2(t - T) - a1 >x A(21il - ICI, t - T) + aTE1(2IL’1 - kl, t - T)x (~a(& T) +c2 x

4ala2Coai - a:
x

t
J0

t[l +

)

2a$][&(O, T) - G(O, T)]) dT]}, (21)exp[alKl - air + a2(k2 - K2) - ai(t - T)l;(B2(k2 - A-2, t - 7) + E2(k2 - K2, t - 7))i

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTJ3 TRANSPORT73

a2( 1 + +)A(K1, T) + yQq(K,, T)x-(5 + 2UlK1 +4a:r)E,(Kl, T)]>-ul(~B,(K,, T)t $(3t 2alJclt 4&)
xEl(Kl, T) - up-A(KI, 7))x(

A@;! - h-2, t - T) t 3B2(k2 - I<27 t - T)--E2(k2 - K2, t - T)]

>I

dr,

(22)

where

B;(k, t) = exp(uft - u;k) erfc (“;y,(24)

E;(k, t) = exp(& $ u;k) erfc (ky$i).

(25)

The integrals appearing in these solutions were evaluated with the help of GaussChebyshev quadrature (Carnahan et al., 1969) whereas the subroutine EXF (vanGenuchten and Alves, 1982) was used to obtain the product of exponential andcomplementary error functions. Note that the two-layer solutions diverge when
the properties of both layers are similar (at = ~2); the simpler homogeneous(one-layer) solution should then be used.
Effluent breakthrough curves are more appropriately described using fluxaveraged rather than volume-averaged concentrations (Kreft and Zuber, 1978;
van Genuchten and Parker, 1984). The well known transformation between fluxaveraged and resident concentrations is most conveniently carried out in the Laplacedomain. The solution for cp given by (10) is transformed according to+-~~_~!!g~ (; + 92>@,2where q2 is defined by (1le). This leads to the following approximate expression

for CT:

c2F M t 2UlCOJa; - a: 0

exp[utli’r - u:r t az(k2 - I(2)74

F. J. LEIJ AND M. TH. VAN GENUCHTEN(-u;<t - T)] .Ly, exp[-(k2s --a;K2)Js7
i)x L,’ ((I; - q&)exp>;;;F))[(exp( -Kr&)+L,’ (c&z; - a&z)(s - a:)(& + al)(+cF_, fi(exp[+2 - 1{2>44s - aiia2(s - a~&> exp(-h&jx Jc,’()I(S--a:)@+a1)

dr

)i

*

(27)

The solution in the real time domain is again obtained with a table of Laplacetransforms:

J

UlCO
t exp[ulKi’l - afr + u2(k2 - 112)
G zui-a: 0-u;<t - T)]X([(B2(k2 - K2, t - T) + E2(k2 - K2, t - 7))

UT - g + a;( 1 + 2u:r)

1

A(lil, T)-[a: + alai(2 + alKl + 24r)]Er(Kr, 7)>+(a~(1 + 2a:r)A(Kr, r) - a42 + u,Kr + 24r)Et(Kr, 7))
x (2A(k2 - K2, t - T) + u2[&(k2 - A-2, t - T)
-E2(k2 - Ir,, t - T)])] dr.

(28)

Unless specified otherwise, all concentrations in this paper are of the resident
ty pe.
3.2. THIN TOP LAYER
Several other approximations may be developed for simplifying conditions (Carslaw
and Jaeger, 1959; Brenner, 1962). As an example, consider the solution for the concentration in the second layer for the case where the first layer is very thin. Thissituation is encountered when the soil surface is ponded (81 = l), when somematerial is placed on top of the soil surface to promote an even water and solute
application, or when a soil crust has developed in the soil surface. The hyperbolic

75

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

functions in Equation (10) can be expanded in power series in the usual manner, toobtainexp[alK1 + a2(h - Ii,) - X2(x - L)]X

7

T/1(92+;) [ r+qq:+... 1 +(qf+iq2) [2a~Ii-*q,+~qf+... 1
where X 1 L = 2at Iilql . Equation (29) is written in powers of at lit, a term similar
to the Peclet number for the first layer (Pt = ~1 LI /Dl), as follows:(30)
Approximate solutions can be obtained for small values of utKt by neglecting
higher order terms. If only the zero-order term in the denominator of (25) is
considered, the following solution can be readily obtained:c;! Lz 2uzCa exp[utKt + a~(& - XT) - ~$1
x tA(k2 - K2, t) + $f?2@2 - Ii,, q
2
1-&[I + 2U&- K2)+ ‘k&]&(k2 - 1<2,

t)} .

(31)

This expression is very similar to the solution for a semi-infinite homogeneous
medium with a third-type inlet condition (Lindstrom et al., 1967). We also consider
the case where, in addition to the zero-order in the denominator of (31), the term
containing at Kl is included in the analysis. The denominator can be factoredaccording to

2UlCO exp[ut1it + u2(k2 - K2) - ~$1
c2 % Ii-1

exp[-(k2 - K2)Js7
xc-’ ( (s - 4)(4X+ a2)(&+ b)1 ’

(32)

76

F. J. LEIJ AND M. TH. VAN GENUCHTEN

where
b=

2 + ;(u: - a;>.

(33)

This leads to the first-order approximation:xE2(IC2 - K2, T) {A(O, t - 7) - bEb(O, t - T)} dr.
>(34)

4. Results and Discussion
4.1. TABULAR RESULTS
The correctness of the previous solutions and the applicability of the series approximations were assessed through comparison with results obtained by numericallyinverting (9) and (10) with an algorithm according to De Hoog et al. (1982). Results
for the numerical inversion of the test function, consisting of exponential and complementary error functions, were identical to the exact (analytical) inversion for atleast five significant digits.
Consider a medium consisting of two layers with contrasting transport properties: values for v were 25 or 40, and those for D were 20 or 40 (any consistent
set of units can be used for these parameters, e.g., cm and days). Table I compares
values for the concentration in both layers according to (21) and (22) with results of
the numerical inversion of (9) and (10) for three different parameter combinations.
Notice that the analytical concentrations are very similar to the values obtainedwith the inversion algorithm for cases 1 and 2. The approximation was found to be
the least accurate when ai and a2 approached the same value, or when concentrations close to the interface were evaluated. As was pointed out earlier, the binomial
approximation is most appropriate for small times when s is relatively large. Notethat large time behavior can often be inferred from the steady-state solution. We
found that inversion of an additional term in the Laplace domain solutions, i.e.,n = 1 in (15) and (16), only slightly improved the comparison with the (exact)
results from numerical inversion.
Case 3 in Table I provides an example where the prediction may become
inaccurate. The parameter values for this case were assigned such that at (= 2.795)
and a2 (= 2.828) approached the same value. The approximate results show thebiggest discrepancies with the concentrations obtained by numerical inversion atAPPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT77

TABLE I. Solute concentration in a two-layer porous medium with v = 25 or 40, D = 20 or

50, and L = 10 as predicted according to (21) and (22), C, or by numerically inverting (9) and(10), as denoted by .L-’
X

C

c-1
tc0.2

C

c-’
tz0.4

C

tz0.6

IF’

c

IT
tz0.a

Case 1: VI = 25, DI = 50, u2 = 40, and Dz = 20 (al = 1.768 and a2 = 4.472)
0
0.884
0.884
0.963
0.963
0.986
0.987
0.995
0.995
2
0.742
0.742
0.915
0.915
0.968
0.969
0.988
0.988
4
0.561
0.561
0.841
0.841
0.940
0.940
0.977
0.977
6
0.375
0.375
0.746
0.746
0.901
0.901
0.961
0.962
8
0.222
0.222
0.644
0.645
0.857
0.858
0.944
0.945
10
0.141
0.142
0.579
0.579
0.828
0.829
0.932
0.933
10
12
14
16
18
20

0.142
0.063
0.021
0.005
0.001
0.000

0.142
0.063
0.021
0.005
0.001
0.000

0.579
0.480
0.372
0.264
0.168
0.094

0.579
0.480
0.372
0.264
0.168
0.094

0.828
0.780
0.722
0.650
0.567
0.473

0.829
0.781
0.722
0.65 1
0.567
0.473

0.932
0.913
0.888
0.857
0.818
0.769

0.933
0.914
0.889
0.858
0.819
0.770

Case 2: V, = 40, DI = 20, ~2 = 25, and D2 = 50 (al = 4.472 and a2 = 1.768)
0
0.999
0.999
1 .0001 .ooo1 .ooo1.000
1 .000
1.000
2
0.988
0.988
1 .ooo1 .000
1 .000
1.000
1.000
1.000
4
0.928
0.928
0.999
0.999
1 .0001 .ooo1 .000
1 .000
0.764
0.764
0.995
6
0.995
1 .000
1 .0001 .ooo1 .ooo
8
0.496
0.496
0.976
0.976
0.998
0.998
0.999
0.999
10
0.152
0.152
0.780
0.780
0.940
0.940
0.979
0.979
10
12
14
16
18
20

0.152
0.049
0.013
0.003
0.000
0.000

0.152
0.049
0.013
0.003
0.000
0.000

0.780
0.600
0.418
0.262
0.148
0.075

0.780
0.600
0.417
0.262
0.148
0.075

0.940
0.870
0.773
0.653
0.522
0.393

0.940
0.870
0.773
0.653
0.522
0.393

0.979
0.952
0.911
0.851
0.774
0.681

0.979
0.952
0.911
0.851
0.774
0.681

Case 3: VI = 25, Di = 20, ~2 = 40, and D2 = 50 (al = 2.795 and a2 = 2.828)
0
0.978
0.978
0.998
0.998
1 .0001 .ooo1 .ooo1 .ooo
2
0.868
0.868
0.984
0.984
0.998
0.998
1 .ooo1 .ooo
4
0.633
0.634
0.942
0.942
0.991
0.99 1
0.998
0.999
6
0.345
0.345
0.848
0.849
0.972
0.972
0.995
0.995
8
0.131
0.131
0.693
0.693
0.930
0.930
0.986
0.986
10
0.033
0.033
0.501
0.496
0.857
0.853
0.968
0.966
10
12
14
16
18
20

0.034
0.011
0.003
0.001
0.000
0.000

0.033
0.011
0.003
0.001
0.000
0.000

0.504
0.372
0.258
0.167
0.100
0.055

0.496
0.370
0.257
0.166
0.098
0.054

0.867
0.789
0.703
0.607
0.506
0.405

0.853
0.784
0.699
0.601
0.498
0.395

0.983
0.95 1
0.920
0.881
0.832
0.771

0.966
0.944
0.913
0.871
0.817
0.75 1TABLE II. Solute concentration in the second layer as predicted for different values of al and a2 with: (i) numerical inversion of (10). (ii) analytical
approximation (22), (iii) small layer approximation (31), and (iv) small layer
approximation (34)
X

(10)

(22)

(31)

(34)

(10)

(22)

(31)

(34)

C a s e 1 w = 5, DI = 10, v2 = 2.5, D2 = l0, L = 0.5 (a,~iPl = 0.125)
t = 1.6
t = 0.8
0.5
1.0

1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0

0.554
0.502
0.450
0.399
0.351
0.305
0.263
0.223
0.188
0.156

0.583
0.530
0.477
0.425
0.374
0.326
0.281
0.240
0.202
0.168

0.658
0.599
0.540
0.482
0.426
0.373
0.323
0.276
0.234
0.196

0.559
0.507
0.455
0.405
0.356
0.310
0.267
0.228
0.192
0.160

0.703
0.667
0.630
0.592
0.554
0.516
0.478
0.441
0.404
0.369

0.716
0.681
0.645
0.608
0.571
0.533
0.495
0.458
0.42 1
0.385

0.816
0.776
0.734
0.692
0.649
0.606
0.563
0.521
0.479
0.438

0.705
0.669
0.632
0.594
0.556
0.518
0.48 1
0.443
0.407
0.372

Case2 VI = 5, DI = 10, v2 = 2.5, D2 = 10, L = 2.5 (u,K, = 0.625)
t= 1.6t = 0.8
2.5
3.0
3.5
4.0
4.5
5.0

0.432
0.378
0.328
0.280
0.237
0.198

0.436
0.382
0.330
0.283
0.239
0.200

1.085
0.987
0.890
0.794
0.702
0.614

0.577
0.518
0.461
0.406
0.355
0.306

0.631
0.591
0.551
0.511
0.472
0.433

0.638
0.598
0.558
0.517
0.477
0.438

1.345
1.279
1.211
1.141
1.070
0.999

0.763
0.721
0.678
0.635
0.592
0.549

Case3 VI = 25, DI = 10, 24 = 20, D2 = 10, L = 0.5 (alIt’ = 0.625)
t =0.lt = 0.2
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0

0.891
0.790
0.657
0.508
0.360
0.233
0.137
0.073
0.035
0.015

0.892
0.791
0.659
0.509
0.362
0.234
0.137
0.073
0.035
0.015

1.762
1.625
1.424
1.172
0.895
0.628
0.402
0.234
0.122
0.058

1.043
0.940
0.801
0.638
0.470
0.318
0.196
0.110
0.055
0.025

0.978
0.955
0.918
0.866
0.796
0.712
0.615
0.511
0.408
0.311

0.977
0.954
0.918
0.866
0.797
0.712
0.615
0.512
0.408
0.312

1.847
1.817
1.768
1.694
1.590
1.455
1.292
1.109
0.916
0.727

1.127
1.103
1.065
1.010
0.937
0.846
0.740
0.624
0.506
0.394

Case4 VI = 25, DI = 10, ~1 = 20, D2 = 10, L = 2.5 (u,K, = 3.125)
t =0.lt = 0.2
2.5
3.0
3.5
4.0
4.5
5.0

0.450
0.298
0.179
0.096
0.047
0.020

0.452
0.300
0.180
0.097
0.047
0.020

21.47
19.80
17.35
14.28
10.90
7.65

4.740
4.179
3.478
2.702
1.944
1.283

0.869
0.792
0.697
0.591
0.480
0.373

0.876
0.798
0.702
0.595
0.483
0.375

22.50
22.14
21.54
20.64
19.37
17.73

5.34
5.19
4.97
4.67
4.28
3.81

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

79

or near the interface when the solute reaches the second layer. Inverting a second
term for the series expression of C2 did not yield significantly different results.
Table II compares values for the concentration in the second layer predicted with(i) numerical inversion of (l0), (ii) the binomial expansion according to (23), (iii)the zero-order approximation for small L, and (iv) the first-order approximation
for small L. The usefulness of the approximations was again assessed by using thenumerical inversion technique as the ‘benchmark’ method. Solute profiles are given
for two values oft assuming a medium with a length typical of soil cores and using
four different parameter sets. Case 1 involves a relatively thin first layer with fairly
small values for at, ~22, and L. The first-order thin-layer approximate solution (34)
is most suitable for this case while approximation (22) performs relatively poor
because of a low Peclet number, pt. Notice that zero-order approximation (31) ismuch less accurate. The same transport parameters as for case 1 were also used forcase 2, except that both layers were now of equal length. Because of the increasein L, the short-layer approximations were found to be less accurate; conversely,
the solutions according to (22) became more reliable. Case 3 of Table II involves
a higher value for at and again a smaller L such that al Kr is the same as forcase 2. The results in Table II show that the relative errors associated-with thefirst-order short-layer approximation are similar for cases 2 and 3 for both shortlayer approximations. Finally, case 4 illustrates that the thin-layer approximations
become progressively worse for larger al Ir’r .

4.2. GRAPHICAL RESULTS
Thus far little attention has been paid to the use of alternative boundary conditions.
The choice of boundary conditions, albeit a somewhat esoteric issue, deservescareful attention when modeling transport in porous media. Figure 1 shows solute
profiles for three different types of interface conditions: (i) a first- or concentration-type condition for layer 2 and an infinite outlet condition for layer 1 (cf. Equation(7)), with a first-type condition at the inlet of layer 1, (ii) a third- or flux-typecondition for layer 2 and an infinite outlet condition for layer 1, with a third-typecondition at the inlet of layer 1, and (iii) a combined first- and third-type condition
as formulated by (5) and (6), also with a third-type condition at the inlet of layer
1. The remaining conditions are defined by (3), (4), and (7). A further discussionof the interface conditions is given by Leij et al. (1991). Figures la and lb present
the three profiles for at = 0.79 and a2 = 3.35 at t = 0.5 and 1, respectively, while
Figures lc and Id show the solute distributions when the layers are switched.
A first-type condition at the interface overpredicts the amount of solute in bothlayers whereas the third-type interface condition predicts a discontinuity in thesolute concentration at the interface. The use of conditions (6) and (7) leads to a
macroscopically continuous concentration at the interface and a mass-conservative
solution (Leij et al., 1991).

80

F. J. LEIJ AND M. TH. VAN GENUCHTEN

0.8

0.8 0
v

0.6

-

0
0

0.6 -

0

0.4

-

0

0.4 0.2 -_ _ _ _ _ _ _ _ _ _ __--__--__

first
lhird
first & lhird_ _ _ _ _ _ _ _ _ _ _ firs*_--__--_.
third
first & third

Fig. 1. Solute profiles in a two-layer medium predicted using (i) a first-type, (ii) a third-type,and (iii) a first- and third-type interface condition with VI = 10, DI = 40, ZQ = 15, and
DZ = 5 at (a) t = 0.5, and(b) t = 1, and with VI = 15, DI = 5, vz = 10, and D2 = 40 at(c) t = 0.5, and (d) t = 1.

Knowledge of the time-dependent behavior of solutes in the subsurface is ofinterest for many practical problems where the concentration is observed or needs
to be predicted at fixed positions. Figure 2 shows breakthrough curves, resulting
from a solute pulse with duration 0.5, just before (z = 9) and after (z = 11)
the interface, and at the outlet of the second layer (z = 20-), for the resident
concentration calculated with (21) and (22). Also included are results for the fluxaveraged concentration at the outlet of the second layer (z = 20+) calculated with
(28). The same parameters are used as for the example shown in Figure 1. Figure 2
may be used to examine how the ordering of the two layers affects a breakthroughcurve. As expected, the curves at 2 = 9 and 11 are different when the layers are
switched. However, notice that the curves in Figures 2c and 2d are also different at
z = 20 after the solute has traversed an equal distance in both layers. In contrast,
for a medium consisting of two ‘semi-infinite’ layers it can be shown explicitly thatthe order of the semi-infinite layers does not affect the outlet concentration (Shamirand Harleman, 1967; Leij et al., 1991). This last finding for the mathematicallysemi-infinite first layer is a result of the assumption that properties of layer 2 will
not influence the concentration of layer 1.

81

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

0.6

1.0

X.11

/

-

I

.

I

1.0

-

b

I

x=20*

0.8

0”

0.6

v

0.4

2

t

’

I

.

I

’

d

0.2
0.0

Fig. 2. Breakthrough curves for a pulse application of solute with duration to = 0.5 to a
two-layer medium with properties VI = 10, DI = 40, ~2 = 15, and D2 = 5 (solid line) or
VI = 15, D1 = 5, ~2 = 10, and D2 = 40 (dashed line): (a) C?(t) at z = 9, (b) C:(t) at
z = 11, (c) C2R(t) at 1: = 20-, and (6) CT(t) at z = 20t.

Two-layer solutions may also be useful for analyzing the effect of dispersion
in the influent or effluent zones of homogeneous porous media typically used incolumn displacement experiments. The ‘one-layer’ analytical solutions normallyused for describing the experimentally determined breakthrough curve will yield
effective transport parameters that are blurred by transport phenomena outside theporous medium. James and Rubin (1972) investigated apparatus-induced dispersion
by independently determining v and D for the porous medium and the outlet
reservoir of a column displacement apparatus using the two-layer solution by
Shamir and Harleman (1967). The parameters reported by James and Rubin for
their experiment A, will be used here to predict flux-averaged concentrations withthe two-layer model according to Equation (28). This solution is compared withthe traditional one-layer ADE solution (case Al by van Genuchten and Alves,
1982) for transport in the soil followed by advection only in the effluent reservoir.
Figure 3 shows the calculated breakthrough curves obtained with the following
transport parameters: VI = 7.55 cm/d and DI = 0.864 cm2/d for the 8.9 cm longsoil and ~2 = 2.56 cm/d and II2 = 0.527 for a bottom compartment of thickness
0.24 cm. The breakthrough curve was calculated at z = 9.14 cm. Notice that the

F. J. LEIJ AND M. TH. VAN GENUCHTEN

82
1.0

0.8
0
0

0.6

0.8 -

4 0.6 -

0‘

0.4 0.2 - l=ld
0.0

0

2

6

4

a

10

X

___________

D=0.527cm2/d

Lko

Fig. 3. Breakthroughcurves, CF(2), at the bottom of adisplacementapparatus (z = 9.14cm)
predicted with and without dispersion in the effluent reservoir assuming VI = 7.55 cm/d,
DI = 0.846 cm*/d, vz = 2.56 cm/d, and, optionally, D2 = 0.527 cm2/d [James and Rubin,1972].

predicted breakthrough curve will be too steep if effluent dispersion is present
but ignored. The results could become even worse if the effluent reservoir, withadditional advective transport, is also neglected.
Similar considerations pertain to dispersion in an influent reservoir. Figure 4
shows an example for such a situation. The resident concentration profiles shownfor t = 0.5 and 1 were calculated with the two-layer solution according to (21)and (22) and the traditional one layer solution for advective-dispersive transport in
the second layer (case A2 by van Genuchten and Alves, 1982) assuming advective
transport in the influent compartment from x = 0-0.24. A somewhat faster andsteeper breakthrough is predicted when influent dispersion is neglected.

83

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

1.0
0.8

_

----

Da=0

4 0.6 -

2

0.4 0.2 0.0

0

1
t

2

Fig. 4. Solute concentration profiles, C”(z), predicted with and without dispersion in theinfluent reservoir, assuming VI = 2.56 cm/d, D1 = 0.527 cm’/d (optional), “2 = 7.55 cm/d,and D2 = 0.846 cm’/d at (a) t = 0.5 d and (b) t = 1 d.

5. Summary and Conclusions
An (approximate) analytical solution was derived for the ADE describing solutetransport during steady one-dimensional flow in a porous medium made up of
two homogeneous layers whose interface was perpendicular to the flow direction.
The solution was obtained with Laplace transforms and the binomial theorem.
The conditions at the interface presume that both the solute concentration andthe solute flux are continuous. Similar analytical solutions have previously been
useful for the early time description of heat conduction or the transport of solutes
in homogeneous finite media. Calculated resident concentrations predicted with
the newly developed analytical solution, compared favorably with those obtained
through numerical inversion of the Laplace transform, except for small values of
the Peclet number or in cases where v~D~/vZD~ is close to unity. Expressions forthe flux-averaged concentration for the second layer were also developed.
Alternative approximate solutions for the concentration in the second layer were
derived by using zero- and first-order expansions of the Peclet number for the firstlayer. The first-order approximation may complement the solution for cases witha thin first layer and a low Pt ; the zero-order approximation was shown to beinaccurate.
Several examples of transport in two-layer media were presented. The solutions
predict continuous concentrations across the interface and conserve mass within the transport system. In contrast, two previously used approaches that use amathematically semi-infinite first layer and invoke either a first- or third-type inlet
condition for both layers lead to mass balance errors or a discontinuous concentration at the interface, respectively. The interface conditions for which the current

84

F. J. LEIJ AND M.TH.VAN GENUCHTEN

tration at the interface, respectively. The interface conditions for which the current
solutions were derived imply that the ordering of the layers will affect the breakthrough curve at the outlet of the medium. The utility of the two-layer solution foranalyzing transport in the influent or effluent regions of laboratory columns was
illustrated using experimental data pertaining to apparatus-induced dispersion.
Acknowledgments
We would like to thank John Knight for providing a copy of the numerical inversion
program according to the De Hoog method, as well as Ed Sudicky for supplying
us a version modified by Chris Neville.
References
Al-Niami, A. N. S. and Rushton, K. R., 1979, Dispersion in stratified porous media: analyticalsolutions, Water Resour. Res. 15, 1044-1048.
Barry, D. A. and Parker, J. C., 1987, Approximations for solute transport through porous media withflow transverse to layering, Transport in Porous Media 2,65-82.
Brenner, H., 1962, The diffusion model of longitudinal mixing in beds of finite length. Numericalvalues, Chem. Eng. Sci. 17,229-243.
Camahan, B., Luther, H. A., and Wilkes, J. O., 1969, Applied Numerical Models, John Wiley, NewYork.
Carslaw, H. S. and Jaeger, J. C., 1959, Conduction of Heat in Solids, Clarendon Press, Oxford.
Dagan, G. and Bresler, E., 1985, Comment on ‘Flux-averaged and volume-averaged concentrations
in continuum approaches to solute transport’, by J. C. Parker and M. Th. van Genuchten, Water
Resour. Res. 21, 1299-1300.
Danckwerts, P. V., 1953, Continuous flow systems, Chem. Eng. Sci. 2, 1-13.
De Hoog, F. R., Knight, J. H., and Stokes, A. N., 1982, An improved method for numerical inversionof Laplace transforms, SIAM J. Sci. Stat. Comput. 3,357-366.
Jacobsen, 0. H., Leij, F. J., and van Genuchten, M. Th., 1992, Lysimeter study of anion transport
during steady flow through layered coarse-textured soil profiles, Soil Sci 154, 196-205.
James, R. V. and Rubin, J., 1972, Accounting for apparatus-induceddispersion in analysis of miscible
displacement experiments, Water Resour. Res. 8,717-721.
Jury, W. A. and Utermann, J., 1992, Solute transport through layered soil profiles: Zero and perfect
travel time correlation models, Transport in Porous Media 8,277-297.
Kreft, A., 1981a, On the residence time distribution in systems with axial dispersed flow, Bull. AC.
Pol. Sci. Tech. 29,509-519.
Kreft, A., 1981b. On the boundary conditions of flow through porous media and conversion ofchemical flow reactors, Bull. AC. Pol. Sci. Tech. 29,521-529.
Kreft, A. and Zuber, A., 1978, On the physical meaning of the dispersion equation and its solutions
for different initial and boundary conditions, Chem. Eng. Sci. 33, 1471-1480.
Leij, F. J., Dane, J. H., and van Genuchten, M. Th., 1991, Mathematical analysis of one-dimensional
solute transport in a layered soil profile, Soil Sci. Soc. Am. J. 55,944-953.
Lindstrom, F. T., Haque, R., Freed, V. H., and Boersma, L., 1967, Theory on the movement of someherbicides in soils: Linear diffusion and convection of chemicals in soils, J. Envir. Sci. Technol.

1,561-565.
Luikov, A, V., 1968. Analytical heat diffusion, Academic Press, New York.
Mikhailov, M. D. andijzigik, M. N., 1984, Unified Analysis and Solutions of Heat and Mass Diffusion,
John Wiley, New York.
Nakayama, S., Takagi, I., Nakai, K., and Higashi. K., 1984, Migration of radionuclide throughtwo-layered geologic media, J. Nucl. Sci. Technol. 21, 139-147.

APPROXIMATE ANALYTICAL SOLUTIONS FOR SOLUTE TRANSPORT

85

Novakowski, K. S., 1992a, An evaluation of boundary conditions for one-dimensional solute transport. 1. Mathematical development, Water Resour. Res. 28,2399-2410.
Novakowski, K. S., 1992b. An evaluation of boundary conditions for one-dimensional solute transport. 2. Column experiments, Water Resour. Res. 28,241 l-2423.
Gzigik, M. N., 1980, Heat Conduction, John Wiley, New York.
Parker, J. C. and van Genuchten, M. Th., 1985, Reply to G. Dagan and E. Bresler, Water Resour: Res.
21, 1301-1302.
Selim, H. M., Davidson, J. M., and Rao, P S. C., 1977, Transport of reactive solutes throughmultilayered soils, Soil Sci. Soc. Am. J. 41,3-10.
Shamir, U. Y. and Harleman, D. R. F., 1967, Dispersion in layered porous media, J. Hydraul. Div.
Proc. ASCE 93(HY5), 237-260.
Taylor, G. I., 1953, Dispersion of soluble matter in solvent flowing through a tube, Proc. Roy. Soc.
London. A 219,186-203.
van der Laan, E. T., 1958, Notes on the diffusion-type model for longitudinal mixing in flow, Chem.
Eng. Sci. 7, 187-191.
van Genuchten, M. Th. and Alves, W. J., 1982, Analytical solutions of the one-dimensional
convective-dispersive solute transport equation, U.S.D.A. Technical Bulletin No. 1661.
van Genuchten, M. Th. and Parker, J. C., 1984, Boundary conditions for displacement experiments
through short laboratory soil columns, Soil Sci. Soc. Am. J. 48,703-708.
Wehner, J. F. and Wilhelm, R. H., 1956, Boundary conditions of flow reactor, Chem. Eng. Sci. 6,
89-93.