Boundary Conditions for Displacement Experiments through Short Laboratory Soil
Columns1
M. TH. VAN GENUCHTEN AND J. C. PARKER2ABSTRACT

Additional Index Words: dispersion coefficient, retardation factor.

lation to column displacement experiments is stilllacking. In this paper we will present a detailed analysis of several analytical solutions that previously havebeen applied to column tracer studies. Using basic
principles of mass conservation, we will show that atleast two of these solutions are inappropriate for analyzing column effluent data.
Recently, a number of papers in both the petroleum
and chemical engineering literature have discussed invery specific terms the physical and mathematical significance of various boundary conditions (Brigham,
1974; Baker, 1977;KreftandZuber, 1978, 1979;Kreft,
1981; Nauman, 1981a,b; among others). In particular,
these papers have pointed out the need to distinguishbetween volume-averaged or resident concentrations
and flux-averaged or flowing concentrations. As will
be discussed further in this paper, an understandingof the fundamental difference between these two concentration modes is pertinent to the selection of properboundary conditions and associated analytical solutions. Recommendations will be given as to which setof boundary conditions are most appropriate for analyzing column effluent data.

Van Genuchten, M. Th., and J.C. Parker. 1984. Boundary conditions for displacement experiments through short laboratory soilcolumns. Soil Sci. Soc. Am. J. 48:703-708.

The following partial differential equation is generally used

This paper presents a discussion of the physical and mathematicalsignificance of various boundary conditions applicable to onedimensional solute transport through relatively short laboratory soilcolumns. Based on mass balance considerations, it is shown that afirst-type or concentration-type condition at the inlet boundary incorrectly predicts the volume-averaged or resident concentration inside both semi-infinite and finite systems. A third-type or flux-type
inlet boundary condition preserves mass in semi-infinite systems,
but underpredicts effluent concentrations from finite columns unlessa local transformation is used to convert volume-averaged concentrations into flux-averaged concentrations. This transformation leads
to an expression for the effluent concentration that is identical to
the solution for the semi-infinite system using a concentration-typeboundary condition. For column Peclet numbers greater than about
five, the resulting analytical expression for the effluent curve is shown
to be nearly identical to the analytical solution for a finite system
based on a flux-type inlet boundary condition and a zero-concentration gradient at the exit boundary. Both solutions correctly preserve
mass in the system; other solutions of the convective-dispersive
transport equation are shown to be inappropriate for analyzing column effluent data.

THEORETICAL
to describe one-dimensional solute transport through a homogeneous medium during steady-state flow
R(dCr/dt) = D(d2Cr/dx2) - v(dCr/dx)[1]
where Cr is the volume-averaged (resident) solution concentration, D is the dispersion coefficient, y is the average porewater velocity, x is distance, and t is time. The retardation
factor R is included to account for linear and reversible equilibrium adsorption. Equation [1] will be applied to finite andsemi-infinite systems that are initially free of any solute:

ROPER formulation of boundary conditions for theP
analysis of column displacement experiments inthe laboratory is critically important to the interpretation of observed data, as well as for subsequent extrapolation of the experimental results to transport
problems in the field. Unfortunately, the most appropriate boundary conditions for a given displacement
experiment are not always intuitively apparent and inmany instances vastly different solutions have been
applied to purportedly identical physical problems. In
reducing the physical complexity of flow and transport
near or across a given boundary, certain simplifying
assumptions are necessary. The exact nature of theseassumptions is often not immediately apparent or appreciated. As pointed out by Bear (1979), intuitively
reasonable boundary conditions may lead to incongruent results, especially for relatively short soil columns. Although a few papers in the soil science literature have dealt with the effects of various boundary
conditions on solute transport (e.g., van Genuchten
and Wierenga, 1974; Parlange and Starr, 1975, 1978),a systematic analysis of the formulation and physicalinterpretation of assumed boundary conditions in re-

CfoO) = 0.

[2]

Analytical solution of Eq. [1] and [2] for various boundary
conditions will be discussed now.

