WATER RESOURCES RESEARCH, VOL. 29, NO. 12, PAGES 41314138, DECEMBER 1993

An Efficient Eulerian-Lagrangian Method for Solving Solute
Transport Problems in Steady and Transient Flow Fields
R. Z H A N G , K. HUANG, AND M. T. VAN G E N U C H T E N
U.S. Salinity Laboratory, Agricultural Research Service, U.S. Departmentof Agriculture, Riverside, California
A computationally efficient, yet relatively simple Eulerian-Lagrangian method is proposed forsolving the one-dimensional convection-dispersion solute transport equation assuming a steady ortransient velocity field. The method uses a modified single-step reverse particle tracking (MSRPT)
technique to handle steep concentration fronts. The scheme utilizes two weighting factors to control
the movement of particles during a backward tracking step. One weighting factor greater than unity isused in the upstream region of the convection front, while another weighting factor less than unity istaken in the downstream region. The two factors were related empirically to the grid Peclet and
Courant numbers. The MSRPT technique is carried out only within the concentration plume at eachtime step. For transient flow fields, the weighting factors were determined using an automaticallyadjustable procedure based on mass balance errors. The MSRPT method maintains the advantages of
the traditional single-step reverse particle tracking (SRPT) procedure, i.e., producing efficient andoscillation-free calculations, but circumvents numerical dispersion introduced by SRPT. A large
number of tests against analytical solutions for one-dimensional transport in uniform flow fields
indicate that the proposed method can handle the entire range of Peclet numbers from zero to infinity.
Numerical tests also show that the MSRPT method is a relatively accurate, efficient and massconservative algorithm for solute transport in transient flow fields. The Courant number at present
cannot exceed 1. The MSRPT approach was found especially useful for convection-dominatedproblems; in fact, an exact numerical solution may be obtained with MSRPT for pure convection.

Convection-dispersion type equations are being widelyused to model solute transport in soil and groundwater.
Owing to the particular combination of hyperbolic andparabolic terms, serious difficulties are often encountered in
obtaining accurate numerical solutions of these equations. A
variety of numerical schemes have been developed to deal
with these difficulties, including an extensive number ofEulerian methods using fixed grid systems, and Lagrangian
approaches involving moving coordinates. There is now a
growing trend to combine the Eulerian and Lagrangian
methods [Garder et al., 1964; Leith, 1965; Heinrich et al.,
1977; Konikow and Bredehoeft, 1978; Varoglu and Finn,
1978; Cheng et al., 1984; Khaleel and Reddell, 1985; Fujinawa, 1986; Dimou and Adams, 1991; Leonard, 1991;
Galeati et al., 1992; Huang et al., 1992, among others].

Eulerian-Lagrangian methods treat the convection part of
the transport problem using a Lagrangian formulation in a
fixed Eulerian grid, while the dispersion term is treated byeither an Eulerian or a Lagrangian formulation. As a typicalexample, Neuman and Sorek [1982] decomposed the concentration field into two parts, one controlled by pure
convection, the other by dispersion. The convection partwas handled with the help of moving particles, while thedispersion problem was solved using finite elements on a
fixed grid. Neuman [ 1984] subsequently presented an adaptive Eulerian-Lagrangian finite element method for solving
the convection-dispersion equation. The method is capableof dealing with the entire range of Peclet numbers from zero
to infinity. However, the scheme is quite troublesome to
implement and time consuming to execute because of theThis paper is not subject to U.S. copyright. Published in 1993 by
the American Geophysical Union.
Paper number 93WR01674.
4131

