WATER RESOURCES RESEARCH, VOL. 30, NO. 2, PAGES 499-507, FEBRUARY 1994

An Eulerian-Lagrangian approach with an adaptively corrected
method of characteristics to simulate
variably saturated water flow
K. Huang, R. Zhang,l and M. T. van Genuchten
U.S. Salinity Laboratory, Agricultural Research Service, USDA, Riverside, California

Abstract. A relatively simple method of characteristics is developed to simulate onedimensional variably saturated water flow. The method uses the Eulerian-Lagrangian
approach to separate the governing flow equation into “convection” and “diffusion”
parts, which are solved with the method of characteristics and the conventional finite
element method, respectively. The method of characteristics combines a single-stepreverse particle tracking technique with a correction strategy to ensure accurate massbalances. The correction process is implemented by weighing the calculated convective
contribution to the pressure head at each node with the pressure head values of two
upstream nodes, using an adaptive weighing factor A. The value of A is automaticallyadjusted by considering the global mass balance at each time step. Numericalexperiments for ponded infiltration are presented to illustrate the scheme’s performancefor situations involving highly nonlinear soil hydraulic properties and extremely dry
initial conditions. Results indicate that the proposed method is mass-conservative,
virtually oscillation-free, and computationally quite efficient. The method is especiallyeffective for simulating highly nonlinear flow scenarios for which traditional finite
difference and finite element numerical methods often fail to converge.

Introduction
Water flow in variably saturated media is an important
topic in several branches of hydrology, soil science, andagricultural engineering dealing with subsurface flow andchemical transport processes. Many authors have used avariety of analytical techniques to solve the governing
equations for infiltration and water flow in soils [e.g., Philip,
1957; Parlange, 1972; Broadbridge and White, 1988; Warrick
et al., 1991]. Analytical solutions provide useful tools forstudying relatively simple unsaturated flow problems such as
infiltration in homogeneous soil profiles. However, they areusually inappropriate for more complicated flow situations
because of their reliance on relatively simple initial andboundary conditions and the need to adopt simplified functions for the unsaturated soil hydraulic properties. Suchsituations generally require the flexibility of a numericalapproach.
Accurate implementation of numerical techniques often
encounters difficulties when nonlinear variably saturated
flow problems are simulated. Numerical problems are usually manifested by a lack of convergence and/or the presenceof undesired oscillations in computed results when sharp
moisture fronts are present. Such situations arise, for example, when calculating ponded infiltration in initially very dry
soils, particularly when the soil hydraulic properties arehighly nonlinear. Several methods have been proposed inefforts to reduce the computational effort and/or provide‘Now at Department of Plant, Soil, and Insect Sciences, University of Wyoming, Laramie.
This paper is not subject to U.S. copyright. Published in 1994 by theAmerican Geophysical Union.
Paper number 93 WR02881.

more stable numerical solutions. For example, Abriola[1986] used hierarchical basis functions to increase theorder, and hence accuracy, of the numerical scheme inelements near the wetting front. Dane and Mathis [1981]
presented an adaptive finite difference scheme that allows afine grid system to move with the wetting front. Morerecently, Gottardi and Venutelli [1992] developed a relatively sophisticated moving finite element model which uses
fewer nodes as compared to fixed-grid methods, and hence
can save CPU time when a relatively large flow domain is
involved. However, initial results indicate that this method
may yield larger mass balance errors than the conventional
finite difference and finite element methods. Gamliel and
Abriola [1992] proposed similar schemes which refine thenumerical mesh in areas where the soil water pressure headchanges rapidly but expand the mesh in less transientregions.
Numerical simulations based on the Richards’ equationfor variably saturated water flow generally use the pressure
head as the dependent variable. Many studies have observed
mass balance problems using this approach, presumably
because of nonlinearities in the soil water capacity term.
Milly [ 1985] proposed a mass-conservative numerical
scheme which greatly reduces the mass balance errors. He
evaluated the storage coefficient as an average value overeach element during each time step. Ross [1990] concluded
that relatively large time steps are possible when the massconservative mixed-form of the Richards’ equation is used.
Relatively large spatial increments were also possible whena hyperbolic sine transform was implemented. Hills et al.
[1989] described a water content based algorithm, whichappears to be superior to the traditional pressure-based
formulation for simulating unsaturated flow in layered soils.
Their scheme, however, could not be applied to conditions

499

Purchased by USDA for Official Use

500

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW

