SOIL SCIENCE SOCIETY OF AMERICAJOURNAL

VOL.

JULY—AUGUST

40

1976

No.

4

DIVISION S-l—SOIL PHYSICS
Mass Transfer Studies in Sorbing Porous Media I. Analytical Solutions1
M. TH. VAN GENUCHTEN 2 AND P. J. WlERENGA 3

symmetrical concentration distributions (Coats and Smith,
1964; Gershon and Nir, 1969). However, several experimental studies have shown strongly nonsigmoid and asymmetrical distributions, either vs. depth (Aylmore andKarim, 1970) or as effluent curves (Kay and Elrick, 1967;
Green etal., 1972). Following others (cf. Giddings, 1963),
we will use furtheron the word "tailing" to indicate the appearance of such highly asymmetrical or nonsigmoid concentration profiles. Tailing has been observed under several
experimental conditions:
1) Unsaturated Flow—Nielsen and Biggar (1961) noted
considerable tailing with decreasing water content at approximately the same flow velocity. They argued that underunsaturated conditions the larger pores are eliminated fortransport and that the proportion of water which does notreadily move within the soil is increased. This water has
been identified as dead, stagnant or immobile water (Turner, 1958; Coats and Smith, 1964; Deans, 1963). A decrease in water content increases the fraction of air-filled
macro-pores, resulting in the creation of additional deadend or blind pores which rely on diffusion processes to attain equilibrium with the displacing solution (Gupta et al.,1973).
2) Aggregated Media—Soils are composed of slowly
and rapidly conducting pore sequences. Aggregates have
many micropores in which displacement is dependent upon
diffusion, since convection in these smaller pores is usually
negligible. This results in slow and incomplete mixing andhence tailing, even under saturated conditions. In large
aggregates the amount of immobile water increases while
the diffusion pathway becomes longer, resulting in

ABSTRACTAn analytical solution is presented for the movement of chemicals
through a sorbing porous medium with lateral or intra-aggregate diffusion. The liquid phase in the porous medium is divided into mobile
and immobile regions. Diffusional transfer between the two liquid
regions is assumed to be proportional to the concentration difference
between the mobile and immobile liquids. Sorption processes in both
the dynamic and stagnant regions of the medium are assumed to be instantaneous and the adsorption isotherm is assumed to be linear. The
analytical model derived here describes the extensive tailing observed
during flow through an unsaturated, aggregated sorbing medium andexplains the often observed early breakthrough of chemicals in theeffluent.
Additional Index Words: intra-aggregate diffusion, adsorption, analytical solution, immobile water, and miscible displacement.

ISPLACEMENT STUDIES have become important tools insoil physics, particularly for predicting the movementDof pesticides, nitrates, heavy metals, and other solutes
through soil. These studies can provide valuable information about both the behavior of a chemical in a porous medium (diffusion, dispersion, anion exclusion, adsorption, or
exchange processes), and about the medium itself (pore size
distribution, aggregation). Much of the work which has
been done to date can be found in reviews by Biggar andNielsen (1967), Fried and Combarnous (1971), and by Bear(1972). Attempts to formulate displacement processes
quantitatively have been based generally on the convectivedispersive equation (Lapidus and Amundson, 1952)
dC

d2C

dC

[1]

'Journal article no. 564, Agric. Exp. Stn., New Mexico State Univ., Las
Cruces, NM 88003. The work upon which this report was based was supported in part by funds obtained from the U.S. Department of the Interior,
Office of Water Resour. Res., as authorized under the Water Resources Act
of 21964. Received 7 July 1975. Approved 15 Mar. 1976.
Formerly Graduate Student, Dep. of Agronomy, New Mexico State
Univ.; presently Research Associate, Dep. of Civil Engineering, Princeton
Univ., Princeton, NJ 08540.
'Associate Professor, Dep. of Agronomy, New Mexico State Univ., Las
Cruces, NM 88003.

where C is the concentration (meq/cm3), D is the dispersion
coefficient (cmVday), v0 is the pore-water velocity (fluxdivided by the volumetric water content 0; cm/day), z is the
distance (cm), and t is time (days). Except for small valuesof vg/D, all solutions of Eq. [1] predict nearly sigmoid or
473