Semi-Infinite Systems
Suppose we add water from an entrance reservoir having
concentration C0 to a semi-infinite profile (0 < x < oo) at
a rate equal to the volumetric flux q = 6v, where 6 is thevolumetric water content of the profile. Mass conservation
across the inlet boundary requires[- D(dCr/dx)
= vC0
[31
+
where x = 0 indicates evaluation at the inlet boundary justinside the profile. Equation [3] is correct for a system inwhich the entrance reservoir physically is not connected tothe column, for example, when the applied solution is slowly
trickled onto the surface of the column (compare with constant rainfall). The equation is also valid for systems where
the column is connected directly to the entrance reservoir,
provided the assumption is made that molecular diffusionand possibly some dispersion in the fore section (x < 0) canbe ignored. This means that the fore section is a perfectly
mixed reservoir. In that case, concentrations across the inlet

1
Joint contribution from the U.S. Salinity Lab., 4500 Glenwoqd
Drive, Riverside, CA 92501, and the Dep. of Agronomy, Virginia
Polytechnic Inst. and State Univ., Blackburg, VA 24061. Received
102 Nov. 1983. Approved 2 Mar. 1984.
Research Soil Scientist, U.S. Salinity Lab., and Assistant Professor, Dep. of Agronomy, Virginia Polytechnic Inst. and State Univ.,
respectively.

703

704

SOIL SCI. SOC. AM. J., VOL. 48, 1984

boundary are macroscopically discontinuous (Danckwerts,
1953; Wehner and Wilhelm, 1956). Equation [3] specifies a
Robin, third-type, or flux-type boundary condition.
For semi-infinite systems, we need an additional boundary condition specifying the behavior of Cj(x,t) when x —>oo. For our discussion it is sufficient to require that(dCJdx) (oo,0 = 0 .
[4]
The analytical solution of Eq. [1] through [4] is (e.g., see
Lindstrom et al, 1967)

Equation [10] can be derived also by considering the actual
solute flux density Js into the column at x = 0:
Js(0,t) = [qC0-9D(dCr/dx)]\x=0. .
[12]This equation indicates that in addition to the convective
mass flux qCm a dispersive component is operative at theboundary. This dispersion term is responsible for the observed mass balance deviations, as can be shown by defining
the relative error E in terms of J~
£=

Co

2

Rx-vt
[2(DRty/2\

\vDR

[13]
X exp [~x =0

J ~ 2 I 1 + ~D + ~DR }

4DRt

Xexp

erfc

[5]

2(DRt) l/2

Equation [5] correctly evaluates volume-averaged or resident concentrations in semi-infinite field profiles. One may
verify that [5] satisfies the mass balance requirement
Too\CjL = R */oI C£x,t) dx.

Substituting Eq. [8] into Eq. [13] also leads to Eq. [10]. Figure 1 shows that the relative error in the mass balance canbe extremely large for small values of the dimensionless parameter f, i.e., at early times, for systems where dispersive
transport is large compared to convective transport, and forsystems where the solute is strongly adsorbed by the solid
phase (large R). Hence, concentration-type boundary condition Eq. [2] fails to satisfy conservation of mass and thediscrepancy is not always negligible.

[6]

In other words, whatever is added at the surface (term on
the left) must be found inside the profile (term on the right).
Instead of Eq. [3], a first-type or concentration-type input
boundary condition has also been used frequently:
Q0,r) = C0 .
[7]This equation assumes that the concentration itself can be
specified at the soil surface, a situation that usually is notpossible in practice. In essence, Eq. [7] indicates that theconcentration is continuous across the inlet boundary. While
the stipulation of such, a concentration continuity at x = 0

Finite Systems
The challenge here is to define a boundary condition that
describes as realistically as possible the physical process offlow and transport across the exit boundary of a finite system. One requirement that always must be satisfied is continuity of the solute velocity across the lower boundary:where L~ indicates evaluation just inside the column, andwhere C, = Cj(t) is the exit concentration. As with Eq. [3]for the inlet boundary, Eq. [14] assumes that diffusion-

may have an intuitive appeal, it is obtained at the cost ofdispersion phenomena in the exit reservoir are negligible, or

maintaining a solute flux continuity across the inlet position.
Consequently, Eq. [7] is incorrect for displacement experiments where the tracer solution is added at a specified rate.
This becomes evident when considering the analytical solution of Eq. [1], [2], [4] and [7] (Lapidus and Amundson,1952):

that the exit reservoir is not in direct physical contact withthe liquid phase inside the column. This last situation may
occur when the effluent solution drips freely from the exit
boundary (e.g., see the fraction collector system describedby Wierenga et al., 1975).
Because of the addition of an extra unknown (Ce), Eq.
[14] leads to an indeterminate system of equations. Hence,an additional equation is needed to fully describe the system. One such equation is based on the intuitive assumptionthat the solute concentration should be continuous acrossthe lower boundary:Rx - vt
vt
l/2

[8]

Substitution of this equation into Eq. [6] shows that themass balance identity is not satisfied. If we define a relative
error E as (see also Eq. [6])CX*,0 dx - vC0t]

E = (1/vCO[9]) = ex/).
--

[15]

2.0
1.6

3 '-2

dx -I

00

then, upon integration, it is easily shown that for Eq. [8]

tn

CO

1
UJ

0.8
0.4

>

[10]
DIMENSIONLESS TIME, f=v a t/DR

where

f = v2t/DR.

[11]

Fig. 1—Plot of the relative mass balance error, E, vs. dimensionlesstime f for analytical solution [8].

VAN GENUCHTEN & PARKER: BOUNDARY CONDITIONS FOR COLUMN DISPLACEMENT EXPERIMENTS

Substitution of this equation into [14] leads to the frequently
used boundary condition[16](dCJdx) (L,f) = 0 .
It is important to realize that Eq. [16] is based on an assumption, namely that the concentration macroscopicallyshould be continuous at x = L (Danckwerts, 1953; Wehnerand Wilhelm, 1956, Pearson, 1959). It is difficult to reason
a priori whether this assumption is reasonable in view of
the uncertainty of the exact physical processes at or near the
lower boundary. However, the requirement of a continuous
distribution at x = L seems inconsistent with Eq. [3] forthe inlet boundary which not only results in a discontinuous
distribution at x = 0, but also in a nonzero concentration
gradient at the column entrance.
Nevertheless, if we accept Eq. [16] as an appropriate
boundary condition for finite columns, then the physical system is completely defined and can be solved. The analyticalsolution of Eq. [1], [2], [3] and [16] is (Brenner, 1962)
C,(x,t)705

the presence of an outflow boundary or effluent collection
system, thus considering the column to be part of an effectively semi-infinite system. The requisite additional stipulation that dC£<x>,t)/dx = 0 may then be invoked, which inturn suggests that Eq. [5] adequately describes the concentration distribution within the column. However, solution[5] does not immediately predict the effluent concentration
since we have no longer imposed Brenner's assumption that
C/0 = C£L,t). An expression for the effluent concentration
can nevertheless be obtained by substitution of Eq. [5] directly into Eq. [14] which yields
RL- vtC.

