Mathematical Analysis of One-Dimensional Solute Transport in a Layered Soil Profile
Feike J. Leij, J. H. Dane,* and M. Th. van Genuchten
ABSTRACTPast research on transport in layered soils has focused mainly on the effects of layering on solute breakthrough curves. Experimental results have beenreported by Shamir and Harleman (1966), Selim et al.
(1977), and Moranville et al. (1977). Theoretical analyses of solute breakthrough curves, using transfer functions and time moments, were carried out by Kreft(198 la) and Barry and Parker (1987). In many cases,
however, it is also important to know the actual solute
distribution in a layered profile. Analytical solutions
for solute transport have been reported by Shamir andHarleman (1967), and Al-Niami and Rushton (1979).
Both studies used a first-type (constant concentration)
inlet condition. Shamir and Harleman ( 1 967) assumed
that the profile consisted of semi-infinite layers, while
Al-Niami and Rushton (1979) used finite layers. Since
only the concentration was continuous at the inlet andinterfaces (Dirichlet conditions), neither solution will
satisfy the principle of mass conservation at interfaces
and, as a result, the applied solute mass is not equal
to the mass accumulated in the soil profile. Al-Niami
and Rushton (1979) also imposed the physically unrealistic assumption that the concentration gradient iszero at the interfaces of the layers.
The objective of this study was to simulate onedimensional solute transport in a layered soil using
both numerical and analytical methods. As with mostprevious studies, the layered profile was approximatedby a series of homogeneous layers, each having its
unique physical and chemical characteristics. Analytical solutions were obtained for a two-layer medium
using three different combinations of inlet and interface conditions. Since the boundary and interface conditions play an important role in the solution of the
ADE for layered media, the applicable boundary conditions are briefly reviewed. The effect of layering on
calculated solute concentrations in a layered soil profile is illustrated qualitatively, using both analyticaland numerical solutions.
THEORY
Analytical Solution for a Two-Layer Medium
Figure 1 schematically shows a layered soil, consisting of
two homogeneous layers, during steady water flow perpendicular to the interface of the two layers. The first layer islocated between x = 0 and x = L, and the second layer
between x = L and infinity. Transport of a linearly exchanging solute in both layers is described with the ADE:
dCk££* t> 0 [1]2

Solute transport studies involving layered media are important forinvestigating how restricting layers affect rates of solute migrationin the soil profile and, more generally, for examining the influenceof soil heterogeneity on solute transport. Analytical solutions of the
one-dimensional advection-dispersion equation (ADE) were obtainedwith the help of Laplace transforms for transport in a two-layeredsoil profile. Assuming that the layers are, in effect, semi-infinite,solutions were obtained for first-type (constant concentration) andthird-type (constant flux) conditions at both the inlet boundary andthe interface of the two layers. Concentration profiles were also obtained for a finite first layer via numerical inversion of the Laplacetransform solution, using a third-type condition at the inlet, and,simultaneously, a first- and third-type condition at the interface.
Volume-averaged or resident-type concentrations were used in allcases. First-type conditions did not meet our criterion of mass conservation, whereas third-type conditions caused discontinuities in theconcentration at the interfaces of layers with differing transport parameters. The concentration at the interface was found to be continuous, and no mass-balance error occurred, when first- and thirdtype conditions were imposed simultaneously at the interface. Several example calculations show the effect of soil layering on solutetransport in a one-dimensional soil profile.

NTEREST IN THE TRANSPORT of chemicals in soils is
Ichemicals
motivated by the potential of agricultural and other(fertilizers, pesticides, industrial wastes) tomove from the soil surface through the unsaturated
zone toward the groundwater table. The leaching ofsolutes in the unsaturated zone may be affected greatly
by the presence of soil layers. Soil stratification is anatural phenomenon that is common to many soils.
Also, artificial barriers (e.g., clay liners) are often used
to slow down or prevent the movement of certain
chemicals. Finally, the macroscopic impact of soil layering is important in more general studies of the effectof soil heterogeneity and anisotropy on solute transport. For example, a heterogeneous soil profile is often
approximated by a series of homogeneous soil layers.
Solute transport in soils is usually described deterministically with the ADE, although alternative
stochastic approaches also exist (Jury, 1982). Exact
and approximate analytical solutions for transport inlayered soils are now available for a limited numberof situations. In most cases, however, numerical methods need to be employed. Unfortunately, comparedwith analytical solutions, numerical solutions frequently do not provide as much insight into the effectsof physical and chemical processes on solute transport,
and often are also susceptible to computational errors.
Hence, analytical solutions remain useful for a variety
of applications, including verification of numerical solutions.

dx

dx

