Transient seepage model for saturatedunsaturated soil systems: a geotechnical engineering approach
L. LAM,D. G. FREDLUND, S. L. BARBOUR AND
Department of Civil Engineering, University of Saskatchewan, Saskatoon, Sask., Canada S7N OW0
Received September 1 1, 1986 Accepted May 26, 1987
A twodimensional finite element model is proposed to simulate transient seepage for complex groundwater flow systems. The complete soil system is treated as a continuum encompassing flow in both saturated and uns
Transient seepage model for saturatedunsaturated soil systems: a geotechnicalengineering approach
L.
LAM,
D.
G.
FREDLUND,
ND
S. L. BARBOUR
Department of Civil Engineering, University of Saskatchewan, Saskatoon, Sask., Canada
S7N
OW0
Received September
1
1, 1986Accepted May
26,
1987
A
twodimensional finite element model is proposed to simulate transient seepage for complex groundwater flow systems.The complete soil system is treated as a continuum encompassing flow in both saturated and unsaturated zones. In the unsaturated zone, the air phase is assumed to be continuous and open to atmospheric pressure. The coefficient of permeability of theunsaturated soil is assumed to be a function of porewater pressure.The governing differential equation is derived within a framework familiar to geotechnical engineers. The stress statevariables and the constitutive relationships for an unsaturated soil are used in the derivation. The finite element solution to thegoverning differential equation is based on the Galerkin weightedresidual method. The nonlinearity of the equation is solvedby iterative procedures.The finite element formulation is implemented into a computer model named TRASEE. The model can be applied to a widevariety of problems involving complex boundary conditions and geometries with arbitrary degrees of heterogeneity andanisotropy. Example problems are presented to demonstrate the capabilities of the model. The results indicate that the quantityof water flow in the unsaturated zone may be substantial, and that the phreatic line is not a flow line. It has been found that thetraditional saturatedonly flownet technique can be approximated as a special case to the proposed saturated unsaturatedmodel.
Key words:
unsaturated flow, finite element model, phreatic line, permeability function, transient seepage.L'on propose un modble
2
deux dimensions en ClCments finis pour simuler l'infiltration transitoire pour des systbmescomplexes d'Ccoulement d'eau souterraine. Le systbme sol complet est trait6 comme un milieu continu d'Ccoulement dans lazone tant saturCe que non saturie. Dans la zone saturke, l'on suppose que la phase gazeuse est continue et ouverte
2
la pressionatmosphCrique. Le coefficient de permCabilitC du sol non saturC est considCrC comme Ctant une fonction de la pressioninterstitielle.L'Cquation diffkrentielle de base est familibre aux ingCnieurs gCotechniciens. Les variables de 1'Ctat de contrainte et les Cquations de comportement pour un sol non saturC sont utilisCes dans la dkrivation. La solution par ClCments finis de 1'Cquationdifferentielle de base est fondCe sur la mCthode Galerkin de rCsidu pondCrC. La non1inCaritC de 1'Cquation est risolue par desprocCdCs itCratifs.La formulation des ClCments finis est exCcutCe dans un modble d'ordinateur nommC TRASEE. Le modble peut &treapplique
B
une grande vari6tC de probbmes impliquant des conditions aux limites complexes et des gComCtries avec des degrCs arbitraires d'hCtCrogtnCitC et d'anisotropie. Des problbmes sont prCsentCs en exemples pour dCmontrer le potentiel de ce modble.Les rCsultats indiquent que la quantitC dlCcoulement d'eau dans la zone non saturCe peut &tresubstantielle, et que la lignephriatique n'est pas une ligne d'Ccoulement. I1 a CtC dCmontrC que la technique du rCseau d'Ccoulement saturC traditionnel peut&treconsidCrC approximativement comme un cas spCcial du modble proposC saturC