need to continuously track the front using numerous particles
at each time step. In a related approach, Yeh [1990] used theEulerian-Lagrangian method with a zoomable (or adjustable)
hidden fine-mesh approach to resolve numerical difficulties.
While this scheme reduces or virtually eliminates numericaldispersion and oscillations, the process of zooming and refining
the elements at each time step is not straightforward in terms of
its practical implementation. Also, the large number of elements needed for this approach requires excessive amounts ofcomputer memory and execution time.
In the adaptive Eulerian-Lagrangian finite element method
used by Neuman [ 1984], the convective components of steep
concentration fronts were tracked forward with moving
particles around each front. Away from the fronts, however,the convective problem was solved with a single-step reverse particle tracking (SRPT) method. The SRPT techniquerequires much less computer time and storage than continuous forward particle tracking. Unfortunately, one drawback of the SRPT method is that it introduces some numerical dispersion near sharp concentration fronts.
The objective of this paper is to present a simple andefficient particle tracking technique to solve the convectiondispersion equation. A modified single-step reverse particle
tracking (MSRPT) method is used to deal with convectiondominated transport problems. We shall show that MSRPT
maintains the advantages of the traditional single-step reverse particle tracking procedure, but eliminates numericaldispersion introduced by SRPT. Simulation results based on
the scheme will be compared with analytical solutions aswell as numerical results obtained with other methods.
GENERAL FORMULATION
Transport of miscible components in a one-dimensionaltransient flow field is described with the convectiondispersion solute transport equation

Purchased by USDA for Official Use

4132

ZHANG ET AL.: EULERIAN-LAGRANGIAN METHOD FOR SOLVING SOLUTE T RANSPORT P ROBLEMS

i3Rg=;

eoz

i

1

-0,;

(1)

where c is the solution concentration, I is time, x is distance,
D is the dispersion coefficient, u = q/O is the pore watervelocity, q is the Darcian fluid flux, 0 is the volumetric watercontent, and R is a retardation factor accounting for linear
equilibrium sorption; i.e.,
R = 1 + pk/O

(2)