2(DRt)l/2[23,
or in dimensionless formCe(T)/C0= l

- T)]

Cn

[24]

2vL

2

D

«2

m=\

Pm

2D

D

Note that Eq. [23] is exactly the same as Eq. [8] when evaluated at x = L. Equation [23] can be derived also by considering Ce to be the unknown closing term of mass balance
Eq. [22]. Differentiating that equation with respect to timeand solving for the effluent concentration gives

2D

vx

[25]
where the eigenvalues ftm are the positive roots of
0m cotCS J - G&ZVvZ.) + (vL/4D) = 0. [18]If applied to effluent curves (x = L) and using dimensionlessvariables, the solution simplifies to
fl2'
>

Ce(T)

PT"

Pm

CJt) = —= 1- m=\
2
[19]
[20]

T = vt/L

P = vL/D[21]

where T is the number of pore volumes leached through thecolumn, and P is the column Peclet number. Computational
programs for evaluating Brenner's series solution are readily
available (van Genuchten and Alves, 1982). The series solution converges only for relatively small values of P. For
large P-values, computationally efficient and very accurate
approximate solutions have also been obtained.
Analogous to Eq. [6] for semi-infinite systems, the massbalance requirement for finite columns isv('[C0-Ce(T)]dT=R(LCr(x,t)dx.
JQ

In other words, the effluent concentration equals the input
concentration minus the change in the amount of solute
stored in the column. Substituting Eq. [5] into Eq. [25] andcarrying out the indicated integration leads directly to Eq.
[23]. A third equivalent expression for Cj(t) is given by

[22]

JQ

This equation states that whatever is added to the column
minus what is leaving the column (left side) must be stored
in that column (right side). Upon substitution, it is readily
shown that Brenner's solution satisfies [22].
We have emphasized that Brenner's solution for Ce satisfies Eq. [3] for the inlet boundary, Eq. [14] for the outlet
boundary, and that the requisite
system of equations is completed by assuming that Q(0 = Cj(L,i). An alternative wayof completing the system of equations is to assume that solute distributions inside the finite column are unaffected by

dx

[26]

which states that the effluent concentration equals the time
derivative of the amount of solute that would be stored in
the region L < x < oo if the column were to be part of asemi-infinite system. Equation [26] can be obtained by taking the time derivative of Eq. [6] and using the resulting
expression to eliminate C0 from Eq. [25].
The analysis above indicates that the flux-type inlet
boundary condition in Eq. [3] properly conserves mass whenused to calculate volume-averaged or resident concentration
distributions in semi-infinite systems. The analytical solution for this problem (Eq. 5) can be used also to derive an
expression for the breakthrough curve, Cj(t), from a finite
column; the expression turns out to be identical to the analytical solution for a first-type input boundary condition
and a semi-infinite system (Eq. 23). This relationship between the two semi-infinite solutions has been described previously (Brigham, 1974; Kreft and Zuber, 1978; Parker andvan Genuchten, 1984; among others). It is a direct consequence of a transformation similar to [14] but applicable atany given distance, x. This transformation will be discussed
next.
Flux-Averaged vs. Volume-Averaged Concentrations
Equation [1] is written in terms of the in-situ, volumeaveraged, or resident concentration, Cn which is the massof solute per unit volume of fluid contained in an elementaryvolume of the system at a given instance. One can also define a flux-averaged or flowing concentration, Cfi which is
the mass of solute per unit volume of fluid passing through706

SOIL SCI. SOC. AM. J., VOL. 48, 1984

a given cross-section during an elementary time interval(Kreft and Zuber, 1978):

Cf=Cr- (D/v) (dCr/dx)

[27]
vx

Substitution of Eq. [27] into Eq. [1] leads to exactly the sametransport equation, except that the resident concentration C,
is replaced by the flux concentration, C/.

R (dCJdt) = D (d2Cj/dx2) - v (dCjdx).

[28]

Hence, spatial and temporal distributions of Cr and C}are
described by the same governing equation. Boundary condition [4] also holds for both C, and C/, this follows immediately by noting that the second spatial derivative of [5]vanishes when x —» oo. Finally, application of Eq. [27] toboundary condition [3] leads to
Cj(0,t) = C0.
[29]
Hence, transport equations that describe resident concentrations in semi-infinite systems subject to first-type inlet
boundary conditions are mathematically identical to transport equations that describe flux concentrations in semiinfinite systems subject to first-type inlet boundary conditions. This observation shows that the solution of Lapidus
and Amundson, given by Eq. [8], predicts flux concentrations in a semi-infinite system, and if extended to a finite
column, properly describes the effluent concentration.
Equation [27] can be applied also to Brenner's solution
for a finite system. Substituting Eq. [17] into Eq. [27] yields
the following expression for the flux concentration Cf
C
VX

V2t

2D

4DR

2D

D

s

L2R

4DR

L2R

m=\

ID

where the am are the positive roots ofaffl cot(am) + (vL/2D) = 0 .

[32]

Note that we ostensibly assumed that Cleary and Adrian's
solution represents flux concentrations rather than resident
concentrations. Because of boundary condition [16], C^L,t)= Cj(L,t) = Cj(t) for Cleary and Adrian's solution, as was
the case for Brenner's solution.
Figure 2 shows that the solution of Cleary and Adrian(1973) deviates substantially from the flux concentration (Eq.
[30]) as derived from the solution of Brenner (1962). Thisdeviation is caused by the fact that boundary condition [16]is not invariant by the transformation from C, to C/ (thesecond spatial derivative at x = L is not zero). Thus, if
Brenner's solution is taken as the correct expression for the
resident concentration, then Cleary and Adrian's solution
leads to an incorrect description of the flux concentration.
Alternatively, one may also argue that if Cleary and Adrian's
solution were to represent flux concentrations, then Brenner's solution cannot represent resident concentrations. This
inconsistency gives additional support to our contention thatEq. [16] is based on a convenient but purely intuitive assumption.
DISCUSSION
To appreciate the differences between various ana-