non saturC.
Mots clb
:
Ccoulement non saturi, modble en ClCments finis, ligne phrkatique, fonction de permCabilitC, infiltration transitoire.[Traduit par la revue]
Can.
Geotech.
J.
24,
565
580 (1987)
Introduction
In geotechnical engineering, seepage analysis may be ofinterest with respect to slope stability analysis, groundwatercontamination control, and the design of hydraulic structures.Since Casagrande's (1937) classic paper Seepage throughdams, seepage problems in geotechnical engineering havebeen generally solved by sketching flow nets. This method isbased on the assumption that water flows only in the saturatedzone. This method of solution is practical for simple steadystate problems where the boundary of the flow region is clearlydefined and the soil conditions are not too complex. However,many seepage problems of practical interest are complex, anda flownet solution is not feasible.With the development of highspeed digital computers,numerical methods are increasingly used in solving seepageproblems. Some of the first attempts at modelling seepageusing the finite element method also considered water flowonly in the saturated zone. An example is the finite elementflow model proposed by Taylor and Brown (1967). In this saturatedonly approach, the phreatic line was assumed tobe the upper boundary of the flow region for unconfined flowproblems. In order to obtain the location of the phreatic line, atrial and error procedure is used to search for a surface that hasboth zero water pressure along it and zero flow across it. Foreach trial, the location of the phreatic line is readjusted, and anew mesh is constructed. This method is tedious, and it is alsobased on the erroneous assumption that the phreatic line is theuppermost flow boundary.
A
saturated unsaturated model would appear to be superiorto the saturatedonly model. By considering water flow in boththe saturated and unsaturated zones, the actual ground surface(rather than the phreatic line) can be used as the uppermostboundary. The phreatic line is obtained by joining points ofzero water pressure in the flow region. The finite differencemodel proposed by Freeze (1971~)
s
an example of asaturated unsaturated flow model. This model allows thecalculation of flow and hydraulic heads in the unsaturatedzone. It also allows the modelling of situations with a surface
Printed in Canada
1
Imprirnf
au Canada
CAN.
GEOTECH.
J. VOL.
24,
1987
FIG.
1.
A
saturatedunsaturated soil profile under hydrostaticconditions.flux. As well, there is no need to adjust the mesh during thesolution of the problem.Much of the research concerning unsaturated flow has beenundertaken by soil scientists. Geotechnical engineers haveoften accepted formulations from researchers with diversebackgrounds without evaluating the assumptions and objectives of their formulation. This paper presents the developmentof a transient finite element seepage model using an approachconsistent with that used by geotechnical engineers. The concepts pertinent to the understanding of unsaturated flow arereviewed. The general governing differential equation fortransient seepage is derived and the finite element solution tothe equation is presented.
Review of unsaturated flow concepts
Unsaturated flow is an example of multiphase flow through aporous media. Two phases, air and water, coexist in the porechannels. Two partial differential equations are required torigorously describe the flow of air and water in unsaturatedsoil. Problems of this kind are referred to as twophase flowproblems (Fredlund 1981). However, for most seepage problems encountered in geotechnical engineering, only flow in thewater phase of the unsaturated zone is of practical interest.This singlephase flow approach is used in this paper. Thefundamental assumption is that the air phase is continuous andthe poreair pressure is equal to atmospheric pressure. This isgenerally the case when the degree of saturation is less thanapproximately 85
%
.
In the case where the degree of saturationis larger than 85%, the air phase is occluded. However, theoccluded air pressure will still be essentially atmospheric andwill not introduce significant error into the singlephase flow approach. Freeze and Cheny (1979) stated that thesinglephase approach to unsaturated flow leads to techniquesof analysis that are accurate enough for almost all practicalpurposes.
Soils are porous materials, with water being free to flowthrough the interconnected pores between the solid particles.Regardless of the degree of saturation of the soil, water flowsunder the influence of a hydraulic gradient. The drivinggradient can be identified by Bernoulli's theorem. Since seepage velocities in soils are normally small, the velocity head canbe neglected, and the driving force is the sum of porewaterpressure head and elevation head:
u
0
I
4
I
I I I
*
10 20 30 40 50 60 70 80
SUCTION
(
kPa)
FIG.
2.
Moisture retention curve for fine sand, silt,
and
Reginaclay.where
h
=
total head,
u,
=
porewater pressure,
z
=
elevationhead above an arbitrary datum, p
=
density of water, and
g
=
acceleration due to gravity.The flow of water through a soil mass can be described usingDarcy's law:where
q
=
discharge per unit area,
i
=
potential gradient, and
k
=
coefficient of permeability (or hydraulic conductivity).Darcy's law, which was srcinally derived for saturated soil,can equally be applied to the flow of water through an unsaturated soil (Richards 1931
;
Childs and CollisGeorge 1950).The only difference is that for flow through an unsaturatedsoil, the coefficient of permeability is no longer a constant butis a function of the matric suction of the soil.Figure 1 illustrates the porewater pressure and moisture distribution within a soil profile under hydrostatic conditions. Thewater table (i.e., phreatic line) is an imaginary surface wherethe porewater pressure is zero. Water may flow up and downacross the water table if there is a hydraulic gradient betweenthe saturated and unsaturated zones. Below the water table, theporewater pressure is positive and increases linearly withdepth. Above the water table, the porewater pressure is negative and decreases with height. The amount of water retainedin the unsaturated zone of a particular soil is a function ofmatric suction
(u,

u,)
of the soil (Fredlund and Morgenstern1977). Since the poreair pressure is atmospheric, the matricsuction is also equal to the absolute value of the negative porewater pressure.Figure
2
shows the laboratory retention curve for fine sand,silt, and Regina clay obtained by Ho (1979). These retentioncurves were obtained using Tempe pressure cells. The operating principles and the testing procedures can be found in the Operating instructions for Tempe pressure cell by SoilMoisture Equipment Corporation, Santa Barbara, California.In order to consider flow in both the saturated and unsaturatedzones of a flow regime, the permeability of the soil at all pointsin both the saturated and unsaturated zones must be known.The saturated coefficient of permeability can be obtained usingstandard laboratory and field procedures. However, the directmeasurement of the unsaturated permeability is often expen
LAM
ET
AL.
567
16f2Lill
6
10
20
3040
50
60 70 80 90 100
SUCTION
(kPa)
FIG.
3.
Unsaturated permeability function of fine sand, silt, andRegina clay.sive and difficult to conduct. Considerable research has beendone on evaluating the unsaturated permeability of a soil basedupon its soil moisture retention curve, which is relativelysimple to obtain (Childs and CollisGeorge 1950; Marshall1958; Green and Corey 1971; Elzeftawy and Cartwright1981). It has been concluded that this method generally givesreliable predictions of experimentally measured permeabilityvalues.The moisture retention curve of a soil is of primary importance in understanding the transient flow of water in theunsaturated zone. The slope of the curve represents the storagecharacteristics of the soil (i.e., the
qW
alue in [9]). The slopeindicates the amount of water taken on or released by the soilas a result of a change in the porewater pressure. The moistureretention curve of a soil also provides a quantitative means ofcalculating the permeability function of the soil usingprocedures developed by Green and Corey (1971) and otherresearchers.The retention curves and the calculated permeability functions as presented in Figs.
2
and
3
are typical for materials ofdifferent grain sizes. A coarser material(e.g., sand) willdesaturate faster than a finer material (e.g., clay). Consequently, as anticipated, a sandy material will have a steeperpermeability function than a clayey material.
General governing equation
Terzaghi's (1943) onedimensional consolidation theory forsaturated soils introduced an approach to seepage problemsthat was particularly suitable to geotechnical engineering. Thegoverning partial differential equation was derived by equatingthe net flow quantity from a soil element to the time derivativeof the constitutive equation for the soil.Fredlund and Morgenstern (1977) proposed the use of twoindependent stress state variables, (a

