Published online November 20, 2007

Analytical Advection–Dispersion Model forTransport and Plant Uptake of Contaminants in
the Root Zone
ORIGINAL RESEARCH

T. H. Skaggs,* N. J. Jarvis, E. M. Pontedeiro, M. Th. van Genuchten, and R. M. Cotta

In regulatory and risk management analyses of environmental contaminants, the vadose zone may be treated as a subcomponent
within a larger environmental modeling framework. For the complexity of the larger system model to remain at manageablelevels, it is desirable that subcomponent models be relatively simple and require few input parameters. In this work, we develop
an advective–dispersive solute transport equation that includes plant uptake of water and solute and present an analytical solution. Assumptions underlying the transport model include linear solute sorption, first-order uptake, and a uniform soil watercontent. We examine the latter assumption in detail and demonstrate the effects of rooting depth, soil texture, and leaching
fraction on the uniformity of the root-zone water content. The new analytical advection–dispersion model should be useful for
estimating the transport and uptake of strongly sorbing and persistent contaminants, where the timescale relevant for assessing
environmental impacts is long (decades) and short-term fluctuations caused by, for example, precipitation can be averaged. As
an illustration, model predictions are made for the uptake of cadmium (Cd) by wheat (Triticum aestivum L.) grown in sludgeamended soil. The predictions are compared with those of a “one-compartment” model that has been proposed previously forrisk analysis and regulatory studies. The comparison shows that the one-compartment model overestimates the long-term,
steady-state Cd concentration in harvested wheat grain. The analytical advection–dispersion model is recommended as a toolfor environmental risk assessment of strongly sorbing, persistent contaminants.

ABBREVIATIONS: ADE, advection–dispersion equation; GITT, generalized integral transform technique.

T

he movement of solutes in the root zone involves complex
interactions of physical, chemical, and biological processes.
Understanding and modeling these processes, as well as theirimpacts on plant growth and water and soil quality, remains a challenging area of research for scientists and engineers. Numericalcomputer simulation models such as HYDRUS (Šimůnek etal., 2005), UNSATCHEM (Šimůnek et al., 1996), MACRO(Larsbo et al., 2005), and RZWQM (Ahuja et al., 2000) permit comprehensive analyses of root zone processes by allowing
for a wide range of soil properties, chemical reactions, and root
functions. Although these models provide great flexibility, theircomplexity has a downside. Using the models requires specifying values for a large number of parameters, many of which areoften unknown or unobtainable. Long-term simulations, such as
T.H. Skaggs and M.Th. van Genuchten, USDA-ARS, U.S. Salinity Lab., 450 W. Big Springs Rd., Riverside, CA 92507; N.J. Jarvis,
Swedish Univ. of Agricultural Sciences, Uppsala, Sweden; E.M.
Pontedeiro, Brazilian Nuclear Energy Commission, Rio de Janeiro,
Brazil; R.M. Cotta, Federal Univ. of Rio de Janeiro, Rio de Janeiro,
Brazil. Received 5 July 2007. *Corresponding author (tskaggs@ussl.
ars.usda.gov).
Vadose Zone J. 6:890–898
doi:10.2136/vzj2007.0124© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmitted in any form or by any means, electronic or mechanical,
including photocopying, recording, or any information storage andretrieval system, without permission in writing from the publisher.
www.vadosezonejournal.org · Vol. 6, No. 4, November 2007

890

may be required when evaluating the transport and plant uptakeof strongly sorbing (i.e., slow-moving) contaminants such as
trace metals and radionuclides, are difficult because of the needto specify time-varying boundary conditions over many years,
and because of lengthy execution times. Also, the complexity ofnumerical codes is such that many of the details of the models
are often not apparent to users other than the model developers,
thus making it difficult for both technical experts and regulatory
agencies to interpret simulation results.
Alternatively, solute transport and plant uptake in the rootzone can be studied using less-complex models, which typicallytake the form of a partial differential equation that can be solved
analytically. While analytical models have the advantage of being
relatively transparent with respect to model inputs and outputs,
this simplicity often comes at the expense of significant assumptions and approximations, such as requiring uniform soil conditions and/or representing nonlinear chemical reactions withlinear models. In some instances, however, these simplifications
may not be too restrictive. For example, in long-term simulations(i.e., over tens or hundreds of years), short-term boundary fluctuations caused by precipitation can be averaged and an analyticalmodel with time-averaged boundary conditions and root-zone
average soil and transport properties may closely approximate
a more complex, transient model (e.g., Destouni, 1991). Suchan approximation may be poor in the case of quickly degrading solutes (e.g., pesticides) but should work well for strongly
sorbing and persistent contaminants (e.g., trace metals, radionuclides) when the relevant timescale for assessing environmental∂( θC + ρ b S ) ∂j s+
+ rs = 0∂t∂zimpacts is measured in decades or centuries rather than seasons
and years.
Among analytically solvable solute transport models, abundant literature exists on the classical advection–dispersion equation (ADE). While analytical solutions to the ADE are widely
available (e.g., van Genuchten and Alves, 1982; Skaggs andLeij, 2002), most do not explicitly incorporate water and solute uptake by roots, thus limiting their usefulness in analyzing
root zone processes. Raats (1975) used the method of characteristics to obtain an analytical transport model that included root
water uptake but neglected diffusion and dispersion processes.
Generalizations and extensions of the Raats (1975) model havebeen provided by Ginn and Murphy (1997) and Schoups andHopmans (2002). The latter study included plant uptake of bothwater and solute. Hoffman and van Genuchten (1983) reviewed
various steady-state transport solutions that include root water
uptake while neglecting diffusion and dispersion, whereas Raats(1981) considered the effect of a constant solute diffusion coefficient on the shape of steady solute profiles. To date, none of
the analytical modeling efforts have included a general (i.e., possibly nonconstant) description of diffusion and dispersion. In
this work, we develop a solute transport model that accounts foradvection and dispersion diffusion, as well as for the uptake ofwater and solutes by roots, and we present an analytical solution.
Additionally, we discuss general characteristics of the proposedtransport model and consider an example application involving
the uptake of cadmium by wheat (Triticum aestivum L.) growing
in sludge-amended soil.

[5]