SOIL SO. SOC. AM. J., VOL. 40, 1976

474

breakthrough curves with more and longer tailing. (Biggar
and Nielsen, 1962; Green et al., 1972; McMahon andThomas, 1974).
3) Pore-water Velocity—Several experimental studies
indicate that tailing becomes more pronounced as the porewater velocity decreases (Biggar and Nielsen, 1962; Villerrrtaux and van Swaay, 1969; Skopp and Warrick, 1974).
From the discussions above it is evident that solutions of
Eq [1] may not always give an accurate description of thephysical processes which occur in aggregated and/or unsaturated soils. Several attempts have been made recently to
account for the observed tailing. Deans (1963) modified eq.
[1] to include transfer by diffusion from dynamic flowing
regions to stagnant zones, but neglected longitudinal dispersion. He used a first-order mass transfer equation to describe this diffusional exchange. Similar concepts were
proposed by Turner (1958), Gottschlich (1963), and more
recently by Skopp and Warrick (1974), who each describedthe diffusional transfer of material from dynamic to stagnant
fluid by the diffusion equation. Of immediate interest in ourstudy is the model described by Coats and Smith (1964),
who expanded Deans' (1963) model to include longitudinal
dispersion. Their model uses the following equations:

dt

dt

dCt,

dz

where &m and 9tm are the fractions of the soil filled withmobile and stagnant water, respectively (cm3/cm3); Cm andCim are the concentrations (meq/cm3) in both the mobile andstagnant regions; vm is the average pore-water velocity inthe mobile liquid; and a is a mass transfer coefficient(day"1). Equation [2] is the same as Eq. [2] of Passioura(1971). Solutions of Eq. [2] and [3] for a semi-infinite column were obtained by Coats and Smith (1964) and byVillermaux and Swaay (1969), the latter authors using a
tracer pulse input, approximated by the Dirac 8-function.
Bennet and Goodridge (1970) obtained the solutions for a
finite column. The last three studies mentioned here eachobtained excellent descriptions of experimental concentration-time curves.
When a sorbing porous medium is considered, an additional term must be added to the left-hand side of Eq. [1] toaccount for the interaction between the chemical and the
medium. If the sorbed concentration is denoted by S (/Ag/gsoil), the general transport equation can be written in thefollowing form:

p dS

e dt

S =

[5]

where K and N are constants. Substitution of Eq. [5] into
Eq. [4] yields a transport equation in terms of one dependent variable, C:

dC

dC

[6]

where the retardation factor R is defined as (Hashimoto etal, 1964)= ! +pKNCN~l/e.

[7]

Instead of equilibrium, several authors have assumed a kinetic approach towards equilibrium. A widely quoted model
stems from Lapidus and Amundson (1952), who proposed
the following rate expression:

[8]

[2]
[3]

= <*(Cm —rate of adsorption (dS/dt). When equilibrium is assumed
between adsorbed and solution concentration, the Freundlich isotherm may be used. This equation is as follows:

Solutions of Eq. [4] and [8] for a semi-infinite medium wereobtained by Lapidus and Amundson, using a step input inconcentration, and by Lindstrom and Narasimhan (1973)
for an initially distributed chemical. Note that Eq. [4] and[8] are mathematically similar to Eq. [2] and [3].
Kinetic models derive their rationale from the assumptionthat the rate of adsorption is slow compared to the rate at
which the chemical moves through the soil. In an earlier
study (van Genuchten et al, 1974), experimental data werecompared with predictions based on Eq. [4], using both theequilibrium and the kinetic models (Eq. [5] and [8], respec-

B
dC

dC

dz

[4]

where p is the bulk density (g/cm3). Equation [4] is an extension of Eq. [1] in that only a sorption term is added.
Because of this, the limitations of Eq. [4] in describing the
movement of a chemical through an unsaturated, aggregated sorbing medium should be similar to those foundfor Eq. [1] in a nonsorbing medium.
Solutions of Eq. [4] can be obtained by specifying the

DYNAMIC SOIL REGION!^^^^9s/Ny\^X
MOBILE WATER

Fig. 1—Schematic diagram of unsaturated aggregated porous medium. (A) Actual model. (B) Simplified model. The shading patternsin A and B represent the same regions.

VAN GENUCHTEN & WIERENGA: MASS TRANSFER STUDIES IN SORBING POROUS MEDIA: I.

tively). Both models were found to describe the effluent
data reasonably well at low pore-water velocities, but were
inadequate in describing the data at higher velocities. Because of the poor predictions, we suggested that only a certain fraction of the sorption sites, i.e. those located aroundthe larger pores, participated actively in the adsorption
process. This would have the same effect as a reduction in
the adsorption constant, leading to a relatively faster movement of the chemical and resulting in an earlier breakthrough. At the same time a considerable portion of thechemical must diffuse to the remaining sites, located insideaggregates and along dead pores. This slow diffusion
process would continuously remove material from the larger
water-saturated pores, resulting in extensive tailing. Basically the same reasoning was used in other studies to account for the deviations between the predictions based on
Eq. [4] and experimental data (Kay and Elrick, 1967; Greenet al, 1972). It is the purpose of this paper to further expand
the model proposed by van Genuchten et al. (1974).
Transfer equations are derived which include diffusion controlled adsorption, and an analytical solution is presented
for the case when the adsorption isotherm is linear.
THEORETICAL
Consider the movement of a chemical through an unsaturated,aggregated sorbing porous medium. Figure 1 represents schematically such a medium. Five different regions can be identified:
1) Air spaces.
2) Mobile (or dynamic) water, located inside the larger (interaggregate) pores. The flow of fluid in the medium is assumed tooccur in this region only. Solute transfer occurs by both convection and longitudinal diffusion.
3) Immobile (dead or stagnant) water, located inside aggregates
and at the contact points of aggregates and/or particles. In saturated media this region is mainly confined to intra-aggregate pores.
Note that air-bubbles and unsaturated conditions may increase theproportion of dead water by creating more dead-end pores.
4) A dynamic soil region, located sufficiently close to themobile water phase for equilibrium (assumed) between the solutein the mobile liquid and that sorbed by this part of the soil mass.
5) A stagnant soil region, that part of the soil matrix where
sorption is diffusion limited. This part of the soil is located mainly
around the micro-pores inside the aggregates, or along dead-end
water pockets. Sorption occurs here only after the chemical has
diffused through the liquid barrier of the immobile liquid phase.
We will characterize the division of the soil matrix into the tworegions, dynamic and stagnant, by the parameter/. Equivalently/
may be defined as that fraction of the sorption sites which is in
direct contact with the mobile liquid. On the basis of the equationof continuity we may now state the following general transport
equation:

