Journal of Hydrology, 151 (1993) 193-228

193

0022-1694/93/$06.00 0 1993 ~ Elsevier Science Publishers B.V. All rights reserved[2]

Analytical solutions for non-equilibrium solutetransport in three-dimensional porous media
Feike J. Leij*, Nobuo Toride, Martinus Th. van Genuchten
US Salinity Laboratory, US Department of Agriculture, Agricultural Research Service, Riverside,
CA 92501, USA(Received 13 January 1993; revision accepted 23 March 1993)

Abstract
The movement of water and chemicals in soils is generally better described with multidimensional nonequilibrium models than with more commonly used one-dimensional and/or equilibrium models. This
paper presents analytical solutions for non-equilibrium solute transport in semi-infinite porous media
during steady unidirectional flow. The solutions can be used to model transport in porous media where
the liquid phase consists of a mobile and an immobile region (physical non-equilibrium) or where solute
sorption is governed by either an equilibrium or a first-order rate process (chemical non-equilibrium). The
transport equation incorporates terms accounting for advection, dispersion, zero-order production, andfirst-order decay. General solutions were derived for the boundary, initial, and production value problems
with the help of Laplace and Fourier transforms. A comprehensive set of specific solutions is presented
using Dirac functions for the input and initial distribution, and/or Heaviside or exponential functions forthe input, initial, and production profiles. A rectangular or circular inflow area was specified for theboundary value problem while for the initial and production value problems the respective initial andproduction profiles were located in parallelepipedal, cylindrical, or spherical regions of the soil. Solutions
are given for both the volume-averaged or resident concentration as well as the flux-averaged or flowing
concentration. Examples of concentration profiles versus time and position are presented for selected
problems. Results show that the effects of non-equilibrium on three-dimensional transport are very similar
to those for one-dimensional transport.

Introduction
Many experimental and theoretical studies have been undertaken to
improve the understanding, management, and prediction of the movementof dissolved substances in soils. These investigations are primarily motivatedby concerns about possible contamination of the subsurface environment.
Although several uncertainties still exist, and probably will remain for sometime, regarding the correct description of water and solute transport processes* Corresponding author.

194

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

in field soils, subsurface solute transport is generally described with theadvection-dispersion equation (ADE). In the deterministic approach,explicit closed-form solutions for the transport problem can often be derived
if the model parameters are constant with respect to time and position. Onedimensional solutions of the ADE have been widely reported for equilibriumtransport (e.g. Van Genuchten and Alves, 1982) and to some extent for nonequilibrium transport (Lassey, 1988; Van Genuchten and Wagenet, 1989). In
spite of the fact that many transport problems are multidimensional, relatively
few analytical solutions have been reported for two- and three-dimensionaltransport (e.g. Wexler, 1989), particularly for non-equilibrium transport(Goltz and Roberts, 1986). This is surprising since such analytical solutions
can be deduced quite easily from one-dimensional solutions.
Goltz and Roberts (1986) provided analytical solutions for two-region-type
mobileeimmobile models applicable to transport from an instantaneous point
source in an infinite porous medium made up of rectangular, cylindrical, or
spherical aggregates. Transport between the two regions was governed by afirst-order rate expression. Analytical solutions for the concentration distribution inside the soil aggregate were obtained by solving Fick’s law for
diffusion assuming different aggregate geometries. The solutions by Goltz
and Roberts (1986) are applicable to transport involving horizontal flow in
infinite aquifers. An alternative, and equally common, pollution scenario
involves downward flow from the soil surface. The semi-infinite nature ofsuch a transport system makes the solution slightly more complicated, butalso more general than the infinite case treated by Goltz and Roberts. Since
the size and geometry of soil aggregates is poorly defined for field problems,
no attempt is made to solve the more complicated problem for solute transport in a porous medium consisting of idealized aggregates.
The problem to be investigated here is of a very similar nature to, but morecomprehensive than, the equilibrium case outlined by Leij et al. (1991). The
solution procedure is based on previous work for one-dimensional nonequilibrium transport assuming a first-order rate exchange (Van Genuchten
and Wagenet, 1989; Toride et al., 1993) and multidimensional transport (Leijet al., 1991). The objective of this paper is to combine these two approaches
and to provide a relatively concise set of solutions of the three-dimensional
non-equilibrium ADE for a wide variety of boundary, initial, and production
value problems.
For clarity of presentation both the ‘chemical’ and ‘physical’ nonequilibrium models will be formulated. These models will be cast in thesame dimensionless form for subsequent mathematical analysis. Afteroutlining the general solution procedure, analytical solutions are presented
for the equilibrium and non-equilibrium concentrations as volume- and flux-

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