where C is the solute concentration (M L−3), θ is the volumetricwater content (L3 L−3), ρb is the soil bulk density (M L−3), S is
the sorbed solute concentration (M M−1), js is the solute massflux (M L−2 T−1), and rs is a solute sink term (M L−3 T−1)
that accounts for the uptake of solute by roots. The solute flux isassumed to follow the standard advection-dispersion model,
j s = q( z )C − θD( z )

∂C
∂z

[6]

where D is the effective diffusion–dispersion coefficient(L2 T−1), which may vary with depth due to the dependence ofdispersion on solute velocity. In this study, we assume that solute
uptake is a first-order process such thatrs = γrwC[7]

where γ is the uptake coefficient (γ ≥ 0). A mechanistic interpretation of γ in terms of transport across root membranes isprovided by Dalton et al. (1975); additional discussion of thefirst-order uptake model in given below in “First-Order Plant
Uptake.” We also assume that the soil water content is approximately uniform with depth [θ(z) ≈ θ] and that sorption is linear with a partitioning coefficient Kd. With these assumptions,introducing Eq. [2], [3], [6], and [7] into Eq. [5] leads to thefollowing transport equation:
R

Theory

∂C
∂ ⎡⎢
∂C ⎤⎥
∂C
=
−[v 0 − v T B( z )]D( z )
⎢
⎥
∂t
∂z ⎣
∂z ⎦
∂z

[8]

+[(1 − γ )v T b( z )]C

Steady Water Flow with Plant Uptake

We are interested only in the case of q0 > T, such that water flow
is vertically downward (q > 0).

where R is the retardation coefficient, v0 ≡ q0/θ, and vT ≡ T/θ.
Note that because q0 > T (see “Steady Water Flow with Plant
Uptake”), it follows that v0 > vT. Formally, R is defined by R =1 + Kdρb/θ, although in practice R is commonly treated as aneffective, empirical parameter. The assumption that the watercontent is uniform with depth is important to obtaining Eq. [8];we examine that approximation in detail in the next section.
The left-hand side of Eq. [8] and the first term on the righthand side are familiar terms common to many advection–dispersion transport models (e.g., Skaggs and Leij, 2002). The second
term on the right-hand side is an advection term that contains a
depth-dependent transport velocity, which decreases with depth
from a maximum value of v0 at the soil surface to a minimum
value of v0 − vT at the bottom of the root zone. The final termin Eq. [8] is a combined sink–source term that accounts for theconcentrating effect of a decreasing water flux and the countering
effect of solute removal. When γ = 0, no solute is taken up andthe concentrating effect in the root zone is maximal. When solute is taken up passively with water (0 < γ ≤ 1), the concentrating
effect is reduced, becoming nonexistent when γ = 1. Specifying
γ > 1 corresponds to active solute uptake (Schoups and Hopmans,
2002), such that solute concentrations are reduced as a result ofplant uptake.

Advective–Dispersive Solute Transport and Plant Uptake

Uniformity of Soil Water Content

Conservation of mass for one-dimensional solute transport
can be expressed as

The derivation of the solute transport equation in the previous section involved the assumption that the soil water content

Conservation of mass for one-dimensional steady waterflow can be written asdq+ rw ( z ) = 0[1]dz
where q is the water flux (L3 L−2 T−1 ), rw(z) is a sink term specifying the rate of water extraction by roots (L3 L−3 T−1 ), and z
is the vertical coordinate (L), oriented positive downward withthe soil surface located at z = 0. The sink term can be expressedas (e.g., Feddes et al., 1978; Skaggs et al., 2006)
rw ( z ) = Tb( z )[2]where T is the transpiration rate (L3 L−2 T−1) and b(z) is thenormalized uptake density (L−1). Integration of Eq. [1] gives thesteady-state water flux profileq( z ) = q 0 − TB( z )[3]where q0 ≡ q(0) and B(z) is the cumulative uptake distribution,z

B ( z ) = ∫ b ( z ′ )dz ′0

www.vadosezonejournal.org · Vol. 6, No. 4, November 2007

[4]

891

⎡
αδ ⎤
αΔh = ln ⎢1 − (1 − 1/ L F )⎥1 + αδ ⎦⎥⎣⎢is approximately uniform with depth. Previous theoretical analyses of water flow with root uptake have found that the pressure
head depth-profile is generally very uniform during steady flow
(Raats, 1974) or quasi-steady flow such as may be achieved withhigh frequency irrigation (Rawlins, 1973). Steady-state watercontent profiles computed by Yuan and Lu (2005, their Fig. 4)
similarly showed a high degree of uniformity in the root zone.
To further justify the uniform water content approximation
invoked above, we present in this section an analysis of steadyflow that extends the work of Raats (1974).
Consider the specific case of a uniform soil where the steadywater flux is specified by the Buckingham-Darcy law, in whichcase Eq. [1] becomesin which
LF ≡

[10]

where Ks (L T−1) is the saturated hydraulic conductivity and α(L−1) is the Gardner parameter that ranges approximately from
0.1 cm−1 for coarse-textured soils to 0.01 cm−1 for fine-textured soils (Raats, 1974). Using the root water uptake term ofRaats (1974), rw(z) is given byrw ( z ) = (T / δ )e−z / δ[11]Equation [11] prescribes an exponential uptake distribution,with the distribution parameter δ (L) corresponding to the depthof an equivalent uniform root system having the same uptake
and transpiration rates (Raats, 1974). Although Eq. [11] accommodates an infinite rooting depth, 99.8% of the uptake occurs
within a depth of 6δ.
For the soil surface boundary condition,q 0 = −K (0)dh(0) + K (0)dzq0 − Tq0[16]

