SOIL SCIENCE SOCIETY OF AMERICAPROCEEDINGS
VOL. 39

NOVEMBER-DECEMBER 1975

No. 6

DIVISION S-l—SOIL PHYSICS
Predicting Anion Movement in Disturbed and Undisturbed Soils1

D. K. CASSEL, M. TH. VAN GENUCHTEN, AND P. J. WiERENGA2
ABSTRACT
Prediction of anion displacement in soils has proven difficult
owing to anion exclusion and the occurrence of anion sources orsinks within the soil system. In addition, variations in bulk density and texture with depth cause mean pore water velocity
and soil water content to be functions of depth. In this paper a
numerical solution of a modified convective diffusion equation,
which accounts for apparent anion exclusion and variations insoil water content and mean pore water velocity as functions ofdepth, is discussed. Moreover, a zero-order rate term is incorporated to account for nitrification occurring in the upper portion of the soil profile.
The model was tested using previously published NO3 and36
C1 displacement data for disturbed and undisturbed cores ofAberdeen loam. The model predicted 36C1 displacement well.
Incorporation of the source term improved prediction of thenitrate breakthrough curve, but agreement of the predictedcurve with observed data was only fair. When average values
for water content and bulk density were used throughout theentire column, instead of depth dependent values, the description of anion movement was not as good. It is suggested, however, that the improvement in accuracy of prediction of anion
displacement may be less than the variation in solute movement
observed among replicate soil cores.
Additional Index Words: miscible displacement, anion exclusion, disturbed soils.

ANY MODELS have been proposed for predicting or deM
scribing solute movement in soils. Nielsen and Biggar(1962) discussed several analytical models for predicting
the concentration of solute in the effluent when water containing solute displaces water in a porous medium initially
1
Contribution from the Department of Soils, North Dakota
State University, Fargo, N.D. 58102. Journal Paper 566. Received 20 May 1974. Accepted 28 July 1975.
2
Associate Professor, North Carolina State Univ. and former
Graduate Research Assistant and Associate Professor, New

Mexico State University, respectively. The senior author is presently with the Department of Soil Science, North Carolina State
Univ., Raleigh, N. C. 27607.

1015

free of solute. They included the diffusion equation and theconvective-diffusion equation. Experimental results showthat these models seldom, if ever, predict the observed
solute behavior accurately, partly because the soils are assumed to be nonreactive porous media. Most agriculturallyimportant soils, however, react with the solute in a variety
of ways.
Adsorption and exchange processes associated withcations and some pesticides tend to delay the transport ofthese solutes through soils [Biggar and Nielsen (1963), Lai
and Jurinak (1971)]. Furthermore, some solute may begained or lost from the soil system due to microbiologicalor chemical processes [Mansell et al. (1968), Corey et al.
(1967), Erh et al. (1967), Misra et al. (1974), Starr et al.
(1974) ]. For negatively charged ions, anion exclusion
phenomena may cause the solute to appear in the effluent
earlier than predicted [Thomas and Swoboda (1970) ].
Solutes such as pesticides and fertilizers which are applied
to soils react with the soil in one or more of the manners
described above.
Corey et al. (1963) used an analytical solution of the
diffusion equation to describe pulse or slug flow of 36C1
and tritiated water through sandstone. They observed that
tritium appeared in the effluent later than predicted while
the reverse was true for 36C1. The delayed appearance oftritium in the effluent was attributed to the surface activityof the sandstone, while the 36C1 behavior was associatedwith anion exclusion. Experiments of this nature indicatethat the models commonly used, whether solved analyticallyor by numerical techniques, rarely predict anion movement
adequately. This is true even for simple homogeneous systems, e.g., Krupp et al. (1972); Thomas and Swoboda(1970); Bresler (1973). One of the more important limitations imposed by the present models is the requirementthat the porous medium or soil be uniform. This requirement leads to further assumptions that bulk density, initial
microbial populations, etc., are also uniform with depth.
Soils very rarely exist in this manner.
The purposes of this study were three-fold: The first was

1016

SOIL SCI. SOC. AMER. PROC., VOL. 39, 1975

20

4O

60

MUTER-CONTENT PORE WATERVELOCITY