also involving saturated flow (i.e., positive pressure heads).
Celia and Pinder [1986], and more recently Celia et al.
[1990], proposed schemes which implement more accurate
and mass-conservative solutions of the mixed-form of the
Richards’ equation. In a different approach, Hansen et al.
[1992] combined the adaptive grid refinement method with
an operator splitting technique to solve the two-phase immiscible flow equations.
Most of the studies above used relatively standard finite
element (FE) or finite difference (FD) methods to simulate
the variably saturated flow process. While the studies have
led to improvements, conventional FE and FD methods
always suffer to some degree from inherent numerical difficulties, such as numerical oscillations, excessive numericaldispersion, and mass balance errors. Besides the presence ofhighly nonlinear hydraulic properties, numerical problems
are often also caused by the hyperbolic nature of thegoverning equation when the gravitational term in the Richards’ equation is relatively important [Sorek and Braester,
1986; Huang, 1988]. Using this last attribute, Sorek andBraester [ 1986] extended the adaptive Eulerian-Lagrangian
approach, initially introduced by Neuman [ 1984] for solving
convection-dispersion transport problems, to unsaturated
flow problems. While their method appears to improve
convergence, no successful simulation example was reported in the paper. Moreover, their method requires a
relatively sophisticated scheme to track the cluster of particles in a nonlinear flow field. After similarly separating thegoverning flow equation into convection and diffusion parts,
Huang [1988] used a different approach by implementing asingle-step reverse particle tracking (SRPT) technique tosolve the convection problem. Preliminary results indicatedseveral advantages of this approach, such as increased
computational efficiency and fewer numerical oscillations.
However, the scheme still exhibited relatively large massbalance errors, as well as some minor numerical oscillations,especially for highly nonlinear soil hydraulic properties.
The purpose of this paper is to extend and generalize theEulerian-Lagrangian solution scheme of Huang [1988]. A
numerically robust and computationally efficient method will
be presented which combines the modified method of characteristics, SRPT, with the conventional finite element approach. The scheme is further improved by incorporating aself-adaptive weighing procedure which virtually eliminatesmass balance errors.

Governing Equations
The governing equation for one-dimensional, vertical water flow in a variably saturated rigid porous medium is givenbyC;=; Kg-K -S(z, t)(1)i1where h is the pressure head, C = d0ldh is the soil watercapacity in which 0 is the volumetric water content, K is thehydraulic conductivity, S( Z, t) represents a source/sink term(e.g., root water uptake rate), z is soil depth assumed to bepositive downward, and t is time. The flow equation issolved subject to the general initial condition
h(z, 0) = hi(z)(2)

and either a first- or second-type boundary condition at thesoil surface (z, = 0), that is,h(O, t) = h,(t)

(3)

or
(4)
where h,(t) and q,,(t) are the prescribed pressure head andnet fluid flux, respectively. The auxiliary condition used here
for the lower boundary (Z = 1) is that of a free draining soil
profile as follows:
ah
-

az I

= 0

(5)

z=l