lytical solutions that are potentially applicable to column displacement experiments, Fig. 3 shows calcu-

m=l

[30]
where the eigenvalues f)m are the same as before (Eq. 18).
Note that Eq. [30] at the exit boundary x = L is the sameas Eq. [19]. This follows also immediately from Eq. [27] andboundary condition [16]. Hence, flux and resident concentrations as calculated with Brenner's solution are identicalat the exit boundary. However, they can be quite different
inside the column. This is shown in Fig. 2 where the twoconcentration modes are plotted versus reduced distance(x/L) using a value of 5 for P. Figure 2 also shows the solution of Cleary and Adrian (1973) for boundary conditions[7] and [16], i.e. (see also van Genuchten and Alves, 1982),

o
•z.
o
cc
h-

z
LJ

O

z
o
o

0° 1.0&LJ

z 0.8<£ 0.6

LJ

a:

0.4

o

s
C r - BRENNER-..
C, - BRENNER
C, - CLEARY AND ADRIAN

o
0.2

0.8

0.4

0.6

0.8

1.0

REDUCED DISTANCE, x/L

Fig. 2—Distributions vs. reduced distance of relative volume-averaged and flux-averaged concentrations predicted with Brenner's
solution, and relative flux-averaged concentrations predicted withthe solution of Cleary and Adrian.