is termed the leaching fraction (i.e., the fraction of q0 that passes
below the root zone).
Figure 1 shows a plot of αΔh versus the dimensionlessparameter αδ for various values of LF. The limit αδ → 0 corresponds to a very fine soil and/or a shallow root distribution,
whereas αδ → ∞ implies a coarse-textured soil and/or very deep
rooting (Raats, 1974). Figure 1 indicates that the difference
between the pressure head at the soil surface and at the effective
base of the root zone is smaller (i.e., the profile becomes more
uniform) when LF is large and/or αδ is small. This means thathigh leaching fractions (large LF), fine-textured soils (small α),
and shallow rooting depths (small δ) promote uniformity in thepressure head profile.
To make these generalities more quantitative, consider thatthe maximum rooting depth, RD, for common agricultural crops
ranges from about 100 to 300 cm (Borg and Grimes, 1986),
with the global-average value being about 210 cm (Canadell etal., 1996). Given that 99.8% of uptake occurs within a depthof 6δ, we may characterize shallow-, average-, and deep-rooted
crops as having δ = R D /6 ≈ 20, 35, and 50 cm, respectively.
Among noncrop plant species, comparable global-average values
are RD = 700 cm for trees (δ ≈ 120 cm), RD = 510 cm for shrubs(δ ≈ 85 cm), and RD = 260 cm for herbaceous plants (δ ≈ 45cm) (Canadell et al., 1996).
Figure 2 shows Δh plotted versus RD (= 6δ) for fine-textured (α = 0.01 cm−1) and coarse-textured (α = 0.1 cm−1) soils
and various values of LF. For α = 0.1 cm−1 (dashed lines), Δhchanges very little as RD increases from 100 (shallow rooted crop)
to 700 cm (tree), indicating that rooting depth is of lesser importance in determining pressure head uniformity in the coarse soil.
The effect of the leaching fraction is more significant, with Δhincreasing approximately from 7 to 22 cm when LF is decreased

⎤
d ⎡⎢
dh
−K
+ K ⎥ + rw ( z ) = 0
[9]
⎥
dz ⎢⎣
dz
⎦
where h is the pressure head (L) and K is the hydraulic conductivity (L T−1). The hydraulic conductivity is assumed to be of
the form (Gardner, 1958)
K ( h ) = K seαh

[15]

[12]

and the lower boundary condition at a depth z = L,dh(z = L ) = 0
dz

[13]

the solution of Eq. [9] for L → ∞ may be expressed as⎡qT ⎛αδ −z / δ ⎞⎟⎤⎥eαh( z ) = ln ⎢ 0 − ⎜⎜⎜1 −⎟⎢ K s K s ⎝ 1 + αδ
⎠⎟⎦⎥
⎣

[14]

Raats (1974, Eq. [21]) presented an equivalent result
expressed in terms of the matric flux potential instead of thepressure head. Equation [14] indicates that the pressure head
will decrease downward, from a high value at the soil surface
down to some constant value as z → ∞ (or, more practically, as
z → ≈6δ). Profile uniformity can be evaluated by calculating thedifference between the pressure head at the soil surface and at
depth, Δh ≡ h(0) − h(∞). From Eq. [14] it follows that this pressure head difference is given bywww.vadosezonejournal.org · Vol. 6, No. 4, November 2007

FIG. 1. Plot of αΔh as a function of αδ for various values of the leaching fraction, LF (α = Gardner conductivity parameter [L−1]; δ = Raats
uptake distribution parameter [L]; and Δh = change in pressure head
between top and bottom of root zone [L]).
892

variations shown in Fig. 3 correspond to the greatest- and leastuniform pressure head profiles plotted in Fig. 2. As expected,
Fig. 3 demonstrates that in addition to the aforementioned
effects of rooting depth and leaching fraction, the uniformity ofthe soil water content profile is affected by the pressure head (orwater content) at the soil surface. In evaluating Fig. 3, our experience suggests that when ΔS is less than about 0.1, the watercontent variability can be safely treated as negligible since thatlevel of variability is similar to the precision that is obtainable
with, for example, electromagnetic methods of measuring watercontents under field conditions. While it is not possible to state
precisely the general conditions for which an effective, uniform
water content approximation is acceptable, we conclude that theapproximation is reasonable whenever the leaching fraction is
sufficiently large (e.g., ≥0.25). For coarse-textured soils wheredθ/dh is large near the wet end of the retention curve, there is anadditional requirement that the soil surface be drier than where
the steep part of the curve is located (for the coarse soil depicted
in Fig. 2, h(z = 0) must be less than approximately −75 cm).
Thus far, this section has been concerned with the internalconsistency of the proposed transport model. That is, we haveaddressed the question of whether the water content can be reasonably approximated as uniform given the model assumptionsof uniform soil properties and steady flow. Once that questionis resolved, a related but separate question that may be asked iswhether a model that specifies a uniform water content can beused to assess transport in fields where the water content is clearlynonuniform (e.g., in layered soil profiles). Evidence exists (e.g.,
Wierenga, 1977; Streck and Piehler, 1998) that solute distributions modeled using an effective uniform water content are notsignificantly different from those modeled using an explicit representation of the water content variability. Thus, we conclude
that Eq. [8] may be useful in the case of nonuniform soil profiles
if an appropriate effective water content can be determined.

FIG. 2. The change in pressure head between the top and bottom of
the root zone (Δh) plotted as a function of rooting depth. Plots areshown for two soil types (α = Gardner conductivity parameter [L−1])
and various values of the leaching fraction (LF).

from 0.5 to 0.1 (across rooting depths of interest, RD > 100cm). For the fine-textured soil (α = 0.01 cm−1, solid lines), therooting depth is far more important. For example, when LF =0.3, Δh increases from 29 to 81 cm when RD is increased from100 to 700 cm. The leaching fraction is also significant for thefine-textured soil: when RD = 300 cm (deep-rooted crop), Δhincreases from 29 to 139 cm when LF is decreased from 0.5 to
0.1.
The largest calculated Δh values were approximately 180
cm for the fine-textured soil (LF = 0.1, RD = 700 cm) and 22cm for the coarse-textured soil (LF = 0.1, RD > 100 cm). These
results correspond to low leaching and deep rooting scenarios,
such as might be found in semiarid and arid climates. Even forthese “worst-case” scenarios, however, the variation in the watercontent may not be excessive, depending on the nonlinearity of
the retention function, θ(h), and the magnitude of dθ/dh. The
nonlinearity in θ(h) means that the change in water content for agiven Δh will depend on the pressure head at the soil surface, h(z= 0). For fine-textured soils where dθ/dh is typically very modestly sized for all h values, a small Δh will correspond to a smallchange in water content, whereas for coarse-textured soils, dθ/dh
may become large such that a small Δh can produce a relatively
large change in water content if h(z = 0) is in the vicinity of thesteep part of the retention curve.
The water content variation can be assessed quantitatively
using the Russo (1988) model to describe the saturation:
2/( 2+m )S( h ) = ⎢⎡e−0.5α h (1 + 0.5α h )⎥⎤⎣⎦First-Order Plant Uptake