u,) and (u,

u,), todescribe the state of stress for an unsaturated soil. Constitutiveequations relating the volume change in the soil structure andfluid phases to stress state variable changes were also proposed(Fredlund and Morgenstern 1976).The net flow through a twodimensional element of unsaturated soil can be expressed asThe constitutive equation for the water phase of an isotropicunsaturated soil is[4]dB,
=
myd(u

u,)
+
m;d(u,

u,)where my
=
slope of the (u

u,) versus 8, plot when d(u,

u,) is zero, m,
=
slope of the (u,

u,) versus 8, plot whend(u

u,) is zero, a
=
total stress in thex and (or) ydirection,u,
=
poreair pressure, and u,
=
porewater pressure.Since my and mcan be assumed to be constant for aparticular time step during the transient process, the timederivative of the constitutive equation can be expressed asCombining [3] and [5] givesIf it is assumed that no external loads are added to the soilmass during the transient process, and that the air phase is continuous in the unsaturated zone (i.e., aulat
=
0 and au,lat
=
O), [6] can be simplified as follows:Expressing the porewater pressure term in [7] in terms oftotal head, the governing differential equation for transientseepage can therefore be written as follows:After incorporating anisotropic soil conditions where thedirection of the major coefficient of permeability is inclined atan arbitrary angle to the xaxis, [8] becomeswhere k
=
k, cos2
a
+
k2 sin2
a,
kyy
=
kl sin2
a
+
k2 cos2
a,
k
=
k
=
(k,

k2) sin
a
cos
a,
k,
=
major coefficient of
YX
permeability
,
k2
=
minor coefficient of permeability, and
a
=
inclined angle between kl and the xaxis.The above governing partial differential flow equation (i.e.,[9] has been developed in a manner most familiar to geotechnical engineers. The equation has been derived on the basis ofunsaturated flow theory (Fredlund 1981) and can be used todescribe continuous flow in saturatedunsaturated soil systems. The termm, represents the rate at which a soil willabsorb or release water when there is a change in matric suc
CAN.
GEOTECH.
J.
VOL.
24.
1987
VOLUMETRICWATERCONTENT
8
POROSITY
with
Li
qual to the area coordinates of the element;
X
=
PgmW
q
=
flow across the perimeter of the element
A
=
area of the element
S
=
perimeter of the element
t
=
time
d8=
my
d(uu,)
.'
do
=
my
~(WU.)+ :d(u.~.)
After numerical integration,
[lo]
can be simplified as
m
u
follows:
SATURATED
'+UNSATURATED
0.1..
I
[l
]
[Dl[hn]
+
[El
[R]
=
[F]
ti
where
+
va
WATER PRESSURE

vs
WATER PRESSURE
(
SUCTION
1
is the stiffness matrix;
FIG.
4.
Moisture retention curve and
m,
or saturatedunsaturatedsoil.tion. Figure
4
illustrates the variation in the value of the term
mW
as the soil changes from saturated to unsaturated condi
is
the
capacitance matrix;tions. The value of
mW
is equal to the slope of the water retention curve. Under the assumptions of a constant total stress anda poreair pressure equal to the porewater pressure,
mW
becomes equal to
my
in the saturated zone. The
mp
value isequivalent io the coeificient of volume change,
m,,
cbmmon tois the flux vector reflecting the boundary conditions;saturated soil mechanics.
Finite element formulation
With the advance of highspeed digital computers, numerical methods provide practical solutions to complex seepageproblems. The finite element method is the most powerful andwidely used numerical method for twodimensional problems.The following section presents the finite element formulationfor the governing equation.The Galerkin solution to
[9]
s given by the following integrals over the area and boundary surface of a triangular element (Lam
1983):
where
Y2

37Y3

Y1,Yl

2
1
[B]
=

A
6.
X27X1

X3,X2

X1
the element;
I
ith
x,,
y,
equal to the Cartesian coordinates of the nodes ofwith
k,, k,, k,,, kyy
equal to the components of the permeability tensor for the element;with
hi
equal to the total head at the element nodes;


is the time derivative of nodal total head.For transient seepage, the time derivative of
[l
.]
can beapproximated by a finite difference procedure. Consequently,the relationship between the nodal heads of an element in twosuccessive time steps can be expressed by the following equations:Equation
[12]
s derived using the central difference approximation, and
[13]
is derived using the backward differenceapproximation. Both of the above time approximations areconsidered to be unconditionally stable. Generally, solutionsbased on the central difference approximation are found to bemore accurate than those of the backward difference approximation. However, the backward difference approximation isfound to be more effective in dampening the numerical oscillations frequently encountered in highly nonlinear flow systems(Neuman and Witherspoon
197
1
;
Neuman
1973).
After the matrices for each element are formed, the algebraicequations for the whole system can be constructed and solvedfor nodal total head. However, due to the nonlinearity of thegeneral seepage equation, an iterative procedure is required toobtain the correct nodal total heads. The iterative procedureinvolves a series of successive approximations. An initial estimate of the coefficient of permeability of an element isrequired in order to calculate a first approximation of the nodaltotal head. The computed nodal total head allows the calculation of the average pressure head of an element. Using thecalculated average pressure heads and the input permeabilityfunction, an improved permeability value can be obtained forthe element. The improved permeability value is used to com