Fig. 1—Variation in (a) bulk density, Db; (b) soil water content, 6; (c) average pore water velocity ,v0 and (d) theanion source or sink coefficient, a with depth in a hypotheticalsoil profile.
to study the effects on calculated breakthrough curves ofapproximating depth-dependent values of water content,pore water velocity, and either anion source or sink termsby average values. A second purpose was to incorporate
apparent anion exclusion into the convective-diffusion
equation to enable a more accurate description of Cl andNO3 displacement through soils. Apparent anion exclusion
as used in this paper makes no distinction between (i) exclusion of anions from water associated with clay surfaces(double-layer phenomenon) and, (ii) an "apparent" exlcusion derived from incomplete mixing of solute with deadwater inside aggregates, i.e., caused by the slow process of
intra-aggregate diffusion. Hence we assume that the exclusion process is related to negatively charged clay surfaces.
The third purpose of this paper was to provide a workablefield model, which can be used for a natural soil bodyand which can be applied with a minimal amount of soils
information. This position is opposite to that of the more
complete models of Misra et al. (1974) or Starr et al.
(1974), whereby rate constants for all constituent reactions
entering into the nitrogen transformations must be measured or estimated before they can be used.

Solving Eq. [1] by numerical methods allows several improvements over the analytical solution. In order to clarify these differences consider the generalized soil property data plotted asa function of soil depth for a nonuniform soil (Fig. 1). Rarely
does bulk density (Db) remain constant with depth, even for
artificially, "uniformly packed" soil systems. Variations in bulk
density in the Ap soil horizon occur throughout the year as aresult of tillage practices. In addition, £>b varies with depth ina given profile due to variations in texture, structure, and horizonation. Figure la shows one possible variation of Db withdepth. For a tilled, agricultural soil, the greater density around
20 cm might result from a plow pan. As a consequence of thisDb-depth variation, the water content 6 is also expected to varywith depth (Fig. Ib). Other possible causes of variation of 9 withdepth may be attributed to changes in soil texture or structure
as they affect the soil water characteristics. The presence ofzones or horizons that restrict water movement may also giverise to temporary changes in soil water content with depth. In
general, the increase or decrease of 9 with depth will not bemonotonic. The broken vertical line in Fig. Ib represents themean soil water content (0) for the entire soil profile; also
shown in Fig. Ib is 0ex, the apparent anion exclusion volume orthat volume fraction of water which apparently repels anions.
9ex may vary with depth due to changes in soil texture, mineralogy, and organic matter content [Bower and Goertzen (1955)].
If we assume a steady-state soil water flux q (cm/day) throughthe soil, the average pore water velocity (v0, Fig. Ic) varies withdepth in a manner inversely proportional to the soil water content 9 at that depth and is given by[3]

If 0ex can be evaluated, the apparent pore water velocity foranion transport v'0(x) is given by= v0(x)6(x)/e'(x)[4]where 6'(x) is defined asl. [5]
THEORETICAL
The one dimensional convective-diffusion equation at variable
water content may be written as:d(6v0Q
"~~dx

dt

[1]

where 9 is the volumetric water content (cmVcm3), C the soluteconcentration in the soil (meq/cm3 soil solution) at soil depthx (cm) and time t (days), D the dispersion coefficient (cmVday)and v0 the pore water velocity (cm/day). Boundary and initialconditions for solving Eq. [1] for a pulse of solute entering thesurface of a semi-infinite profile are:[2a]C(x, 0) = 0

lim

, dec
dx

v0ec).' = Voec0=0o<t<t0f>fo

In general, /?(*) varies with depth and must be measured experimentally.
Variations of 9' and v'0 in a hypothetical soil are shown in
Fig. Ib and c, respectively. Physically, 9'(x) is that part of thesoil water through which the solute is transported. One would
expect a priori that the dispersion coefficient would be larger
when the exclusion phenomenon is taken into account because
the same spreading of the solute must be obtained in a smaller
volume of water and in a shorter time period. Therefore, we
may write:

D' = D/p[5a]V'0 = V0/^[5b]f = pt[5c]o' = pe[5d]C' = C/p

[5e]

[2b]

where
C0 is the concentration of the incoming solution (meq/cm3) and ta is the pulse period.
Under conditions of constant water content (8 = ~6), 8 can beeliminated from Eq. [1] and [2], v0 assumed constant, and theequations solved analytically [Lindstrom et al. (1967)].

Many reports in the literature indicate that the dispersion coefficient D (cmVday) is a function of mean pore water velocity,
v0. We contend, furthermore, than an anion dispersion coefficientis also a function of the apparent exclusion volume 9eJi and henceof 0(x) as shown in Eq. [5]. Hence the equation £>' = £>//?(*)
describes the apparent dispersion coefficient D'. In reality, D'CASSEL ET AL.: PREDICTING ANION MOVEMENT IN SOILS

is a function of x since ft and D are both functions of x; D is afunction of depth because v0 and 9 are depth-dependent.
Using these relations, Eq. [1] in terms of 9', v'0, D', and C'
and /' becomes now, including a zero-order reaction term a(meq/cm3/day):