The proposed transport model assumes that solute uptake
is a first-order process (Eq. [7]). The plant uptake of solutes is[17]

where S is effective saturation (0 ≤ S ≤ 1) and m is Mualem’s(1976) tortuosity parameter, which we assume has the value m= 0.5. The difference in saturation between the surface and thebase of the root zone is ΔS ≡ S[h(0)] − S[h(∞)] = S[h(0)] − S[h(0)− Δh]. Figure 3 shows plots of ΔS as a function of h(0). Theseplots were made using the leaching fraction and root depth combinations that produced the largest and smallest values of Δhobserved in Fig. 2 (where RD ≥ 100 cm). Thus, the water content
www.vadosezonejournal.org · Vol. 6, No. 4, November 2007

FIG. 3. Changes in soil water saturation between the top and bottomof the root zone (ΔS) plotted as a function of the soil surface pressure head for two soil types. Plots are shown for a low-leaching,
deep-rooting scenario and a high-leaching, shallow-rooting scenario.
(α = Gardner conductivity parameter; LF = leaching fraction; RD =maximum rooting depth; Δh = difference in pressure head between
the top and bottom of root zone; h = pressure head; and z = depth.)893

C ( z , 0) = C I

a highly complex process that depends on plant species, solute
species, and solution composition (e.g., the uptake of cadmium
is affected by the presence of other divalent cations) (Tinker
and Nye, 2000). Detailed, mechanistic descriptions of uptake–absorption processes are possible, but these defy a simple mathematical description. Measurements of uptake kinetics commonly
show that uptake increases linearly with solution concentration
at low concentrations and then levels off and becomes constantat high concentrations. Consistent with these observations,uptake is often modeled assuming Michaelis-Menten kinetics(Tinker and Nye, 2000):
rs ∝

C
KM +C

where C0 is the concentration of the infiltrating solution and CI
is the initial concentration of the soil water. An analytical solution may be obtained using the generalized integral transform
technique, or GITT (Ozisik and Murray, 1974; Cotta, 1993;
Liu et al., 2000). Liu et al. (2000) reported the GITT solutionof a generic transport equation of which Eq. [8] is a special case.
For the present problem, the solution given by Liu et al. (2000)
may be written asM

C ( z , t ) = C 0 + ∑ N n−1/ 2 ϕn ( z )Tn ( t )where the infinite summation required in the formal solution
has been truncated at M terms and where the eigenfunction ϕnand its norm Nn are, respectively,

where KM is the Michaelis–Menten constant. Alternatively,
uptake may be modeled using two intersecting lines (Tinker
and Nye, 2000), one representing a linear increase in uptake
below some critical or threshold concentration, and one representing a constant maximum uptake rate above the criticalconcentration (Fig. 4).
When the concentration remains below the threshold concentration, the first-order uptake model is consistent with thetwo-line model, with the uptake coefficient γ being proportional
to the slope of the first line segment (Fig. 4). A connection with
the Michaelis–Menten representation also exists at low concentrations (KM >> C), where the Michaelis–Menten model may
be approximated as rs ∝ C/KM. In this case, the uptake coefficient is γ ∝ 1/KM. Note that if γ is related to KM in this way,then uptake calculated for a given concentration with Eq. [7]will always be greater than or equal to that specified by the corresponding Michaelis–Menten model (Fig. 4).

[23]

⎞
1⎛
v D
N n = ⎜⎜⎜ L + 2 02 0 2 ⎟⎟⎟2 ⎝⎜βn D0 + v 0 ⎠⎟[24]v 0 cos( βn L ) − D0 βn sin( βn L ) = 0[25]

Because of the specific form of Eq. [8] and the boundary conditions imposed, it is possible to simplify the expression presented
by Liu et al. (2000) for Tn(t). Letting T(t) be a M-length vector with nth element Tn(t), values for Tn(t) in Eq. [22] can beobtained by computing (cf. Liu et al., 2000, Eq. [13])T( t ) = exp(−A −1Bt )[ T(0) − B−1G ] + B−1G

[26]
where elements of the M × M matrices A and B and the Mlength vectors G and T(0) are given by

We define D0 ≡ D(z = 0) and consider the solution of Eq.
[8] subject to the following boundary and initial conditions∂C(L, t ) = 0∂zϕn ( z ) = cos[ βn ( L − z )]The eigenvalues βn are the first M zeros of the equation

Analytical Solution

z =0

[22]

n =1

[18]

⎛
∂C ⎞⎟
⎟
v 0C 0 = ⎜⎜v 0C − D0⎜⎝∂z ⎠⎟⎟[21]

L

Anr = N n−1/ 2 N r−1/ 2 ∫ R ϕn ( z )ϕr ( z )dz

[27]

0

[19]

⎧
⎪ L
dϕr ( z )Bnr = N n−1/ 2 N r−1/ 2 ⎪dz⎨∫0 [ v 0 − v T B( z )]ϕn ( z )
⎪
dz
⎪
⎩
L
dϕ ( z ) dϕn ( z )+ ∫ D( z ) r
dz
0
dz
dz

[20]

[28]

}

L

−∫ (1 − γ )v T b( z )ϕn ( z )ϕr ( z )dz + v 0 ϕr (0)ϕn (0)0L

Gn = C 0 N n−1/ 2 ∫ (1 − γ )v T b( z )ϕn ( z )dz[29]

Tn (0) = (C I − C 0 )N n−1/ 2 sin( L βn )/ βn

[30]

0

A computer program implementing the presented solution isavailable from the senior author. For simplicity we consideredin this study a constant inlet concentration, C0(t) = C0, and aconstant initial concentration, CI(z) = CI. A solution for the caseof a time-varying inlet and depth-varying initial concentration
can also be obtained, although in that case the computation of
Tn is more complicated due to the time dependence of Gn (seeLiu et al., 2000).