dC

475

[10]

We further assume that adsorption is instantaneous and that therelation between adsorption and solution concentration can be
described by the Freundlich relationship (Eq. [5]). Hence

Substitution of Eq. [11] in Eq. [10] results in(0m +fpKNCmN-1)^

01..

+ [0lm + (1
.. d

Z/~^
<-m

n

%/^
OC f f[12]

A relationship between the concentrations of solute in the mobile and immobile liquids has to be established in order to complete the description of the transport mechanism. Diffusionaltransfer between these two liquid regions is defined here by a relation similar to Eq. [3]. Since a sorbing medium is considered andsince adsorption takes place from the immobile liquid to the immobile part of the porous medium a sorption term must be added to
the left-hand side of Eq. [3]. The resulting transfer equation is
m

at

m m .

at

[13]

Substitution of Eq. [11] in Eq. [13] leads finally to1]= «(Cm - C,M).

[14]

Linear adsorption
For linear adsorption (N = 1) Eq. [12] and [14] reduce to:[15]= «(Cm - Clm).

[16]

For a semi-infinite column and a pulse input of solute, the boundary and initial conditions are
Ost<tlv[17]

Hrn[C m (z,t)] = 0[18]

Cm(z,0) = Clm(z,0) = 0.

[19]

To- obtain a more amenable system of equations, the following
dimensionless variables are introduced:
where the subscripts m and t'm refer to mobile and immobile
regions-, respectively. This study assumes a constant water flux(q); hence D, Om, 0,-OT,/and vm (= q/9^) are constants and Eq. [9]simplifies to
dC,-m

dSi

T=Vot/L=-vmt<t>/L[20]

x = z/L[21]= 0m/0 = 0J(om + eim)

[22]
[23]

476

SOIL SCI. SOC. AM. J., VOL. 40,

P = v m L/D[24]

a =aL/(vmem) = aL/q.

[25]

[39]

Equations [15] and [16] can now be written in the following form:
<>£,„
dx

[26]

1976

c2(x ,T) = «r G(X,T

G(x,T) = erfc{(P/4/3RT)ll2(/3Rx

[27]

[40]

Jo

- T)}

- Px + PT/pR)e\p(Px)erfc{(P/4(3RT)l/2(pRxwhere Rm andRim are retardation factors based on the dynamic andstagnant regions, respectively (Fig. 1):+ T)}- T)2/4/3RT}