where R is the retardation factor, C is the resident concentration3 expressed in mass of solute per volume of solution[M L- ], t is time [T], D is the dispersion coefficient[L2 T-'],
1v is the mean pore-water velocity [L T" ], x is distance inthe direction of flow [L], and the subscript k refers to thesoil layer (k = 1,2). We emphasize that, in this study, C
represents volume-averaged or resident concentrations (vanGenuchten and Parker, 1984) rather than flux-averaged concentrations as employed by Barry and Parker (1987). Volume- and flux-averaged concentrations are related according
to (e.g., Kreft and Zuber, 1978)

F.J. Leij and M.Th. van Genuchten, U.S. Salinity Lab., 4500 Glenwood Drive, Riverside, CA 92501; and J.H. Dane, Dep. of Agronomy and Soils, Auburn University, Auburn, AL 36849. Alabama
Agric. Exp. Stn. Journal no. 3-892244P. Received 23 Apr. 1990.
'Corresponding author.
Published in Soil Sci. Soc. Am. J. 55:944-953 (1991).

944

LEU ET AL: ONE-DIMENSIONAL SOLUTE TRANSPORTvCF = vC - D —dx[2]where the superscript F denotes a flux-averaged value.
The dispersion coefficient D in Eq. [1] reflects two mechanisms responsible for solute spreading, molecular diffusion
and mechanical dispersion. A common expression for thedispersion coefficient is

D = \D0 + «|v|

[3]

where D0 is theionic or molecular diffusion coefficient infree water [L2T-1], X is a tortuosity factor (0 < X < 1), anda is the dispersivity [L]. Equation [3] assumes that, macroscopically, the active process of diffusion can be combined with the more passive process of mechanicaldispersion. It is unlikely that this assumption correctly describes the microscopic diffusion-dispersion process in thevicinity of interfaces. Hence, the ADE as given by Eq. [1]can, at best, approximate solute transport in the neighborhood of interfaces.
Two different approaches were used to solve Eq. [1] analytically. In one approach, we assumed that each homogeneous layer is part of an effectively semi-infinite system.
This approach assumes that concentrations in the first layer
are not affected by the presence of a second layer and, hence,that analytical solutions for homogeneous media can be applied immediately to the first layer. Barry and Parker (1987)
derived conditions that indicate that this approach is justified in most cases. Transport in the second layer was solved
for a first- or third-type condition at the interface of the twosoil layers using the solute concentration of the first layer as
the inlet condition for the second layer. In concept, this approach is relatively straightforward and can be extended toany desired number of layers (although the solution process
becomes more tedious when the number of layers increases).
In the second approach, we assumed that the first layer
was finite. The interface conditions were then formulated sothat solute fluxes and solute concentrations became continuous. This approach caused the concentration of the firstlayer to become dependent on the transport properties ofthe second layer. Both approaches are discussed separately
below.

8, v,

Fig. 1. Schematic of transport in a two-layer medium (6 = volumetric water content, v = pore-water velocity, Cin = concentrationof the incoming solution, D — dispersion coefficient, R = retar-

dation factor, and Cout = concentration of the effluent).

is the more appropriate one for solute displacement experiments (van Genuchten and Parker, 1984). Because solutions
obtained with a first-type condition are useful for situations
where the concentration is measured as a flux-averaged value(e.g., for column effluent curves), we also considered the firsttype boundary condition. We noted that adoption of the thirdtype condition leads to a macroscopic discontinuity of the(volume-averaged) resident concentration at the inlet, eventhough, microscopically, the concentration distribution may
remain continuous for all times t > 0 (see discussions by
Dagan and Bresler, 1985; Parker and van Genuchten, 1985).
This discontinuity results from the simplifying assumptions
that (i) the incoming solution is well mixed and has a constantconcentration, and (ii) the dispersion coefficient at or immediately below the soil surface is the same as that of the soil
profile farther away from the inlet boundary.
If we assume that the concentration of the first layer isunaffected by the second layer, then one can readily invokeboundary conditions for homogeneous media at the outletof that first layer. A commonly used condition assumes thatthe concentration is macroscopically continuous across theoutlet boundary, at x = L:dx

0<x <Ldx[4]where g, is assumed to be a constant. For the inlet boundary
at x = 0, we imposed step functions of duration /„ using
either a first-type condition of the form

= c0

0 < t < t0
t> *„

[5]

or a third-type condition^jv.C;

0<*<*0t> tn[6]where x[Q, sometimes denoted also as x—.0% implies that x= 0 is approached from inside the soil profile. The physicalbasis of Eq. [5] and [6] has been discussed at length in theliterature (e.g., Danckwerts, 1953; Wehner and Wilhelm,
1956; Pearson, 1958; Kreft and Zuber, 1978; van Genuchten
and Parker, 1984). Adoption of a first-type condition always
leads to continuity in the resident concentration, whereas athird-type condition always leads to continuity in the solute
flux. Because of these properties, the first-type condition is
also known as a concentration-type condition, and the thirdtype condition as a flux-type condition. The latter condition

= 0

t> 0

[7]

However, it is difficult to imagine how dispersive transport
between two layers will vanish during transient advective
transport across the outlet boundary. Alternatively, the presence of an effectively semi-infinite medium could be invokedby imposing the condition

Semi-infinite First Layer
The initial condition imposed on Eq. [1] is

C, (x,0) = gl

945

= 0

t> 0

[8]

Because of considerations outlined in detail by van Genuchten and Parker (1984), and because of possible mathematical simplifications later, we decided to use Eq. [8] forour analytical solutions. We emphasize that, at least for homogeneous soils, differences between solutions for finite andsemi-infinite media are usually rather small (van Genuchten
and Alves, 1982).
The analytical solution of Eq. [1] subject to Eq. [4] and[8] for a first- or third-type condition at the inlet is:(x,t)

Si + (Q - ft) A(x,t)+ (C0 - ft) A(x,t)

0<t<t0t> t0 [9]where, for a first-type condition (Lapidus and Amundson,1952),

ertc
CTfc+ 2exp A
and, for a third-type condition (Lindstrom et al., 1967),

946

SOIL SCI. SOC AM. J., VOL. 55, JULY-AUGUST 1991

where T is an integration variable. Similarly, the solution of
Eq. [1], [12], [14], and [15] for third-type conditions is given
by Eq. [16] with

exp

v.Q
4DlR1t\ f erfc [ *t - v.r |J [[ (4R.D, Tr)l/2

A(x,t) =_ ± 1 1 + M + v '*
2I

vx

X exp

V

A

DtR

+ exp 1 -

erfc

X

1

exp -

X

Analytical solutions for the second layer were obtainedwith the following boundary and initial conditions:
[12]
x> L
= g2
[13]
(first type)or(4R2D2)1'2[R2(x - L) - v2(f4/?2Z)2 (t - T)

exp. ^ - L >(x - L) - v2t]2 }4R2D2t}l/2(third type)

Ax

[14]

dC2
[15]
=0
dx
where 6 is the volumetric water content and g2 is assumed
to be constant. Interface Condition [13] is used in conjunction with Eq. [9] and [10] while Condition [14] is appliedwith Eq. [9] and [11]. Notice that Eq. [13] implies that theconcentration is continuous across the interface, while Eq.
[14] implies that the solute flux is continuous.

The analytical solution of Eq. [1] with k = 2, subject to
Conditions [12], [13], and [15] for the first-type inlet andinterface conditions, was obtained by Laplace transformation (for details, see Leij and Dane, 1989). The solution may
be written in the form

(QC2 (x,0g2(C0 - g,M(*,0 - C<A(x, t -0<t<t0t > t0 [16](* ~ £) ~D2R2 (x - L) + v2t

exp

[18b]

The integrals appearing in these solutions were evaluatedwith the help of the Gauss-Chebyshev formula (see Appendix A). Equations [17] and [18] show that transport properties of both layers influence the solute concentrations of
the second layer.
The above solutions hold for a two-layer soil profile. Laplace domain solutions are easily obtained for profiles consisting of more than two layers. Consider step displacement(ta _> oo) in a medium containing n homogeneous layers withuniform initial concentration, g, in the entire profile. For
first-type conditions at the interfaces of two consecutive layers, the Laplace solution for layer k is
k - 1

where
A(x,t) -

R2D2

I V,
x-L

erfc

R{L exp [\k (x - Lk .

<n)

[19]

exp [\k (x - Lk . ,)] +

[20]

J

and, for third-type conditions,

+ exp

exp

(x,s) =

C0~gs (v* k - 1

j- i

where s is the Laplace transform variable, / is the length oflayer j, Lk_i is the distance from the inlet to the interface
between layers k — 1 and k, and
V- \2

2D,

2Djj

1/2

[21]

LEU ET AL: ONE-DIMENSIONAL SOLUTE TRANSPORT

Solutions in the real time domain are complicated, but
can be formally obtained by making use of the convolution
theorem. The important points illustrated by Eq. [19] and[20] are that (i) the order of the soil layers from layer 1 tok — 1 is irrelevant when calculating the solute concentrationof layer k, and (ii) the concentration of any layer k is notaffected by any layer beyond k. Similar conclusions wereobtained by Shamir and Harleman (1967), Selim et al.,(1977), and Barry and Parker (1987).
Finite First Layer
If we deal with resident concentrations, mass conservation
will require that the solute flux be continuous across theinterface. In addition, one could reason that the concentration across the interface must be continuous as well, at least
microscopically, because of the impact of ionic or molecular
diffusion. The transport problem where both the flux andthe concentration are continuous resembles the problem ofheat flow in composite media with no contact resistance at
the interface (Carslaw and Jaeger, 1959), and was applied
previously to steady-state solute transport by Wehner andWilhelm (1956). Rather than assuming that the first layer ispart of an effectively semi-infinite system, the two interface
conditions now serve as an outlet condition for the first layer
and an inlet condition for the second layer (Kreft, 1981b).
Consequently, Eq. [1] was solved subject to the following
conditions

C, (x,0) = ft

0 < x <L

[22]

C2 (x,0) = g2

x>L

[23]

x i 0

v,C0
0

[24]

t>t
t> 0

(

[25]

\

6lvlCl - 0,-D, —dx ] x I L- 622D2
dC2

dx

xlL

= 0

t> 0

t>0

[26]

[27]

Equations [25] and [26] stipulate continuity in concentrationand flux across the interface. Using Eq. [25], the third-type
interface condition can also be written as a second-type(Neumann) condition:^ 2= n82D2 —t > 0 (second type)C, (x,s) =947

sC0 - g, - C0- ft) exp (-a,.

- + q, (0, + 02) exp (w-L)- <72) exp

X

exp

[29]

C2(x,s) = 52gi[Co - g, - C0 exp (-<„£)] exp (v,L/2g,)- «z) exP (~

exp

[30]

where

2D

dd"2
^J

[31]

<,,. - ^JVJHVJ[32]Note that the concentration profile in the first layer now
depends on the transport properties of the second layer. Attempts to analytically invert Eq. [29] and [30] were not successful. Instead, we obtained concentration profiles using
numerical inversion of the Laplace transform according to
the method of Crump (1976). The solution in the Laplacedomain can also be readily used for time moment analysis(e.g., Kreft, 1981a). Finally, we note that the two-layer solution given by Eq. [29] and [30] provides a method for
examining concentration profiles at the inlet and outlet
boundaries of a homogeneous medium by adopting appropriately small values for Z>, and D2, respectively (see also
Wehner and Wilhelm, 1956; Pearson, 1958).
Numerical Solution for a Layered Medium
Analytical solutions are often not available, or cannot bederived (e.g., for nonlinear sorption or time-dependent
boundary conditions), in which case, numerical methods
must be employed. In this section, we briefly outline themathematical problem and its solution. For a soil column
consisting of n layers, Eq. [1] was solved for the following
conditions:

0<x<LnCt (x,0) = &( - A f + v,C,)

x i

[33]
t>0

[34]

[28]

The analytical solution of the problem, subject to Eq. [22]to [27], is considerably more difficult than the previous solutions for either the first-type or third-type inlet and interface conditions. The solution in the Laplace domain isdx

x 1 Lk(second type)dx

f c = 1,2.. n

[35]

948

SOIL SCI. SOC. AM. J., VOL. 55, JULY-AUGUST 1991

or
(SvC)k-{,D-[36](third type)

[37]

t>0

=0

dx

Note that only one condition is needed at each interface,i.e., either Eq. [35] or [36]. Equation [1], subject to Eq. [33]through [37], was solved numerically using an implicit
Crank-Nicolsen-type finite-difference method. Details of thenumerical solution are given by Leij and Dane (1989).
RESULTS AND DISCUSSION
Several examples are presented to illustrate the effects of layering and the use of different boundary conditions on computed solute distributions. First, we
assessed the magnitude of mass-balance errors associated with the three analytical solutions that were derived for the two-layer medium. A second example
shows concentration profiles obtained analytically fora two-layer medium having layers with different advective and dispersive properties. A third example
considers two different scenarios for reducing contaminant transport in a soil profile. Concentration profiles
for this last example were obtained numerically. Table
1 lists numerical values of the transport parameters
used for these examples.
Analytical Solutions for a Two-Layer Medium
Appendix B describes how the mass-balance error
is denned and how errors are obtained for the three
Table 1. List of transport parameters for the simulations (v = porewater velocity, D = dispersion coefficient, 6 = volumetric watercontent, and R = retardation factor).tFig.

Layer (fc)

V

D

#

analytical solutions. These errors were obtained bycomparing the amount of solute accumulated in thesoil profile with the cumulative solute flux entering
the soil (van Genuchten and Parker, 1984). The solutions are distinguished according to the invoked interface condition, i.e., a first-type solution, given by
Eq. [9], [10], [16], and [17]; a third-type solution given
by Eq. [9], [11], [16], and [18]; and a combined firstand third-type solution given by Eq. [29] and [30]. Fora third-type condition, the relative error was exactly
zero, whereas, for a combined first- and third-type condition, the error was found to be almost zero (<0.1%)
due to approximations involving the numerical inversion of the Laplace transform. The error associatedwith the first-type boundary and interface conditionsin terms of dimensionless parameters is (see Appendix
B):p2
1
2 expi - £W,)"^exp(7>)erfc(^i I - (1 - P + f,)

erfc

-fi

^-Jlexp^erfcfli^l]erfc

P-Si-y

[38]

R

cm/d2
3a
1
1.0
40
0.4
1.0
0.4
2
5
1.0
3b
1
0.4
5
1.0
0.4
2
40
1
1.0
4a
0.2
50
20
20
0.5
2
1.0
20
1
0.5
1.0
4b
20
20
0.2
2
1.0
20
50
1
4.25
5a
0.4
100
52
0.4
2
4.25
52
100
0.4
3
100
52
4.25
4.25
0.4
1
7
Sb
10
0.5
14.0
2
18
8
0.4
7
4.25
3
10
1
0.4
4.25
5c
10
7
0.5
14.0
2
18
8
0.4
4.25
7
10
3
6a
1
0.2
7.5
7
10
0.2
7.5
2
7
10
7.5
0.2
3
10
7
1
0.2
4
7.5
1
6b
1.0
0.4
2
0.5
2.5
0.2
7.5
4
3
1
6c
1
7.5
4
0.2
1
0.4
1.0
2
2.5
0.5
1
0.2
7.5
3
4
t £*/C0 — 0.05, te, is the initial concentration of layer k and C, is the concentration of the incoming solute pulse).

cm/d
10
10
10
10

Fig. 2. Relative mass-balance error as a function of the dimensionlesstimes f i and f2 for the analytical solution of transport in a twolayer medium using a first-type condition.

LEU ET AL: ONE-DIMENSIONAL SOLUTE TRANSPORT

layer, whereas the opposite is true for Fig. 3b. This
causes the concentration profile in Fig. 3a to become
steeper in the second than in the first layer, whereas
in Fig. 3b the reverse can be observed. Previous studies(e.g., van Genuchten and Parker, 1984) and Eq. [38]showed that a first-type condition at the inlet boundary
introduces a mass-balance error that depends on thevalues for D and v. This error is proportional to D and
can be evaluated qualitatively by comparing the firstwith the third-type solution. Figures 3a and b showthat the largest errors occur in the first and second
layer, respectively.
A discontinuity in concentration develops when athird-type interface condition is used. This situation
is similar to the inlet boundary where the flux-averaged concentration is continuous, and the resident
concentration discontinuous. From Eq. [2] and [14],it follows that flux-averaged concentrations are continuous at the interface. Finally, the solution for a
combined first- and third-type condition at the interface exhibits both mass conservation and continuityof the (resident) solute concentration. The differences
between the solute profiles in the first layer for thirdand combined first- and third-type conditions show
the influence of the second layer on transport in thefirst layer (Fig. 3a). The differences in concentration
are less pronounced for lower values of D in the firstlayer (Fig. 3b).
The effect of different values for v on concentration
profiles during steady flow are shown in Fig. 4. In Fig.

where the dimensionless time ffc = vfe/(RkDk) and thePeclet number P = v,.L/.D,. For the one-layer case
(L —> °°), Et reduces to Eq. [10] of van Genuchten andParker (1984). We note that, if the solute influx were
determined just inside the porous medium as the sumof the advective and dispersive fluxes, the error resulting from the inlet boundary would be equal to zero
because of the mass-conserving nature of the ADE(first-type conditions at the interface of soil layers
would still lead to errors). Although the use of suchan error criterion seems attractive for mathematicalreasons, the solute influx just inside the medium may
provide an incorrect estimate of the actual (physical)
solute flux just outside the medium, since no flux continuity is mandated for a first-type inlet condition. Figure 2 shows Er (%) vs. the dimensionless parameter f t
for both layers. The error is largest for small values off! and f2> with f j having a slightly larger effect on theerror than f2. If f i = £2, we obtain Fig. 1 of van Genuchten and Parker (1984).
Figures 3 and 4 show calculated solute distributions
for transport of a nonreactive solute (R = 1) in a twolayer medium at various times for different values of
D and v. Relatively large differences in D and v wereused to illustrate the effects of soil layering. The typesof inlet and interface conditions are specified in the
figures (note that a third-type soil surface condition
was always used in conjunction with a combined firstand third-type interface condition). The value of D forthe first layer of Fig. 3a is larger than that of the second
Time=

1.0

0.4

949

Time = O.8 dd. _ . _ . _ . Ist-type————— 3rd-type
I - - - - — 1st- and 3rd-type

0.8

1st-type
3rd-type
1st- and 3rd-type

o 0.6

o

0.4
layer Iv«IO cm/d

0.2

layer II
v=IO cm/d
D=5 cmz/dlayer I
v=IO cm/d0=40 cm2/d

D-4Ocm2/d

10

12

10

140

x (cm)14

x (cm)

Fig. 3a. Relative concentration (C/C0) as a function of distance (x) for the analytical solution of the advection-dispersion equation for a twolayer medium with v, = v2 (v = pore-water velocity) and />, > D2 (D = dispersion coefficient) at two different times using three different
interface conditions.
Time* 0.4 d
I.U

x\

0.8

Time= 0.8 d~ J ---^Ji .-.-.'-. ist-iype^ijj.
———— 3rd-type\sN.----— 1st- and 3rd -type. _ . _ . _ . Ist-type————— 3rd -type- —— — 1st- and 3rd -type"^•'5X '
\\
\\

0.6
0°

\l
>l

•

'\
'v

X

\\

0.4
0.2

layer 1v=K> cm/dz

0=5 cm /dlayer IIV^A

v= 10 cm/d

^ __D* 40 cm2/d .

^•*-& '**"'•—.

O^^=^"S".~x (cm)'.

rlayer 1v»IO cm/d
0*5 cmVd

|

x

^^,

X

"-

^^^^
^'~
^^i^..
|
layer II^^~i
VBlO cm/d| D'40 cmVdx (cm)

Fig. 3b. Relative concentration (C/C0) as a function of distance (x) for the analytical solution of the advection-dispersion equation for a twolayer medium with v, = v2 (v = pore-water velocity) and Dt < D2 (D = dispersion coefficient) at two different times using three different
interface conditions.

950

SOIL SCI. SOC. AM. J., VOL. 55, JULY-AUGUST 1991

Time= O.I d
I.U

^-*>^
%\

'

Time = 0.2 d

j .-.-.-. »-rype————— 3rd-type'\ \ ————— 3rd- type

0.8
l\\

\

o 0.6
o

\ \

V
Vf .
\ \

0.4

layer 1v .50 cm/d
0.2

D-2Ocmz/d_C)layer 1v=50cm/d
D=20cmVd

layer II
"20 cm/d
D'20cm2/d

\

\
\ \
\^.

layer II .
v=20cm/d
D=20cm2/d

^i'x

i

.

2

4

6

.

1

8

K>

\Z

\ 4C)

2

4

i

.

i

6

8

10

"" ~-V12l<x (cm)x (cm)

Fig. 4a. Relative concentration (C/C0) as a function of distance (x) for the analytical solution of the advection-dispersion equation for a twolayer medium with Z), = D2 (D = dispersion coefficient) and v, > v2 (v = pore-water velocity) at two different times using three different
interface conditions.
Time* 0.2 d
1.0

7'^.

0.8

'
\

'

Time= 0.4 d——v^<^'.-.-.'-. IM-type.—— .'— Ist-type^^<- i^l - -----SS

1st - and 3rd - type.

N w

^ '^*^ ^ *s^
o 0.6
0
~v

0

0.4
0.2

°C)

^*5^^

\\ '
\\ 1
N>Xk
"•»•
"
soemM
N ^«sj^n.on._lM<.~-. ^D'ZO cmVdi™- i
V 20 cm/d
D.20crrfVd

T*
2

x%t

4

6

"•

8

10

•
12

x (cm)

1 * C)

layer 1v= 20 cm/d
D=20cm^d
2

4layer II
v=50cm/d
D«20cnr/d
6

8

10

12

1'

x (cm)

Fig. 4b. Relative concentration (C/C0) as a function of distance (x) for the analytical solution of the advection-dispersion equation for a twolayer medium with Dl = D2(D = dispersion coefficient) and v, < v2 (v = pore-water velocity) at two different times using three different
interface conditions.

4a, the pore-water velocity is largest in the first layer,
while in Fig. 4b the second layer has the larger velocity.
The solute front is steeper at the interface in the lowvelocity layer. The mass-balance error associated withthe first-type condition can again be evaluated by comparison with the third-type condition. This error isinversely proportional to v. For a third-type condition,
the concentration at the interface is lowest in the lowvelocity layer. The values of the advective and dispersive fluxes are, respectively, lower and higher inthe low-velocity layer than in the high-velocity layer
because of the differences in 6. We note again that, fora third-type condition, the flux-averaged concentration is continuous at the interface. The differences between the third- and the combined first- and thirdtype conditions are minor except in the vicinity of theinterface. The analytical solution for a third-type condition can, therefore, be used in many instances, rather
than the more complicated solution for a combined
first- and third-type condition, which requires numerical inversion of the Laplace transform.
Numerical Solutions for a Three-Layer Medium
The next two examples are meant to illustrate theeffect of a "barrier" on contaminant transport. Concentration profiles were obtained numerically using athird-type condition at the inlet and a second- or thirdtype condition at the interface.
Consider a homogeneous sandy soil with a low value
for R, which is saturated or close to fluid saturation.
Figure 5a shows calculated solute distributions for thisprofile, using the parameter values listed in Table 1.
Transport of the solute can be slowed by incorporating
a fine-textured layer that reduces the solute flux bydecreasing the water flow rate and increasing adsorption. Figure 5 shows the effect on the solute-concentration profiles of embedding a 2-cm-thick clay layer
in the otherwise homogeneous profile. Results are forthird-type (Fig. 5b) and second-type (Fig. 5c) conditions at the interface of the soil layers. The presenceof the clay layer was assumed to result in a tenfold
reduction of v (Table 1). The value of D was
calculated
according to Eq. [3], using \D0 = 2 cm2/d, while a
was assumed equal to 0.5 cm for the sand and 2 cm
for the clay.
Figure 5b shows that the concentration in the clay
layer is considerably lower than in the sand layer atboth interfaces. This is caused by the use of two different parameter values in the ADE at the interface(i.e., different values for D/R and v/R). One could argue that, at least microscopically, the concentration in
the liquid phase should be continuous at the interface.
The use of a second-type condition is, therefore, preferred, although the profiles in Fig. 5b and c are very

951

LEU ET AL: ONE-DIMENSIONAL SOLUTE TRANSPORT
Sand

1.0

Ctay

0.8

0.6

l-10d

0.4

0.2
0.0

Ctay-sand-ctay

1.0

l-lOOd
0.8
0.6
0.4
0.2
0.0

1.0

Clay-sand-clay

0.8
0.6
0.4
0.2
0.0
10

20

X(cm)

_ i-20d

10

20

30

X(cm)

Fig. 5. Relative concentration (C/C0) as a function of distance (x)for different times (t) in (a) a hypothetical sand without a clay
layer, (b) with a hypothetical clay layer and a third-type condition,and (c) with a clay layer and a second-type condition.

Fig. 6. Relative concentration (C/C0) as a function of distance (x)for different times (t) in (a) a hypothetical clay without a sandlayer, (b) with a hypothetical sand layer and a third-type condition,and (c) with a sand layer and a second-type condition.

similar, except for the clay layer. Locations of the front
for zero dispersion (shown by the five vertical dashedlines in Fig. 5b) were included to show the effectivenessof the added adsorption capacity of the clay layer. The
resulting step fronts correspond to times of 2, 4.25, 6,
7.75, and 10 d. The solute moves more slowly throughthe sand when no dispersion occurs. Because of dispersion, the adsorption capacity of the clay layer is not
fully utilized. This example shows that the magnitudeof the dispersion coefficient needs to be known for
accurately predicting solute movement under these
conditions.
Figure 5 showed results for situations where a finetextured layer may retard solute transport in an otherwise coarse-textured, relatively wet soil. Solute transport may also be slowed by including a coarse-textured
layer in a relatively dry porous medium (Winograd,
1974). Figure 6 compares concentration profiles for ahomogeneous clay (Fig. 6a) with those for a clay-sandclay system using either a third-type (Fig. 6b) or asecond-type (Fig. 6c) condition at the interface of thesoil layers; a third-type condition was always used at
the inlet of the soil. The inclusion of the sand layer
was assumed to result in a tenfold reduction of thewater flow rate (Table 1) because of the low hydraulicconductivity of the sand. Since adsorption is now absent, the rate of solute movement is reduced to a lesser

extent. Because of the lower value of P for the threelayer system, considerable time is needed for the inlet
concentration to become equal to the input concentration. Notice that the discontinuity at the interface
behaves in an opposite manner than for the clay layer
in the previous example.
CONCLUSIONS
Analytical solutions of the one-dimensional ADE
were obtained for transport in a two-layer medium
assuming both a first- and a third-type condition at
the inlet boundary and the soil interface. Both layers
were assumed to be part of a semi-infinite system. In
addition, concentration profiles were obtained for a
combined first- and third-type condition at the interface, and a third-type condition at the inlet, using numerical inversion of the Laplace transform. In thiscase, the first layer was taken to be finite. The ADE
was also solved with a numerical finite-difference
scheme, assuming a third-type condition at the inlet
and a second- or third-type condition at the interface.
The solution of the ADE subject to concentrationtype (first-type) inlet and interface conditions may leadto substantial mass-balance errors if the solute profile
is subjected to a known solute flux. These errors are
avoided by using flux- or third-type conditions, lead-952

SOIL SCI. SOC. AM. J., VOL. 55, JULY-AUGUST 1991

ing to a discontinuity at the soil surface and at theinterface (the flux-averaged concentration, in that case,would remain continuous across the interface). A very
realistic solution is obtained by requiring that both theflux and the concentration be continuous. This implies
that, for the analytical solution, a first- and a thirdtype condition are imposed simultaneously at the interface, whereas, for the numerical solution, a secondtype condition must be used at the interface. Analytical solutions for those conditions can be closely approximated by solutions for a third-type condition.
Two numerical examples demonstrated that soil stratification (i.e., soil layers with different physical orchemical properties) can have major impacts on thetransport of solutes in soils.
Appendix A. Numerical Integration with
the Gauss-Chebyshev Formula
Numerical integration of the convolution integrals appearing in Eq. [ 1 7a] and [ 1 8a] was performed with the GaussChebyshev method, which is generally more accurate at theend of the integration interval than Gaussian quadrature(Vandergraft, 1983). Also Gauss-Chebyshev offers greater
flexibility in choosing the number of integration points.
The integration variable T, varying between 0 and t, needs
to be transformed into a variable ?/ varying between — 1 and+ 1. This was accomplished by the linear transformation T= t(i) + l)/2. For a function f(r) the integration changes asfollows:

mixed, and situations involving relatively high pore-watervelocities. The total solute influx, IN, across a unit crosssectional area at x = 0 is then denned as
IN =(9iV,C0 dT[Bl]

which should be equal to the amount of solute accumulated
in the profile, AC, as predicted by one of the analytical solutions
AC = f fl,0,C dx+ I R262C dx = AC, + AC2[B2]

A relative mass-balance error can be determined according
to (see van Genuchten and Parker, 1984):[B3]

Equation [B2] is most easily integrated in the Laplace domain, and then inverted back into the real-time domain. The
Laplace-domain solutions for first-type inlet and interface
conditions, for which the solutions in the real-time domain
were given by Eq. [9], [10], [16] and [17], are:C, (x,s) =exp (X,x)

C2 (x,s) = -^ exp [X,L + X2 (x - L)]

[B4a]
[B4b]

O

f(r) dr = o

[Al]

f[t(n
-i

The Gauss-Chebyshev formula states that, for a function

Integration in the Laplace domain, followed by inversion of
the Laplace transform, yields:
AC, = 0,C0(A*,)1/2 {exp(-a?0+[A2]

f (1 -

J

where Ej is the truncation error, and t\ and a are the integration points and weighting coefficients, respectively:= COS

1 - exp-k -(y - IK2n

j = 1,2,... ,n4t+ 2att - { — + ait I erfc (a,/1/2)

[A3]

fl

'"n
where « is the number ofintegration points.
Using g(r,) = (1 - i/2)"2 ffo), Eq. [Al] can then be evaluated according to
1V2,

[B5a]

k - 2a, (t - T)2(t - r)'/2[A4]

Appendix B. Mass Balance Calculations for
the Analytical Solutions
The validity of the imposed boundary conditions and associated solutions can be evaluated by formulating a massbalance equation in the Laplace domain, followed by analytical inversion of that equation. For gl = g2 = 0 and ta—> °°, assume that the solute influx is equal to the solute flux
just outside the inlet boundary of the soil, and solely determined by advective transport with flowing water. Ignoring
diffusion and dispersion (back diffusion) into the supply solution is realistic for most practical cases, including drip
irrigation, leaching from solution reservoirs that are well

+ exp (2atk) erfc

fc + 2a, (f - T)2(t - r)1/2[B5b]

where ak = vk (4R*A)-1/2 and k = LCR,/A)1/2. Equation [38]can be readily derived from Eq. [B3] and [B5].
For third-type conditions at the inlet and interface, thesolutions in the Laplace domain are:C, (x,s)— exp(X.x)
5

[B6a]

LEIJ ET AL: ONE-DIMENSIONAL SOLUTE TRANSPORT

C2 (X,S) =2th

a2 + (al + s)1/2 s

X exp [X,L + \2 (x - L)][B6b]
Substitution of Eq. [B6] into Eq. [B2] and inverting gives
AC = 0,v, C0t[B7]
A similar approach was followed for the problem involving a finite first layer and the simultaneous use of first- andthird-type conditions at the interface. However, the Laplace
inversion must now be obtained numerically, since no
closed-form expression in the real-time domain is available.

953