FIG. 4. Illustration of the Michaelis–Menten and two-line uptake
models. Also shown is a dilute solution approximation of the Michaelis–Menten model.
www.vadosezonejournal.org · Vol. 6, No. 4, November 2007894

It is also of interest to consider the steady-state (∂C/∂t = 0)
solution of the transport equation (Eq. [8]) for the special caseof no dispersion (D = 0). In this case, a solution can be obtainedby separation-of-variables, with the result being

Figure 5 shows concentration profiles computed for successive times. At very early times (5 d, in this example), the solute profile is qualitatively similar to that predicted with standard advection–dispersion models. The similarity ends quickly,
however, as the concentration behind the front becomes larger
than C0 while the leading edge retains the familiar dispersive
shape. This solute build-up is due to the concentrating effects ofa water flux that decreases with depth. At later times the solute
profile reaches a steady state, with solute concentrations increasing from a minimal value at the soil surface to a maximal value
at and below the base of the root zone. The maximum concentration that will be achieved is a function of the amount ofwater and solute being taken up, as determined by the leaching
fraction (LF) and the solute uptake coefficient (γ). When γ ≠ 0,
the maximum concentration is also affected by the value of thedispersion coefficient.
Figure 6 illustrates the effects of the uptake coefficient (γ)
and dispersivity (αL) on the solute profile. The profiles shown as
solid lines are for a large time value such that solute concentrations within the root zone (z ≤ 6δ) have obtained approximately

⎡
⎤1−γ
1
⎥
C ( z )/ C 0 = ⎢[31]⎢ 1 − (1 − L F )B( z ) ⎥⎣⎦where the leaching fraction is defined as before, that is, LF = (v0− vT)/v0 = (q0 − T)/q0.

Results and Discussion
General Model Characteristics

Figures 5 and 6 contain plots that illustrate general characteristics of the advection–dispersion model with uptake of waterand solute. The scenarios depicted in these figures involve, beginning at time t = 0 d, the application of water with solute concentration C0 to a soil that has an initially uniform concentration of
CI = 0.25 C0. The solid lines in the plots are concentration vs.
depth profiles computed with the analytical solution presented
in the previous section (Eq. [22]). The dashed lines are steadystate profiles computed with Eq. [31] for D = 0. The simulations
assume an exponential water uptake distribution (Eq. [11]) anda dispersion coefficient that was specified asD( z ) = α L v ( z ) = α L [ v 0 − v T B( z )][32]

where αL is the dispersivity (L). Model and algorithm parametervalues are given in the figure captions. Note that the computational domain length L was much larger than the plotted depth
such that the plotted results are for an effectively semi-infinite
profile. Also, for simplicity we assumed a value of R = 1 for theretardation factor; using a larger R value merely shifts the plots
to later times.

FIG. 6. Plots of solute concentration vs. depth computed for various
values of the transport parameters γ and αL. Except where indicatedotherwise, the concentration profiles were computed using the following parameter values: t = 1000 d, R = 1, δ = 20 cm, v0 = 1 cm d−1,
vT = 0.75 cm d−1, γ = 0.1, αL = 10 cm, CI = 0.25 C0, L = 1000 cm,
and M = 80. Dashed lines are steady-state results computed for D= 0. (t = time; R = retardation factor; δ = root distribution parameter;q0 = water flux at soil surface; T = transpiration rate; θ = volumetricwater content; v0 ≡ q0/θ; vT ≡ T/θ; γ = first-order solute uptake coefficient; αL = dispersivity; D = dispersion coefficient; CI = initial soil
solution concentration; C0 = infiltrating solution concentration; L =depth of computational domain; M = number of terms summed inanalytical solution.)

FIG. 5. Plots of solute concentration vs. depth computed for successive times with the advection–dispersion transport model. The
concentration profiles were computed using the following parametervalues: R = 1, δ = 20 cm, v0 = 1 cm d−1, vT = 0.75 cm d−1, γ = 0.1,
αL = 10 cm, CI = 0.25 C0, L = 1000 cm, and M = 80. Dashed lines are
steady-state results computed for D = 0. (R = retardation factor; δ =root distribution parameter; q0 = water flux at soil surface; T = transpiration rate; θ = volumetric water content; v0 ≡ q0/θ; vT ≡ T/θ; γ =first-order solute uptake coefficient; αL = dispersivity; D = dispersion
coefficient; CI = initial soil solution concentration; C0 = infiltrating solution concentration; L = depth of computational domain; M = numberof terms summed in analytical solution.)
www.vadosezonejournal.org · Vol. 6, No. 4, November 2007895

their steady values, whereas concentrations a short distance
below the root zone are still increasing. Figure 6a shows thatwhen no solute is excluded from the uptake process (γ = 1), rootzone solute concentrations do not increase above C0; the steadystate profile is simply C(z) = C0. For γ = 0, no solute is taken
up by the roots and the build-up of solute is maximized. Also
shown is the result for the intermediate case γ = 0.5. For γ = 1.5,
active uptake of solute occurs. In this particular case the steadysolute concentration profile decreases with depth, approaching a
constant value that is less than C0 but greater than CI.
Figure 6b illustrates the effects of dispersion on solute profile development. The profiles were computed for three valuesof the dispersivity: 2, 10, and 30 cm. As expected, the results
show that increasing the value of the dispersivity produces more
solute spreading at the front, which increases the time required
for the profile to reach steady state: the αL = 2 cm profile at t =1000 d has very nearly reached steady state, whereas the αL = 30
cm profile is still increasing in the bottom of the root zone. Also
as expected, the computed solute profile approaches the steadystate, D = 0 solution (dashed line, Eq. [31]) as the dispersivity
approaches zero.
The computed profiles in Fig. 6b also illustrate a possible
limitation of the advection–dispersion model as applied to transport with simultaneous water and solute uptake. The problem isthat when solute builds up in the root zone, a positive concentration gradient develops (i.e., directed back toward the soil surface),
which, according to the advection–dispersion model, produces
an accompanying diffusive–dispersive flux back toward the soil
surface. The magnitude of this flux in reality should reflect only
diffusive processes rather than field-scale dispersion processes.
The use of a diffusion–dispersion coefficient value that reflects
field-scale dispersion processes hence causes an excessively large
upward component of the solute flux, leading to near-surface
root zone concentrations that are too large. The higher the dispersivity, the more the surface concentration will exceed C0. The
field-scale dispersivity typically has a value in the range of 5 to
20 cm (Jury et al., 1991), although this parameter often carries a high degree of uncertainty. Given this uncertainty and thepossible limitations of the advection–dispersion model, it would
seem prudent to err on the side of choosing a relatively smallvalue of αL (or D) when using this model. This limitation is not
exclusive to the model developed here; the same considerations
should apply to any advection–dispersion model that featuresroot uptake, such as the comprehensive numerical models notedin the introduction.