PORE VOLUMES, T
Fig. 3—Relative effluent concentration profiles for analytical solutions CA (Cleary and Adrian, 1973), LA (Lapidus and Amundson,1952), LB (Lindstrom et al., 1967) and BR (Brenner, 1962). The
effluent curves are plotted for three values of the column Peclet
number, P.

VAN GENUCHTEN & PARKER: BOUNDARY CONDITIONS FOR COLUMN DISPLACEMENT EXPERIMENTS

lated effluent curves for three values of the column
Peclet number, P. Plotted are the solutions of Clearyand Adrian (1972), Lapidus and Amundson (1952),
Lindstrom et al. (1967), and Brenner (1962), referred
to respectively by CA, LA, LB, and BR and given inthis same order by Eq. [31], [8], [5], and [17]. Notethat the four solutions deviate drastically when P =1, but slowly converge to each other as P increases.
Analytical solutions LA and BR are very similar whenP = 5, while essentially no differences between these
two solutions are present when P reaches a value of20.
An important attribute of a breakthrough curve froma finite column is the area above each curve. This area,
sometimes referred to as the holdup (H), should represent the amount of solute that can be stored in thecolumn. Mathematically, H is given by[1 - CJ,T)/C0] dT.

[33]

One may verify that H = R for the two curves designated LA and BR in Fig. 3. For the curve CA (Eq.
[31]), His given by
H = R[l - (l/P) - (e~p/P)][34]
while for curve LB (Eq. [5])H = R[l + (l/P)].
[35]
Figure 4 gives a plot of the relative holdup (H/R)
vs. P for the four analytical solutions for CJ(T). Notethat H/R deviates substantially from 1.0 for the solutions of Cleary and Adrian (1972) and Lindstrom etal. (1967). These deviations illustrate that highly inaccurate estimates for R can be obtained when thesolutions CA or LB are fitted to observed effluent data,
especially when P is small. For example, analyticalsolution LB overestimates R by no less than 50% whenP = 2.
The solutions LA and BR always lead to correct
estimates for R, independently of the value of the Peclet number. However, the two solutions may yieldsomewhat different estimates of P, especially when P
becomes less than about 5 (Fig. 3). As discussed earlier, LA and BR each are based on different assumptions regarding the physics of flow and transport at ornear the lower boundary. Because of the apparent inconsistency in the stipulation of finite and zero concentration gradients at the inlet and outlet boundaries
for Brenner's solution, and because of its numericallymore tedious form compared to Eq. [8], we recomOL