195

averaged types. General solutions will be given for: (1) the boundary value
problem (BVP) with an arbitrary solute input, which depends on time andtransverse spatial coordinates; (2) the initial value problem (IVP) using an
arbitrary initial distribution, which depends on all three spatial coordinates;(3) the production value problem (PVP) based on an arbitrary production
profile, which depends on time and the three spatial coordinates. Although thegeneral solutions are relatively inconvenient for further use, they can begreatly simplified for some practical input, initial, and production profiles
described by Dirac, Heaviside, or exponential functions.

Model description
Solute movement in an isotropic soil during one-dimensional downward
flow and with three-dimensional dispersion is modeled analytically by assuming that both the velocity and the dispersion coefficients are constant withrespect to time and space. These assumptions are approximate but may beadequate for many applications. Furthermore, solutes are assumed to be
subject to linear retardation (i.e. equilibrium sorption between solution andsorbed phases can be described with a linear isotherm), and zero- and firstorder production or degradation processes. Non-equilibrium is attributed todifferences in exchange sites (assuming equilibrium sorption for ‘type-l’ andnon-equilibrium sorption for ‘type-2’ sites), or to mobile and immobile fluid
regions (Van Genuchten and Wagenet, 1989).
In the two-site model, non-equilibrium sorption for type-2 sites is typicallydescribed with a first-order rate law (Selim et al., 1976; Cameron and Klute,1977). The governing equations for this ‘chemical’ non-equilibrium model
may be written as

-

ask

at =