Legislative standards for allowable or critical metal loadings
to arable land vary by country. One emerging method for risk
assessment involves simple mass balance calculations wherein
the soil is treated as a single compartment subject to first-order
uptake into harvested plant parts, linear sorption partitioning
in the soil, and advective leaching below the root zone (Pačes,
1998; Tiktak et al., 1998; Blombäck et al., 2000; Bergkvist et al.,2005). Such an approach implicitly assumes uniform well-mixed
conditions within the soil compartment and ignores depth variations of transport and uptake of both water and trace metals
in the root zone. Given these assumptions, and using the same
notation as above, the change in stored amount A (M L−2) isgiven bydA
= q 0C 0 − q zC − γTC[33]dt
where qz is the water flow at the base of the root zone and A = (θ+ ρbKd)RdC. Noting that qz = q0 – T leads todC= I − kC
[34]
dt
where I = (v0C0)/(RRD) and k = [v0 + vT(γ − 1)]/(RRD). When
C(0) = CI, integrating Eq. [34] from t = 0 to t gives⎛I ⎞C ( t ) = ⎜⎜⎜ ⎟⎟⎟(1 − e−kt )+ C Ie−kt⎝k ⎠In the one-compartment model, the concentration of harvested
parts of annual plants Cp (M m−1) is given by
γTC
[36]
Bp
where Bp is the harvested biomass yield (M L−2 T−1). This
expression is appropriate only for plants whose growth cycle isshort relative to the timescale of solute transport in the root zone
(e.g., annual crops). Setting dC/dt = 0 in the one-compartment
model gives the steady-state soil solution concentration as
C p(t ) =C0
1 + (γ − 1)(1 − L F )and the steady-state plant concentration as
C=

CP =

γTC 0
Bp [1 + (γ − 1)(1 − L F )]

[37]

[38]

In our example application, we compare the analyticaladvection–dispersion and one-compartment models for the caseof cadmium uptake into wheat grain, an important staple foodin many parts of the world. The advection–dispersion model inthis example was implemented using a constant dispersion coefficient, D(z) = D, and a linear uptake distribution,

Example Application

Cadmium is one of the most mobile and bioavailable trace
metals and constitutes a potential human health risk throughdietary intake (Buchet et al., 1990). In many locations the cadmium content of arable soils has significantly increased over
the years, primarily from application of cadmium-containing
phosphorous fertilizers and from atmospheric deposition. The
cadmium content in Swedish arable soils, for example, hasapparently increased by about 30% since the beginning of thetwentieth century (Andersson, 1992). The current soil content
of cadmium can result in concentrations in harvested plant parts
that are close to acceptable limits set for human consumption in
many countries (e.g., Eriksson et al., 2000).
www.vadosezonejournal.org · Vol. 6, No. 4, November 2007[35]⎧⎪ 0.8 1.6⎪⎪−zb( z ) = ⎨ R D R D2
⎪⎪
0
⎪⎪⎩

z ≤ RD[39]

z > RD

This uptake distribution conforms to the “40:30:20:10 rule”(e.g., Raats, 1974; Hoffman and van Genuchten, 1983), suchthat 40% of the uptake is in the top quarter of the root zone,
30% in the second quarter, 20% in the third quarter, and 10%
in the bottom quarter. Omitting details, an expression equiva896

Conclusions

lent to Eq. [36] may be developed for the advection–dispersion
model to compute the harvested plant part concentration:
C P(t ) =

1 RD
1 RD
rsdz = ∫ γTb( z )C ( z , t )dz
∫
BP 0
BP 0

Although a variety of comprehensive numerical simulation
models exist for analyzing solute transport processes in the root
zone, regulatory and risk management analyses often require
simpler modeling tools with fewer input parameters. Simpler
analytical models may be especially useful for assessing environmental impacts of strongly sorbing and persistent solutes suchas trace metals. This is because the relevant timescales for suchassessments are measured in decades or centuries, in which casetransport and plant uptake can be modeled effectively using
annual and root-zone average parameters.
In this work we developed an analytical advection–dispersion model that includes plant uptake of water and solute.
Although considerably simpler than many popular comprehensive numerical models, the analytical advection–dispersion
model—with its explicit accounting of advection, dispersion,and uptake processes within the root zone—is much more realistic than one-compartment, well-mixed models that are currently
being used or considered for use in risk assessments. As an example, we compared predictions made with the advection–dispersion and one-compartment models for the uptake of cadmiumby wheat. The comparison showed that the one-compartment
model significantly overestimated the long-term, steady-state
Cd concentration in harvested wheat grain. For this reason, we
recommend the analytical advection–dispersion model as a preferred tool for environmental risk assessment of strongly sorbing,persistent solutes.
Possible limitations of the presented model arise due toassumptions of linear solute sorption and first-order plant
uptake. We thus suggest that model applications be restricted
to cases of slight to moderate soil contamination where theseassumptions are likely to be reasonable. Additionally, the model[40]