2.0

I
Q."

1.5

Q
_l

o

1.0

Ul
>

1.5
UJ

tr

0

4

8
PECLET

12
NUMBER,

16

20

P

Fig. 4—Relative holdup (H/R) vs. column Peclet number, P, for analytical solutions CA, LA, LB and BR.

707

mend that only the solution of Lapidus and Amundson (1952) be used to calculate effluent curves fromfinite columns. Of course, this same solution (Eq. [8])should also be used to calculate solute fluxes in thefield. Similarly, we recommend that only the solutionof Lindstrom et al. (1967) be used to evaluate in situ,
volume-averaged concentrations.
In some cases, breakthrough curves are obtained
within laboratory columns by means of suction cups
or other extraction devices. In that case, it is difficult
to reason which concentration mode will be observed.
Because the extracted soil solution is obtained at a
fixed location in the soil, it is unlikely that the observed data are flux concentrations. However, becauseof the transient behavior of displacement experiments
and the uncertainty of how exactly flow lines are disrupted by the installation and performance of suction
cups, the observed concentration mode is probably
not exactly that of a resident concentration either. This
problem is pertinent for both laboratory and in situ
field measurements.
It should be noted that transformation [27] fromresident to flux concentrations has been shown to hold
only for Eq. [1] and its analytical solutions for a semiinfinite system (LB and LA). The transformation may
not necessarily hold for more complicated transportmodels (e.g., for certain kinetic adsorption models).
In that case, flux concentrations—and hence breakthrough curves—are best calculated by first deriving
the analytical solution for the resident concentrationand subsequently applying Eq. [14] to obtain an
expression for the effluent curve.
SUMMARY AND CONCLUSIONSThis study shows that solute transport in a semiinfinite field profile is described properly with analytical solution [5] for a flux-type soil surface input
boundary condition. The analytical solution for thisproblem can be used also to describe resident (volume-averaged) solute distributions inside finite soilcolumns. However, the solution needs to be modified
if applied to effluent curves from finite columns. Basedon mass balance calculations, and assuming that thecolumn exit boundary does not alter flow lines upstream of the boundary, an equation for the breakthrough curve (Eq. [23]) can be derived that is mathematically similar to the solution for a concentrationtype soil surface inlet boundary condition. This similarity is not accidential but results from a transformation from volume-averaged (Cr) to flux-averaged(C/) concentrations. A similar transformation fails for
two available solutions applicable to finite systems.
Of the two solutions for a finite system, only thesolution of Brenner (1962) for the flux-type inlet
boundary condition conserves mass in the system. Thissolution leads to an expression for the effluent curvethat is nearly identical to the derived expression [23]for the breakthrough curve from a finite profile, provided the column Peclet number is not much less thanabout 5.
Because of some inconsistencies in the stipulationof concentration continuities at the inlet and outlet
boundaries for Brenner's solution and because of its
numerically more tedious form compared to the so-708

SOIL SCI. SOC. AM. J., VOL. 48, 1984

lution of Lapidus and Amundson (1952), we recommend that the latter solution always be used to evaluate flux-averaged concentrations, whether they
pertain to finite systems (effluent curves) or semiinfinite field profiles. We similarly recommend that
only the solution of Lindstrom et al. (1967) be usedfor volume-averaged concentrations.