1017

ABERDEEN(DISTURBED)O

O 0.25- NUMERICAL

D'- 5.6 cmVdo,(D

X - 8 7 cm

tO 0.20t' 0.815

V Hi)

O

[6]

0.15

UJ

with the following initial and boundary conditions for a semiinfinite profile:|_ 0.10[7a]C(x, 0) = 0

UJ 0.05cc

0.5

x-H>0->
0<t'<t'0= 0

f > t'0

1.0

PORE
[7b]

The a term in Eq. [6] may represent a source (+) or a sink(—) for the anion under consideration; a may, or may not, be afunction of depth in a given soil profile. In this study, a wasassumed to be zero for Cl-movement and positive for NO 3 ~,varying with depth in a manner similar to that shown in Fig. 1 d.
In soils experiencing net denitrification, a is negative. Ideally,an independent measurement of a as a function of soil depth isdesirable.
In the desire to obtain a model which has some degree of
utility an attempt was made to keep that portion of the numerical model dealing with the source term « as simple as possible.
Experimental data for displacement of nitrate and chloride
through the soils involved were analyzed and a value of aobtained by curve fitting. A zero-order term was found to obtain an acceptable fit. Zero-order reactions are applicable totransformations when substrate is plentiful; due to the 5 to 6
weight percent organic matter in the A horizon, there apparently
was sufficient substrate for nitrification of organic matter initially present. It is not the purpose of this paper to dwell on theproper order of rate reaction: this will depend upon the soil inquestion, its water content distribution, the soil water flux, andmany other factors. Kirda et al. (1974) concluded that thesetransformations in a field situation may fall between zero- andfirst-order. McLaren (1969a, b) and Starr et al. (1974) treat thissubject in detail.
Equation [6] subject to conditions [7] for a semi-infinite profile was solved numerically using S/360 CSMP computer language. The program used is a modification of the model discussed extensively by van Genuchten and Wierenga (1974,
example 3). The modification included proper formulation of
the upper boundary condition, inclusion of the anion exclusion
phenomena as discussed above, and allowing v, a, and e to varywith depth.
EXPERIMENTAL
Aberdeen loam, a soil belonging to the fine, montmorillonite
family of the Glossic Udic Natriborolls, was used. The 22-cm
deep Ap horizon is underlain by a clay loam B2 having strong
medium prismatic structure. Undisturbed Aberdeen soil cores
30 cm in diameter and 80 cm long were collected from the field(Cassel et al., 1974). Disturbed Aberdeen soil cores, 87 cm long,
were packed horizon by horizon to reconstruct the original
morphological horizonation. After prewetting and preleaching
the core, a 2.54 cm pulse of water containing 71.4 meq/liter
NO3-N and tagged with 36C1 was applied to the soil surface. The
relative concentrations of NO3 and 36C1 in the effluent wereanalyzed and plotted versus pore volumes of column effluent.
The number of pore volumes (T) displaced through the soil was
calculated by dividing the volume of effluent by P0. The volumetric water capacity is defined by

1.5

2.0

VOLUMES EFFLUENT

Fig. 2—Relative3636C1 concentration versus pore volume for displacement of C1 through disturbed Aberdeen loam. Eachdatum represents the mean of three cores; the vertical lines
represent the range. The curve is the numericalsolution of
Eq. [6] with /J = 0.815 and D' — 5.6 cm2/day.

•=r

BAdx

[8]

where A is the crosssectional area of the soil core and L thecolumn height.
Apparent anion exclusion volumes for the 0-22, 22-56, and56-87 cm depth intervals of Aberdeen loam were measured asfollows: Columns of uniformly packed soil 15 cm long and 7.6
cm_diameter were wetted and drained to a mean water content
of J = 0.383 cm3/cm3. 36C1 solution, introduced into the soil
through fritted glass end plates, was used to displace the soil
solution from the column. 36C1 concentration in the effluent
was plotted versus T. The volume of soil water which apparently
excluded anions, Vex (cm3) was calculated by

V
~

QV

=P
*

n

[9]