The uptake coefficient γ is an important parameter since it
determines the extent to which cadmium in the soil solution is
excluded from the plant. Measurements of cadmium uptake haveindicated that uptake is linear over a range of concentrations
relevant for sludge applications (Hamon et al., 1999; McGrathet al., 2000). Measurements of cadmium concentrations in soil
solutions and wheat grain (Bergkvist et al., 2005), together withtypical values of the water use efficiency T/Bp for wheat (see Eq.
[36]), suggest a value for γ of approximately 0.05, which we usedin the calculations presented below. The remaining parametervalues are shown in the caption of Fig. 7, which compares predictions of the two models. The assumed cadmium loading rate
(1.6 mg m−2 yr−1) reflects all sources of cadmium. In the wheatgrowing areas of Europe, diffuse aerial deposition, together with
impurities in lime, commercial fertilizer, and feeds, would notsupply more than approximately 0.3 to 0.5 mg Cd m−2 yr−1.
Thus, for the scenario depicted in Fig. 7, we implicitly assume
a significant additional local cadmium source such as industrial
emissions or sludge applications. The assumed loading rate is
not unrealistic in that the current allowable limit in the EU for
cadmium supplied with sludge is set at 15 mg m−2 yr−1 (EUdirective, 86/278/EEC).
Figure 7 shows that for the above parameterization, theanalytical advection–dispersion model predicts a final steadystate Cd concentration in wheat grain of approximately 0.2 mg
kg−1, reached after about 500 yr, which equals the allowable
concentration under EU directives. The root depth in the onecompartment model is usually set to the plow depth (?0.25 m)
since this best approximates the assumption of uniform mixing (e.g., Tiktak et al., 1998). This assumption ignores subsoil
uptake of cadmium, which can be considerable (Johnsson et al.,2002). Figure 7 shows that with RD = 0.25 m, the one-compartment model substantially overestimates the grain cadmium
concentration, although the pattern of increase is similar to theanalytical advection–dispersion model, with a steady-state value
approached after approximately 500 yr. Increasing RD to 1 m, tomatch the maximum depth of roots in the advection–dispersion
model, produces good agreement with the advection–dispersion
model at early times (up to 300 yr). However, the final steadystate value of the one-compartment model is independent of RD(see Eq. [38]), being about 75% larger than the value predictedby the advection–dispersion model. In terms of risk assessment,such a difference in model results is highly significant. We conclude that the simplifications in the one-compartment model
can lead to erroneous decisions, although the predicted deviations in the present example would become significant only aftera few hundred years depending on the value of RD used in the
one-compartment model. Finally, we note that the analyticaladvection–dispersion model requires only a few extra parameters, which are easily estimated. The advection–dispersion model
therefore should be recommended as the preferred tool for environmental risk assessment of persistent contaminants.

www.vadosezonejournal.org · Vol. 6, No. 4, November 2007

FIG. 7. A comparison of predictions of the analytical advection-dispersion and one-compartment models for cadmium uptake into wheat
grain using the following parameter values: q0 = 0.7 m yr−1, T = 0.5m yr−1, γ = 0.05, Bp = 0.5 kg m−2 yr−1, θ = 0.4 m3 m−3, R = 376, D= 0.015 m2 yr−1, CI = 1 mg m−3 and C0 = 2.286 mg m−3, equivalent
to a cadmium loading of 1.6 mg m−2 yr−1. The root depth, RD, in theadvection–dispersion equation model was set to 1 m, and two valuesof RD are shown for the one-compartment model (0.25 and 1 m). (q0= water flux at soil surface; T = transpiration rate; γ = first-order solute uptake coefficient; Bp = harvested biomass yield; θ = volumetricwater content; R = retardation factor; D = dispersion coefficient; CI =initial soil solution concentration; C0 = infiltrating solution concentration; RD = maximum rooting depth.)897

Taylor, W.R. Jordan, and T.R. Sinclair (ed.) Limitations to efficient water
use in crop production. ASA, CSSA, and SSSA, Madison, WI.
Johnsson, L., I. Öborn, G. Jansson, and D. Berggren. 2002. Evidence for spring
wheat (Triticum aestivum) uptake from subsurface soils. In Y. Zhu and Y.
Tong (ed.) Proc. of the Int. Workshop of Soil-Plant Interactions, Beijing,
China. 27–30 Oct. 2002. Kluwer Academic, Dordrecht, the Netherlands.
Jury, W.A., W.R. Gardner, and W.H. Gardner. 1991. Soil physics. John Wiley& Sons, New York.
Larsbo, M., S. Roulier, F. Stenemo, R. Kasteel, and N.J. Jarvis. 2005. An
improved dual-permeability model of water flow and solute transport in
the vadose zone. Vadose Zone J. 4:398–406.
Liu, C., J.E. Szecsody, J.M. Zachara, and W.P. Ball. 2000. Use of the generalized
integral transform method for solving equations of solute transport inporous media. Adv. Water Resour. 23:483–492.
McGrath, S.P., F.J. Zhao, S.J. Dunham, A.R. Crosland, and K. Coleman. 2000.
Long-term changes in the extractability and bioavailability of zinc andcadmium after sludge application. J. Environ. Qual. 29:875–882.
Mualem, Y. 1976. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour. Res. 12:513–522.
Ozisik, M.N., and R.L. Murray. 1974. On the solution of linear diffusion
problems with variable boundary condition parameters. J. Heat Transfer
96:48–51.
Pačes, T. 1998. Critical loads of trace metals in soils: A method of calculation.
Water Air Soil Pollut. 105:451–458.
Raats, P.A.C. 1974. Steady flows of water and salt in uniform soil profiles withplant roots. Soil Sci. Soc. Am. Proc. 38:717–722.
Raats, P.A.C. 1975. Distributions of salts in the root zone. J. Hydrol. 27:237–
248.
Raats, P.A.C. 1981. Residence times of water and solutes within and below theroot zone. Agric. Water Manage. 4:63–82.
Rawlins, S.L. 1973. Principles of managing high frequency irrigation. Soil Sci.
Soc. Am. Proc. 37:626–629.
Russo, D. 1988. Determining soil hydraulic properties by parameter estimation:
On the selection of a model for the hydraulic properties. Water Resour.
Res. 24:453–459.
Schoups, G., and J.W. Hopmans. 2002. Analytical model for vadose zone solutetransport with root water and solute uptake. Vadose Zone J. 1:158–171.
Šimůnek, J., D.L. Suarez, and M. Šejna. 1996. The UNSATCHEM software
package for simulating one-dimensional variably saturated water flow, heattransport, carbon dioxide production and transport, and multicomponent
solute transport with major ion equilibrium and kinetic chemistry, Version
2.0. Res. Rep. no. 141. U.S. Salinity Lab., Riverside, CA.
Šimůnek, J., M.Th. van Genuchten, and M. Šejna. 2005. The HYDRUS1D software package for simulating the one-dimensional movementof water, heat and multiple solutes in variably-saturated media, Version
3.0. HYDRUS Software Ser. no. 1. Dep. of Environ. Sci., Univ. of Calif.,
Riverside.
Skaggs, T.H., and F.J. Leij. 2002. Solute transport: Theoretical background. p.
1353–1380. In J. Dane and C. Topp (ed.) Methods of soil analysis. Part 4.
ASA and SSSA, Madison, WI.
Skaggs, T.H., M.Th. van Genuchten, P.J. Shouse, and J.A. Poss. 2006.
Macroscopic approaches to root water uptake as a function of water andsalinity stress. Agric. Water Manage. 86:140–149.
Streck, T., and H. Piehler. 1998. On field-scale dispersion of strongly sorbing
solutes in soils. Water Resour. Res. 34:2769–2773.
Tiktak, A., R. Alkemade, H. van Grinsven, and K. Makaske. 1998. Modelingcadmium accumulation at a regional scale in the Netherlands. Nutr.
Cycling Agroecosyst. 50:209–222.
Tinker, P.B., and P.H. Nye. 2000. Solute movement in the rhizosphere. Oxford
Univ. Press, New York.
van Genuchten, M.Th., and W.J. Alves. 1982. Analytical solutions of the onedimensional convective-dispersive solute transport equation. USDA Tech.
Bull. 1661. U.S. Gov. Print. Office, Washington, DC.
Wierenga, P.J. 1977. Solute distribution profiles computed with steady-state andtransient water movement models. Soil Sci. Soc. Am. J. 41:1050–1055.
Yuan, F., and Z. Lu. 2005. Analytical solutions for vertical flow in unsaturated,rooted soils with variable surface fluxes. Vadose Zone J. 4:1210–1218.