We further assume that the soil hydraulic properties, B(h)and K(B), can be described by the parametric functions of
van Genuchten [ 19801

e(h) = ‘r +(1 + (ah,“)!71

K(S,) = K,S;‘*[l - (1 - s;‘m)m]2

(6)
(7)

in which S, is the effective saturatione - Brs, = ~
es - or(8)

In these equations, or and es are the residual and saturated
water contents, respectively; K, is the saturated hydraulicconductivity; (Y and n are shape parameters; and m = 1 l/n.

Numerical Implementation
As was pointed out by Sorek and Braester [1986] andHuang [1988], flow equation (1) has a parabolic-hyperbolic
nature similar to the convection-dispersion solute transport
equation. Numerical solutions of the flow and transport
equations often exhibit numerical oscillations and/or artificial dispersion. These difficulties are most often encounteredwhen the convection term in the transport equation, or
equivalently, the gravitational term in the flow equation,
dominates [Neuman, 1984; Casulli, 1988]. Huang [1988]
extended the Eulerian-Lagrangian method of Neuman [1984]
to the solution of the variably saturated flow equation by
applying a simple, efficient SRPT solution scheme to thegravity term in (1). While successful for some tests, themethod suffered from mass balance errors and numericaloscillations when the flow problem was highly nonlinear
and/or extremely dry initial soil profiles were considered.
We shall show below that the accuracy and efficiency of themethod of Huang [1988] can be improved considerably byimplementing an adaptively corrected SRPT method.
Reformulation of the Richards’ equation

Following Huang [1988], (1) can be rewritten in the form-$g-S(z, t)

(9)

501

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW

which has a similar structure as the convection-dispersion
solute transport equation. Therefore we can follow theEulerian-Lagrangian approach which has proved to be very
effective for solving the convection-dispersion problem[Neuman, 1984]. For this purpose, (9) is first rearranged asfollows

where N is the number of nodes, and pi are chapeau
functions which vary linearly between neighboring nodes,
being one at node i and zero at all other nodes. If thepressure head hf at time t is known, hf+‘* at time level t +At can be solved using steps outlined in the next four
sections.
Single-step reverse particle tracking (SRPT)

Using the substantial or Lagrangian derivative
Dh

ah

Dt

dt+“*az

-=

ah
(11)

(10) can be formulated in Lagrangian form as( 1CD;=; Kg -S(z, t)(12)

The convective component K can be obtained using the
SRPT technique described by Neuman [1984] and Galeati etal. [1992]. The method considers each node i to be a
fictitious particle which moves from location Zi at time t
along its characteristic path to location zi, the coordinate ofnode i, at time t + At. The location of the fictitious moving
particle i at time t can be tracked backward according to(13), that is,zi = q -(13)

dZldt = v*

N

h<O

(14)

h>O

We emphasize that u* is different from the pore watervelocity, u = q/O, in which q is the Darcian fluid flux
density. The transformed Richards’ equation (12) holds now
for a predominantly parabolic type diffusion problem that iswell suited for solution using the FE method. Numericalsolution of (12) will lead to a symmetric matrix equation.
Equation (12) subject to the initial and boundary conditions can be solved on a fixed grid using the EulerianLagrangian approach by adopting the following FD approximation of the Lagrangian derivative [Neuman, 1984]:

ot=

(18)

j=l

Formally, the computed convective component can bedirectly substituted into (15), after which a conventional FE
or FD method can be used to solve (12) for hf+Ar. However,similar to numerical solutions for convection-dominatedsolute transport [Neuman, 1984], the SRPT technique tends
to produce numerical dispersion near the wetting front whenwater flow is dominated by the gravity component aKlaz in(1). In addition, numerical oscillations and mass balanceerrors may appear when the hydraulic properties are highly
nonlinear. The following adaptive correction procedure wasdeveloped to overcome these numerical difficulties.
As before, let pi be the convective component obtained by
the SRPT. A modified convective component h: may be
calculated at each time step by weighing pi with the pressure
heads of the first two nodes upstream of node i, that is,

,tj

At

(15)

where At is the time increment and 6 is the ‘convective
component’, that is, the contribution of the gravitationalterm to the pressure head. This contribution can be obtainedby implementing an adaptively corrected method of characteristics to be described later.
Assume that the soil profile is discretized into N - 1
elements of variable size, Azi (i = 1, 2, * * * , N - 1). The
dependent variable h can then be approximated using linear
finite elements as follows
N

h ~ ~ L C hi(t)~i(z)i=lpi = C hf+j(Zii,
Adaptive correction

1 aK aK~*~~~~_c ah
ae

ht+A’_(17)

Since the velocity u* is a function of z and t, an iterative
procedure is required to accomplish the above integration.
Once location zi is computed, the corresponding convective
component of node i can be obtained by interpolation
between nodal values as follows:

In the above equations, Z represents the characteristic pathof a moving particle, and v* is the velocity of the fluid
particle along the characteristic path:

Dh

u* dt

t

in which h no longer represents the pressure head at aparticular point in space and time but rather the pressure
head of a fluid particle moving along a characteristic pathdescribed by the equationv* = 0I

r+Ar

F; = ; ([(l - h)hf_, + Ah;_,] + &)(19)

where 0 5 A I 1 is a correction factor. The value of A is
automatically adjusted during each time step, according to
the global mass balance error 6, defined ass = (Wq + wi - wOut) - wStOr x
W star

looy

0

(20)

where W,, Wi, Wout, and W,,,, are the volumes of water
infiltrated, initially present, drained, and currently stored in
the soil profile respectively.
When the mass balance error at time t is negative and its
absolute value is larger than the error at t - At, A is
increased at time t + At by some factor (e.g., 0.05). When

502

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW

the mass balance error at time t is positive and greater thanthe error at t - At, A is decreased. However, A is kept
unchanged when the absolute mass balance error remains
the same or decreases. This approach assumes that the statevariable h at a node is affected primarily by the pressure
heads of the first two upstream nodes.
Galerkin finite element representation

We used the Galerkin finite element method to solve thetransformed flow equation (12). The method is similar to thescheme described in detail by van Genuchten [1982] andKool and van Genuchten [1991]. Substituting (16) into the
Lagrangian governing equation (12), performing the Galerkinorthogonalization, and applying the L2 mass-lumping
scheme used by van Genuchten [1982] and Milly [1985],
leads to the matrix equation

D{h}[A]{h} + [ B ] z = {F}(21)[AijI =

(22a)

where

[Bij] = Ci

’ +i dz(22b)

I0{Fi} = (-4 + K)$ilb -I1
0

S4i dZ(22c)

in which the subscripts i and j represent nodal indices. The
first term on the right-hand side of (22c) is not needed forprescribed pressure head boundary conditions. Analogous to
the pressure head, the nonlinear coefficients K and C are
interpolated linearly across each element as follows:N

K e C Ki+i(Z)

(23a)

i=l

N

C s 2 Ci+i(Z)

(23b)

i=l

where Ki and Ci are nodal values of the respective coefficients. A similar interpolation was also used for the root
uptake term S.
According to the Eulerian-Lagrangian approach, the nodal
Lagrangian derivative in (21) may be discretized by means ofbackward finite differences using the approximation
Dh.I h!+“‘h’
-z I
Dt
At

(24)

where the modified convective component ET is calculatedusing (19). This approximation is consistent with the notionthat node i is regarded as a particle reaching Zi at tktl
[Neuman, 1984]. Substituting (24) into (21) leads to the final
matrix equation([a] t-h [B,){h)““’ = & [B](E*] + {F}(25)

Time-stepping procedure and evaluationof nonlinear coefficients

Since the hydraulic conductivity K and the soil watercapacity C in (1) are nonlinear functions of the pressure headh, an iterative procedure is necessary at every time step toobtain the correct solution. We used the Picard iteration
method with under-relaxation [van Genuchten, 1982; Kool
and van Genuchten, 19911. Convergence was assumed when
the following condition was satisfied:&!‘A’ - h;+“l % 811hf+A'I + 6,

(26)

where 12 !+At and h ftAt are, respectively, the projected andcomputed values of the pressure head at the ith node at the
new time level, and 6i and a2 are relative and absoluteconvergence criteria. The value of At was automatically
updated for each time step depending on how quickly thesolution converged at the previous time step. For thispurpose we used the adaptive time-stepping scheme implemented previously by Kool and van Genuchten [1991].
Convergence and accuracy of the flow solution arestrongly influenced by how the nonlinear hydraulic conductivity K and the soil water capacity C are evaluated in the
numerical scheme. Different ways of evaluating the hydraulic conductivity, and especially the soil water capacity, canhave pronounced effects on the mass balance error, particularly when steep moisture fronts are present. For suchsituations, alternative solution schemes using water content
as the dependent variable have been proposed [e.g., Hills etal., 1989]. The advantage of O-based schemes is that they are
inherently mass conservative. However, as pointed out by
Milly [1985], a mass-conservative solution may also be
obtained with the head-based formulation if the capacity
term is evaluated in an appropriate manner. Mass conservation may be achieved by defining the soil water capacity sothat at the element level the following equality holds:

I

C,(httAt - h’) de =(OrtAr - 0’) de

(27)

Ie

e

in which C, is the effective element soil water capacity givenby (23b) and where integration is performed on an element
by element basis. According to the L2 mass-lumping scheme[Milly, 1985] the resulting expressions for the nodal values of
C at time t are given by(28)
When changes in the nodal pressure heads are very small,the denominator of (28) approaches zero, in which case thesoil water capacity is best computed analytically as thederivative of the soil water retention curve evaluated at thehalf-time level hft1’2A’ = 1/2(h’ + hffAt), that is,de (h)
ci = dhh ,+112&r(29)

In our simulations, nodal values of hydraulic conductivity
were always evaluated at the half-time level using (7), that is,
K.I = K(h!+“*A’).
I

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW

503

0 FE

Figure 1.

Calculated water content (0) profiles versus depth (z) for ponded infiltration in soils with (Y =0.01 cm-’ and (a) n = 1.3, (b) n = 2.0, and (c) n = 5.0. Solid lines represent characteristic finite element(CFE) solutions; circles give finite elements (FE) results.

Numerical Verification
The proposed characteristic finite element scheme (CFE)
was tested by simulating ponded infiltration into a 200-cmdeep soil profile. Results obtained with the CFE method will
be compared with solutions of the mass-conservative finite
element (FE) method of Kool and van Genuchten [1991]. To
facilitate comparisons between CFE and FE, all simulations
used L2 mass lumping, employed the same adaptive timestepping scheme, and implemented the same spatial discretization, using uniform spatial increments AZ of 2 cm. The
simulations were carried out on an IBM 486/33 personal
computer. CPU time was used to compare the computational
efficiency of the different methods, while the numericalaccuracy was evaluated in terms of the mass balance error 6
given by (20). Simulations were carried out for different
hydraulic parameters and different initial pressure heads hi.
The residual and saturated water contents in (6) were assumed to be 0.05 and 0.45, respectively, in all simulations.
Results will be presented for ponded infiltration assuming azero pressure head (ha = 0) at the soil surface.
Figure 1 compares water content profiles calculated using
the CFE method (solid lines) with those obtained using
standard FE (open circles) assuming hi = - 15,000 cm, (Y =0.01 cm-‘, K, = 50 cm d-‘, and three values of n; that is,n = 1.3, typical of a relatively fine-textured soil, n = 2 fora medium- to coarse-textured medium, and n = 5 for acoarse-textured soil [Camel and Parrish, 1988; Yates et al.,1992]. The CFE and FE simulation results are in close
agreement except for those in Figure la where the FE
solution appears to be slightly ahead of the CFE solution.
Both methods yielded mass balance errors of less than 0.05%
in this example. For n = 1.3, the CFE simulation requiredslightly less CPU time (13.0 min) than the FE scheme (15.4
min). The improved efficiency of the CFE method is more
apparent when n = 5 (Figure lc), for which the CFE and FE
schemes used 9.1 and 22.6 min CPU time, respectively. By
comparison the CFE and FE schemes for n = 2 (Figure lb)
required 9.8 and 9.2 min, respectively.
Figure 2 shows pressure head distributions corresponding

to the water content profiles in Figure 1. It is evident from
Figures 1 and 2 that the pressure head wetting fronts aremuch steeper than those for the water content. All CFE
simulations were found to be stable (no oscillations wereobserved). The FE solutions, on the other hand, exhibitedsome minor oscillations when n = 1.3 in the near-saturated
zones between the wetting fronts and the soil surface. The
FE scheme for n = 1.3 also produced slightly positive
pressure head values (0.5 ~ 3.0 cm) in this zone. These
results, not further noticeable in Figure 2a, are inconsistent
with the simulated ponded infiltration problem which invoked a zero pressure head boundary condition at the soil
surface. The simulated positive pressure heads probably also
caused the slightly steeper wetting fronts in Figures la and2a for the FE scheme as compared to the CFE method.
The FE results in Figures 1 and 2 were obtained with theL, mass-lumping scheme of Milly [1985]. Implementing the
mass-conservative mixed-form solution of the Richards’
equation [Celia et al., 1990] did not significantly improve theresults for n = 1.3. This scheme also exhibited oscillatory
solutions and, like the standard FE method, produced positive head values between the wetting fronts and the soil
surface. By comparison, the CFE never produced oscillations in computed pressure head distributions. The numerical oscillations and positive pressure values in the FE results
were found to increase at increasingly lower II values, andlikely were caused by the shape of the conductivity function
near saturation. As shown elsewhere [van Genuchten andNielsen, 1985; van Genuchten et al., 1991], relatively small nvalues (e.g., 1.0 < n < 1.3) in (6) and (7) will lead to
extremely nonlinear K(h) functions near saturation, sometimes with important repercussions for accurately simulating
variably saturated flow [Vogel et al., 1992]. Adoption ofalternative K(h) functions, such as those by Brooks andCorey [1964], or modifications of (6) and (7) as proposed by
Vogel and Cislerova [1988], may significantly decrease thenonlinearities in the soil hydraulic functions and lead to morestable flow (but dramatically different) simulation results[Vogel et al., 1992].

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW504

a-0.01. n-2.0h,=-IWO0 em
t-0.2 d

(‘4

(a)
50

(c)

t=0.2 d

0.6
0.4

0.8
1.0

OFE
- CFE

0’

-150000.

. -h (cm)h (cm)h (cm).

Calculated pressure head (h) distributions versus depth (z), corresponding to the results in
Figure 1. (a) n = 1.3; (b) IZ = 2.0; and (c) IZ = 5.0.

Figure 2.

Figure 3 examines the effect of the parameter (Y on
numerical simulations assuming n = 2, K, = 50 cm d-l, an
initial pressure head hi of - 10,000 cm, and three values ofthe parameter (Y: 0.1, 0.01, and 0.001 cm-‘. The CFE
solutions show relatively good agreement with the FE results for all three values of (Y. The mass balance error forboth methods was always less than 0.02%. The results in
Figure 3 are indicative of the fact that simulations withhigher (Y values (Figure 3a) generally produce steeper wetting fronts than simulations with smaller (Y values (Figure3c), with the CFE scheme requiring somewhat less computer
time at the higher cr as compared to the FE method. The
CPU times for the CFE and FE schemes were, respectively,
3.3 and 7.8 mins for a! = 0.1, 8.3 and 7.4 min for (Y = 0.01,
and 4.5 and 4.2 min for (Y = 0.001.
Next, we examined the performance of the CFE scheme
for infiltration in soils having relatively large II values. Figure
4 demonstrates the capability of the CFE method to solve

the flow problem with an extremely high value of 10 for II.
Soils of this type are characterized by very narrow pore andparticle size distributions (usually coarse-textured soils) andtypically have very nonlinear 0 (h) as well as K(h) hydraulicfunctions [van Genuchten and Nielsen, 1985]. The simulations again pertained to ponded infiltration into a soil profile
having an initial pressure head hi of -10,000 cm and asaturated hydraulic conductivity K, of 500 cm d -’ andassuming (Y = 0.1 cm-‘. The CFE computed water content
and pressure head distributions are shown in Figures 4a and4b, respectively. CPU time for this example was 6.0 min,while the mass balance error was 0.53%. The standard finite
element (FE) method assuming the same finite element grid
system as used for CFE (constant AZ of 2 cm) failed toconverge for this infiltration problem. The mixed-form
method of Celia et al. [ 1990] also failed to converge when thestandard 2 cm for the grid spacing A Z was used. Smaller
nodal spacings also did not lead to convergent FE solutions.

a=0.01. n-2.0h,=-10+X0 cm

(b)

0

h

l-0.2 d

OFE
0.0

I
0.1

.

- CFE
I ’ I * I ’
0.2 0.3
0.4(
e-I
0.2 0.3 0.4
0.5eSimulated water content (0) distributions versus depth (z) for ponded infiltration in soils withn = 2.0 and (a) (Y = 0.1, (b) cy = 0.01, and (c) (Y = 0.001 cm-‘.

Figure 3.

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW
0

50.5

0

.

(a)

t=0.02 d

5050

0.04

N

2
s
N

0.08

150-

100

0.1
a-0.1. n-10.0

150h,--10000 cm
KS-500 cm/d
200

,

,

0.00 0.10.,.,.,
0.20 0.30 0.40““‘-shod”” 0h (cm)

0

Figure 4. Calculated (a) water content, 8, and (b) pressure
head, h, distributions versus depth, z, obtained with the

CFE method for ponded infiltration in a very coarse-textured(n = 10) soil profile.

Accurate and stable simulations of infiltration into very
dry soils are often difficult to obtain with standard numericalmethods [Zaidel and Russo, 1992] because of the development of steep wetting fronts. Figure 5 shows one set ofresults obtained with the CFE and FE schemes for ponded
infiltration into a soil having a very low initial pressure headhi of -50,000 cm. The simulation concerned a medium- tocoarse-textured soil assuming (Y = 0.01 cm -I, n = 2, and K,= 50 cm d-i. Simulation of the O.l-day infiltration event
required 10.7 and 12.0 min for the FE and CFE methods,
respectively. Notice the close agreement between the two
methods in Figure 5.
All calculations thus far were obtained with 2-cm ele-

1

-I

I.50

(b)

-50000-25000h (cm)

Figure 5. Calculated (a) water content (0) and (b) pressure
head (h) distributions versus depth, Z, for ponded infiltration
into an initially very dry (hi = -50,000 cm) soil profile.

a=O.124. n=2.28
e,=o.41. e,=o.osi
hi=-10000 cm---- Ar=4.0 cm
0 Az=2.0 cm- Az=l.O cm-10000

200
c

)

I
0.10

'
0.20

I

,

0.30

7

I
0.40

e

Figure 6. Effect of spatial discretization, A Z, on simulated
water content (0) distributions versus depth (z) using the
CFE method.

ments. We also examined the effects of different spatial
increments on the results, again for ponded infiltration.
Figure 6 compares calculated water content distributions
assuming element sizes AZ of 1, 2, and 4 cm. The corresponding CPU times were 7.5, 2.6, and 0.98 min, and theintroduced mass balance errors were 0.0397%, 0.055%, and0.0125% for AZ = 1, 2, and 4 cm, respectively. The results
indicate that the solutions with A Z = 1 cm and AZ = 2 cm
are very close, while some relatively minor deviations
appear near the wetting front when A Z = 4 cm. This showsthat the CFE method can produce accurate solutions using
relatively large spatial increments, thus further improving
the computational efficiency as compared to the FE method.
The above results, as well as those of many other simulations, indicate that the proposed CFE method is computationally somewhat more efficient than the standard FE forsituations leading to relatively steep pressure head wetting
fronts. Such situations most often arise for soils having verylarge or very small n values, and/or relatively large (Y values.
However, the FE method is equally or slightly more efficientthan the CFE scheme for more moderate nonlinear infiltration problems involving intermediate n and relatively small (Yvalues. This feature results from the fact that the CFE
method requires additional computations for solving (17)that are associated with the gravity term in the flow equation. Intermediate values of n in (6) and (7) do not lead tostrongly nonlinear soil hydraulic properties and steep wetting fronts, and hence there appears little advantage toemploying the particle tracking method. These conclusions
are consistent with studies by van Genuchten and Nielsen[1985] and van Genuchten et al. [1991] who show that theunsaturated hydraulic properties, especially K(h), become
very nonlinear for soils having relatively small (e.g., 1 .O < n
5 1.3) or large (n 2 3) values of n. Numerical simulationsof infiltration in such soils using standard FE or FD methods
are generally the most difficult, especially for initially very

506

HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW

dry soil profiles. The most extreme cases arise when thegravity term of (1) is relatively dominant, in which case theflow equation acquires hyperbolic properties which are besttreated with the proposed CFE method.
The results in this paper pertain to ponded infiltration. We
also carried out several simulations for constant flux infiltration at the soil surface. The proposed CFE method wasfound to work efficiently and accurately and found to be
oscillation-free for most numerical tests. However, someminor oscillations were observed when the applied surface
boundary flux q0 was close to K, for infiltration into
fine-textured soils having extremely small n values (such asn = 1.05).

Summary and Conclusions
The rigorous splitting approach based on the EulerianLagrangian concept, initially introduced for solving convection-dispersion type solute transport problems, was extended to variably saturated water flow. The approachpermits one to separate the original flow problem into
“convection” and “diffusion” parts, which are then solved
with the method of characteristics and the Galerkin finite
element method, respectively. A self-adaptive correction
method was developed to modify the “convection contribution” obtained by the single-step reverse particle tracking
technique. The correction procedure consists of weighting
the convection contribution with the pressure heads of two
upstream nodes. The weighting factor, 0 5 A 5 1, wascontinually updated in the program so that the global massbalance error could be minimized.
Numerical simulations of ponded infiltration for a wide
range of soils showed that the Eulerian-Lagrangian approachworks well for variably saturated flow problems describedwith the Richards’ equation, including situations involvingvery dry initial conditions. While maintaining the simplicityof the single-step reverse particle tracking (SRPT) method,
the proposed adaptive correction technique significantly
improved the mass balance accuracy of previous characteristic finite element (CFE) methods. The proposed CFE
method yielded results which were also virtually free ofoscillations. The only situations where some minor oscillations were observed were for flux infiltration into soils withvery small n values when the infiltration rate q. approached
the value of the saturated hydraulic conductivity, K,. The
CFE scheme appears especially attractive for highly nonlinear flow problems, such as for fine-textured soils having
small n values (n < 1.3) and coarse-textured soils having
large values of Iz and (Y, particularly when very dry initial soil
profiles are present. Such nonlinear problems are very
difficult to simulate with most standard finite difference or
finite element methods.
For simulating water flow in soils having intermediate nvalues (1.5 < n < 3.0) the CFE method required approximately the same or slightly more CPU time than the standard
finite element (FE) method, mostly because the CFE requires additional calculations associated with the particle
tracking scheme. However, the CFE method is more efficient and numerically far more robust than the.standard FE
method for problems having highly nonlinear hydraulicproperties and/or very dry initial soil profiles. As such, we
believe that the CFE method may become an attractive
alternative to existing methods for simulating variably saturated water flow.

References
Abriola, L. M., Finite element solution of the unsaturated flow
equation using hierarchic basis functions, in Proceedings of the
6th International Conference on Finite Elements in Water Resources, edited by A. Sa Costa, A. Melo Baptista, W. G. Gray,

C. A. Brebbia, and G. F. Pinder, pp. 125-133, Springer-Verlag,
New York, 1986.
Broadbridge, P., and I. White, Constant rate rainfall infiltration: A
versatile nonlinear model, 1, Analytical solution, Water Resour.
Res., 24, 145-154, 1988.
Brooks, R. H., and A. T. Corey, Hydraulic properties of porous
media, Hydrol. Pap. 3, 27 pp., Colo. State Univ., Fort Collins,
1964.
Carsel, R. F., and R. S. Parrish, Developing joint probability
distributions of soil water retention characteristics, Water Resour. Res., 24, 7.55-769, 1988.
Casulli, V., Eulerian-Lagrangian methods for hyperbolic and convection dominated parabolic problems, in Modeling Surface andSubsurface Flow, Developments in Water Science Series, vol. 35,
edited by M. A. Celia, L.-A. Ferrand, C. A. Brebbia, W. G. Gray,and G. F. Pinder. p p . 183-191. Elsevier. New York. 1988.
Celia, M. A., and G. F. Pinder, An alternating-direction collocation
solution for the unsaturated flow equation, in Proceedings of the
6th International Conference on Finite Elements in Water Resources, edited by A. Sa Costa, A. Melo Baptista, W. G. Gray,
C. A. Brebbia, and G. F. Pinder, pp. 395-410, Springer-Verlag,
New York, 1986.
Celia, M. A., E. T. Bouloutas, and R. L. Zorba, A generalmass-conservative numerical solution for the unsaturated flow
equation, Water Resour. Res., 26, 1483-1496, 1990.
Dane, J. H., and F. H. Mathis, An adaptive finite difference scheme
for the one-dimensional water flow equation, Soil Sci. Soc. Am.
J., 45, 1048-1054, 1981.
Galeati, G., G. Gambolati, and S. P. Neuman, Coupled and partiallycoupled Eulerian-Lagrangian model of freshwater-seawater mixing, Water Resour. Res., 28, 149-165, 1992.
Gamliel, A., and L. M. Abriola, A one-dimensional moving grid
solution for the coupled non-linear equations governing multiphase flow in porous media, 1, Model development, Znt. J.
Numer. Methods Fluids, 14, 25-45, 1992.
Gottardi, G., and M. Venutelli, Moving finite element model for
one-dimensional infiltration in unsaturated soil, Water Resow.
Res., 28, 3259-3267, 1992.
Hansen, R., M. S. Espedal, and 0. Nygaard, An operator splitting
technique for two-phase immiscible flow dominated by gravity
and capillary forces, in Computational Methods in Water Resources IX, vol. 1, Numerical Methods in Water Resources,
edited by T. F. Russell, R. E. Ewing, C. A. Brebbia, W. G. Gray,and G. F. Pinder. DD. 59-67. Elsevier Science. New York. 1992.
Hills, R. G., I. Porro, D. B. Hudson, and P. J. Wierenga, Modeling
one-dimensional infiltration into very dry soils, 1, Model development and evaluation, Water Resour. Res., 25, 1259-1269, 1989.
Huang, K., Solution of saturated-unsaturated flow by finite element
or finite difference methods combined with characteristic technique, in Proceedings of the 7th Computational Methods in Water
Resources, vol. 2, edited by M. A. Celia, L. A. Ferrand, C. A.
Brebbia, W. G. Gray, and G. F. Pinder, pp. 241-246, Elsevier,
New York, 1988.
Kool, J. B., and M. T. van Genuchten, HYDRUS: One-dimensionalvariably saturated flow and transport model, including hysteresis
and root water uptake, Res. Rep. 124, U.S. Salinity Lab., Agric.
Res. Serv., U.S. Dep. of Agric., Riverside, Calif., 1991.
Milly, P. C. D., A mass-conservative procedure for time-stepping inmodels of unsaturated flow, Adv. Water Resour., 8, 32-36, 1985.
Neuman, S. P., Adaptive Eulerian-Lagrangian finite element
method for advection-dispersion, Znt. J. Numer. Methods Eng.,20, 321-337, 1984.
Parlange, J., Theory of water movement in soils, 6, Effect of water
depth over soil, Soil Sci., 113, 308-312, 1972.
Philip, J. R., The theory of infiltration, 1, The infiltration equationand its solution, Soil Sci., 83, 345-357, 1957.
Ross, P. J., Efficient numerical methods for infiltration using Richards’ equation, Water Resour. Res., 26, 279-290, 1990.
Sorek, S., and C. Braester, An adaptive Eulerian-Lagrangian approach for the numerical simulation of unsaturated flow, inProceedings of the 6th International Conference on Finite Ele--HUANG ET AL.: EULERIAN-LAGRANGIAN APPROACH TO SIMULATE WATER FLOW
ments in Water Resources, edited by A. S a Costa, A. Melo

Baptista, W. G. Gray, C. A. Brebbia, and G. Pinder, pp. 87-100,
Springer Verlag, New York, 1986.
van Genuchten, M. T., A closed-form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J.,
44, 892-898, 1980.

van Genuchten, M. T., A comparison of numerical solutions of the
one-dimensional unsaturated-saturated flow and transport equations, Adv. Water Resow., 5, 47-55, 1982.
van Genuchten, M. T., and D. R. Nielsen, On describing andpredicting the hydraulic properties of unsaturated soils, Ann.
Geophys., 3, 615-628, 1985.
van Genuchten, M. T., F. J. Leij, and S. R. Yates, The RETC codefor quantifying the hydraulic properties of unsaturated soils,
EPA/600/2-91/065, 85 pp., R. S. Kerr Environ. Res. Lab., U.S.
Environ. Prot. Agency, Ada, Okla., 1991.
Vogel, T. N., and M. Cislerova, On the reliability of unsaturated
hydraulic conductivity calculated from the moisture retention
curve, Transp. Porous Media, 3, l-5, 1988.
Vogel, T. N., M. T. van Genuchten, R. Zhang, and M. Cislerova,
Effect of the description of soil hydraulic properties near satura-

507

tion on numerical simulation of soil water movement, Agron.
Abstr., 84, 230, 1992.

Warrick, A. W., A. Islas, and D. 0. Lomen, An analytical solution
to Richards’ equation for time-varying infiltration, Water Resow.
Res., 27, 763-766, 1991.

Yates, S. R., M. T. van Genuchten, and F. J. Leij, Analysis ofpredicted hydraulic conductivities using RETC, in Proceedings ofthe International Workshop on Indirect Methods for Estimating
the Hydraulic Properties of Unsaturated Soils, edited by M. Th.
van Genuchten, F. J. Leij and L. J. Lund, pp. 273-283, Universityof California, Riverside, Calif., 1992.
Zaidel. J.. and D. Russo. Estimation of finite difference interblock
conductivities for simulation of infiltration into initially dry soils,
Water Resour. Res., 28, 2285-2295, 1992.

K. Huang, M. T. van Genuchten, U.S. Salinity Laboratory, 4500
Glenwood Drive, Riverside, CA 92501.
R. Zhang, Department of Plant, Soil, and Insect Sciences, University of Wyoming, P.O. Box 3354, Laramie, WY 82071-3354.
(Received March 11, 1993; revised August 30, 1993;
accepted October 6, 1993.)