INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 1623—1637 (1997)

A THIRD-ORDER NUMERICAL SCHEME WITH UPWIND
WEIGHTING FOR SOLVING THE SOLUTE
TRANSPORT EQUATIONKANGLE HUANG, JIR[ Ï S[ IMUs NEK AND M. TH. VAN GENUCHTENº.S. Salinity ¸aboratory, ºSDA, ARS Riverside, CA, º.S.A.

SUMMARY
Solute transport in the subsurface is generally described quantitatively with the convection—dispersion
transport equation. Accurate numerical solutions of this equation are important to ensure physicallyrealistic predictions of contaminant transport in a variety of applications. An accurate third-order in time
numerical approximation of the solute transport equation was derived. The approach leads to corrections
for both the dispersion coefficient and the convective velocity when used in numerical solutions of thetransport equation. The developed algorithm is an extension of previous work to solute transport conditions
involving transient variably saturated fluid flow and non-linear adsorption. The third-order algorithm isshown to yield very accurately solutions near sharp concentration fronts, thereby showing its ability toeliminate numerical dispersion. However, the scheme does suffer from numerical oscillations. The oscillations could be avoided by employing upwind weighting techniques in the numerical scheme. Solutions
obtained with the proposed method were free of numerical oscillations and exhibited negligible numericaldispersion. Results for several examples, including those involving highly non-linear sorption and infiltration into initially dry soils, were found to be very accurate when compared to other solutions. ( 1997 by
John Wiley & Sons, Ltd.
KEY WORDS: third-order; numerical method; transport; upwind weighting

INTRODUCTION
The quality of groundwater depends upon the quality of water infiltrating the vadose zone andthe ensuing physical, chemical and biological reactions and transformations involving this waterand its dissolved solutes as they move through the unsaturated zone toward the groundwatertable. Solute transport in variably saturated soils and groundwater is generally described using
convection—dispersion-type transport equations. Much attention has been devoted to the development of accurate numerical solutions of these equations. Especially challenging has been theformulation of numerical solutions that are free of oscillations and artificial dispersion. The mostpopular solution techniques remain those based upon finite differences and finite elements.1,2 Of
these two approaches, finite element methods have become the more popular ones because of
their relative ease of considering irregular flow domains and associated boundary conditions.
Whether based upon finite differences or finite elements, typical numerical schemes are theimplicit and Crank—Nicolson solution algorithms. The implicit scheme gives a first-order,
backward-in-time approximation of the time derivative, and as such may introduce considerable
numerical dispersion in computed concentration distributions, especially when relatively sharp
fronts are to be simulated. By contrast, the time-centred, second-order accurate Crank—Nicolson
CCC 0029—5981/97/091623—15$17.50( 1997 by John Wiley & Sons, Ltd.

Received 5 February 1996
Revised 9 June 1996

1624

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

scheme limits or may even eliminate numerical dispersion, but often at the cost of significant
phase errors1 which can lead to physically unrealistic oscillations in the numerical solutions. Van
Genuchten and Gray3 previously developed a dispersion-corrected linear finite element scheme
for solving the transport equation. Their scheme is essentially a time-centred, Crank—Nicolsontype approximation, yet provides a third-order approximation of the time derivative through theuse of adjusted dispersion coefficients. As compared to the standard Crank—Nicolson scheme, thedispersion-corrected scheme suffers from less numerical oscillations while at the same timeyielding a much more precise description of sharp concentration fronts (i.e., producing lessnumerical dispersion).
A parallel and increasingly popular approach for limiting or avoiding numerical oscillations
has been the use of upstream weighting methods.4~13 By introducing corrections to the convective velocities and consequently adding an artificial dispersion to increase the stability, upstream
weighting ensures that the calculated concentrations remain physically realistic and free ofnumerical oscillations (i.e. do not overshoot or undershoot the correct solutions), regardless ofgrid size and time step.10, 14 Calculated non-physical values (e.g. negative concentrations) sometimes can cause the computations to abort, for example, when non-linear or multicomponent
solute reactions are involved.15, 16 While upstream weighting may be effective in eliminating
oscillatory behaviour, the resulting scheme is also known to lead to numerical dispersion for mostpractical discretizations (mesh sizes).
The objective of this paper is to develop a numerical solution method for the variably saturated
solute transport equation which is not only free of numerical oscillations, but also produces no or
limited numerical dispersion. The objective is achieved by first extending the dispersion-correctednumerical scheme of van Genuchten and Gray3 to transport under transient variably saturated
flow conditions and non-linear adsorption. Subsequently, the derived third-order approximationof the transport equation will be solved using the upwind Petrov/Galerkin finite element method.
FORMULATIONThe one-dimensional transport of miscible components in a variably saturated medium isgenerally described with the convection-dispersion equation as follows:

A B

Lhc Los
L
Lc
Lqc
# "
hD
! !k hc!k os#c h#c o!Sc
8
4
8
4
4
Lt
Lt
Lx
Lx
Lx

(1)

where c is the solute concentration in solution [ML~3], s is the adsorbed concentration[MM~1], h is the volumetric water content [L3L~3], o is the soil bulk density [ML~3], D is thedispersion coefficient [L2T~1], q is the Darcian fluid flux density [LT~1], k and k are
8
4
first-order decay coefficients [T~1] for the liquid and solid phases, respectively; c and c are
84
zero-order rate constants for the liquid [ML~3 T~1] and solid [T~1] phases, respectively; S isa sink term in the water flow equation (e.g. to account for root water uptake) [T~1], c is the4
solute concentration of the sink term [ML~3], x is distance from the soil surface downward [L],and t is time [T]. The Darcian fluid flux density, q, in (1) is defined as
q"!K

Lh
#K
Lx

(2)

in which K is the hydraulic conductivity [LT~1]. For transport in the unsaturated zone, q, K and
h are generally functions of the soil—water pressure head, h [L], which is obtained by solving the
INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)( 1997 by John Wiley & Sons, Ltd.

1625

A THIRD-ORDER NUMERICAL SCHEME

variably saturated flow equation:

A

B

Lh
L
Lh
"
K !K !S
Lt Lx
Lx

(3)

The influence of sorption on contaminant transport is potentially important since sorption willaffect the contaminant transport velocity in the subsurface. For relatively simple systems, therelationship between the adsorbed concentration, s, and the solution concentration, c, can be
described with a sorption isotherm. In this study we assume that the adsorbed and solution
concentrations are always in local equilibrium. A common expression for the equilibrium
adsorption isotherm is the non-linear Freundlich equation given by
s"kcg

(4)

in which k and g are empirical coefficients. Note that when the exponent g equals 1·0 theFreundlich equation is linear, in which case the parameter k reduces to the familiar solid-liquid
phase partitioning coefficient, k [L3M~1]. When the exponent g is different from 1·0, the$Freundlich isotherm, and hence also the solute transport equation itself, become non-linear. In
that case an iterative solution scheme, similar to the one applied to the variably saturated flow
equation, must be used to solve the transport equation.
Substitution of (4) into (1) leads to the final form of the solute transport equation considered inthis paper:

A B

LhRc
L
Lc
Lqc
"
hD
! !khc#!
Lt
Lx
Lx
Lx

(5)

okcg~1
k"k #
k
8
4
h

(6)

okcg~1
R"1#
h

(7)

where

and
!"c h#c o!Sc
(8)
8
4
4
Note that the retardation factor R as defined by (7) is different from its conventional definition(e.g. Reference 17) in which the Freundlich exponent g appears as an additional multiplier in thelast term of (7). Both k and R are non-linear coefficients of the concentration c, unless g"1.
THIRD-ORDER APPROXIMATION IN TIME
The dispersion-corrected numerical solution of (5) is based on the following third-order accurate
finite difference equation in time:3

C

D

C

D

f t`*t!f t 1 Lf *t L2 f t`*t 1 Lf *t L2 f t
#
"
!
#
2 Lt
2 Lt
6 Lt2
6 Lt2
*t

(9)

where the superscripts t and t#*t indicate evaluation at the old and new time levels in the
numerical solution, and where *t represents the time increment used in the calculations. The( 1997 by John Wiley & Sons, Ltd.

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

1626

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

above equation holds for any third-degree polynomial, f, in time. Using (hRc) for f in equation (9)
gives

C

D

C

D

(hRc)t`*t!(hRc)t 1 L(hRc) *t L2(hRc) t`*t 1 L(hRc) *t L2 (hRc) t
#
"
!
#
Lt
2
Lt
*t
2
6 Lt2
6 Lt2

(10)

The left-hand side of this equation will be associated with the first-order time derivative of thegoverning transport equation (4). An approximation for the second-order derivative in (10) may
be obtained using the rule of differentiationCL2 (hRc) L L (hRc)
"
Lt2
Lt
Lt

D

(11)

To simplify (11), we assume momentarily that dispersion and other reaction terms in thetransport equation can be neglected as high-order effects, i.e.,L(hRc)
L(qc)
+!
Lt
Lx

(12)

Incorporating (12) into (11) leads to

C D

L2(hRc)
L L(qc)
+!
Lt2
Lt Lx

(13)

Reversing the order of differentiation of the right-hand side of (13) and expanding the resulting
time derivative givesA

L2 (hRc)
L
Lc
Lq
+!
q #c
Lt2
Lx
Lt
Lt

B

(14)

The second term on the right-hand side of (14) may be neglected when the water flow regime issteady-state, or when temporal changes in the flux, q, are relatively small. In general, however, thesecond term of the right-hand side of (14) should be considered, as is done in our study.
Expanding both sides of (12) and incorporating (7) leads to
hR

A

B A

Lc
Lh
Lcg~1
Lc
Lq
#c
#ok
+! q #c
Lt
Lt
Lt
Lx
Lx

B

(15)

Next, consider the continuity equation for water flow, i.e.,
Lh
Lq
"! !S
Lt
Lx

(16)

in which S is the source/sink term as defined earlier. Substituting (16) into (15) gives
hR

Lc
Lcg~1
Lc
!cS#cok
+!q
Lt
Lt
Lx

(17)

Differentiating the third term of (17), and rearranging the resulting equation, leads to
q

A

q
Lc
Lc
"
cS!q
Lx
Lt hR/

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

B(18)( 1997 by John Wiley & Sons, Ltd.

A THIRD-ORDER NUMERICAL SCHEME

1627

ok (g!1) cg~1
/"1#
hR

(19)

where

or with (7)A

1/"1#(g!1) 1!
R

B

(20)

Incorporating (18) into (14) leads to the following approximation for the second-order time
derivative:

C A

B

L2 (hRc)
q
L
Lc
Lq
+!
cS!q
#c
Lt2
Lx hR/
Lx
Lt

D

(21)

Finally, substituting equations (5) and (21) into (10) yields the following third-order finite
difference approximation of the transport equation

C A

B

D

Lc
t`*t
(hRc)t`*t!(hRc)t 1 L"
hD~ !q~c !khc#!
Lx
*t
2 Lx
(22)

C A

B

D

Lc
t
1 L
#
hD` !q` c !khc#!
Lx
2 Lx
where
q2*t
D~"D!
,
6Rh2/

q2*t
D`"D#
6Rh2/

(23)

and
*tqS *tLq*tqS *tLq
q~"q!
!
, q`"q#
#
6Rh/
6 Lt
6Rh/
6 Lt

(24)

One may verify that the dispersion corrections given by (23) reduce to those derived by van
Genuchten and Gray3 when /"1, i.e., for the case of linear sorption (g"1) or no adsorption atall (R"1). Equation (24) shows that now also corrections are implemented for the Darcian fluid
flux density, q. This result is different from, and more accurate than, the approach originally usedby van Genuchten18 for solute transport under transient variably saturated water flow condition.
The derivation above holds for Freundlich-type adsorption. One may obtain similar higherorder numerical schemes for any arbitrary non-linear isotherm of the general forms"F (c)(25)

The retardation factor to be used in (22)—(24) becomes theno F (c)R"1#h c( 1997 by John Wiley & Sons, Ltd.

(26)

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

1628

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

The corrected equations for D and q remain the same as (23) and (24) above, provided that / isnow defined asCo LF F (c)
!
/"1#
c
hR Lc

D

(27)

The above derivations are also easily extended to multidimensional transport problems
described by
LhRc"+ · (hD+c!qc)!khc#!
Lt
1

(28)

where q is the water flux density vector, and the dispersion coefficient tensor:
D

xx

D" D
yx
1
D
zx

D
xy
D
yy
D
zy

D
xz
D
yz
D
zz

(29)

The corresponding three-dimensional variably saturated flow equation is given byLh"!$ · q!S
Lt

(30)

The approximation for the second-order time derivative in (10) is now of the formCq
L2(hRc)Lq"!+ ·(cS!q · $c)#c
hR/
Lt2
Lt

D

(31)

Substituting (28) and (31) into (10) results in the following finite difference approximation of thethree-dimensional solute transport equation:(hRc)t`*t!(hRc)t 1" [+ · (hD~+c!q~c)!khc#!]t`*t*t2
11# [$ · (hD`+c!q`c)!khc#!]t
2
1

(32)

where
D

xx

D " D
yx
1
D
zx
$

D
D
q2 q q q qx
x y x zxyxz*tqq
q2 q q
D
D $
y
y z
yy
yz
6h2R/ y x
qq qq
q2
D
D
z
z x
z y
zy
zz

(33)

and

G HA

GH

q
q
x
x
*t
*t L
q " q
1$
S $
(34)
q
y
6hR/
6 Lt yqqzz
in which the coefficient / is the same as before, i.e., equation (20) for a Freundlich isotherm, or (27)for an arbitrary isotherm F (c).
$INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

B( 1997 by John Wiley & Sons, Ltd.

1629

A THIRD-ORDER NUMERICAL SCHEME

UPWIND/PETROV—GALERKIN FORMULATION
To avoid numerical oscillation, the third-order finite difference approximation of the transport
equation (22) is discretized by using the upwind-weighting Galerkin finite element method as
follows:

P AC A

B
B

D
DB
A BD

Lc
l (hRc) t`*t!(hRc) t
l1
L
t`*t
/ (x) dx"hD~
!khc#!
i
*t
Lx
2 Lx
0
0
Lc
t
1 L
hD`
!khc#!
/ (x) dx
#
i
Lx
2 Lx

P

C A
PCA B

!

l

(35)

1 Lq~c t`*t 1 Lq`c t#w (x) dx#*t
i
2 Lx
2 Lx

0
where l is the length of the flow domain and is assumed to be discretized into (N!1) elementswith increment *x "(x !x ), / are the standard linear basis functions, and w are asymmeti
i
i~1 i
i
ric weighting functions defined byw (x)"/ (m)!3a/ (m) / (m)
1
1
1
2

(36)
w (x)"/ (m)#3a/ (m) / (m)21
12
in which the subscripts 1 and 2 represent the end nodes of each finite element, and m is the local
co-ordinate defined in terms of the global co-ordinate, x, as2(x!x )
1 !1, x3[x , x ]
m"
1 2
x !x
2
1
The weighting factor a in (36) is optimally chosen as4
a"coth

A B

Pe
2
!
2
Pe

(37)

(38)

where Pe is the grid Peclet number. After substituting the usual finite element approximation forc in (35) and performing the integrations, the following final matrix equation is obtained:Act`*t"Bct#F(39)

where A and B are both tridiagonal matrices, and c and F are concentration and flux vectors,
respectively. The below-diagonal (b), diagonal (d), and above-diagonal (a) entries of A are given
by
1
a
1
b "!
(h D~ #h D~ )! (2q~ #q~ )! q~ , i"2, N
i
i i
i~1
i
i~1 i~1
12
4 i~1
4*x
i~1
1
1
d"
(h D~ #h D~ )#(h D~#h D~ )i i
i`1 i`1i 4*x
i~1 i~1
4*x i ii
i~1*x*x
# i~1 (h º #2h º )# i (2h º #h º )i~1 i~1i i
i ii`1 i`1
6*t
6*t

(40)

(41)

a
1
# (q~ !q~ )# q~
i~1
i`1
2 i
12
( 1997 by John Wiley & Sons, Ltd.

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

1630

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

1
1
a
a "!
(h D~#h D~ )# (q~#2q~ )! q~ , i"1, N!1
i
i
i
i`1
i`1
i
i`1
4*x
12
4 i`1
i

(42)

k *t
º "R # i
i
i
2

(43)

where

The entries of matrix B are the same as those of A except that D~ is replaced by !D`, q~ by!q` a by !a, and k by !k. Finally, the entries of vector F are given by*x
f " 1 (2! #! )#J
1
1
2
0
6
1
f " [*x (2! #! )#*x (2! #! )] ,
i 6
i~1
i
i~1
i
i
i`1
*x
f " N~1 (2! #!
)!J
N
N
N~1
l
6

(44)

2)i)N!1

(45)

(46)

in which J and J are solute fluxes across the inlet (x"0) and outlet (x"l) boundaries,0lrespectively. The resulting system of equations is solved using the Thomas’ algorithm, yielding at
each time step the desired concentration vector ct`*t .
NUMERICAL TESTSThis section presents results of several numerical tests to illustrate the performance of the
proposed numerical algorithm and its relative accuracy. The algorithm was implemented in thecomputer code HYDRUS (version 5·0)19 which is a revised and updated version of the original
HYDRUS code developed by Kool and van Genuchten.20 The updated code simulates onedimensional water flow, solute transport and heat movement in a variably saturated porous
medium using linear finite elements.
All numerical examples were run assuming a solute flux boundary condition at the soil surface(x"0) and a zero concentration gradient condition at the bottom boundary (x"l). For
simplicity, each simulation assumed a uniform initial concentration profile. The proposed
numerical algorithm was first verified for an idealized solute transport scenario for which ananalytical solution is available. The algorithm was then applied to transport problems involving
transient variably saturated water flow and non-linear adsorption. Solutions obtained with theproposed method will be compared with the commonly used first-order implicit schemes. The
solution domain in all cases was assumed to be 200 cm in length and, except where indicated,discretized using uniform spatial increments, *x, of 2 cm. Molecular diffusion was assumed tobe negligible during the calculation, giving the following approximation for the dispersion
coefficient:
D"e

DqD
h

(47)

where e is the dispersivity (L) of the medium.
INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)( 1997 by John Wiley & Sons, Ltd.

A THIRD-ORDER NUMERICAL SCHEME

1631

Steady-state flow
The first test problem deals with solute transport in a medium having a constant convective
velocity (l"q/h) of 4 cm/d, and a relatively low dispersivity (e) of 0·04 cm. No solute sorption,i.e., s"0, was considered in this example. A tracer solution having a relative concentration, C , of0
1· 0 was applied to an initially solute-free medium for a period of 20 d, after which solute-freewater was again added. The grid Peclet number, defined as Pe"l *x/D"*x/e, was equal to 25,
which represents a case of convection-dominated transport. A constant time step *t"0·1 day
was used.
Figure 1 shows that the solutions obtained with the proposed method (Central-UW) using*x"1 cm is completely free of numerical oscillations. The implicit methods, with and without
upwind weighting, also produced stable solutions. The implicit algorithm without upwind
weighting suffered from the least numerical dispersion for this example, while the proposed
method caused less numerical dispersion than the upwind weighted implicit scheme. Notice thatthe dispersion-corrected algorithm without upwind weighting predicted the concentration fronts
accurately, albeit with some serious oscillations. When the Peclet number was increased to
100 (*x"2 cm and e"0·02 cm), the implicit method without upwind weighting became unstable (Figure 2), while the proposed method still produced stable solutions (see Figure 2).

Figure 1. Concentration distributions calculated with the proposed method (Central-UW), the third-order algorithm(Central), the implicit algorithm (Implicit), the upwind weighted algorithm (Implicit-UW), and the exact solution (steadystate flow; Pe"25)( 1997 by John Wiley & Sons, Ltd.

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

1632

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

Figure 2. Comparison of calculated concentration distributions obtained with the proposed method and several othernumerical algorithms (steady-state flow; Pe"100)Moreover, the numerical dispersion induced by the proposed scheme is again somewhat smaller
than that generated with the upwind-weighted implicit algorithm.
Steady-state flow and nonlinear adsorption
The second example deals with the transport of a solute subject to non-linear adsorption. We
used for this example the Freundlich isotherm (4) with coefficient values k"0·5 and g"1·5, thus
representing a rather strongly non-linear case. The flow and transport conditions were exactly thesame as above example. Because of the presence of a non-linear isotherm, the scheme required
numerical iterations. Convergence was assumed to be reached when the absolute error in
calculated concentrations between two consecutive iterations became less than 0·0001. Varying
time steps were employed in this and the following examples. The simulation started by using an
initial time step, *t "10~5 d. The ith time step *t 3 [*t , *t ] was automatically set
0
i
.*/
.!9
according to the iteration convergency history at the previous time step, that is, if the iteration for
INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)( 1997 by John Wiley & Sons, Ltd.

A THIRD-ORDER NUMERICAL SCHEME

1633

Figure 3. Concentration distributions obtained with the proposed method and the implicit scheme (steady-state flow andnon-linear adsorption; Pe"20).

the previous time step converged slowly, the current time step is reduced to be smaller than theprevious one; otherwise the current time step is stretched to be larger than the previous. The
minimum time step *t "10~6 d and the maximum one *t "0·1 day were arbitrarily.*/.!9chosen in our tests. Shown in Figure 3 are calculated concentration distributions for Pe"20(*x"2 cm and e"0·1 cm). Numerical solutions obtained with the proposed (Central-UW)
scheme were completely free of numerical oscillations, even near the concentration peaks
which generally are the most vulnerable to non-physical oscillations. By contrast, solutions
obtained with the widely-used backward implicit finite element scheme suffered from serious
numerical oscillations. Considerable overshoot and especially undershoot were produced with
the implicit algorithm for Pe"100 (e"0·02 cm), while the proposed method still worked very
well (Figure 4).
Transient flow and non-linear adsorption
The above comparisons pertain to transport during steady-state water flow. Simulations werealso carried out for solute transport under highly transient flow conditions. We considered( 1997 by John Wiley & Sons, Ltd.

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

1634

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

Figure 4. Simulated concentration distributions obtained with the proposed method (Central-UW) and the implicit
algorithm (steady-state flow and non-linear adsorption; Pe"100)Figure 5. Calculated pressure head distributions for ponded infiltration in a loamy soil
INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)( 1997 by John Wiley & Sons, Ltd.

A THIRD-ORDER NUMERICAL SCHEME

1635

ponded infiltration into an initially relatively dry (!103 cm of water pressure head) loamy
soil. The soil-hydraulic properties were described with the closed-form functions of van
Genuchten21 assuming hydraulic parameters a"0·015 (1/cm), n"2·0, h "0·45, h "0·05,
4
3
and K "2 cm/d. An initially square concentration distribution (Figure 6a, t"0) was assumed
4in this example. The solute was again subjected to Freundlich-type non-linear equilibrium
adsorption (k"0·5 and g"1·5). The assumed transport experiment started by applying
solute-free irrigation water (ponding) to the soil surface for a period of 20 days. This leachingprocess was followed by adding water having a solute concentration of 1·0 to the soil profile.
Correct simulation of this example by means of standard numerical methods can be very difficult
because of the highly non-linear nature of both the flow and transport process. For each time
step, accurate solutions for the pressure head, h, and the water content, h, were first obtainedby solving the mixed-form Richards equation using the modified Picard iteration algorithm22
coupled with a newly formulated convergence criterion.23 Procedures of numerical iteration
for the non-linear transport solution are the same as the previous example. Figure 5 shows
calculated pressure head distributions during the leaching process. Notice the very sharp
distributions of the moisture front. Calculated flux density q by Darcy’s law (2) from the h-profilesin Figure 5 were subsequently used in numerical transport algorithm to obtain solutions of soluteconcentration.

Figure 6. Calculated concentration distributions obtained with the proposed method (Central-UW) for solute transport
during transient infiltration and involving non-linear adsorption (Pe"20)( 1997 by John Wiley & Sons, Ltd.

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)

1636

K. HUANG, J. S[ IMUs NEK AND M. TH. VAN GENUCHTEN

Figure 7. Calculated concentration profiles obtained with the proposed method for transport during transient flow andinvolving nonlinear adsorption (Pe"100)Figure 6 presents predicted concentration distributions for Pe"20 (*x"1 cm, e"0·05 cm).
The results show that the proposed Central-UW method produced very stable solutions, i.e.,being oscillation-free and having very little if any numerical dispersion. Notice that the concentration fronts are very steep, which is typical for cases where the numerical solutions are mostvulnerable to both numerical oscillations and artificial dispersion. This example was purposely
designed so as to test the performance of the proposed method in situations where multiple solutepeaks are present. A relatively poor performance of a numerical scheme in such situations often
leads to peak clipping (yielding lower concentration peaks) and valley elevation (producing
relatively higher minimum concentrations) as compared to the correct solutions.24 Figure 7
shows that numerical solutions for Pe"100 (e"0·01 cm) also produced excellent results.
Finally, we note that for all of the simulations discussed above we obtained extremely good massconservation, with mass balance errors always being less than 0·1 per cent. Although a good massbalance per se does not necessarily mean that the predicted solutions are accurate, a correct massbalance is at least a requisite for an accurate numerical solution.
CONCLUSIONS
An accurate third-order approximation of the solute transport equation was derived. The scheme
leads to corrections for both the dispersion coefficient and the convective velocity in the solutetransport equation. The algorithm is an extension of previous work by van Genuchten and Gray3
and van Genuchten18 to conditions involving transient fluid flow and non-linear adsorption. The
third-order algorithm accurately describes concentration distributions near sharp concentration
fronts, thereby showing its ability to eliminate artificial dispersion. Numerical oscillations in thealgorithm could be avoided by employing upwind weighting. Solutions obtained with theproposed method were completely free of numerical oscillations and exhibited negligible numerical dispersion. Results for several examples involving non-linear sorption and infiltration intoa very dry soil, were found to be very accurate when compared to other solution methods.
REFERENCES
1. G. F. Pinder and W. G. Gray, Finite Element Simulation in Surface and Subsurface Hydrology, Academic Press,
New York, 1977.
2. L. Lapidus and G. F. Pinder, Numerical Solution of Partial Differential Equations in Science and Engineering, John
Wiley & Sons, New York, 1982.
INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)( 1997 by John Wiley & Sons, Ltd.

A THIRD-ORDER NUMERICAL SCHEME

1637

3. M. Th. van Genuchten and W. G. Gray, ‘Analysis of some dispersion corrected numerical schemes for solution of thetransport equation’, Int. j. numer. methods eng., 12, 387—404 (1978).
4. J. C. Heinrich, P. S. Huyakorn, O. C. Zienkiewicz and A. R. Mitchell, ‘An ‘upwind’ finite element scheme fortwo-dimensional convective-transport equations’, Int. j. numer. methods eng., 11, 131—143 (1977).
5. B. P. Leonard, ‘A survey of finite differences of opinion on numerical modeling of the incomprehensible defective
confusion equation’, Finite Element Methods for Convection Dominated Flows, Annual Meeting of the ASCE,New York, 1979, pp. 1—17.
6. P. S. Huyakorn and K. Nilkuha, ‘Solution of transient transport equation using an upstream finite element scheme’,Appl. Math. Modelling, 1, 187—195 (1979).
7. A. N. Brooks and T. J. R. Hughes, ‘Streamline upwind/Petrov-Galerkin formulations for convection dominated flow
with particular emphasis on the impressible Navier—Stokes equations’, Comput. Methods Appl. Mech. Eng., 32,
199—259 (1982).
8. N.-Z. Sun and W. W.-G. Yeh, ‘A proposed upstream weight numerical method for simulating pollutant transport ingroundwater’, ¼ater Resource Res., 19, 1489—1500 (1983).
9. A. Mizukami and T. J. R. Hughes, ‘A Petrov—Galerkin finite element method for convection dominated flows: anaccurate upwinding technique for satisfying the maximum principle’, Comput. Methods Appl. Mech. Eng., 50 (1985).
10. P. H. Sammon, ‘An analysis of upstream weighting’, Soc. Petrol. Eng. J. Res. Eng., 3, 1053—1056 (1988).
11. J. J. Westerink and D. Shea, ‘Consistent higher degree Petrov—Galerkin methods for the solution of the transientconvection—diffusion equation’, Int. j. numer. methods eng., 28, 1077—1101 (1989).
12. M. Putti, W. W.-G. Yeh and W. A. Mulder, ‘A triangular finite volume approach with high-resolution upwind terms
for the solution of groundwater transport equations’, ¼ater Resource Res., 26, 2865—2880 (1990).
13. C. T. Miller and F. H. Cornew, ‘A Petrov—Galerkin method for resolving advection-dominated transport’, in T. F.
Russell, R. E. Ewing, C. A. Brebbia, W. G. Gray, G. F. Pinder (eds.), Numerical Methods in ¼ater Resources, Vol. 1,
Proc. Computational Methods in ¼ater Resources IX, 1992, pp. 157—164.
14. F. W. Letniowski and P. A. Forsyth, ‘A control volume finite element method for three-dimensional NAPL
groundwater contamination’, Int. j. numer. methods fluids, 13, 955—970 (1991).
15. P. A. Forsyth, ‘Three-dimensional modeling of steam flush for DNAPL site remediation’, Int. j. numer. methods fluids,
19, 1055—1081 (1994).
16. J. S[ imuns ek and D. L. Suarez, ‘Modeling of carbon dioxide transport and production in soil, 1. Model development’,¼ater Resource Res., 29, 487—497 (1993).
17. M. Th. van Genuchten and R. W. Cleary, ‘Movement of solutes in soil: Computer-simulated and laboratory results’,
in G. H. Bolt (ed.), Soil Chemistry, B. Physico-Chemical Models, Elsevier, Amsterdam, 1979, pp. 349—386.
18. M. Th. van Genuchten, ‘Mass transport in saturated-unsaturated media: one-dimensional solutions’, Research Report
78-WR-11, Water Resources Program Princeton University, 1978.
19. T. Vogel, K. Huang, R. Zhang and M. Th. van Genuchten, ‘The HYDRUS (version 5·0) code for simulating
one-dimensional water flow, solute transport, and heat movement in variably-saturated media’, Research Report No.
140, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 1996.
20. J. B. Kool and M. Th. van Genuchten, ‘HYDRUS (Version 3.3): One-dimensional variably saturated flow andtransport model, including hysteresis and root water uptake’, Research Report No. 124, U.S. Salinity Laboratory,
USDA, ARS, Riverside, CA, 1991.
21. M. Th. van Genuchten, ‘A closed-form equation for predicting the hydraulic conductivity of unsaturated soils’, Soil
Sci. Soc. Am. J., 44, 892—898 (1980).
22. M. A. Celia, E. T. Bouloutas and R. L. Zarba, ‘A general mass-conservative numerical solution for the unsaturated
flow equation’, ¼ater Resource Res., 26, 1483—1496 (1990).
23. K. Huang, B. P. Mohanty and M. Th. van Genuchten, ‘A new convergence criterion for the modified Picard iteration
method to solve the variably-saturated flow equation’, J. Hydrol, 178, 69—96 (1996).
24. G. T. Yeh, J. R. Chang, H. P. Cheng and C. H. Sung, ‘An adaptive local grid refinement based on the exact peakcapture and oscillation free scheme to solve transport equations’, Comput. Fluids, 24, 293—332 (1995).

( 1997 by John Wiley & Sons, Ltd.

INT. J. NUMER. METHODS ENG., VOL. 40: 1623—1637 (1997)