was derived based on an assumption of uniform soil water content. Among other factors, water content uniformity is affectedby soil texture, rooting depth, and leaching fraction, with finetexture, shallow rooting, and high leaching all promoting uniformity. In general, a leaching fraction of greater than about 0.25
corresponds to a fairly uniform profile when steady-state waterflow occurs, except perhaps in very coarse-textured soils where
substantial nonuniformity could arise if the surface soil were to
be maintained at a high degree of saturation. Where nonuniform
soil conditions prevail, the model may be used if an appropriate,
effective water content can be identified.
ACKNOWLEDGMENTS

This research was supported in part by the Swedish Environmental
Protection Agency funded project “Cadmium leaching from arable land
and uptake into crops: calculation of critical loads and developmentof a simulation tool for dynamic long-term predictions”; a fellowshipfrom CNPq (Brazil) to E. M. Pontedeiro; and the SAHRA Science andTechnology Center as part of NSF grant EAR-9876800.

References
Ahuja, L.R., K.W. Rojas, J.D. Hanson, M.J. Shaffer, and L. Ma (ed.). 2000. Rootzone water quality model: Modeling management effects on water quality
and crop production. Water Resources Publications, LLC, Highlands
Ranch, CO.
Andersson, A. 1992. Trace elements in agricultural soils: Fluxes, balances, andbackground values. Report 4077. Swedish Environmental Protection
Agency, Solna.
Bergkvist, P., N.J. Jarvis, J. Eriksson, and L. Rapp. 2005. Critical load of cadmium
on arable soils in Sweden. Emergo no. 4, Studies in the BiogeophysicalEnvironment. Dep. of Soil Sci., Swedish Univ. of Agricultural Sciences,
Uppsala.
Blombäck, K., M. Kalvi, J. Eriksson, and I. Öborn. 2000. Cadmium accumulationin agricultural soils and uptake by plants as an effect of cadmium inphosphate fertilizers: Model estimations. p. 11–31. In Assessment of risks
to health and the environment in Sweden from cadmium in fertilizers.
Report 4/00. National Chemicals Inspectorate, Stockholm.
Borg, H., and D.W. Grimes. 1986. Depth development of roots with time: An
empirical description. Trans. ASAE 29:194–197.
Buchet, J.P., R. Lauwerys, H. Roels, A. Bernard, P. Bruaux, F. Claeys, G.
Ducoffre, P. de Plaen, J. Staessen, A. Amery, P. Lijnen, L. Thijs, D. Rondia,
F. Sartor, A. Saint Remy, and L. Nick. 1990. Renal effects of cadmium
body burden of the general population. Lancet 336:699–702.
Canadell, J., R.B. Jackson, J.R. Ehleringer, H.A. Mooney, O.E. Sala, and E.-D.
Schulze. 1996. Maximum rooting depth of vegetation types at the global
scale. Oecologia 108:583–595.
Cotta, R.M. 1993. Integral transforms in computational heat and fluid flow.
CRC Press, Boca Raton, FL.
Dalton, F.N., P.A.C. Raats, and W.R. Gardner. 1975. Simultaneous uptake ofwater and solutes by plant roots. Agron. J. 67:334–339.
Destouni, G. 1991. Applicability of the steady state flow assumption for solute
advection in field soils. Water Resour. Res. 27:2129–2140.
Eriksson, J., A. Andersson, B. Stenberg, and R. Andersson. 2000. Current statusof Swedish arable soils and cereals. Report 5062. (Summary in English.)
Swedish Environmental Protection Agency, Stockholm.
Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water use
and crop yield. PUDOC, Wageningen, the Netherlands.
Gardner, W.R. 1958. Some steady-state solutions of the unsaturated moistureflow equation with application to evaporation from a water table. Soil Sci.
85:228–232.
Ginn, T.R., and E.M. Murphy. 1997. A transient flux model for convective
infiltration: Forward and inverse solutions for chloride mass balance
studies. Water Resour. Res. 33:2065–2079.
Hamon, R.E., P.E. Holm, S.E. Lorenz, S.P. McGrath, and T.H. Christensen.
1999. Metal uptake by plants from sludge-amended soils: Caution isrequired in the plateau interpretation. Plant Soil 216:53–64.
Hoffman, G.J., and M.Th. van Genuchten. 1983. Soil properties and efficientwater use: Water management for salinity control. p. 73–85. In H.M.
www.vadosezonejournal.org · Vol. 6, No. 4, November 2007

898