where P(cm3) is the total volume of effluent collected at the
moment when the relative 36C1 concentration in the effluent isequal to 0.50. The apparent anion exclusion volume 0ex is thencalculated by[10]
RESULTS AND DISCUSSION
Disturbed Aberdeen Loam
Figure 2 presents data for the displacement of a pulse of
C1 through disturbed Aberdeen soil (Cassel et al., 1974).
Three separate disturbed soil cores were used for this experiment. The vertical lines indicate the range of 36C1
measured for the three replicates. The seemingly large range
in relative 36C1 concentrations may be explained by thevariation in the total amount of water held in each core:P0 ranged from 19.8 to 23.3 liters even though gamma radiation attenuation data showed that the three cores had nearly
identical bulk density distributions (Fig. 3). In Fig. 2, theabcissa refers to the mean pore volume for the three cores.
Chloride appeared earlier in the effluent of the core having
the smallest pore volume since the soil water flux was identical for all cores, being 0.383 cm/day.
The solid curve in Fig. 2 was obtained by solution of Eq.
[6] with /? = 0.815, a = 0, and an apparent dispersion co36

1018

SOIL SCI. SOC. AMER. PROC., VOL. 39, 1975
BULK

1.0

DENSITY

1.20

1.40

(gin/cm9)

0.20
DISTURBED
0.0007 meq /cm3/
week
D'=5j6cm*/day

1.60

V f(x)
9= f(x)-I20

Db- 1.47 g/cm3
6 • 0.346

0

140

1.0

2.0

PORE VOLUMES EFFLUENT

a.

Fig. 5—Mean and range of relative concentration versus porevolume for displacement of nitrate through disturbed Aberdeen loam. The curve is thenumerical solution of Eq. [6]using « = 0.0001 meq/cm3 per day.
O 60

0.20

§
o
80
O.K)
0.30

WATER CONTENT

0.40

3