in which p is the porous medium bulk density and k is anempirical distribution constant. The equation will be solved
here subject to the initial and boundary conditionsc(x, 0) = c(J

(3)

8C

-/3D-+v(c-CJ=O

ax

ac
ax

-zz 0

inflow boundary

outflow boundary

(4)

(5)

where Co and C, are prescribed constant concentrations,while p controls the type of boundary condition imposed at
the inlet position (p = 0 for a prescribed concentration
condition and /3 = 1 for a prescribed flux condition).
By using the Lagrangian derivative

where L is the solution domain. Neuman and Narasimhan[1977], among others, showed that parabolic problems suchas (10) are more effectively solved by implementing a mass
lumping scheme in which the first term of (10) is approximated bydci
f3R;ip;dx=;j;
IL

dt=

At

(12)

(13)

where
,Dzzdx

Eij=

(14)

IL
Hi = -qi

inflow boundary

Hi = 0outflow boundary(15)

wjj= 6ij

ORq; dx(16)

JL

where c no longer represents the concentration at a point inspace and time, but rather the concentration of a fluid
particle moving along the characteristic path described by
the equation

N

ck+’
- ci
I

ck+r} = [G] + F (~1(7)

The Lagrangian formulation of the governing equation eliminates the convection term so that the equation takes on apurely parabolic appearance which can be solved efficiently
with a finite element method.
Let us define the finite element approximation of thesolute concentration at location x and time t as

dc.I

where Fi is the “convective concentration” of node i to be
discussed in detail in the following section, and c,!‘+’ is theconcentration of node i at tk+, . Following the procedure of
Neuman [1984], we substitute (11) and (12) into (10) and use
the linear basis functions to eventually obtain a system of
linear equations for the concentration at time step tk+l =tk + At aS fOllOWS:

where v* = v/R, (1) can be rewritten in the Lagrangian form
as

(8)

(11)

The Lagrangian derivative may be approximated by [Neuman, 19841(6)

dxldt = v*

BR~i dxIL

Gi = -qiCt

inflow boundary(17)

Gi = 0outflow boundary

The global matrices Hi and Gi in these equations were
derived assuming /3 = I (i.e., for prescribed flux conditions).
When p = 0, c is known at the inflow boundary.
The relative accuracy of the different schemes will be
evaluated, in part, by using the mass balance error (ME) a s
follows:
ME = lOO(C, - C,)IC,(%)(18)

Cm = Cini + Ci, - Co”,

(19)

with

c(X, t) = 2 ci(t)9i(X)

(9)

i=l

where N is the number of nodes in the solution domain,vi(X) represents the linear shape function, and Ci is thec o n c e n t r a t i o n a t n o d e i (i = 1, * * * , N). Applying theGalerkin procedure to (7) and incorporating (9) leads toI, [,R ~-~ (eD 3]~i dX=O(10)

i= 1, 2, ***, N

where C,, Cini, Gin, and C,,, are the amounts of solute
stored in the profile, initially present in the profile, cumulative input, and cumulative output, respectively.

M ODIFIED S INGLE-STEP R E V E R S E
P ARTICLE T RACKING
The convective concentration, c in (12), may be computed
using a single-step reverse particle tracking (SRPT) tech-

4133

ZHANG ET AL .: EULERIAN - LAGRANGIAN M ETHOD FOR S OLVING S OLUTE T RANSPORT P ROBLEMS

nique. That is, a fictitious particle from each node i is sent
backward to the point
fk+lx; = xi -u* dt- cu = 0.5
A-A Cu = 0.25
H Cu = 0.125

(20)

I tk
during each time step. This means that a particle leaving xi
at tk will reach grid point location Xi exactly at tk+, . The
convective concentration is subsequently computed by thefinite element approximation; i.e.,Nc(Xi, tk) = C(X:, tk) z c c(Xj, tk)‘Pj(X:)

(21)

j=l

l/Pe

This technique works well in regions where the concentration gradients are small; however, the method introduces
artificial dispersion near sharp concentration fronts.
To eliminate the numerical problem, a modified single-stepreverse particle tracking (MSRFT) technique is formulated
to compute the convective concentration. First, the singlestep particle tracking method is carried out within a certain
region determined by
xt = x0 + x, + xd

(22)

where x0 is the initial center location (at t = 0) of a sharp
front, and x, and xd are the travel distances related toconvection and dispersion, respectively; i.e.,
t

x, =

0.0

u* dt

(23)

I0

x;=xi_

x; = xi -fkil

Relationships between (top) w1 and l/Pe and (bottom) w2

w1 = 1 + (l/Cu - 1) exp [-2.63(1/PeC~)‘.~‘~]
lIPeCu > 0.05

(26b)w2 = 1 - (1 - Cu) exp [-2.9(C~lPe)~.‘~]

(27)

where

w,u* dt

Pe = u*AxlD
upstream region(25a)downstream region(25b)w2v* dt

I ft
x0 + x, 5 x 5 xt

1 .o

1

I It
x < x0 + x,

0.8

(24)

The term x0 + x, in (22) indicates the center location of the
sharp front at time t, while (24) is an empirical expression
based on numerical experiments. For elements within xt, thelocation at tk of a paI?iCk which reaches Xi at tk+l iS
computed by
fkfl

0.6

and l/Pe for varying values of Cu.

l/2
0
I

0.4

l/Pe
Fig. 1.

’ D dtxd = 40.2

where w , and w2 are weighting factors for the upstream anddownstream regions, respectively. Next the convective concentration is computed by (21) using the calculated xi.
The weighting factors wt and w2 were found to changewith different flow velocities, dispersion coefficients, andspace and time discretizations. They were related empirically to the grid Peclet number (Pe) and Courant number(Cu) as follows:
w1 = 1 + (l/Cu - 1) exp [-4l.O2(1/PeC~)‘.‘~~]
lIPeCu IO.05(26a)

Cu = v*AtlAx(28)

Equations (26) and (27) were obtained by carrying out a large
number of simulations for a steady flow field using different
values for Pe and CU. For each pair of Pe and Cu values, w,
and w2 were adjusted by minimizing the mass balance error,
and by minimizing deviations between the numerical and theavailable analytical solutions. Results of (26) and (27) arecurrently limited to C u s 1. Since the weighting factors
depend only on the grid Peclet and Courant numbers and arecomputed for each grid, the MSRPT method can be used forproblems with both uniform and varying element sizes.
Figure 1 shows families of curves of wt and w2 versus
l/Pe for different values of Cu. The weighting factor w, isgreater than unity, and increases as Pe increases and/or C u
decreases. On the other hand, the weighting factor w2 is lessthan unity, and increases as Pe decreases and/or Cu increases. Figure 2 presents similar relationships between w,
and w2 and 1/Pe for Cu = 0.25. Notice that when Pe + 00,
Wl = l/Cu, and xd = 0. This means that in the region x 5

ZHANG ET AL.: EULERIAN-LAGRANGIAN METHOD FOR SOLVING SOLUTE TRANSPORT PROBLEMS

4134

51

where iV is the number of nodes, CT is the simulated
concentration, and Ci is the presumed correct solution of theconcentration at node i. For steady state flow, Ci will be
evaluated using analytical solutions. No useful analyticalsolutions exist for problems of solute transport in transient
flow fields. For these problems we used as our correct
solution simulated results obtained with a finite element
scheme which assumed very fine grid and time discretizations.

- W,
- PYI

EXAMPLES
Problem 1 deals with the solution of an advancing concentration front. The analytical solution of the transport equation in a uniform flow field, i.e., forl/Pe
Fig. 2. Dependency of wr and w2 on l/Pe for Cu = 0.25.

(33)
subject toxi, i.e., within x0 + xc, the location xi of a moving particle
beginning from node i is given exactly by
X: = xi - AX = Xi-1

(29)

In other words, the moving particles coincide perfectly with
the nodes of the fixed grid system, and hence~(Xi, tk) = C(X:, tk) = C(Xi-I, tk)(30)

Since there is now no interpolation error for c, exact results
can be obtained for the pure convection problem. For
dispersion-dominated problems, both wr and w2 approach
unity (Figure 2), and subsequently the MSRPT method
reduces to the traditional SRPT technique which provides
sufficiently accurate solutions for such problems [Neuman,
1984]. We note that the SRPT method was only used tocompute the convective concentration of node i outside theregion of x f.
For transient problems the flow velocity and dispersion
coefficient change with time and space, in which case thePeclet and Courant numbers become functions of time andspatially dependent. We introduced a correction factor a,with the weighting factors wr and w2 in (25) to account for
the transient effects; i.e.,w\(t) = UfWl

Wi(t) = urw2

(31)

where w 1 and w2 are computed with (26) and (27) as before.
The correcting factor a, in the numerical scheme is determined using an adaptive procedure which depends on thevalue of the mass balance error. The procedure consists of
two steps. First, the initial value a0 of a, at t = 0 is set equal
to unity. Next, at each time step, the mass balance error(ME) is computed. If ME is within a preset tolerance (e.g.,
IME] < 0.l%), the value of a, is kept unchanged; i.e., u,+~= a,. Otherwise, a,,, is set equal to 0.95 a, if ME < 0, andto 1.05 a, if ME > 0.
In addition to the global mass balance error, the relative
accuracy of the different schemes are also compared in termsof the sum-squared error E defined as
E ( t ) = i [Ci(t) - CT(t)12
i=l

(32)

c(x, 0) = 0

OSX<~

(34)

c(0, t) = 1

t>O

(35)

E (00, t) = 0

t>O

(36)

is given by [van Genuchten and Alves, 1982]
1 x-ut 1
1 I
c(x, r) = 2 [erfc L2(DRt) 1121

+ exp (:I erfc [,~h,:,:lJ(37)

which, for the limiting case when D -+ 0, reduces toc(x, t) = 1

x < vtc(x, t) = 0

otherwise(38)

Three simulations were carried out with the MSRPT
scheme. For case 1, the following parameters were used:
Ax = 200, At = 100, t = 9600, u = 0.5, R = 1, and D =0 (Pe = a), and assuming a domain of 0 5 x 5 12800 (anyconsistent set of units may be used for these parameters).
Cases 2 and 3 were the same as for case 1 except for D which
was set at 2 (Pe = 50) and 50 (Pe = 2), respectively. Figure
3 displays the calculated concentration profiles for the threecases. Results obtained with the SRPT and MSRPT methods
are compared with the analytical solution. The calculations
show that over the whole range of 0 I Pe < 00, the MSRPT
eliminates both numerical oscillation and artificial dispersion. The SRPT method on the other hand produces significant numerical dispersion, except for relatively small Pe
values. Figure 4 further compares results obtained withMSRPT, the finite element method (FEM), and the upstream
weighted method (UWM) [Heinrich et al., 1977] against theanalytical solution. The upstream weighted method uses
nonlinear weighting functions so that greater weight isplaced on variables associated with the upstream node of anelement. The parameters used for this example were thesame as those of case 1 except for D = 1 (Pe = 100). Noticethe excellent results obtained with the MSRPT method. By

ZHANG ET AL .: EU L E R I A N - LAGRANGIAN M ETHOD FOR S OLVING S OLUTE T RANSPORT

comparison, the FEM exhibits numerical oscillation, while
the computations using UWM show extensive numericaldispersion.
Another class of problems amenable to indirect solution
by the MSRPT is the initial value problem with a “stepfunction” type initial condition. Problem 2 again concerns
one-dimensional transport in a steady flow field, but now
with the initial and boundary conditions given byc(x, 0) = 1

1.2

PROBLEMS- Analyticalo MSRPT
A-A UpStWOm
i--s F E M

0.8
o

4135

0.6

OlX<X,c(x, 0) = 0x 2x1c(0, t) = 0

t>O

(39)
i
1

(40)

4000

I

.

8000

’

I

12000

X

1

Fig. 4. Comparison of concentration profiles computed with themodified single-step reverse particle tracking method, the finite
element method, the upstream weighted method, and the analyticalsolution for v = 0.5 and D = 1 (Pe = 100).

The exact solution for a semi-infinite profile is given by [vanGenuchten and Alves, 1982]
c(x, t) =f {erfc[‘:~~~);~‘]+ exp (vxlD) erfc [(~~~)~~t]}

0.8
0

0.6
D=2
Pe = 50

0

4000

8000

12000

X

1.2
1 .o

- Analyticalo MSRPT---- SRPT

0.8
v

0.6
0.4

Fig. 3. Concentration profiles computed with the modified single-step reverse particle tracking method, the single-step reverse
particle tracking method, and the analytical solution for u = 0.5,
R = 1, and (top) D = 0 (Pe = m), (middle) D = 2 (Pe = 50), and(bottom) D = 50 (Pe = 2).

(41)

The initial condition (39) represents a unit step function near
the inlet boundary. The center of the plume at time t is
located at x, = 0.5~~ + ut. MSRPT was used to only
compute F for nodes in the half region of the plume from thecenter to the extended area, i.e., between 0.5~~ + vt andx, + vt + 4(Dt)“2. The convective concentrations of thenodes were subsequently mirrored to the other half region,i.e., between vt - 4(Dt)‘12 and 0.5~~ + vt. Based on the
mirrored concentration values, the convective concentrationof the nodes in this region was calculated by linear interpolation. Using Ax = 200, At = 100, t = 9600, v = 0.5, D =0.2, R = 1, and x1 = 1200, we simulated solute movementin the same domain as problem 1, utilizing both SRPT andMSRPT. The results are compared in Figure 5 with theanalytical solution. Again, MSRPT matches the analyticalsolution almost exactly, whereas SRPT leads to serious
numerical dispersion.
Problem 3 involves the infiltration of water and a dissolved
solute in an initially solute-free soil profile. The problem wasdesigned to show the performance of MSRPT when simulating solute transport during transient unsaturated water flow.
For this example we used MSRPT to solve (1) subject to thesame initial and boundary conditions as Problem 1. Becauseof the transient flow conditions, 0 and v, and hence also D
and R in (1) are functions of time and space. At each time
step, f3 and v were obtained by numerically solving the
Richards equation; i.e.,(42)
where C is the soil water capacity, K is the hydraulicconductivity, and h is the pressure head. Equation (42) was4136

ZHANG ET AL.: EULERIAN-LAGRANGIAN METHOD FOR SOLVING SOLUTE TRANSPORT PROBLEMSe MSRPT

1.21- Analytical-.-. S R P T

o MSRPT: Ax = 1- FIX: Ax = 0 . 2
II FEM: Ax = 1

19

o

0.6

0

4000

8000

12000

X
X

Fig. 5. Concentration distributions computed with the modified
single-step reverse particle tracking method, the single-step reverse
particle tracking method, and the analytical solution for v = 0.5 andD = 0.2 (Pe = 500).

solved using a mass-lumped linear finite element scheme as
discussed by Kool and van Genuchten [1991]. The soil
hydraulic properties in this equation were given by [vanGenuchten, 1980]
0 - or[1 + (~lhl)Tm(43)

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

(44)

Se=-=

e _ e
s
I

where 0, and Bs are the residual and saturated watercontents, respectively, K, is the saturated hydraulic conductivity, (Y and n are model parameters determined by experimental data, and m = 1 - l/n. The dispersion coefficient in(1) was computed as [Bear, 1972]
D=A]v] +007

(45)

where A is the dispersivity of the medium, Do is the ionic ormolecular diffusion coefficient in free water, and T is atortuosity factor. The tortuosity factor was evaluated as afunction of the water content [Millington and Quirk, 1961] asfollows:$- = 0 10’3/(9,2Fig. 6. Solute transport simulations at times t = 2, 4, 8, and 12days for a transient flow field using the modified single-step reverse
particle tracking method (Ax = I), the finite element method with
Ax = 1, and the FEM method assuming a more refined grid
(Ax = 0.2).

element size Ax of 1 cm. However, a much larger element
size is possible with this scheme; this feature may further
increase the computational efficiency of MSRPT relative to
FEM. Figure 7 shows that some numerical dispersion occurred when the element size was increased. This numericaldispersion may have been caused by the process of choosing
the weighting factors according to the mass balance error.
The mass stored in the profile was calculated by summationof the element contributions. For smaller elements, the massbalance should be more accurate at each time step, which
would result in more appropriate values for the weighting
factors. Additionally, changing the element size likely affected the accuracy of the water flow simulation, which inturn should have impacted the accuracy of the solute transport simulation.
Transport simulations for a variably saturated flow field
were also carried out using a solute flux boundary conditionat the soil surface. The solute dispersivity was now taken as
0.01 cm. We used the simulated FEM results with the refinedelements of Ax = 0.05 cm (Pe = 5) as a standard for(46)

Calculations were obtained for a 300-cm-deep soil profile,an infiltration rate of 8 cm day-‘, and the following hydraulicparameters: n = 1.8, (Y = 0.1 cm-t, 0, = 0.05, 8, = 0.4, Do= 0, and K, = 10 cm day-‘. The concentration of the input
solution was set at unity. We selected an element size Ax of
1 cm and a relatively small value of 0.04 cm for thedispersivity A, such that Pe = Ax/h = 25. Numerical results
obtained with the MSRPT and FEM methods are shown in
Figure 6 at times t = 2, 4, 8, and 12 days. To eliminate orminimize oscillations and numerical dispersion in the FEM
results, we had to further refine the element size. Figure 6
also displays the FEM concentration profiles using therefined grid (Ax = 0.2 cm, Pe = 5). The profiles obtainedwith FEM using the smaller elements (Ax = 0.2 cm) and
MSRPT (Ax = 1 cm) are quite comparable; however, theFEM method took 25 times more computer time thanMSRPT (approximately 5 hours versus 12 min on a 486 PC).
The MSRPT results in Figure 6 were obtained with an

1.2- FEM: Ax = 0.2

1

0.0

1
0

4L

I

I

100

200

I

300

X

Fig. 7. Effect of element size on numerical results for a transient
flow field obtained with the modified single-step reverse particle
tracking method.

ZHANG ET AL.: EU L E R I A N - LAGRANGIAN M ETHOD FOR SOLVING SOLUTE TRANSPORT PR O B L E M STABLE 1. Mass Balance Errors (ME) of the Modified SingleStep Reverse Particle Tracking (MSRPT) Method and
Sum-Squared Errors (E) of the MSRPT and
Finite Element Methods (FEM)

Day(MSRPT)

Sum-Squared
Error, E(MSRPT)1
2
4
5
8-0.08970.0012-0.0318-0.0126-0.0343

1.20
0.958
1.45
2.85
12.5

Mass Error,
ME, %Sum-Squared
Error, E
(FEM)
3.15
6.58
9.47
12.0
73.7

4137

1.4
1.2

---. water content- Concentration

1 .0

All parameters used in the simulations were kept the same for thetwo methods.

comparison. The sum-squared errors of MSRPT and FEM
with Ax = 0.5 cm (Pe = 50) are compared in Table 1, which
also lists the mass balance errors of the MSRPT scheme atdifferent times. Notice that they are less than 0.1% for all
simulations. Calculated solute distributions obtained with
the MSRPT, FEM and the refined FEM techniques areshown in Figure 8 at times of 2, 4, and 8 days. The results inTable 1 and Figure 8 clearly demonstrate that MSRPT
provides for this example much more accurate results thanthe oscillation-prone FEM. While MSRPT yields distributions which compare closely with those simulated using therefined FEM, the MSRPT results required about 80 times
less computer time as compared to the refined FEM.
Finally, Figure 9 presents the concentration profiles for apurely convective transport process (h = 0) in a transient
flow field as simulated with the MSRPT method. As for theuniform flow field (Figure 3 (top)), MSRPT gives an exact
description of the sharp front without producing any numerical oscillation and artificial dispersion.
C ONCLUSIONS
A computationally efficient Eulerian-Lagrangian method
is presented for solving one-dimensional convectiondispersion problems in steady and transient flow fields.

1.21- FEM: Ax = 0.05
o MSRPT: Ax = 0.5* FEM: Ax = 0.5

Fig. 9. Concentration profiles for a purely convective transport
process simulated with the modified single-step reverse particle
tracking (MSRPT) method at times of t = 2, 4, and 8 days in atransient flow field.

The relatively simple algorithm is a modification of the
SRPT method. The method employs two weighting
factors to control the particle tracking process. One weighting factor greater than unity is used in the upstream
region, and another less than unity is used in the downstream
region. The two weighting factors approach unity fordispersion-dominated problems, in which case MSRPT becomes the traditional SRPT method. A relatively pragmatic
approach was followed by relating the weighting factors
to the grid Peclet and Courant numbers by means of empirical equations. The weighting factors for transient flow
fields were selected using an automatically adjustable procedure which minimizes mass balance errors at each time
step.
Several problems were used to compare the MSRPT
predictions against available analytical solutions and numerical results obtained with a linear finite element method
assuming a more refined grid system. Results indicate that
MSRPT is effective and accurate for a variety of problems.
The method can handle solute transport problems over theentire range of Peclet numbers form zero to infinity for both
steady and transient flow fields.
The MSRPT scheme was found to eliminate numericaloscillations and to reduce or virtually eliminate numericaldispersion in comparison to the finite element method, thesingle-step reverse particle tracking technique, and the upstream weighted method. MSRPT also dramatically reduced
computational times as compared to the finite element
method with a refined grid system. The MSRPT method
appears especially useful for convection-dominated transport problems. The algorithm in this paper was formulated
for one-dimensional problems and is currently limited to grid
Courant numbers less than or equal to 1. Further research is
being carried out to extend the MSRPT concept to problems
with larger Courant numbers and multidimensional flow
fields.
R EFERENCES

Fig. 8. Solute transport simulations at times t = 2,4, and 8 days
for a transient flow field using the modified single-step reverse
particle tracking method (Ax = 0.5), the finite element method with
Ax = 0.5, and the FEM method using a more refined grid system (Ax= 0.05).

Bear, J., Dynamics of Fluid in Porous Media, Elsevier, New York,1972.
Cheng, R. T., V. Casulli, and S. N. Milford, Eulerian-Lagrangiansolution of the convection-dispersion equation in natural coordinates, Water Resow. Res., 20, 944-952, 1984.

4138

Z HANG ET AL.: EULERIAN-LAGRANGIAN METHOD FOR SOLVING SOLUTE T RANSPORT P ROBLEMS

Dimou, K. N., and E. E. Adams, Representation of sources in a 3-D
Eulerian-Lagrangian mass transport model, in Water Pollution:
Modelling, Measuring and Prediction, edited by L. C. Wrobel
and C. A. Brebbia, pp. 251-264, Elsevier, New York, 1991.
Fujinawa, K. A., “Characteristic” finite element scheme for convective-dispersive transport with non-equilibrium reaction, Int. J.
Numer. Methods Eng., 23, 1161-1178, 1986.

Galeati, G., G. Gambolati, and S. P. Neuman, Coupled and partiallycoupled Eulerian-Lagrangian model of freshwater-seawater mixing, Water Resow. Res., 28(l), 149-165, 1992.
Garder, A. O., D. W. Peaceman, and A. L. Pozzi, Numericalcalculation of multidimensional miscible displacement by themethod of characteristics, Soc. Pet. Eng. J., 4, 26-36, 1964.
Heinrich, J. E., P. S. Huyakom, 0. C. Zienkiewicz, and A. R.
Mitchell, An “upwind” finite element scheme for two dimensional convective transport equation, Znt. J. Numer. MethodsEng., 11, 131-143, 1977.

Huang, K., R. Zhang, and M. T. van Genuchten, A simple particle
tracking technique for solving the convection-dispersion equation, in Computational Methods in Water Resources IX, vol. 1,
Numerical Methods in Water Resources, edited by T. F. Russell,
pp. 87-96, CMP, Boston, Mass., and Elsevier Science, New
York, 1992.
Khaleel, R., and D. L. Reddell, Miscible displacement in porous
media: MOC solution, J. Zrrig. Drain. Eng.,-1ll, 45-64, 1985.
Konikow, L. F.. and J. D. Bredehoeft. Computer model of 2-D
solute transport and dispersion in ground water, in Techniques ofWater Resources Investigations of the United States Geological
Survev. chap. C2. U.S. Government Printing Office, Washington.

D.C., 1978.‘
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, 101 pp., U.S. Salinity Lab.,Riverside, Calif., 1991.
Leith, C. E., Numerical simulation of the earth’s atmosphere, in
Methods on Computational Physics, vol. 4, pp. 1-28, Academic,
San Diego, Calif., 1965.

Leonard, B. P., The ULTIMATE conservative difference scheme
applied to unsteady one-dimensional advection, Comput. Methods Appl. Mech. Eng., 88, 17-74, 1991.
Millington, R. J., and J. M. Quirk, Permeability of porous solids,
Trans. Faraday Soc., 57, 1200-1207, 1961.

Neuman, S. P., Adaptive Eulerian-Lagrangian finite element
method for advection-dispersion, Znt. J. Numer. Methods Eng.,20, 321-337, 1984.
Neuman, S. P., and T. N. Narasimhan, Mixed explicit-implicit
iterative finite element scheme for diffusion type problems, I,
Theory, Znt. J. Numer. Methods Eng., 11, 309-323, 1977.
Neuman, S. P., and S. Sorek, Eulerian-Lagrangian methods foradvection-dispersion, in Finite Element in Water Resources,
edited by K. P. Holz et al., pp. 14.41-14.68, Springer-Verlag,
New York, 1982.
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. and W. J. Alves, Analytical solutions of the
one-dimensional convective-dispersive solute transport equation,
Tech. Bull. 1661, 151 pp., Agric. Res. Serv., U.S. Dep. of Agric.,
Washington, D. C., 1982.
Varoglu, E., and W. D. L. Finn, A finite element method for the
diffusion-convection equation with constant coefficients, Adv.
Water Resour., 1, 337-343, 1978.

Yeh, G. T., A Lagrangian-Eulerian method with zoomable hidden
fine-mesh approach to solving advection-dispersion equations,
Water Resour. Res., 26, 1133-1144, 1990.

K. Huang, M. T. van Genuchten, and R. Zhang, U.S. Salinity
Laboratory, Agricultural Research Service, U.S. Department ofAgriculture, 4500 Glenwood Drive, Riverside, CA 92501.
(Received December 14, 1992;
revised June 4, 1993;
accepted June 21, 1993.)