7 [(l -f)kc - s,J - plc -f+c + X1(x, y, 2, f)a[(1 -f>kc - Sk] - /i,s ksk + c1 -J)&,k(X>_h z, f>

where c is the volume-averaged solute concentration in the liquid phase(MLe3), Sk is the concentration in the sorbed phase for type-2 sites

196

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228(MM-‘); t is time (T); x is position (L) in the direction of flow; y and z arerectangular coordinates perpendicular to the flow direction (L), i.e. parallel tothe soil surface; D, is a dispersion coefficient (L2T-‘) in the x direction (i.e.
representing longitudinal dispersion); D, and D, are transverse dispersion
coefficients in the y and z directions, respectively; 13 is the volumetric watercontent (L3LP3); v is the macroscopically averaged downward pore-watervelocity (LT-‘);fis the fraction of type-l exchange sites, p is the bulk density(MLe3); k is a distribution coefficient for linear sorption (MLP3), LY is a firstorder kinetic rate coefficient (T-l); pI and ps are first-order decay coefficientsfor degradation in the liquid and sorbed phases, respectively (T-l); X1(MLm3TP’) and A, (MM-‘T-l) are time- and position-dependent zero-order
rate coefficients for solute production in the liquid and sorbed phases, respectively; while the subscripts e and k refer to equilibrium and kinetic exchangesites, and 1 and s refer to the liquid and sorbed phases.
The ‘physical’ non-equilibrium model assumes that the liquid region can bepartitioned into mobile (or flowing) and immobile (or stagnant) regions, thatsolute movement occurs by advection and dispersion in the mobile region, andthat solute exchange between the two regions occurs by first-order diffusion(Coats and Smith, 1964; Van Genuchten and Wierenga, 1976). The governing
equations for the two-region model are[ei, + (1 -flpk] 2 = Ck!(C, - Cim) - [4mPl,im + (l -J)PkPs,itnlCim+ eimXl,im(X~Y~z~ l> + (l -J)PXs,im(X,Y7z, l)( 2b)where in this casefrepresents the fraction of sorption sites in equilibrium withthe fluid of the mobile region and a is a first-order mass transfer coefficient(T-l); the subscripts m and im refer to the mobile and immobile liquid regions,
respectively, with 13 = 13~ + Oim.
The two-site and the two-region models can be cast in the same dimensionless form using the parameters listed in Table 1. Notice that several
parameters in Table 1 are expressed in terms of an arbitrary characteristic

197

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

Table 1
Dimensionless parameters for the two-site and two-region transport models
ParameterTwo-site modelTwo-region model

TX,y,zpx,p,,pz

4 +fpk

P
R
W

Cl

c2

B+pk

I+@e
cy(1 - P)RL
u
C
-

1+@

e

CYL
ev
cm

co

co
Sk

Cim
-

(1 -fPco

co

PI

P2

XI

x2

; [(I -f )Pk’4,kl& (eh +fPu
& It1 -f)ds,kl& PiA,im i-(1 -f)P&,im]concentration, co, and characteristic length, L. The following dimensionlessmodel can be obtained (see also Nkedi-Kizza et al., 1984)(1 - @)I?$+ = w(C1 - C2) - P2C2 + X2Wl y, 6 T)

(3b)where ,L? is a partition coefficient, R is a retardation factor, Ct and C2 areequilibrium and non-equilibrium resident concentrations, respectively, T istime; X is the longitudinal coordinate; Y and Z are transverse coordinates;

198

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

P,, Py and P, are the Peclet numbers in the X, Y and Z directions; w is a masstransfer coefficient; p and X are dimensionless first-order decay and zero-orderproduction terms; and the subscripts 1 and 2 refer to the equilibrium and nonequilibrium phases, respectively. The magnitude of several dimensionlessparameters in Eq. (3a,b) depends not only on the dimensional model
parameters but also on the arbitrary constants co and L. Unless specified
otherwise, we assume that the solute concentration is of the volume-averaged
or resident concentration type.
The initial and boundary conditions for Eq. (3) arec, (X, Y, z, 0) = c&K K z, 0) = f(K y, z>$3, Y,Z, T) = 0
Cl (X, fcq z, T) = 0~(x,*cm,Z,T) =o(6a, b)c&Y, Y, *co, T) = 0

g+(X, Y,hqT) =o(7a, b)wherefand g are arbitrary functions which will be specified later, along with
A, to illustrate pertinent transport problems. Because advective-dispersive
transport occurs exclusively in phase 1, the boundary conditions are only
cast in terms of Ct. The linearity of the problem enables us to solve theboundary (BVP), initial (IVP), and production (PVP) value problems
separately. In the BVP it is assumed that Xi and X2 in Eq. (3) and g(Y, Z,
T) in Eq. (5a) are equal to zero. The IVP is solved by dropping the production
terms, X, and by setting f(X, Y, Z) equal to zero. The PVP requires that the
input function, g( Y, Z, T), and the initial distribution, f(X, Y, Z), are set tozero.
Concentrations obtained from the effluent of finite soil columns or fromsolution samplers are usually viewed as representing the equilibrium concentration, Cr. On the other hand, soil coring or non-destructive techniques suchas TDR or electrical conductivity will give total concentrations, i.e. thosedefined as mass of solute per unit volume of soil, with weighted contributions
from both the equilibrium (C,) and non-equilibrium (C,) phases according toC, = ,BRCl + (1 - P) RC,

(8)

199

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

The injection and detection mode for several solute displacement experiments may require the use of flux-averaged or flowing concentrations(Kreft and Zuber, 1978; Parker and Van Genuchten, 1984). They can beobtained by substituting the analytical expression for the resident concentration, C or C’, obtained for a third-type inlet condition, into the well-known
transformation
1 acr

cf=cr---

P, ax
This equation can be readily applied if the expression for C’ is differentiablewith respect to X and if the flow field is uniform, i.e. having a constant Peclet
number.
General solutions
The solution of Eq. (3), subject to Eqs. (4)-(7) was obtained using Laplacetransforms with respect to X and T, and a double Fourier transform withrespect to Y and Z. This approach was discussed previously by Leij et al.
(199 l), among others. Further information regarding the use of integral transforms can be found in Spiegel (1965) and Sneddon (1972). The solution of the
non-equilibrium three-dimensional ADE is tedious but relatively straightforward. Details regarding the derivation of all presented analytical solutions
are available from the senior author upon request. The inverse Laplace transform with respect to Xutilizes the shifting property, the convolution theorem,
and a table of Laplace transforms. Subsequently, the last remaining boundary
condition (Eq. 5b) is used in conjunction with Leibnitz’ rule for differentiationof an integral. The Laplace inversion with respect to T is carried out by using
the property that the iterated Laplace transform of a function is equal to theordinary Laplace transform of the generalized convolution integral of thisfunction (De Smedt and Wierenga, 1979a; Walker, 1987). This method wasapplied by Lindstrom and Narasimhan (1973) to derive the solutions for an
initial value problem involving non-equilibrium sorption. The inverse Fouriertransformation constituted the last step of the solution procedure.
For an arbitrary input function, g( Y, Z, T ) the general solution to the BVP,
CB, isc,” =g(Y-gS,Z--R,T-r)+f’~(M,+~(~,r)

1dfld+dr214

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

useful foundation for deriving specific analytical solutions applicable tosimple boundary, initial, or production conditions. We derived solutions for
the BVP and IVP assuming a Dirac distribution, and for the BVP, IVP, or
PVP assuming Heaviside or exponential distributions in (1) semi-infinite
rectangular, (2) finite rectangular, (3) cylindrical, and (4) spherical regions.
All solutions are based on Fourier transforms and a Cartesian coordinate
system.
The different areas where solute is applied (BVP), initially distributed (IVP),or produced (PVP), are sketched in Fig. 1. The first three cases are in thetransverse plane while the spherical case, which pertains only to the IVP and
PVP, includes the longitudinal coordinate. The magnitude for g needs to be
specified as a function of time for the BVP whilef for the IVP and X for the
PVP should be given as a function of the longitudinal distance for cases I, II
and III, or the radial distance for case IV. The input functions for the BVP, theinitial concentrations for the PVP, and the X profile for the PVP are given
schematically in Fig. 2. We continue by presenting the solutions for thevarious cases of the BVP, IVP and PVP, along with some brief comments
regarding the method of solution.
Table 4(a) lists solutions for the BVP, and specifications of the input functions and the type of transverse auxiliary function on which they are based.
The derivation of these transverse functions, I’, is fairly straightforward; thespecific independent variables for each problem are listed in the left-handcolumn of Table 4(a). The entire solution can be obtained by inserting theauxiliary functions given in Tables 2 and 3. For the Dirac input, i.e. for thecase where a solute mass m is applied instantaneously at position ( YO, Za) at
the surface at time TO, the expression for C, can be readily derived from thegeneral solution. The solution for C2 was obtained with the help of Fubini’s
theorem. In the case of a Heaviside-type input, an arbitrary number of solute
pulses with different constant concentrations is applied. The problem was firstcast in terms of the J-function, after which the expressions could be simplified
by changing the order of integration and integration by parts using some ofthe properties of the J-function listed in appendix A of Van Genuchten (1981).
The integration with respect to the transverse coordinates, Y and Z, was basedon elementary properties of the error function (Gautschi, 1965). Case I concerns solute application to a semi-infinite rectangular area; its finite corner
position ( Yi, ZJ may vary for each pulse. Similarly, cases II and III involve
solute application to finite rectangular and circular areas.
The third input function consists of a term which decreases exponentiallywith time, as specified in Table 4(a). The solute in this case is applied to thesame region as for the Heaviside input function; however, the transverse
dimensions are time invariant. Evaluation of the general solution containing

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

215

go is very similar as for the Heaviside input. Evaluation of the gl function withexponential decrease was done by applying Fubini’s theorem and partial
integration. The resulting expression consisted of the product of an exponential function and the zero-order modified Bessel function. Since no expression
for this integral appears to be available (cf. Luke, 1962) the series representation 9.6.12 of Olver (1965)

z,[$izj] =l+‘T+f$T2+gT3+..
(1!)2
!2
!2

(16)

was used to integrate (2.321.2) of Gradshteyn and Ryzhik (1980). This
procedure led to’ exp (Y&K b-1 dr
I11

qv) =

(

Ty -p
n=l

w (T-P and 33 = YB - (:‘$

k=* yBk(n _ k)!

)uw

In most cases a convergent solution was obtained for a maximum y1 value of25; however, the solution was less accurate for, among other factors, relatively
large values for a and T. Similar problems occur for the series approximationof the J-function by De Smedt and Wierenga (1979b) (cf. Lassey, 1982).
Table 4(b) contains the different types of initial conditions for which the
IVP was solved, along with the definitions of variables in the auxiliary transverse function, and the semi-infinite and finite rectangular, cylindrical andspherical regions. For the Dirac initial condition, i.e. for an instantaneous
release of a solute mass m at (X0, Yo, Z,), the IVP can be readily solved. Notethat evaluation of the term containing f in Ci is impractical for the Diraccondition. Specific solutions for the Heaviside-type initial distribution, consisting of y1 semi-infinite or finite parallelepipeds, cylinders, or spherical shells,
were obtained through integration with respect to the longitudinal and transverse coordinates. The identity

could be established with properties of the error function and the Laplacetransform; the expression for Gf was included in Table 2. The solution for C:
again involved the application of Fubini’s theorem to evaluate one of the

216

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

integrals with respect to time. For the spherical IVP, where the solute concentration is uniformly distributed in shells, the evaluation of the integrals washampered by the presence of variable integration limits. For the exponentiallydecaying distribution with depth (0 < X < oo), the integration with respect tothe transverse coordinates was the same as for the Heaviside distribution withan arbitrary number of steps. Integration in the longitudinal direction was
based on Eq. (18) andO3 exp (-~rW’(K 7) dt = G#, 7)(19)
10
which was derived through integration by parts and using properties of theerror function and expression (3.322.2) of Gradshteyn and Ryzhik (1980).
This approach is an alternative to evaluating integrals using Laplace transforms. Note that the IVP was not solved for an initial profile that decreases
exponentially from the center of a spherical region.
Finally, several specific solutions were obtained for the PVP assuming only
a spatial dependency of Xk. Time dependency can easily be introduced (cf.
Lindstrom and Boersma, 1989), although it is then quite difficult to simplify
the general solution significantly. Table 4(c) contains the production functions
and the expressions for the concentrations and auxiliary functions, I’j, forseven problems. The Xk profiles were described with Heaviside and exponential functions for the same geometries as the PVP. The methodology to obtainthe specific solutions is similar as used previously for the BVP and the IVP.
All integrals in the analytical solutions of Table 4 were evaluated with
Gauss-Chebyshev quadrature (Carnahan et al., 1969). This method allows
greater flexibility in the number of integration points, as opposed to, forexample, Gaussian integration. A disadvantage, particularly for the solutions
involving multiple integration, is the increased amount of computer time. If
this is of major concern, one may prefer to solve the entire problem numerically. Further information regarding the computer program to evaluate thesolutions in Table 4 can be obtained from the senior author.
ExamplesIn this section we provide several illustrations of results obtained with theanalytical solutions from Table 4. Since the three-dimensional non-equilibrium model is based on existing and generally well-known one-dimensional
non-equilibrium and three-dimensional equilibrium models, we will not focus
on the conceptual aspects of the predicted concentration profiles, nor will we
provide a complete overview of results for all different cases. For example,
Dirac functions, spherical geometries, and flux-averaged concentrations willF.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

219

Fig. 4. IVP involving a Heaviside profile of finite rectangular regions with v = 50, R = cO = 1, D, = 20,
Dy = D, = 5, Jo, = kL2 = 0, and L = 100. Two-dimensional solute profiles at z = 0: (a) initial distribution
f0rC1 = c2 with cfl,f2,hfi)=(~, 1, 0.5, O), (XI> x2, x3, x4)=@, 5,10, 15), CY2,.~3)=(~2r 23) =(5, 10); (b) CI
profile at t = 0.5 for p = 0.5 and w = 1; (c) C2 profile at t = 0.5 for p = 0.5 and w = 1; (d) CT profile att = 0.5 for 0 = 0.5 and w = 1; (e) C, profile at r = 0.5 for p = 0.9 and w = 5.

not be treated. The model parameters are dimensional in all examples given
here, except for CL, X and w; no units will be specified for these rate parameters,
and hence the results can be interpreted using any dimensionally consistent setof parameters. The dimensionless concentration is investigated as a functionof dimensional time and position.
The first example involves breakthrough curves stemming from a result of a
Heaviside step input or initial condition for a semi-infinite transverse region.
Figure 3(a) presents the equilibrium resident concentration for the BVP as afunction of t and y for pulse and step input, with intermediate application ofsolute-free water, at the surface region given by y < 0 and z < 0. Duringtransport through the soil the magnitude of the first solute peak, gl = 5,

220

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

has been reduced considerably at the soil outlet. Notice that a steady-state
concentration equal to g2 = 2 is established for t > 3. Figure 3(b) shows thebreakthrough curve for the IVP with a similar initial profile as the input
function for the BVP in Fig. 3(a). The solute peaks are somewhat higher
and a clear dip between solute peaks can be distinguished at about t = 1, in
contrast to the BVP. The transverse solute profiles for both cases are fairly
symmetrical about y = 0.
The second example involves solute distributions in the xy-plane at t = 0.5
and z = 0 resulting from a uniform initial solute distribution in the form of
two parallelepipeds. The initial distribution is sketched in Fig. 4(a). The twodimensional profiles for C,, C2, and Cr are sketched in Figs. 4(b), (c) and (d),respectively, using ,D = 0.5 and w = 1. The concentration in the equilibriumphase is largely determined by advection and dispersion; the initially distinct
solute peaks have disappeared as a result of dispersion (Fig. 4(b)). Because no
advective or dispersive transport can occur in phase 2, the C2 profile still
resembles the initial profile at t = 0.5. The solute in the region x > 25, i.e.
near the toe of the C2 profile, arrived there by movement in phase 1 andtransfer between the two phases according to Eq. (3b). Exchange between
the two phases also causes a tail to develop in the equilibrium profile (Fig.
4b). The resulting total concentration, Cr, is plotted in Fig. 4(d). A similar C,
profile for the ‘equilibrium’ condition (,0 = 0.9, w = 5) is provided in Fig. 4(e).
Notice that the non-equilibrium profile is contained within a smaller region of
the soil and with higher maximum concentrations.
A third example is used to illustrate the effects of changes in the first-orderdecay term, p. Figure 5(a) shows how the solute is initially distributed between
depths x = 5 and 15 in two parallelepipedal areas. The solute will move
downward as a result of the application of solute-free water at the surface.
Consider a hypothetical case where, for instance through bioremediation, one
hopes to decrease the concentration to sufficiently low levels before the solute
reaches an underlying aquifer. Figures 5(b)-(d) show the contours of the total
concentration, Cr, in the xy-plane at t = 1 for an increasing ~1 assuming equal
decay in the equilibrium and non-equilibrium phases. This rather simplistic
example clearly demonstrates decreasing solute concentrations in both directions as a result of decay. Because dispersion and diffusion are relatively
important compared with advection, part of the solute moves upward fromthe initial area toward the soil surface.
Fig. 5. (opposite) IVP involving a Heaviside profile of finite rectangular regions with Y = 5, R = q, = 1,
D, = 40, Dy = D, = 10, ~1 = ~2 = 0, w = 5, /‘J = 0.9. (a) Initial distribution with (fi,f2,h,f4) =(O, 2, l,O),(xl, x2, xl, x4) = (0, 5, 10, 15), (_I+. ~3) = (5, 2.5), and (~2, ~3) = (5, 2.5); (b) CT contours at z = 0 for
fi, = ~2 = 0; (c) CT contours at z = 0 for p, = ~2 = 0.5; (d) CT contours at z = 0 for p, = JL* = 1.0.

225

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

0

5I

I
15

I
10

I

20

I

25

X
Fig. 7. (Continued)

A more detailed picture of the transverse solute profile as a result of nonequilibrium conditions is given in Fig. 6. This example illustrates a BVP
involving the steady application of a solute to a circular region of the soil
surface. Figures 6(a) and (b) give the Ct and C2 profiles, respectively, at
x = 10 and t = 1. The highest concentration in the equilibrium phase in
Fig. 6(a) occurs when non-equilibrium conditions dominate (,D = 0.1 andw = 0.1). Note again that although these concentrations will occur in theeffluent, concentrations in the soil expressed as solute mass per total liquidor bulk volume are considerably lower. Figure 6(b) shows relatively highvalues for C2 during minor non-equilibrium conditions; in this case the solutetransfers relatively easy from the equilibrium to the non-equilibrium phasebecause of the increased value of the transfer coefficient and a large proportion of the equilibrium phase. Figures 6(c) and (d) show transverse solute
profiles for t = 5. The C1 profiles are virtually identical and approach asteady-state profile for all three degrees of non-equilibrium. On the other
hand, significant differences still exist between the C, profiles, with the nonequilibrium profile for ,L3 = 0.1 and w = 0.1 still slowly increasing. The difference in transverse profiles as a result of non-equilibrium is very similar as forthe longitudinal profiles already investigated for one-dimensional transport(e.g. Toride et al., 1993). Eventually, all equilibrium and non-equilibrium
profiles will become identical.

226

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

The last example deals with different Heaviside production profiles for theequilibrium and non-equilibrium phases. Figure 7(a) shows the assumed
production profiles X,(x) and X,(x) in a semi-infinite rectangular area(- cc < y < 0 and - cc < z < 0). The equilibrium and non-equilibrium concentrations at y = z = - 10 for t = 1 and 10 are given in Figs. 7(b) and (c),respectively. The non-equilibrium C, profile resembles the corresponding X2
distribution much more closely than the C, profile follows the X1 distribution.
The explanation for this behavior is similar to that for the distributions in Fig.
4. The phase 2 solute profile is considerably less smooth than the phase 1
curve. Furthermore, the concentration below the production zone decreases
as a result of dispersion. These observations hold true for the initial period(t = 1) as well as for t = 10, when a steady-state situation has been establishedin which the production and advection-dispersion processes balance eachother.

Summary and conclusions
General analytical solutions were derived for three-dimensional non-equilibrium advectiveedispersive transport during one-dimensional steady flow ina semi-infinite soil system. The transport model accounts for two-site chemical
non-equilibrium or two-region physical non-equilibrium, and contains terms
for first-order degradation, and zero-order production. The analytical procedure consisted of taking Laplace transforms with respect to time and the semiinfinite longitudinal spatial coordinate, and Fourier transforms for the infinite
transverse directions. The inverse transformations were carried out according
to procedures previously reported for one-dimensional non-equilibrium andthree-dimensional equilibrium transport. The resulting general solutions of
the BVP, IVP and PVP, apply to flux- and volume-averaged concentrations
for phases 1 and 2. Solutions for a Dirac input or initial distribution could bereadily derived from the general solutions. Other useful expressions were
derived by specifying the transverse geometry where the solute was applied(BVP), initially distributed (IVP), or produced (PVP), namely: (1) a semiinfinite rectangle; (2) a finite rectangle; (3) a circle; (4) a sphere. Time-dependent solute input for the BVP was modeled as a series of step or Heavisidefunctions and an exponential function. The initial and production profile
distribution were also specified vs. depth as Heaviside or exponentialfunctions. For the Heaviside function this led to an initial or production
domain that, corresponding to the previous four transverse geometries,
consisted of a series of: (1) semi-infinite parallelepipeds; (2) finite parallelepipeds; (3) cylinders; (4) concentric shells. All specific analytical solutions wereF.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

227

presented in tabular form. Several graphical illustrations of concentration
profiles and breakthrough curves were given for results obtained withsolutions to the BVP, IVP and PVP. These results are generally straightforward extensions of previous findings for one-dimensional non-equilibriumtransport.
References
Cameron, D.A. and Klute, A., 1977. Convectiveedispersive solute transport with a combined
equilibrium and kinetic adsorption model. Water Resour. Res., 19: 7188724.
Carnahan, B., Luther, H.A. and Wilkes, J.O., 1969. Applied Numerical Methods. John Wiley,
New York, 604 pp.
Coats, K.H. and Smith, B.D., 1964. Dead-end pore volume and dispersion in porous media.
Soc. Petrol. Eng. J., 4: 73-84.
De Smedt, F. and Wierenga, P.J., 1979a. A generalized solution for solute flow in soils withmobile and immobile water. Water Resour. Res., 15: 1137-l 141.
De Smedt, F. and Wierenga, P.J., 1 9 7 9 b . Mass transfer in porous media with immobile water. J.
Hydrol., 41: 59967.
Gautschi, W., 1965. Error function and Fresnel integrals. In: M. Abramowitz and I.A. Stegun(Editors), Handbook of Mathematical Functions. Dover Publishing Co., New York, pp.
295-329.
Goldstein, F.R.S., 1953. On the mathematics of exchange processes in fixed columns. 1. Mathematical solutions and asymptotic expansions. Proc. R. Soc. London., 219: 15 I ~ 185.
Goltz, M.N. and Roberts, P.V., 1986. Three-dimensional solutions for solute transport in an
infinite medium with mobile and immobile zones. Water Resour. Res., 22: 1139- 1148.
Gradshteyn, I.S. and Ryzhik, I.M., 1980. Table of Integrals, Series, and Products. Academic
Press, San Diego, 1160 pp.
Kreft, A. and Zuber, A., 1978. On the physical meaning of the dispersion equation and its
solutions for different initial and boundary conditions. Chem. Eng. Sci., 33: 1471-1480.
Lassey, K.R., 1982. On the computation of certain integrals containing the modified Bessel
function Zs(<). Math. Comput., 39: 6255637.
Lassey, K.R., 1988. Unidimensional solute transport incorporating equilibrium and rate-limited isotherms with first-order loss. 1. Model conceptualizations and analytic solutions.
Water Resour. Res., 3: 343-350.
Leij, F.J., Skaggs, T.H. and van Genuchten, M.Th., 1991. An analytical solutions for solutetransport in three-dimensional semi-infinite porous media. Water Resour. Res., 27: 27192733.
Lindstrom, F.T. and Boersma, L., 1989. Analytical solutions for convectiveedispersive transport in confined aquifers with different initial and boundary conditions. Water Resour. Res.,
25: 241-255.
Lindstrom, F.T. and Narasimhan, M.N.L., 1973. Mathematical theory of a kinetic model for
dispersion of previously distributed chemicals in a sorbing porous medium. SIAM J. Appl.
Math., 24: 4966510.
Luke, Y.L., 1962. Integrals of Bessel Functions. McGraw-Hill, New York.
Nkedi-Kizza, P., Biggar, J.W., Selim, H.M., van Genuchten, M.Th., Wierenga, P.J., Davidson,
J.M. and Nielsen, D.R., 1984. On the equivalence of two conceptual models for describing

228

F.J. Leij et al. / Journal of Hydrology 151 (1993) 193-228

ion exchange during transport through an aggregated oxisol. Water Resour. Res., 20: 11233
1130.
Olver, F.W.J., 1965. Bessel functions of integer order. In: M. Abramowitz and I.A. Stegun(Editors), Handbook of Mathematical Functions. Dover Publishing Co., New York, pp.
3555433.
Parker, J.C. and van Genuchten, M.Th., 1984. Flux-averaged and volume-averaged concentrations in continuum approaches to solute transport. Water Resour. Res., 20: 866-872.
Selim, H.M., Davidson, J.M. and Mansell, R.S., 1976. Evaluation of a two-site adsorptiondesorption model for describing solute transport in soils. In: Proc. Summer Computer
Simulation Conf. Am. Inst. Chem. Eng., Washington, DC., pp. 444-448.
Sneddon, I.H., 1972. The Use of Integral Transforms. McGraw-Hill, New York.
Spiegel, M.R., 1965. Theory and Problems of Laplace Transforms. Schaum’s Outline Series,
McGraw-Hill, New York, 261 pp.
Toride, N., Leij, F.J. and van Genuchten, M.Th., 1993. Analytical solutions for nonequilibrium
solute transport with first-order decay and zero-order production. Water Resour. Res., 29:216772182.
Van Genuchten, M.Th., 198 1. Non-equilibrium transport parameters from miscible displacement experiments. Res. Rep. No. 119, US Salinity Laboratory, USDA, ARS, Riverside,CA, 88 pp.
Van Genuchten, M.Th. and Alves, W.J., 1982. Analytical solutions of the one-dimensional
convectiveedispersive solute transport equation. USDA Tech. Bull. 1661, U.S. Gov. Print
Office, Washington, DC., 149 pp.
Van Genuchten, M.Th. and Wagenet, R.J., 1989. Two-site/two-region models for pesticidetransport and degradation: theoretical development and analytical solutions. Soil Sci. Soc.
Am. J., 53: 1303-1310.
Van Genuchten, M.Th. and Wierenga, P.J., 1976. Mass transfer studies in sorbing porous
media. I. Analytical solutions. Soil Sci. Soc. Am. J., 40: 473-481.
Walker, G.R., 1987. Solution to a class of coupled linear partial differential equations. IMA J.
Appl. Math., 38: 35-48.
Wexler, E.J., 1989. Analytical solutions for one-, two- and three-dimensional solute transport
in ground-water systems with uniform flow. U.S. Geol. Surv., Open-file Rep. 89-56, 250 pp.