(cm'/cm )Fig. 3—Mean and range of bulk density and soil water content
distributions for disturbed cores of Aberdeen loam.
0.25

o

0.15

0:05

UJ
(E

0.5

1.0
VOLUMES

1.5

VOLUMES

2.0
EFFLUENT

Fig. 6—Mean36relative concentration data for displacement of
nitrate and C1 through undisturbed cores of Aberdeen loam.
The curves are numericalsolutions of Eq. [6] with /} =0.815, D' = 7.6 cm2/day.

o.io-

PORE
36

1.0
PORE

a-0.815, D'"5.6cm/day
e-f(x),
v 0 -f(x)J-0.815, D''5.6cm2/day
9=0.383, va> 0383 cm/day
6" IJOO, D -4.5 cm'/doy
9-0.383, v «0.383 cm/day. 0.20bl

bl

0
ABERDEEN (DISTURBED)

o

-

O.5O

2.0

EFFLUENT

Fig. 4—Relative C1 concentration versus pore volume for
disturbed Aberdeen loam. The solid line is the numericalsolution of Eq. [6]. Anion exclusion was ignored for thedashed curve. The dashed-dotted curve is the analytical solution of Eq. [6], assuming a constant soil water content.
efficient D' equal to 5.6 cm2/day. Values of 6 used in thecomputation varied with depth as shown in Fig. 3. The
apparent anion exclusion volume 0ex was measured independently of the above experiment using disturbed, uniformly packed 15-cm long soil columns as described above.
/? was also calculated using this information and was essentially identical for the three layers; hence a value of 0.815
was used for the entire column. The curve in Fig. 2 describes the displacement of 36C1 adequately, considering thevariability in soil water content and bulk density between
the three replicates.
The effect of ignoring anion exclusion and variation ofwater content with depth is shown in Fig. 4. Curve 1 is identical to the curve shown in Fig. 2. Curve 2 was obtainedusing an average water content 9 of 0.383. Curve 3 represents results obtained when ignoring the anion exclusion.
Comparison of Curve 1 and Curve 2 indicates little advantage in using the actual water content distribution as opposed to the mean water content for this soil under the conditions of this experiment. It should be mentioned thatcurves 2 and 3 could have been obtained with an analyticalsolution as mentioned in the theoretical part of this paper.
Curves 2 and 3 could also have been obtained with ananalytical solution of Eq. [6] [Lindstrom et al. (1967),
and van der Ploeg and Benecke, (1974)].
Data similar to that for 36C1 in Fig. 2 are shown for NO3
in Fig. 5. Once again considerable variation between thethree replicate disturbed Abderdeen soil cores was observed. The solid line represents the numerical solution withthe following parameters: a = 0.0007 meq NO3/cm3 per
week, D' = 5.6 cm2/day, v0 = f(x). The source term wasincluded to account for a net increase in NO3 resulting frommicrobial nitrification: a was obtained by curve fitting andwas found to have a value of approximately 10~4 meq/
cm3 per day in the A horizon or the upper 22 cm of thesoil profile; a was assumed to be zero below 22 cm, a choicebased partially upon Dancer and Peterson's (1974) conclusion that nitrification is largely confined to the Ap horizon. Moreover, organic matter content decreases rapidly
at this depth. Although the calculated curve does not accurately describe the observed nitrate concentration distribu-

CASSEL ET AL.: PREDICTING ANION MOVEMENT IN SOILS
tion in the effluent, it better approximates the data than forthe case where a — 0 throughout the profile (see Fig. 4).
A more exact manner of accounting for nitrification within
the soil profile in a steady-state flow system has been proposed by McLaren (1969a, 1969b). McLaren's model
treats each step in the nitrification process and employs
specific rate constants for each reaction. Starr et al. (1974)
have verified experimentally a portion of McLaren's nitrogen transformation model using uniformly packed columnsof Yolo loam. In order to use this model, however, appropriate rate constants for each reaction must be measured
for each soil and possibly each horizon in a given soil.
Hence such a model is limited in use at the field level
whereas the simpler approach used in this paper might serve
to provide a realistic estimate of the nitrate produced.
Undisturbed Aberdeen Loam
Finally, Fig. 6 compares observed and calculated curvesof NO3 and 36C1 from undisturbed cores of Aberdeen
loam. The 36C1 effluent curves are described well; prediction of the effluent is not good but is slightly better thanthe prediction of NO3 for the disturbed Aberdeen loamcores (Fig. 5). Notice that a is greater for the undisturbed
than for the disturbed Aberdeen cores indicating a greater
net nitrification than for the disturbed soil. This behavior is further evidence that the use of disturbed soils inexperiments can lead to results very dissimilar from identical experiments performed using natural soil cores. In addition, the apparent dispersion coefficient is greater for undisturbed Abderdeen being 7.6 cm2/day compared with 5.6
cm2/day for disturbed Aberdeen.
Throughout this paper we have dealt with an apparent
anion exclusion volume which varied with depth. No
attempt was made to separate the "double layer" type anion
exclusion from exclusion arising from incomplete intraaggregate mixing, sometimes referred to as "dead end"
pores. Saxena et al. (1974) have recently questioned theconcept of "dead pores."
Immediately after adding a 2.5-cm pulse of water to thesoil surface, the water moves downward through the soil
and flows at a slower rate as time and depth increase, i.e.,
the process is transient. Nearly all water will have passed
the 5-cm depth, for example, after 1 hour, but it takes 6 to
7 days for this same quantity of water to flow past the lower
end of the core. Hence the deeper in the soil core, the moretime there is available for diffusion of anions out of deadend pores to occur. In the upper portion less time for diffusion is available as water moves relatively rapidly. For thisreason the ratio of "double layer" to "dead end" pore typeanion exclusion, if indeed, both types are manifested, willchange with depth. Therefore, separation of these two processes in a model designed for field application in nonuniform soils is academic. We conclude that a useful anion
displacement model, capable of working in the field situation, needs only an apparent anion exclusion volume orcoefficient.

1019

ERRATA
Predicting Anion Movement in Disturbed andUndisturbed Soils

90' C' 9 Ioir.,dC'\ d , . , , „ , , , ,~^~ ' K \6 ° "to ) ~ T* (° V°C ] + a

D. K. CASSEL, M. TH. VAN GENUCHTEN AND P. J. WIERENGASOIL Sci. Soc. AMER. PROC. 39(6): 1015- 1019. 1975.

with the following initial and boundary conditions for a semiinfinite profile:

Several confusing errors occurred in publication of the above
article in the November-December 1975 issue of the SSSA
Proceedings.
I) The first paragraph of the "Theoretical" section shouldread:
The one dimensional convective-diffusion equation at variable
water content may be written as

C(x,Q) = 0

dec _ a Lsc\

d(0Voc)~ ~ ~ ~ - ~ -rL1J,

where 0 is the volumetric water content (cnvVcm3), C the soluteconcentration (meq/cm3 soil solution) at soil depth x (cm) and timet (days), D the dispersion coefficient (cnvVday), and v,, the porewater velocity (cm/day). Boundary and initial conditions forsolving Eq. [1] for a pulse of solute entering the surface of a semiinfinite profile areCU.O) = 0

x >O

=0

[2a]

f&t.

[2b]

where C0 is the concentration of the incoming solution (meq/cm3
solution) and /„ is the pulse period.
2) The first paragraph of page 1017 should read:Using these relations, Eq. [1] in terms of d' , v0', D' , C' , and /'now becomes, i n c l u d i n g a zero-order reaction term a'(meq/cm3/day) equal to a//3 where P is constant with depth:

x>0

[7a]

\im(-6'D'——+v'08'c}=v'00'C'0Q^t'<t'0*-»° \
Sx
I
=

"

r,n
[6]

' — ' <>

[7b]