[41]
12

ff^r.r) = exp(-« - v) {/„(£) //3 + 71(|)(«/v) ' /(l - /3)} [42]

Rm =

[28]

-f)pK

Rim = I +fl2(T,T) =.exp(-« - v){/ 0 (f)/(l - ffi + /.(^(v/i*)1'*/^}

[43]

I/2

[44]

f = 2(«v)[29]

DISCUSSION
The average retardation factor R for linear adsorption is (Eq. [7]):R=\+(pK/6).

[30]

Define further the dimensionless parameter, ($, as
pfK

(8 =

6+pK

Om Rm

[31]

e R'

With these definitions of R and /3, the transport equations finallybecomeCim __

Cim

cm

-----

.).

, ,
[32]

[33]

The initial and boundary conditions in dimensionless form are
0 5= T < T,

r

hm \c
'

p

[34]

** / ~ IL 0\>m[cm(x,T)] = 0[35]cm(x,0) = cim(x,0) = 0[36]

One technique for solving Eq. [32]-[36] is to apply the Laplacetransform and then to treat the transformed equations as ordinary
differential equations. Lindstrom and Narasimhan (1973) haveapplied this technique to a similar mathematical problem and obtained the solutions of Eq. [4] and [8] for an initially distributed
chemical. With straight-forward, albeit lengthy, transformationsof their solutions, it can be shown that4 the solutions of Eq.
[32]-[36] become (van Genuchten, 1974, Appendix C):
cm(x,T) =[37]\c2(x,T)
ctm(x,T) = •0 < T < T!

[38]

Equations [37]-[44], representing the solute concentrations in both mobile and immobile regions of the soil,
depend on four dimensionless parameters: R, f t , a and P.
Both R and /3 (Eq. [30] and [31], respectively) depend on
the adsorption constant K, and hence are determined byboth the chemical and the porous medium. The masstransfer coefficient a will also depend on the chemical andthe medium, because this parameter involves the diffusion
coefficient of the chemical, the length over which theaverage concentration difference is considered and thecross-sectional area through which the lateral diffusion
process takes place (Coats and Smith, 1964). The column
Peclet number P, however, will depend solely on the medium, provided the contribution of longitudinal diffusion to
the dispersion coefficient is negligibly small (only mechanical dispersion).
In order to evaluate the different parameters and their
sensitivity, effluent curves were calculated for a 30-cm longsoil column (L = 30), leached with 3 pore volumes of achemical (7\ = 3) at a constant flux (q) of 10 cm/day. Ignoring exit effects and assuming that only mobile liquid leaves
the basal plane of the column, the effluent concentrations
can be calculated with Eq. [37], with x = 1. Results areshown in Fig. 2-5.
Figure 2 shows the influence of the fraction mobile water(<f>) on the shape of the effluent concentration distribution.
Values of the different parameters used in the calculations
are given in the figure. When the amount of immobile liquidincreases (i.e. <t> decreases), convective transfer is increasingly confined to a smaller cross- sectional area of the column, resulting in a higher solute velocity in the mobile liquid (vm = q/0m) and hence leading to an earlier breakthrough of the chemical in the effluent. Also, the slope ofthe effluent curve at the higher concentrations is affected bya change in (f>. When <f> = I, no immobile water is present,
i.e. no intra-aggregate diffusion takes place, and the curve
becomes more or less symmetrical (curve not included inFig. 2). Equations [15] and [16] then reduce to Eq. [6] and[7] with N = 1 and p replaced by /p. In that case/ becomes
identical to the variable FREQ ("fraction near equilibri4

M. Th. van Genuchten. 1974, Mass transfer studies in sorbing porous

media. Ph.D. Thesis. New Mexico State Univ, Las Cruces.

477

VAN GENUCHTEN & WIERENGA: MASS TRANSFER STUDIES IN SORBING POROUS MEDIA: I.

q = 10cm/day8=0.400=30 cm 2 /day
p= 1.30 a/cm 3
f = O.4O

a=O.I5

I/day

T,=

Kads=0.50

1

1

2

3

4

5

6

7

PORE VOLUME-T

Fig. 2—Calculated effluent curves for a sorbing medium, as influenced by the fraction mobile liquid. ($).

um") which van Genuchten et al. (1974) used in theircalculations to account for the early arrival of picloram in
the effluent. FREQ is hence that fraction of the sorption
sites which participates actively in the adsorption process,the remaining fraction being inaccessible for the displacing
fluid. However, as long as some immobile liquid is present(4> = 0.999 or <1.) extensive tailing will occur, since approximately 60% of the chemical (f = 0.40) has to diffuse
into the stagnant region of the soil, i.e. to those sorption
sites which are in equilibrium with the immobile liquid.
Approximately the same set of curves as shown in Fig. 2
is obtained when the parameter/is allowed to vary (Fig. 3).
This parameter describes the amount of adsorption taking
place inside the dynamic region of the soil, as a fraction of

the total adsorption. When/increases, i.e. when relatively
more adsorption occurs in the dynamic region and relatively
less inside the aggregates, the total adsorption remaining thesame, the chemical will appear later in the effluent. When/= 1, adsorption takes place only in the dynamic region of
the soil. Some solute may then still diffuse into the stagnant
region; however, no adsorption occurs here and the chemical can only be stored in the immobile liquid.
Both/(fraction adsorption sites in dynamic region) and </>(fraction mobile water) have the same effect on the shape ofthe effluent curves. This is not surprising since both parameters are included in the dimensionless variable /3 (Eq. 31)
which determines the fraction of the retardation occurring in
the dynamic region of the soil.

1.0 -

0= 30

cm /day$=0.65a =0.15

I/day

T,= 3
0.50

0

1

2

3

4
PORE VOLUME-T

5

6

7

Fig. 3—Calculated effluent curves for a sorbing medium, as influenced by the fraction of the adsorption sites in the dynamic region (f).

478

SOIL SCI. SOC. AM. J., VOL. 40,

1976

q B 10cm /day0=0.40

1.30

g/cm'

0,65
f= 0.40D= 30

3cm^/day

4
PORE VOLUME-T

Fig. 4—Calculated effluent curves for a sorbing medium, as influenced by the mass transfer coefficient (a).

Figure 4 shows the influence of the mass transfer coefficient a (day"1). When a = 0, there is no diffusion of material into the immobile liquid, and the solution for the mobile
concentration (cm) reduces toG(x,T)
cm(x,T) =[45]

peak concentration and considerable tailing. With increasing values of the mass transfer coefficient, the rate of exchange between the two regions increases, eventually leading to a new "equilibrium," where the concentrations inboth mobile and immobile liquid are identical: cm = cim.
The linear adsorption model (Eq .[15] and [ 1 6]) then reduces
to a single equation (a—»oo);. G(x,T - T,)with G(x, T) given by Eq. [41]. Note that the effluent curveacquires a more or less symmetrical shape. When a issmall, a slow exchange of material between mobile and immobile regions takes place, causing a significant decrease in

1

2

3

(0 + p K } . = 9mD_ Vm9m[46]

The curve «-»<» represents the solution of Eq. [46] for identical initial and boundary conditions as used before.

4

5

6

7

PORE VOLUME-T

Fig. 5—Calculated effluent curves for a sorbing medium, as influenced by the dispersion coefficient (D).

VAN GENUCHTEN & WIERENGA: MASS TRANSFER STUDIES IN SORBING POROUS MEDIA: I.

INFLUENCE

Kd,
°dS

x'
E
V

'

x'
/

2

! .e s
§
d
K

O.O

q

Cm/d y"°°fi=0.40_^-————\^^—-"""^\^s^^p= 1.30 g/cm3^^_JL^^^^^^

D« 3O cm 2 /doy\\*° °*65
f=0.40^^^^*^\^-—*^"^*^~ VV^s^*^ ^^^*~^^-~~~~~~^~\~^t*^^/^/''^ sS^~^^*~^\\/ o/ /XX^50/ / °;
/ /
7S
III
°* '

479

\\\

___

0=0.15v3 I/day\\AV~~~X\ X\ \\\ \ \ \ \2>0\•: ijjjj y ^^^I / / / // / / /
I / / / /

O

I

/
-°
/
2

2

\ \ a50\ °-75'\\ \ 25 \ \\ °- \ \ \

3

4

5

\

\

\

6

7

PORE VOLUME-T

Fig. 6—Calculated effluent curves for a sorbing medium, as influenced by the adsorption constant (K,Mk).
The influence of the dispersion coefficient D (cm2/day)on the effluent curves is shown in Fig. 5. The influence of D
on the shape of the calculated curves is not very large. Note.,he discontinuities,.
..
»u 11.3
1 and
j /i4.3
T pore volumesiefor^-theat. uboth
limiting case when D becomes zero. From Fig. 5 we may
conclude that for the purpose of predicting effluent concentration distributions an estimated value of the dispersion coefficient may often be adequate.
Figure 6 finally shows the influence of the adsorption
constant K on the shape and position of the effluent concentration distributions. The curve K = 0 represents the solution for a nonsorbing medium, given by Eq. [33]—[40] with
R = 1. Note that the curves are highly asymmetrical, evenfor the nonsorbing medium. When K increases, the chemical appears later in the effluent, while the peak concentration decreases somewhat. Also the tailing becomes more
, ...
.
.
, ,, * . .
,pronounced with increasing values of K. This is to be expected since more and more solute is allowed to diffuse intoand be adsorbed by the stagnant region of the soil. After
passage of the solute pulse the chemical slowly diffusesback into the mobile liquid, resulting in extensive tailing.
In conclusion, the calculated effluent curves in Fig. 2-5
clearly show the very significant effects of adsorption, immobile water and lateral diffusion on the position and shapeof effluent concentration distributions. In two following
papers (van Genuchten and Wierenga, 1976; van Genuchten
et al., 1976) the solutions presented here will be used todescribe observed effluent data from soil columns leachedwith tritium (3H2O) and 2,4,5-T (2,4,5-Trichlorophenoxyacetic acid), respectively.
APPENDIX—NOTATION
C0= influent concentration (/u.g/cm3)
Cj= concentration
breakthrough side effluent curve
Oig/cm3)
Cm,Cjm = concentration 3of mobile and immobile liquid, respectively (/ig/cm )cm,cim = relative concentrations (Cm/C0;dJC0, respectively)= dispersion coefficient (cm2/day)/ ,= fraction adsorption sites in dynamic region='o>'i
Bessel functions£ fr = modified
CQnstants jn Freundlich isotherm: s = KC»L— iengtn Of column (cm)P
= Peclet number of column: P = vmL/Dq= volumetric flow velocity (cm/day)R= average retardation factor: K = 1 + pK/0£"
= retardation factor of dynamic region:/?m = 1 +fpK/Om=*'»^f0" factor °f stagnantreglon: Rim = l + (l ~5m 5jm = adsorpti;n in dynamic and stagnant region> respec.
tively Og/g)t= time (days)h= pulse period (days)T= P°re volume: T = v0t/L = vmt4>IL
TI= dimensionless Pulse Pflod: ^ = V^>L
v0
= average pore-water velocity (cm/day): v0 = q/9Vm= average pore.water velocity in dynamic region(cm/day): vm = q/0m
x
= dimensionless distance: x = z/Lz= distance (cm)«= mass transfer coefficient (day-1)
a= dmiensionless mass transfer coefficient: a = aLlq^
= fra(£on mobile water. ^ = 9J0= bulk density (g/cm3)p
9= water content: 0 = 6m + 0im(cm3/cm3)
Om= mobile water content (cnrYcm3)0^= immobile water content (cm3/cm3)
D

480

SOIL sci. soc. AM. J., VOL. 40, 1976