PLAXIS LE

PLAXIS LE
Groundwater
1D/2D/3D SATURATED / UNSATURATED
FINITE ELEMENT GROUNDWATER MODELING
Theory Manual
Written by:
The Bentley Systems Team
Last Updated: Thursday, December 16, 2021
Bentley Systems Incorporated

COPYRIGHT NOTICE
Copyright © 2021, Plaxis bv, Bentley Systems, Incorporated. All Rights Reserved.
Including software, file formats, and audiovisual displays; may only be used pursuant to applicable software license
agreement; contains confidential and proprietary information of Plaxis bv and/or third parties, including Bentley Systems,
Incorporated, which is protected by copyright and trade secret law and may not be provided or otherwise made available
without proper authorization.
Copyright PLAXIS program by: Plaxis bv P.O. Box 572, 2600 AN DELFT, Netherlands Fax: +31 (0)15 257 3107; Internet site:
www.bentley.com
These manuals may not be reproduced, in whole or in part, by photo-copy or print or any other means, without written
permission from Plaxis bv
TRADEMARK NOTICE
Bentley, "B" Bentley logo, SoilVision logo, and SOILVISION, SVSLOPE, SVOFFICE, SVOFFICE 5/GE, SVOFFICE 5/GT, SVOFFICE
5/WR, SVSOILS, SVFLUX, SVSOLID, SVCHEM, SVAIR, SVHEAT, SVSEISMIC and SVDESIGNER are either registered or
unregistered trademarks or service marks of Bentley Systems, Incorporated. All other marks are the property of their
respective owners.
BENTLEY SYSTEMS Table of Contents Page 3 of 53
1 INTRODUCTION .............................................................................................................................................................. 5
2 FUNDAMENTALS OF SEEPAGE THEORY........................................................................................................................... 6
2.1 CONSERVATION OF MASS ..................................................................................................................................................... 7
2.1.1 Changes in Volume of Stored Water .......................................................................................................................................................................7
2.2 FLOW LAWS ....................................................................................................................................................................... 8
2.2.1 Flow of Liquid Water ..............................................................................................................................................................................................8
2.2.2 Anisotropic Flow of Liquid Water............................................................................................................................................................................9
3 PDE’S FOR SEEPAGE ANALYSIS...................................................................................................................................... 11

3.1 ONE-DIMENSIONAL SEEPAGE ............................................................................................................................................... 11
3.1.1 1D Seepage...........................................................................................................................................................................................................11
3.2 TWO-DIMENSIONAL SEEPAGE .............................................................................................................................................. 11
3.2.1 2D Seepage...........................................................................................................................................................................................................11
3.3 THREE-DIMENSIONAL SEEPAGE ............................................................................................................................................ 12
3.3.1 3D Seepage...........................................................................................................................................................................................................12
3.4 UNSATURATED FLOW EQUATION.......................................................................................................................................... 13
3.4.1 H-Based Formulation ............................................................................................................................................................................................13
4 BOUNDARY CONDITIONS IN SEEPAGE .......................................................................................................................... 14

4.1 NATURAL BOUNDARY CONDITION ......................................................................................................................................... 14
4.1.1 Differences Between Actual and Applied Boundary Flux ......................................................................................................................................14
4.2 ESSENTIAL BOUNDARY CONDITION ....................................................................................................................................... 14
4.2.1 Review Boundary Condition (Drain) ......................................................................................................................................................................15
4.3 GRADIENT OR UNIT GRADIENT BOUNDARY CONDITION ............................................................................................................ 16
4.4 INITIAL WATER TABLE ........................................................................................................................................................ 16
5 FLUX SECTIONS ............................................................................................................................................................. 18
5.1 FLUX SECTION THEORY CALCULATION.................................................................................................................................... 18
6 MATERIAL PROPERTIES ................................................................................................................................................ 21
6.1 SOIL-WATER CHARACTERISTIC CURVE ................................................................................................................................... 21
6.1.1 Fredlund and Xing (1994) Equation ......................................................................................................................................................................21
6.1.2 Fredlund (2000) Bimodal Equation .......................................................................................................................................................................21
6.1.3 Van Genuchten (1980) Equation ...........................................................................................................................................................................22
6.1.4 Van Genuchten (1980) and Mualem (1976) Equation ..........................................................................................................................................22
6.1.5 Gardner (1958) Equation ......................................................................................................................................................................................23
6.1.6 Brooks and Corey (1964) Equation .......................................................................................................................................................................23
6.1.7 Zapata Estimation (2000) .....................................................................................................................................................................................24
6.2 HYDRAULIC CONDUCTIVITY FIT METHODS .............................................................................................................................. 24
6.2.1 Gardner’s (1958) Equation ...................................................................................................................................................................................24
6.3 HYDRAULIC CONDUCTIVITY ESTIMATION ................................................................................................................................ 26
6.3.1 Brooks and Corey Estimation ................................................................................................................................................................................26
6.3.2 Modified Campbell (1996) Estimation ..................................................................................................................................................................27
6.3.3 Fredlund, Xing and Huang (1994) Estimation .......................................................................................................................................................28
6.3.4 Van Genuchten (1980) and Mualem (1976) Estimation .......................................................................................................................................29
6.3.5 Van Genuchten (1980) Estimation ........................................................................................................................................................................29
6.3.6 Leong and Rahardjo (1997) Estimation ................................................................................................................................................................30
7 SOIL ATMOSPHERE MODELING .................................................................................................................................... 32

7.1 ATMOSPHERIC FLUX BALANCE ............................................................................................................................................. 32
7.2 INFILTRATION AND RUNOFF WITHOUT CONSIDERATION OF EVAPORATION ................................................................................... 33
7.2.1 Pressure Head Calculated Runoff Correction ........................................................................................................................................................33
7.3 EVAPORATION .................................................................................................................................................................. 34
7.3.1 Potential Evaporation ...........................................................................................................................................................................................34
7.3.1.1 Penman (1948) Equation ....................................................................................................................................................34
7.3.1.2 Thornthwaite (1948) Equation ...........................................................................................................................................35
7.3.2 Actual Evaporation ...............................................................................................................................................................................................36
7.3.3 Case 1 Solution: “Limiting-Function” for AE/PE Model ........................................................................................................................................36
7.3.3.1 Mass Flow Governing Equation ...........................................................................................................................................36
7.3.3.2 Initial soil water contents ....................................................................................................................................................37
BENTLEY SYSTEMS Table of Contents Page 4 of 53
7.3.3.3 Boundary condition for moisture flow ................................................................................................................................37

7.3.3.4 Actual Evaporation ..............................................................................................................................................................38
7.3.3.5 Determining vapor pressure values throughout each day ..................................................................................................38
7.3.4 Case 2 Solution: Wilson-Penman (1994) Model ...................................................................................................................................................41
7.3.4.1 Governing Equations ...........................................................................................................................................................42
7.3.4.2 Initial soil water content profile ..........................................................................................................................................42
7.3.4.3 Initial soil surface temperature ...........................................................................................................................................42
7.3.4.5 Actual Evaporation ..............................................................................................................................................................43
7.3.4.6 Determining Vapor Pressure Values throughout each Day .................................................................................................43
7.3.5 Case 3 Solution: “Experimental-Based” for AE/PE Model ....................................................................................................................................44
7.3.5.1 Governing Equation .............................................................................................................................................................44
7.3.5.2 Initial water content conditions ..........................................................................................................................................44
7.3.5.4 Determine an empirical expression for the ratio of actual evaporation, AE, to potential evaporation, PE ........................44
8 STREAM TRACES ........................................................................................................................................................... 47

9 PARTICLE TRACKING ..................................................................................................................................................... 48
10 NUMERICAL IMPLEMENTATION ............................................................................................................................... 49
10.1 INTERPOLATION ORDER ...................................................................................................................................................... 49
10.2 TIME STEPPING ................................................................................................................................................................. 49
10.2.1 Time Integration ..............................................................................................................................................................................................49
10.2.2 Convergence Criteria .......................................................................................................................................................................................50
10.2.3 Time-Step Control ............................................................................................................................................................................................50
11 REFERENCES ............................................................................................................................................................. 51
BENTLEY SYSTEMS Introduction 5 of 53
1 INTRODUCTION
Groundwater flow problems in geotechnical and geo-environmental engineering involve the solution of a partial differential
equation referred to as a PDE. The PDE must be solved for all “finite elements” which when combined form a “continuum” (or
the geometry of the problem). The theory of groundwater flow expressed in mathematical form embraces the physical behavior
of the material (e.g., constitutive laws) and the conservative laws of physics (i.e., conservation of energy). The physical behavior
of many materials, (particularly unsaturated soils), is nonlinear and as a consequence, the PDE becomes nonlinear in character.
It is well known that the solution of nonlinear PDEs can present a challenge for numerical modeling.
The purpose of the theory manual is to provide the user with details regarding the theoretical formulation of the PDE as well as
the numerical method used in the solution. The intent of the theory manual is not to provide an exhaustive summary of all
theories associated with groundwater flow. Rather, the intent is to clearly describe details of the theory used in the PLAXIS LE
- Groundwater software.
The generic finite element solver solves the partial differential equation for groundwater flow. The solver algorithm has
implemented cutting-edge numerical solution techniques that can accommodate linear and highly nonlinear PDEs. The solution
technique utilizes adaptive time steps algorithm and automatic mathematically designed mesh generation. The application of
these advanced numerical techniques is particularly valuable in solving highly nonlinear and complex problems. Most commonly
it is the unsaturated soil portion of the soil continuum that brings in nonlinear soil behavior. The advanced solvers make it
possible to obtain converged and accurate solutions for many problems that were previously unsolvable.
The primary attributes of the solution process are as follows:
• Fully automatic mesh generation,

• Integrated climatic calculation of actual evaporation rates
• Handling of seepage face boundary conditions,
• Fully implicit approach in the solver, which provides for a robust solution of difficult models with convergence
issues,
• Fluid mass-balance tracking,
• 3-noded triangle and 4-noded quadrilateral elements for 2D analysis and 4-noded tetrahedron elements for 3D
analysis,
• Adaptive time stepping
• Automatic generation and control of time steps,
• Newton-Raphson convergence iteration schemes,
• The default preconditioner is the diagonal-block inverse matrix, and
BENTLEY SYSTEMS Fundamentals of Seepage Theory 6 of 53
2 FUNDAMENTALS OF SEEPAGE THEORY

The following section presents the fundamentals of the seepage theory implemented in PLAXIS LE - Groundwater. An overview
of the PDE’s used by the Groundwater module for the analysis of saturated/unsaturated seepage is presented. Several
phenomena are considered, including liquid water flow and water vapor flow. The seepage theory is developed based on Darcian
flow law, Fick’s flow law, and the conservation of mass. The theory presented below is a general outline of the theory needed
to solve most seepage problems.
Continuum mechanics principles and partial differential equations (PDEs) have been traditionally used for modeling seepage in
saturated/unsaturated soil systems. The partial differential equations governing seepage may involve transient coupled soil-
atmosphere processes with nonlinear and heterogeneous soil properties along with nonlinear boundary conditions. Relatively
simple steady-state saturated confined flow can also be addressed.
Seepage can be modeled as follows within the context of continuum mechanics principles:
• Identify the physical processes of concern associated with the problem at hand,
• Establish the “continuum variables” acting upon a representative elemental volume (REV) of the medium,
• Develop field equations governing the physical processes of concern by making the assumption that the medium
can be considered as a continuum from a macroscopic standpoint (i.e., considering a REV of soil) while using
measurable soil properties:
o Apply conservation laws,

o Apply verified constitutive laws, and
o Develop a final system of well-posed determinate partial differential equations.
• Establish initial, internal, and boundary conditions for the problem, and
• Provide a mathematical solution for the PDE or system of PDEs.
A series of assumptions form the backdrop for the derivation of the partial differential equations governing seepage. The
following set of assumptions can be considered generally valid:
• Soil phases can be described using a continuum mechanics approach,

• Pore-air and all of its constituents (including water vapor) behave as ideal gases,
• Local thermodynamic equilibrium between the liquid water and water vapor phases exists at all times at any point
in the soil, and
• Atmospheric pressure gradients are negligible.
The four assumptions described above may become inadequate under certain situations. For instance, the compressibility of
water may have a significant effect in an analysis of regional groundwater systems (i.e., large domains). Therefore, the software
provides the capability of taking the compressibility of water into account as part of aquifer storativity (Freeze et al. 1979).
2.1 CONSERVATION OF MASS

The conservation of mass of water in a referential element is used to derive the governing equation for saturated/unsaturated
seepage. A continuum mechanics framework is used, resulting in the derivation of a differential calculus equation to represent
seepage. The assumption is made that the variables involved are continuous and valid from a macroscopic, phenomenological
standpoint.
A differential equation for the conservation of mass of water can be derived by considering a REV of soil (Figure 1). The continuity
equation can be applied by taking into consideration the flow rates in and out of the REV and equating the difference to the rate
of change of mass (or heat) to storage within the REV with time. The following differential equation is obtained by considering
three-dimensional flow conditions using the Cartesian coordinate system:
qxw q w
y qzw 1 M w [1]
− − − =
x y z Vo t
where:
= total water flow rate in the i-direction across a unit area of the soil, kg/m2-
qi w
s; qiw = wviw, kg/m2-s,
w = density of water, 1000 kg/m3,

viw = water and air flow rate in the i-direction across a unit area of the soil, m/s,
Vo = referential volume, Vo = dx dy dz , m3,
Mw = mass of water within the representative elemental volume, kg, and
t = time, s.
The total water flow rate, vw, also known as specific discharge, is a macroscopic measure of the rate of flow through soils. A
measure of the “average actual flow velocity” for a saturated soil can be obtained by dividing vw by the soil porosity (n = Vv
/V). The total water flow rate, vw, can occur as liquid water and/or water vapor flow.
q y
qy + dy
y q z
qz + dz
y z
qx
z q x
dy qx + dx
qz x
dz
x
O qy
qx, qy, qz – rate of flow of mass
dx of air, mass of water, or heat
Figure 1 Soil representative elemental volume and fluxes at the element faces
2.1.1 Changes in Volume of Stored Water

The constitutive relationship for the amount of water stored in the soil pores is usually written in terms of the volume of water
referenced to the overall total volume. The change in volume of water stored in the soil pores can be expressed as a coefficient
of water storage, m2w, as follows:
dVw
= m2w d ( ua − uw ) [2]
Vo
where:
m2w = d (Vw / Vo ) = e dS
d (ua − uw ) 1 + e d (ua − uw )
Vw /Vo = volumetric water content,

e = void ratio,
S = degree of saturation, and

(ua – uw) = matric suction.
The above equation is based on the assumption that changes in the volume of pore-water stored in the soil are a function of
soil suction and are independent of changes in total stress. The soil property, m2w, is obtained by taking the derivative of the
soil-water characteristic curve, swcc, (i.e., the slope), as shown in Figure 2.
Using the derivative of the SWCC provides a smooth transition between saturated and unsaturated conditions, provided that
appropriate coefficients of water storage are used. As the soil saturates, the effects of changes in soil suction and changes in
effective stresses (in a saturated soil) become equal (i.e., m2w = mv). Consequently, for saturated conditions, changes in water
volume can be referenced to changes in void ratio. Usually, the value of mv is considerably lower than the maximum value of
m2w, which occurs as the soil desaturates.
Numerical difficulties can arise from the use of extremely low values of m2w. This is a possible source of convergence problems
when modeling ground surface infiltration problems. In order to alleviate convergent difficulties, the value of mv must be slightly
(but not excessively) raised as illustrated in Figure 2.
Figure 2 Soil-water characteristic curve showing the water storage characterization at low suction values.
2.2 FLOW LAWS

Table 1 presents an overview of flow laws traditionally used for modeling saturated/unsaturated soil flow. The flow laws establish
relationships between measures of flow and driving potentials. Driving potentials can be established based on spatial gradients
of the energy stored per unit volume (Bear, 1972). There are several flow equations that have a similar form, but have distinct
potentials and different material properties. The flow laws presented in Table 1 have well established equations that have been
experimentally verified.
Table 1 Overview of types of flow within an unsaturated soil and the corresponding mechanisms, driving potentials, and flow laws
Flow mechanism Driving Potential Flow Law
(1) (2) (3)
Liquid water Hydraulic head, h (m) Darcy’s law
Mass concentration of vapor per unit
Water vapor diffusion Modified Fick’s law
volume of soil, Cv (kg/m3)
Pore-air and pore-water have both miscible and immiscible mixture characteristics. Water can flow as liquid water or as water
vapor diffusing through the free air-phase. PLAXIS LE - Groundwater takes into account both liquid and vapor flow. Both flow
mechanisms are essential in the modeling of certain water flow conditions. For instance, evaporation requires the consideration
of the phase change from liquid to vapor water and the flow of water vapor (Wilson et al., 1994).
2.2.1 Flow of Liquid Water

The rate of flow of liquid water in a saturated/unsaturated soil can be described using a generalization of Darcy’s Law (Bear,
1972), where the driving mechanism is the total hydraulic head gradient. The hydraulic conductivity is assumed to mainly vary
with matric suction in the soil. The generalized Darcy’s law can be written as follows:
h
v wx = − k wx ( )
h ; vwy = −k wy ( ) ; v wz = − k wz ( )
h [3]
x y z
where:
vwi = liquid pore-water flow rate in the i-direction across a unit area of the soil
due to hydraulic head gradients, m/s,
= hydraulic conductivity in the i-direction, m/s. For unsaturated soil, it is the
kwi()
function of soil suction.
= soil suction, kPa, which is equal to matric suction, ua – uw, plus osmotic

suction, ,
h = hydraulic head, m, which is equal to uw + y

w
uw = pore-water pressure, kPa,

ua = pore-air (gauge) pressure, kPa,
w = unit weight of water, kN/m3, and
y = elevation, m.
The hydraulic conductivity function, kw() provides the relationship between the hydraulic conductivity and the soil suction. The
hydraulic conductivity function can also be written in terms of volumetric water content. As a soil dries, there is less and less
water present in the soil. Since water will flow only where there is water present, the hydraulic conductivity decreases
accordingly as the volumetric water content decreases. This behavior is represented in PLAXIS LE - Groundwater by entering a
hydraulic conductivity function for each soil. In the following context, the kwi() is simplified as kw.
The use of a continuous kw function provides a smooth transition between saturated and unsaturated soil conditions. For
saturated conditions, kwsat is generally considered a constant and equal to the saturated hydraulic conductivity.
The hydraulic conductivity function can be obtained experimentally using laboratory tests or field tests. The hydraulic
conductivity function can also be estimated using the saturated hydraulic conductivity and the soil-water characteristic curve
(Fredlund and Xing, 1994). PLAXIS LE - Groundwater provides several options for estimating the hydraulic conductivity function.
2.2.2 Anisotropic Flow of Liquid Water

Natural deposition of soil layers often results in angled layering. The angled layering results in dominant water flow in a direction
parallel to the direction of the layering. This observed phenomenon can be simulated in PLAXIS LE - Groundwater by specifying
anisotropic angled behavior in either 2D or 3D.
Figure 3 illustrates how angled anisotropy is considered. The angle origin, reference orientation, is in the horizontal direction
and the angle increases counter-clockwise.
A mathematical transformation must be performed to translate rotated anisotropic permeability parameters on to the coordinate
system used to solve the problem. The general form of the conductivity matrix in two-dimensional problems is as follows (Bear,
1972):
 h 
  − k wxx
vwx  − k wxy   x 
 [4]
v  = − k   h 
− k wyy   

 wy 
  wyx

 y 

where:
kwxx = kw1 + kw2 + kw1 − kw2 cos 2

2 2
k +k k −k
kwyy = w1 w2 + w1 w2 cos 2
2 2
kwxy = kw1 − kw2 sin 2

2
kw1 and kw2 are the values of hydraulic conductivity at the principal directions of anisotropy. These values can have constant
values for saturated conditions or can vary according to water content in unsaturated soil conditions. The angle  defines the
principal directions of anisotropy (Figure 3).
The transformation in 3D is slightly more complex than in 2D. Three angles are required to describe the hydraulic conductivity
transformation. The theory for the hydraulic conductivity transformations is not presented in this manual.
Figure 3 Principal directions of anisotropy in two-dimensions

BENTLEY SYSTEMS PDE’s For Seepage Analysis 11 of 53
3 PDE’S FOR SEEPAGE ANALYSIS

The flow law equations (Darcy law and Fick’s law), and a water volume change constitutive equation must be combined with
the continuity of water mass equation to obtain the partial differential equation that governs the conservation of water mass
(i.e., both liquid and vapor). The next sections present the partial differential equations used by PLAXIS LE - Groundwater for
the analysis of saturated-unsaturated soil seepage problems.
3.1 ONE-DIMENSIONAL SEEPAGE

One-dimensional seepage analysis is often used for modeling large planar areas. Cover system design is one of many typical
applications. The number of nodes normally used in one-dimensional analyses is relatively low when compared to the number
of nodes used in two-dimensional analyses. Therefore, computation times are dramatically reduced when adopting a one-
dimensional model.
3.1.1 1D Seepage
Considering the reference volume, Vo, and assuming that water is incompressible, the following equation is obtained for one-
dimensional transient saturated/unsaturated seepage:
  h u w  w h
k wy + kvd  = − w m2 [5]
y  y y  t
where:
y = coordinate in vertical direction, (corresponding to elevation).
The PDE governing the flow and storage of water within a saturated/unsaturated soil is presented using total head, h, as the
primary variable. However, pore-water pressure, uw, can also be used producing identical results provided the geometry
dimension in the y-coordinate is small. It is also possible to use pore-water pressure as the primary variable when studying the
dissipation of excess pore-water pressures.
Three soil property functions can be identified for the transient seepage PDE; namely:
• hydraulic conductivity function, kw,

• vapor conductivity function, kvd, and
• soil-water characteristic curve, whose derivative with respect to soil suction is represented by m2w.
The above-mentioned soil properties functions vary with soil suction. Therefore, the PDE is physically nonlinear.
The partial differential equation for water flow is based on the assumption that the rate of water mass flow across a REV is
continuously distributed in space. Therefore, the spatial distribution of water flow rate can be described using the partial
derivative of water flow in a particular direction. These comments apply to the all the seepage PDEs presented in the next
sections.
For steady-state conditions, the PDE for liquid and vapor water flow reduces to the following equation:
  h u 
k wy + kvd w  = 0 [6]
y  y y 
The hydraulic conductivity can be considered as being constant when solving saturated seepage problems. Therefore, the PDE
for saturated and unsaturated water seepage has the same form.
3.2 TWO-DIMENSIONAL SEEPAGE

Two-dimensional seepage analysis is used for modeling cross-sections passed through elongated geometries, such as an earth
dam.
3.2.1 2D Seepage
Assuming the reference volume, Vo, remains constant and the water is incompressible, the following equation can be written
for transient saturated/unsaturated seepage:
  h u    h u  h
+ kvd w  + k wy + kvd w  = − wm2w [7]
x 
k wx
x x  y  y y  t
where:
x = horizontal direction, and
y = vertical direction, (corresponding to elevation).
The partial differential equation, PDE, is presented for anisotropic properties with the principal direction of anisotropy
corresponding to the x- and y-directions. Anisotropic material properties that do not coincide with the Cartesian coordinate axis
can be considered as presented in the previous chapter.
For steady-state conditions, the water storage portion of the equation is set to zero and the PDE reduces to the following
equation:
  h u    h u 
+ kvd w  + k wy + kvd w  = 0 [8]
x 
k wx
x x  y  y y 
The governing PDE for steady state seepage can be further reduced by assuming vapor flow is negligible and soil is saturated.
The resulting PDE can be considered to be in its simplest form for two-dimensional seepage.
  h    h 
k wx  + k wy  = 0 [9]
x  x  y  y 
3.3 THREE-DIMENSIONAL SEEPAGE

Three-dimensional seepage analysis is often used when the problem geometry cannot be adequately represented by a one-
dimensional column, a two-dimensional, a plan, or an axisymmetric model.
3.3.1 3D Seepage
The PDE used in PLAXIS LE - Groundwater for the solution of three-dimensional transient saturated/unsaturated seepage
problems is as follows:
  h u    h u    h u 
+ kvd w  +  kwz + kvd w  +  kwy + kvd w 
x 
kwx
x x  z  z z  y  y y 
[10]
h
= − w m2w
t
where:
x = first horizontal direction,
z = second horizontal direction, orthogonal to the x-direction, and
y = vertical direction, corresponding to elevation
The partial differential equation for flow in the three main orthogonal directions is equal to the flow along each direction, (i.e.,
the x, y, and z-directions). For steady-state conditions, the PDE for seepage reduces to the following equation:
  h u    h u    h u 
k wx + kvd w  + k wz + kvd w  + k wy + kvd w  = 0 [11]
x  x x  z  z z  y  y y 
Neglecting vapor flow and assuming the soil is saturated, the PDE governing steady state seepage reduces to:
  h    h    h 
k wx + k wz + k wy =0 [12]
x  x  z  z  y  y 
3.4 UNSATURATED FLOW EQUATION

A transient-state seepage problem in which the permeability and volumetric water content variables change in accordance with
soil suction is common to geotechnical engineering. The fluid motion in unsaturated soils is generally assumed to obey the
partial differential equations presented in the previous chapter. These equations are similar to the classical Richards equation
(Hillel, 1980). The PDEs for water seepage in unsaturated soils can be written in several forms. The three most common
formulations of the unsaturated flow equation are identified as the “h-based” form, the “q-based” form, and the “mixed form”.
There are advantages and disadvantages associated with each of the formulations. The advantages and disadvantages become
apparent when solving particular unsaturated seepage problems.
There are several advantages associated with the q-based form. One advantage is that it can be formulated to be perfectly
water mass conservative. It is not commonly used, however, because this form degenerates in fully saturated systems and
material discontinuities produce discontinuous q profiles.
The h-based form of the PDE seepage equation is the most commonly implemented form. Its primary drawback is that it can
suffer from poor water mass balance when solving transient seepage problems. The poor water mass balance problem is
exacerbated in situations where the soil-water characteristic curve for the material is highly nonlinear.
Celia (1990) proposed a “mixed form” of the Richards equation. The “mixed form” was designed to improve the water mass-
balance of the “h-based” formulation.
3.4.1 H-Based Formulation

The most commonly implemented form of the governing partial differential equation for transient seepage is shown below. The
PLAXIS LE - Groundwater formulation is presented here and is based on total head:
 

x  x  y 
(
) y

(k wx + kvd ) h  +   k wy + kvd h − kvd  = − wm2w h
t
[13]

The partial differential equation governing water flow and the storage of water within a saturated/unsaturated soil is formulated
using total head, h, as the primary variable. However, pore-water pressure, uw, (written as a water head, uw/w), can be used
when solving certain seepage problems.
The formulations presented above satisfy the condition that the difference between the water flow entering or leaving a unit
volume is equal to the change in volumetric water content. Under steady-state conditions, the water flux entering and leaving
a unit volume is the same and as a consequence, the storage term (i.e., right-hand side of the equation) becomes zero.
The above formulation makes the assumptions that there is no loading or unloading of the soil mass during the transient
process. Pore-air pressures are assumed to remain constant and at atmospheric pressure. Changes in volumetric water content
are assumed to be strictly dependent on changes in the soil suction state variable.
BENTLEY SYSTEMS Boundary Conditions In Seepage 14 of 53
4 BOUNDARY CONDITIONS IN SEEPAGE

Several types of boundary conditions can be applied when solving seepage problems. The boundary conditions associated with
seepage are as follows:
• Natural (or Neumann) boundary condition: flux,

• Essential (or Dirichlet) boundary conditions: imposed value, and
• Special boundary conditions: combinations of the above conditions.
The two basic boundary conditions intrinsically related to the formulations normally used in the finite element solutions are:
Essential (or Value) boundary condition and Natural boundary conditions. Essential boundary conditions are assigned to nodes
as fixed values. Natural boundary conditions are assigned to the sides of elements (and are defined by a surface integral).
Natural boundary conditions correspond to the surface integral of the term inside the outer derivative of second order partial
differential equations.
4.1 NATURAL BOUNDARY CONDITION

Natural (or Neumann) boundary conditions arise from the integration by parts of second order derivatives. The integration
results in a surface integral that corresponds to a flux quantity. Therefore, natural boundary condition associated with the
seepage PDE corresponds to the total amount of water flow normal to the surface (m3/s).
Natural boundary condition =  ( v  n ) dS

[14]

where:
 = boundary area,
v = total flux vector,
n = outward surface-normal vector, and
dS = an infinitesimal element at the boundary.
Natural boundary conditions are called “flux boundary conditions” in PLAXIS LE - Groundwater.
Natural boundary conditions are appropriate choices for the representation of situations such as simple soil-atmosphere fluxes,
the water uptake inside a well, and the groundwater flow taking place at the bottom of a domain. The absence of a boundary
condition corresponds to a zero flux natural boundary condition.
Natural boundary conditions can be applied as constant values). An expression can be a function of time, space, or any other
meaningful variable. Flux expressions can be used to represent several hypothetical and real world scenarios such as the
increase in water uptake in a well during the course of a day or the increase in groundwater flow with depth.
The natural boundary conditions associated with the seepage PDEs do not make a distinction between the types of flow (i.e.,
whether it is liquid or vapor flow). The determination of the amount of flow taking place in the form of a liquid or a vapor is not
required for the application of a natural boundary condition.
For instance, the imposition of a negative flux at a soil surface can result in both liquid and vapor fluxes at the surface. The
partitioning of the imposed total flux into vapor and liquid flux will depend on the soil properties and pore-water pressures.
Nevertheless, the total amount of flux at the surface will always correspond to the applied boundary conditions. The fractions
of liquid and vapor flow can be determined from the resulting pore-water pressure gradients, temperature gradients and the
soil property functions.
4.1.1 Differences Between Actual and Applied Boundary Flux

It is possible that a natural boundary condition might apply more fluid at the boundary than the model can reasonably accept.
This is particularly true when running numerical models which are saturated. The mass-balance of the modeling domain must
be of primary consideration. Therefore, if the user applies a large boundary flux which cannot be accepted by the modeling
domain, then the natural boundary condition will need to be reduced to represent the flux which can be applied while maintaining
an appropriate mass balance over the modeling domain. The excess flux which is not applied is conceptually similar to “runoff”
but this type of runoff may be present even if the “runoff” specifier in the software is not specified.
4.2 ESSENTIAL BOUNDARY CONDITION

Essential (or Dirichlet) boundary conditions correspond to predetermined hydraulic head values.
Essential Boundary Condition = h [15]
where:
uw
h = specified hydraulic (total) head value = +y
w
uw = pore-water pressure,
w = unit weight of water,
y = elevation (y in a model).
This type of boundary condition is called “head” or “hydraulic head” in PLAXIS LE - Groundwater. Essential boundary conditions
can be used to represent numerous situations such as the head imposed by a water reservoir (Figure 4) or the head at the
bottom of a domain where the water table is relatively constant. Essential boundary conditions are always required in steady-
state problems. Transient problem may or may not present an essential boundary condition.
Forced negative Phreatic

Original Reservoir
pore-water line
Elevation
pressures
Specified head
Free surface
(Pressure = 0)
Reservoir
Elevation
Figure 4 Head upstream boundary conditions in the case of rapid drawdown simulation
Essential boundary conditions can be applied as constant values.Similar to the natural boundary condition, an expression can
be a function of time, space, or any other meaningful variable. Head expressions can be used to represent several hypothetical
and real world scenarios such as the increase in the filling of a reservoir or transient heads imposed by tides.
Essential boundary conditions must be consistent with the initial conditions of a problem. Discontinuity not only misrepresents
the actual problem, but also results in numerical oscillations. For instance, the simulation of a sudden reduction in pore-water
pressure at a surface should always be done using the ramping of head over time, starting with an initial head equal to initial
conditions. An appropriate ramping time interval that represents the actual rate of change should be selected.
The user must be careful when applying a head boundary condition to the upstream side of an earth levee or dam during a
rapid drawdown scenario. The head boundary condition may be inappropriate for some situations as it inherently forces
unsaturated soil conditions above the water table. Professional judgment is needed to determine whether or not the “forcing”
of negative pore-water pressures is appropriate within the context of the problem under consideration.
4.2.1 Review Boundary Condition (Drain)
More complex boundary conditions are required in order to model certain seepage problems. An example is the situation of a
free drainage surfaces with an unknown seepage exit point. This situation commonly occurs on the downstream slope of an
earth fill dam (Figure 5). The downstream face of the dam can be represented as a modification of the natural and essential
boundary conditions.
Figure 5 “Review boundary” condition for a free flow surface
The boundary condition used to represent the free flow condition is shown in Figure 5. The boundary is called a “Review boundary
condition” in PLAXIS LE - Groundwater. The terminology “drain boundary condition” is sometimes used to refer to this type of
boundary condition.
A review boundary condition is applied as follows:
• If pore-water pressures are negative, then the boundary condition has zero flux,
• If pore-water pressure is positive, then the boundary condition is equal to a negative (outward) flux that brings
pore-water pressures on the surface to zero.
If the amount of negative flux is large, changes at the boundary will take place nearly instantaneously. A large negative flux is
equivalent to setting the essential boundary condition equal to the elevation relative to the datum:
h= y [16]
The essential boundary condition creates another degree of nonlinearity in the system that must be solved through the use of
an iterative process. Another type of special boundary condition available in PLAXIS LE - Groundwater is called the “climate
boundary condition”.
Climate boundary conditions are used to model complex soil-atmosphere fluxes. Additional input data is required in order to
use climatic boundary conditions. The following chapter presents the theory for this type of boundary condition.
4.3 GRADIENT OR UNIT GRADIENT BOUNDARY CONDITION

The gradient or unit gradient boundary condition may be applicable in certain cirmustances, where the user wants to control
the flow out of the model by controlling of hydraulic head gradient. The rate of water flow can be described using Darcy’s law
(see section 2.2.1). The rate of water flow can be written as:
dh
q = −kA [17]
dl
where:
q = rate of water flow (volume per unit time)
k = hydraulic conductivity of the medium,
A = cross-sectional area of the column through which the water flows through,
dh/dl = hydraulic gradient, that is, the change in head over the length of interest.
For a unit gradient boundary condition, dh/dl = 1. This type of boundary condition is typically applied to the bottom of a 1D
numerical model when performing cover design as there is reasonable precedence for the unit gradient boundary condition in
the cover design application. The solver implicitly applies dh/dl = 1 and no value is required to input from the user. For a
gradient boundary condition, a value other than 1 can be entered by the user. Since the gradient boundary conditions assume
the flow is always outward of the model, the entered gradient value should not be negative. The recommended value should be
in the range of 0 < dh/dl <= 1.
The current implementation allows control of the gradient in one direction only, i.e., the flow to move outward in the vertical
direction only. Hence, the recommended use of a gradient or unit gradient boundary condition is at the bottom of the model
only.
4.4 INITIAL WATER TABLE

Initial pore-water pressure conditions can be calculated assuming hydrostatic conditions and the following equation:
uw =  w ( h − y ) [18]
where:
uw = pressure, kPa,
w = unit weight of water, kN/m3,
h = total water head, m, and
y = elevation, m.
Therefore,
when h=y uw = 0; represents the phreatic surface
when h>y uw > 0; represents the saturated region

when h<y uw < 0; represents the unsaturated region
Defining a water table as initial conditions is well-suited for some situations requiring a quick solution. For more complex
situations, using hydrostatic conditions may not be realistic and may cause convergence problems. To overcome convergence
problems with complex models it is best to import the initial head conditions from a steady-state run of the model.
BENTLEY SYSTEMS Flux Sections 18 of 53
5 FLUX SECTIONS
In PLAXIS LE - Groundwater, the flux of water across a user-defined section can be determined in either steady state or transient
analysis. This flux section tool is used to calculate water fluxes into or out of a specific line segment, region segment or other
places of interest. Some common applications are to determine the amount of seepage on the downstream face of a dam or
below cut-off walls.
5.1 FLUX SECTION THEORY CALCULATION

The theoretical calculation of flux section is illustrated in this section. Refer to Figure 6 for a diagram illustrating the theory. A
3-noded linear triangular element is used to demonstrate the theory and the flux of water across section A-A is calculated. The
line section A-A crosses the element at points a and b with Cartesian (global) coordinates, (xa, ya) and (xb, yb), respectively.
Point C (x, y) is at the center of the flux section, and water flux passing this point is defined as (Bathe, 1982):

q x 
q =   = − DB( x, y )H [19]

 
q y 
where:
qx = flux in x-direction
qy = flux in y-direction
D = material property matrix
B = matrix for head gradient interpolation, and
H = vector of of total head at the element nodes
The D matrix for an isotropic material is defined as:
k x 0
D=  [20]
 0 k y 
The B matrix is defined as:
 N1 N 2 N 3 
 x 
B =  x x
N 3  [21]
 N1 N 2
 y y y 
where:
N1, N2, N3 = shape functions at nodes 1, 2 and 3, respectively
The H is the total head vector of the element nodes:
 H1 
 
H = H 2  [22]
 H 3 
where:
H1, H2, H3 = the total head at nodes 1, 2 and 3, respectively
The total head at a point (x, y) within the element can be defined using the shape functions and head values at the nodes as:
3
H ( x, y ) = N H
i =1
i i [23]
The above shape functions are usually expressed in local (natural) coordinates and for triangular elements (shown in Figure 7)
they are defined as:
N1 = 1 −  − 
N2 =  [24]
N3 = 
The B matrix can be rewritten as:
−1
 x y   N1 N 2 N3   N1 N 2 N3 
     
       −1     
B= =J [25]
 x x   N1 N1 N1   N1 N1 N1 
             
  
 x y 
   
The matrix J =  is called Jacobian matrix and J–1 is the inverse of J.
 x x 
   

For iso-parameter elements, the same shape functions are used to determine coordinates of a point (x, y) within the element
in global coordinates.
3
x= N x
i =1
i i
[26]
3
y= N y
i =1
i i
Combining the equations [19], [20], [22] and [25], the flux vector q(qx, qy) at the centre point of the flux section is
determined. The flux crossing the flux section is determined as:
q = q x2 + q 2y
l = ( xa − xb ) 2 + ( ya − yb ) 2 [27]

Qx  
q x 0  xa − xb 

 =  
  
Q y  0 qy 
 ya − yb 

The normal flux across the flux section is determined as:
Qn = ( q  n )  l [28]
where:
n = unit normal vector to the flux section shown in Figure 6
l = length of the flux section
= fluxes in x-, y-directions and normal flux across the flux section,
Qx, Qy, Qn
respectively
Figure 6. Flux section calculation for a 3-noded element
Figure 7. 3-noded element in local (natural) coordinates

BENTLEY SYSTEMS Material Properties 21 of 53
6 MATERIAL PROPERTIES
This section will present the theory behind material properties used in the PLAXIS LE - Groundwater modeling software.
6.1 SOIL-WATER CHARACTERISTIC CURVE

The soil-water characteristic curve is central to the application of unsaturated soil mechanics. It defines the nonlinear
relationship between the amount of water in the soil (i.e., water content) and soil suction. Historically, the amount of water in
the soil has been represented using the volumetric water content variable. However, other designations such as the degree of
saturation of the soil prove to be superior designations of the amount of water in the soil when the soil undergoes volume
change as soil suction is increased. Representation of the soil-water characteristic curve is accomplished either through fitting
existing data or estimating the curve from grain-size information.
6.1.1 Fredlund and Xing (1994) Equation

Fredlund and Xing (1994) presented a three-parameter equation with the flexibility to fit a wide range of materials. The equation
also contained a correction variable that provided increased accuracy in the high suction range. The parameters of the equation
were typically found using a least-squares algorithm. The original form of the equation was as follows.
 
     
 In 1 +    
hr    
 w =  s 1 − 
1
[29]
  106      nf

mf 
 In 1 +     
  hr     In exp (1) +    
    af  
 
where:
 = soil suction value (kPa).
w = volumetric water content at soil suction, ,
s = saturated volumetric water content,
= material parameter which is primarily related to the air-entry value of the
af
soil in kPa,
= material parameter which is primarily a function of the rate of water
nf
extraction from the soil once the air-entry value has been exceeded,
mf = material parameter which is primarily a function of the residual water

content,
hr = suction at residual water content (kPa), and
Fitting method: Least squares nonlinear regression

Required input: Drying laboratory data consisting of points on the curve of volumetric water content versus
soil suction.
Applicable material types: All soils
6.1.2 Fredlund (2000) Bimodal Equation

The bimodal equation can be viewed as two superimposed unimodal SWCC curves. The fitting algorithm fits the bimodal equation
by subdividing the overall curve into an upper and lower portion. Each of the two portions is then fit with a nonlinear least
squares regression algorithm and the results are then combined through the use of superposition. The breaking point between
the two curves is designated using the w parameter which must be selected by the analyst.
    
          
       ln 1 +  
 1 1 
  1 −   3000  
w = ws  s   + (1 − s )  [30]
   n fb m fb
 a fb      k fb l fb
 j        1000000   
     ln 1 + 

 ln  exp(1) +      ln  exp(1) +  fb   3000   
                
     
where:
w = gravimetric water content at any soil suction,
ws = gravimetric water content at any soil suction,
 = soil suction, kPa,
afb = fitting parameter,

nfb = fitting parameter,
mfb = fitting parameter,
jfb = fitting parameter, kPa,
kfb = fitting parameter,
lfb = fitting parameter, and
s = Fredlund bimodal split

Required input: Drying laboratory data consisting of points on the curve of volumetric water content versus
soil suction.
6.1.3 Van Genuchten (1980) Equation

Van Genuchten (1980) presented a three-parameter equation with the flexibility to fit a wide range of materials. The parameters
of the equation could be found using a least-squares algorithm.
 
 
+ ( ws − wrvg ) 
1
ww = wrvg  [31]
 1 + ( avg ) vg 
mvg

n
   
where:
ww = gravimetric water content at any soil suction,
wrvg = residual gravimetric water content,
ws = saturated gravimetric water content,
avg = material parameter which is primarily a function of the air-entry value of

the soil in kPa,
nvg = material parameter which is primarily a function of the rate of water

extraction from the soil once the air-entry value has been exceeded, and
mvg = fitting parameter, and
 = soil suction (kPa).

Required input: Drying laboratory data consisting of points on the curve of gravimetric water content versus
soil suction. It should be noted that data points well beyond residual suction conditions may
distort the best-fit analysis.
6.1.4 Van Genuchten (1980) and Mualem (1976) Equation

Van Genuchten (1980) proposed a closed form simplification for solving the Mualem (1976) integral equation. The relationship
between the m and n parameters of the van Genuchten (1980) equation was first prescribed in order to reduce the number of
fitting parameters from three to two. The simplification proposed by Mualem (1976) is shown below.
 
 1 
ww = wrm + ( ws − wrm )   1 
 [32]
 nm 1− n  
 1 + ( am )   m  
where:
wrm = residual gravimetric water content,

am = material parameter which is primarily a function of the air-entry value of

the soil in kPa,
nm = material parameter which is primarily a function of the rate of water

extraction from the soil once the air-entry value has been exceeded, and

soil suction.
6.1.5 Gardner (1958) Equation

Gardner (1958) presented a continuous equation for the first coefficient of permeability function. The form of the equation has
subsequently been the basis for the soil-water characteristic curve as well as many other equations proposed in subsequent
literature. However, it should be noted that the equation was originally proposed as an equation to best-fit measured
permeability data.
 
ww = wrg + ( ws − wrg )  1  [33]
 1 + ag ng  
 
where:
wrg = residual gravimetric water content,
ag = material parameter which is primarily a function of the air-entry value of

the soil in kPa,
ng = material parameter which is primarily a function of the rate of water

extraction from the soil once the air entry value has been exceeded, and

soil suction.
6.1.6 Brooks and Corey (1964) Equation

Brooks and Corey (1964) proposed a power-law relationship for the SWCC. The model represented an attempt to use an
equation to describe the soil-water characteristic curve. The equation can be written as follows:
nc
a 
ww = wr + ( ws − wr )  c  [34]
 
where:
wr = residual gravimetric water content,
ac = bubbling pressure (kPa),
nc = pore size distribution index (dimensionless), and


soil suction.
6.1.7 Zapata Estimation (2000)

Zapata (2000) performed a statistical regression analysis on a sampling of soils from 50 states in the USA. The regression led
to “average” recommended curves based on categorization groupings as fine or coarse-grained soils. The average values for
the SWCCs were further grouped according to plasticity index (PI). This estimation method provides a simple approximate
method for the user to enter simple information about the SWCC for a soil and can be best-fit with the Fredlund and Xing (1994)
fitting parameters.
6.2 HYDRAULIC CONDUCTIVITY FIT METHODS

The PLAXIS LE - Groundwater software implements Gardner’s equation for fitting unsaturated permeability data.
6.2.1 Gardner’s (1958) Equation

Gardner (1958) permeability function for unsaturated soils is expressed as a function of suction:
ks
kw =
n
   [35]
1+ a 
  w g 
where:
kw = hydraulic conductivity or permeability of the water phase,
ks = saturated hydraulic conductivity of the water,
w = density of water,
a = fitting parameter,
n = fitting parameter,
g = acceleration of gravity, and

Required input: Laboratory data consisting of at least three points.
Applicable material types: All soils.
The Gardner (1958) equation provides a flexible permeability function that is defined using two parameters, a and n. The
parameter, n defines the slope of the function, and a is a parameter related to the breaking point of the function. The Gardner
(1958) equation is meant to be obtained from laboratory data.
Figure 8 shows the sensitivity of the parameters, a and n. The permeability function has been quite often used in saturated-
unsaturated flow modeling. The Gardner (1958) equation is sensitive to the air-entry value of the soil and the rate of
desaturation. These features are modeled in a continuous manner.
Figure 8 Sensitivity of Gardner’s (1958) equation for the coefficient of permeability as a function of the matric suction (from
Fredlund and Rahardjo, 1993)
A set of data is presented in Figure 9 to demonstrate the ability of the Gardner (1958) equation to fit laboratory data of the
coefficient of permeability for various soils. The equation is easy to apply for saturated-unsaturated flow modeling.
Figure 9 Comparison between the measured and the predicted coefficient of permeability values for different materials using the
Gardner (1958) equation (from Huang et al., 1994)
6.3 HYDRAULIC CONDUCTIVITY ESTIMATION

Estimation methods (also known as pedo-transfer functions) are provided in the PLAXIS LE - Groundwater software in order to
facilitate the estimation of the unsaturated portion of the hydraulic conductivity curve. Most estimation methods are based on
a description of the soil-water characteristic curve and therefore require a specific fit to be present in the software. Commonly
used estimation methods are provided in the following sections.
6.3.1 Brooks and Corey Estimation

Estimation methods (also known as pedo-transfer functions) are provided in the PLAXIS LE - Groundwater software in order to
facilitate the estimation of the unsaturated portion of the hydraulic conductivity curve. Most estimation methods are based on
a description of the soil-water characteristic curve and the saturated coefficient of permeability. Specific fits for the hydraulic
conductivity are present in the PLAXIS LE - Groundwater software. Commonly used estimation methods are provided in the
following sections.
2 + 3
 
k = ksat  b  for suction,   b
  [36]
k = ksat for suction,   b
where:
k = hydraulic conductivity (or coefficient of permeability) with respect to the
water phase, m/s
ksat = saturated hydraulic conductivity with respect to the water, m/s
b = Brooks and Corey (1964) soil-water characteristic curve fitting parameter,
 = Brooks and Corey (1964) soil-water characteristic curve fitting parameter,
 = soil suction, kPa.
Required input: Saturated hydraulic conductivity and a fit of the soil-water characteristic curve using the
Brooks and Corey (1964) equation.
The Brooks and Corey (1964) equation that fits some soil-water characteristic curve data can be written in the form of a power-
law relationship.

 
 =  b  for suction,   b [37]
 
where:
 = normalized water content (defined in Equation [38]),
b = air-entry value,
 = any suction, and
 = pore-size distribution index.
The normalized volumetric water content, , is defined as follows:

 − r
= [38]
s − r
where:
s = saturated volumetric water content, and
r = residual volumetric water content.
Equation [40] is suitable for fitting laboratory SWCC data for coarse materials that have a low air-entry value.
Brooks and Corey (1964) also suggested a procedure for estimating the residual water content. The Brooks and Corey (1964)
permeability function is based on the model of a porous media developed by Burdine (1953), Kozeny (1927), and Wyllie and
Gardner (1958). The recommended function is shown below:
k w = k s for    b [39]
k w = k s   for    b [40]
where:
kw = coefficient of permeability with respect to the water phase for the soil
saturation (i.e., S = 100%),
 = empirical constant.
The empirical constant,  in turn is related to the pore-size distribution index.
The Brooks and Corey (1964) model is simple to use and appears to be quite reasonable for coarse-grained soils such as sands
and gravels.
6.3.2 Modified Campbell (1996) Estimation

The Modified Campbell (1996) equation is implemented into PLAXIS LE - Groundwater to provide a hydraulic conductivity
equation that levels off at high soil suctions. The shape of the function is consistent with theoretical considerations for hydraulic
conductivity of an unsaturated soil when water flow transitions to the vapor phase. The point of residual suction can be assumed
to be the point at which water movement becomes discontinuous or transitions to vapor flow. The Campbell (1973) equation
was modified to produce an equation that tends to level off at approximately the residual suction for the soil. The modified
equation as implemented into the PLAXIS LE - Groundwater software is presented below (Fredlund, 1996).
p
  
  
     
 ln  1 +    
 hr    
  1
k ( ) = ( k s − kmin )  1 −     + kmin [41]
  10    
6

mf

  ln 1 + 
nf
  
  hr    ln exp (1) +   

   af  
         

where:
k = hydraulic conductivity or permeability of the water phase, m/s
ks = saturated hydraulic conductivity of the water phase, m/s
kmin = calculated minimum hydraulic conductivity, m/s
p = parameter used to control the modified Campbell (1973) estimation of
hydraulic conductivity,
af = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
nf = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
mf = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
hr = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
and
 = soil suction.
Required input: Saturated hydraulic conductivity and a fit of the soil-water characteristic curve by the
Fredlund and Xing (1994) equation.
As a material dries, there is less and less water present in the soil matrix. The hydraulic conductivity then decreases accordingly
as the volumetric water content decreases. The modified Campbell equation reflects this behavior by using the following
equation.
k = ks  p ( ) [42]
where:
k = hydraulic conductivity at any level of suction, m/s
ks = saturated hydraulic conductivity, m/s,
 = soil suction, kPa,
 = normalized volumetric water content or θw/θs represented with any equation
(i.e., van Genuchten, 1980; Fredlund and Xing, 1994), and
p = power factor to adjust the prediction (same as in Equation [41]).
A modification was made to Campbell’s (1973) equation before it was implemented into the PLAXIS LE - Groundwater software.
The modification adjusts the Campbell equation such that the function flattens once a minimum permeability has been reached.
The hydraulic conductivity remains relatively constant once the water phase in the soil becomes discontinuous. Water flow in
the soil is then primarily the result of vapor diffusion through air. The vapor phase flow can be accommodated through use of
the Campbell (1973) equation as shown below:
k ( ) = ( ks − kmin )  p ( ) + kmin [43]

where:
kmin = minimum permeability.
The above equation allows the hydraulic conductivity versus soil suction function to level off after a particular soil suction has
been reached. Initially, it was suggested that the equation could be set to level off once the residual water content conditions
had been reached. However, it was observed for some laboratory data that the permeability function tended to flatten at about
one log cycle of suction higher than the suction corresponding to the residual water content.
The method proposed by Campbell (1973) is implemented in PLAXIS LE - Groundwater. The implemented algorithm uses the
soil-water characteristic curve and the saturated hydraulic conductivity to estimate the hydraulic conductivity of a soil at all
levels of suction.
6.3.3 Fredlund, Xing and Huang (1994) Estimation

Fredlund et al., (1994) presented a modification of the Mualem (1976) integration method for estimating the hydraulic
conductivity of a material as a function of soil suction. The integration procedure is complex and a closed-form solution is not
available. PLAXIS LE - Groundwater performs the integration and computes a series of data points that can be written onto a
graph as x-y data.
Solution method: Integration by Simpsons rule

Required input: Saturated hydraulic conductivity and Fredlund and Xing (1994) fit of soil-water characteristic
curve
Equations available in the literature for predicting the coefficient of permeability use the soil-water characteristic curve data
only for a limited range of suction values (e.g., Brooks and Corey, 1964; Mualem, 1976; van Genuchten, 1980). These equations
require knowledge of the residual water content. The residual water content θr is the water content below which a large increase
in suction is required to remove additional water. Kunze et al., (1968) investigated the effect of using a partial soil-water
characteristic curve for the prediction of coefficient of permeability and concluded that the accuracy of prediction significantly
improved when the complete soil-water characteristic curve was used.
Fredlund et al., (1994) proposed an equation to estimate the coefficient of permeability of a soil over an extended range of soil
suction values. The estimation procedure makes use of the soil-water characteristic curve data for the entire suction range of
0 to 1,000,000 kPa. This equation tends to be more practical for the estimation of the coefficient of permeability over a large
range of suction values. The coefficient of permeability function is of interest at large suction values particularly for structures
such as soil covers, as well as other near-ground-surface structures.
The equation suggested by Fredlund et al., (1994) for predicting the coefficient of permeability is given below :
b
 ( e y ) −  ( )
   ( e y ) dy
ln ( ) ey
kr ( ) = [44]
b
 ( e y ) −  ( aev )
   ( e ) dy
y
ln ( aev )
ey
where:
b = ln(1,000,000),
y = dummy variable of integration representing the logarithm of suction, and
 = soil suction, given a function of volumetric water content, and
aev = air entry value of the soil under consideration.
The Fredlund et al., (1994) permeability equation makes use of the Fredlund and Xing (1994) equation (i.e., equation for fitting
the soil-water characteristic curve data for the entire range of suctions). The Fredlund and Xing (1994) equation has been found
to fit the soil-water characteristic data for essentially all type of soils and over all suction ranges (Benson et al., 1997; Leong
and Rahardjo, 1997). More details are available in Fredlund et al., (1994).
6.3.4 Van Genuchten (1980) and Mualem (1976) Estimation

The equation for calculating the permeability function by the van Genuchten (1980) and Mualem (1976) estimation method is
based on Equation [32]. The van Genuchten and Mualem SWCC fitting equation is as follows.
 2

n −1− n  
 1
 
 1
1 − ( a )  n  1 + ( a ) 
   
n1− 
  

k ( ) = ks    1
[45]
  1−  2
1 + ( a ) n   n  
 
 
where:
k = hydraulic conductivity or permeability of the water phase, m/s,
ks = saturated hydraulic conductivity of the water phase, m/s,
a = van Genuchten and Mualem soil-water characteristic curve fitting

parameter,
n = van Genuchten and Mualem soil-water characteristic curve fitting
parameter, and
Required input: Saturated hydraulic conductivity and van Genuchten and Mualem fit of the soil-water
characteristic curve
6.3.5 Van Genuchten (1980) Estimation

Several investigators such as Brooks and Corey (1964) and Mualem (1976) have proposed closed-form equations for estimating
the coefficient of permeability of unsaturated soils based on Burdine’s theory (1953). Brooks and Corey (1964) equation may
not converge rapidly when used in numerical simulations of seepage in saturated-unsaturated soils. The Mualem (1976)
equation is in an integral form. However, it is possible to derive a closed-form analytical equation provided there is a fixed
relationship between two of the fitting parameters.
The equation proposed for fitting the soil-water characteristic curve by van Genuchten (1980) is flexible, continuous and has a
continuous slope. The closed-form equation proposed for estimating the coefficient of permeability function has been extensively
used for saturated-unsaturated soils flow modeling.
 −m

2

 1 − ( a ) 1 + ( a )  
nm n

  
k ( ) = ks   
 [46]
1 + ( a ) n 
m 2
   
 
where:
ks = saturated hydraulic conductivity of the water phase, m/s,
a = van Genuchten soil-water characteristic curve fitting parameter, kPa,

n = van Genuchten soil-water characteristic curve fitting parameter,
m = van Genuchten soil-water characteristic curve fitting parameter, and
Required input: Saturated hydraulic conductivity and van Genuchten fit of the soil-water characteristic curve
The van-Genuchten’s equation (1980) for fitting the soil-water characteristic curve data is given below:
 = r +
( s − r )
[47]
( )
m
1 + a n 
 
where:
 = volumetric water content,
s = saturated volumetric water content,
r = residual volumetric water content, and
 a, n, m = material constants.
van Genuchten (1980) suggests the use of 1,500 kPa to represent residual conditions for a soil. For many soils, a volumetric
water content corresponding to a residual suction of 1,500 kPa is a reasonable approximation. An analytical procedure has also
been suggested for estimating the residual water content.
Figure 10 provides the comparison between the predicted and measured values of the soil-water characteristic curve along the
drying and wetting paths with respect to suction for Guelph loam (van Genuchten, 1980). Also shown is the variation in the
coefficient of permeability. The equations proposed by van Genuchten (1980) provide excellent fits for many soil types.
Figure 10 Comparison between the predicted (continuous solid lines) and measured values (circles) of the soil-water characteristic
curve along drying and wetting paths and the variation of coefficient of permeability with respect to suction (from van Genuchten,
1980)
6.3.6 Leong and Rahardjo (1997) Estimation

Leong and Rahardjo (1997) proposed a permeability function for predicting the unsaturated coefficient of permeability. The
estimation is based on a fit of the soil-water characteristic curve with the Fredlund and Xing (1994) equation. The equation for
estimation of the unsaturated soil permeability function can be written as follows:
p
 
 
 
 1 
k ( ) = ks   [48]
  
mf
nf

 ln exp(1) +   
  
   af   
       
where:
ks = saturated hydraulic conductivity of the water, m/s,
p = parameter used to control the Leong and Rahardjo (1997) estimation of
hydraulic conductivity,
af = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
kPa,
nf = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
mf = Fredlund and Xing (1994) soil-water characteristic curve fitting parameter,
In Leong and Rahardjo (1997), the best fitted permeability function was used when comparing the predicted and measured
coefficient of permeability for several soil types. The results were tested for both the wetting and drying curves. A good fit
was obtained for a wide range of experimental data. It was found that if the exponent p was known for a given soil, the
coefficient of permeability could be obtained indirectly from the soil-water characteristic curve. Otherwise, p needed to be
determined using a curve fitting process with permeability data. The value of p varied from 4.3 to 52.1 for the soils studied.
BENTLEY SYSTEMS Soil Atmosphere Modeling 32 of 53
7 SOIL ATMOSPHERE MODELING

Soil-atmosphere moisture fluxes can be modeled by using “Climate boundary conditions”. Soil-atmosphere moisture fluxes are
influenced by the following types of processes:
• Infiltration fluxes,
• Runoff, and
• Evaporative fluxes.
Moisture infiltration and runoff can be modeled by considering the amount of precipitation and the infiltration capacity of the
soil. The infiltration capacity is determined by applying a modified version of the “Review boundary condition”.
Evaporative fluxes are modeled by using an appropriate negative flow boundary condition and by considering moisture
movement through vapor flow. Consideration of vapor flow becomes important when modeling soil-atmosphere fluxes. Liquid
flow alone cannot represent the entire moisture migration associated with evaporation at the soil surface. Eventually there will
essentially be a shut-off of both liquid and vapor moisture movement near the ground surface. The vapor flux component theory
has been presented in the previous chapters.
The following section presents a detailed description of the theory of soil-atmosphere modeling.
7.1 ATMOSPHERIC FLUX BALANCE

There is an atmospheric moisture flux balance and a thermal flux balance that must be satisfied at the ground surface when
calculating Actual Evaporation, AE. Basically, water falling on the ground surface either infiltrates the soil (or runs off) or else
rises to the sky through the process called “Actual Evaporation, AE”. The ground surface moisture and thermal flux equation
can be written as follows.
P = AE + NP + Roff [49]
Qn = Qh + Ql + Qg [50]
where:
P = precipitation, m3/m2/day, or m/day
AE = actual evaporation from ground surface, m3/m2/day or m/day
NP = net Percolation or infiltration, m3/m2/day or m/day
Roff = runoff, m3/m2/day or m/day
Qn = net radiation, kJ/m2/day, or equivalently converted into m/day,
= sensible heat transferring from ground surface to air, kJ/m2/day, or
Qh
equivalently converted into m/day,
= latent heat associated with the water phase change including evaporation
Ql
or freezing, kJ/m2/day, or equivalently converted into m/day, and
Qg = ground heat flux, kJ/m2/day, or equivalently converted into m/day.
Precipitation information can be obtained from weather station records and is usually provided on a daily basis. Preferably
precipitation data should be collected on an hourly basis when modeling near-ground-surface phenomena. The mechanics of
net infiltration, NP, can be described by Darcy law. Net radiation, Qn, can also be obtained from weather station records or it
can be approximated using an equation suggested by Penman in (1948). The latent heat component, Ql, can be estimated using
actual evaporation, AE, or the formation of ice near the ground surface during freezing. The sensible heat component, Qh,
reflected from the ground surface to the air is described as follows (Penman, 1948; Gray, 1970; Wilson 1990):
Qh = C f  f (u )(Ts − Ta ) [51]
where:
Qh = sensible heat, m/day,
Cf = conversion factor, (i.e., 1 kPa = 0.0075 mHg),
 = psychometric constant, 0.06733 kPa/ oC at 20 oC,
f(u) = function depending on wind speed, f(u) = 0.35(1+ 0.146 Ww), and
Ww = wind speed, km/hr.
Actual Evaporation, AE, is difficult to measure directly but can be calculated from fundamental thermodynamic considerations.
Equations [49] and [50] are fundamental to describing the coupling of moisture and heat flow processes. Actual evaporation,
AE, depends on the water content and temperature of the soil at ground surface. In addition, the rate of evaporation also
depends on the air temperature and air relative humidity. The air temperature and soil temperature at the ground surface are
generally not the same but are inter-related through net radiation, Qn, latent heat, Ql, and sensible heat, Qh. The available
surface water is controlled by total precipitation, actual evaporation, and runoff. These variables play an important role in
partitioning convective heat flux into sensible heat and latent heat (Wetzel and Boone, 1995).
7.2 INFILTRATION AND RUNOFF WITHOUT CONSIDERATION

OF EVAPORATION
Infiltration, evaporation and runoff are determined simultaneously. If the moisture evaporation can be neglected, the moisture
flux balance at the soil surface is rewritten:
P = NP + Roff [52]
The amount of precipitation and the hydraulic properties of the soil are the main variables required when calculating net
infiltration. Runoff is computed in an iterative manner and the amount of infiltration can be determined by consideration of
previous conditions.
Normal flux boundary conditions can be set to include the effects of runoff.
Depending on the amount of moisture flux applied to the uppermost boundary, the pore-water pressure, uw, may tend to
become higher than zero. The higher the influx, the higher the pore-water pressure, uw becomes. If uw becomes higher than
zero, a condition occurs that corresponds to ponding.
However, if the ground surface is assumed to be well-drained with no ponding, a mechanism must be implemented in order to
limit the amount of infiltration to a lowest possible value. The low value would keep the uw value at the surface equal to zero.
One common way of applying this limiting condition is to switch the boundary condition to an essential boundary condition equal
to zero. This would occur when the pore-water pressure, uw, becomes equal to zero. Another option, switches the boundary
condition to a different natural boundary condition that is equivalent to the essential boundary condition (Gitirana Jr., 2004).
The natural boundary condition is similar to the “Review Boundary Condition” presented in the previous chapter.
The ponding height for a model is by default set to zero. This means that when water is applied to the boundary of a model,
the maximum pore-water pressure will be restricted to a maximum of 0 kPa. Any additional water applied which causes the
upper boundary to exceed 0 kPa and will be re-classified as runoff. If a value greater than zero is specified for ponding then the
maximum pressure allowed at the ground surface is increased to the pond height times w. Runoff conditions will not occur until
the increased maximum pressure is encountered.
The calculation of runoff adds significant complexity to the calculations in a numerical model. It is possible that adding a runoff
calculation might double or triple computational times. A simple model should be set up and solved prior to the addition of a
model implementing a runoff calculation.
7.2.1 Pressure Head Calculated Runoff Correction

This runoff correction option is similar to the equation presented by Gitirana (Gitirana Jr., 2004), but uses the pressure head
variable instead of the hydraulic conductivities and the FAC. This equation, developed by SoilVision Systems Ltd., is able to
achieve the same goal as the Estimated Runoff Correction option, while offering increased stability. It also removes the need
for a user-defined FACTOR input.
if uws  0 then ( pond − hps )


Natural BC =  [53]
 if uws  0 then n p cos 

where:
uws = pore-water pressure at the surface, kPa,
hps = pressure head at the surface, m,
pond = the height of the pond, being 0 if no pond is allowed,
np = is the net precipitation available at the ground surface, and
 = slope of the ground surface with the horizontal direction.
This is the runoff correction option implemented in PLAXIS LE - Groundwater.
7.3 EVAPORATION
The effects of evaporation on a soil near the ground surface depend on the vapor pressure gradient between the soil surface
and the atmosphere. Atmospheric coupling is achieved in PLAXIS LE - Groundwater in the form of an evaporative flux boundary
condition. PLAXIS LE - Groundwater provides a number of methods for defining the evaporative flux in terms of potential
evaporation (PE) or actual evaporation (AE).
When an evaporative boundary condition is being considered in PLAXIS LE - Groundwater, the governing equations must include
vapor pressure gradients.
7.3.1 Potential Evaporation

Potential evaporation, PE, is the amount of evaporation that would occur from a saturated soil surface. The potential evaporation
can be entered: i.) as measured data (i.e., Pan evaporation), ii.) as a constantor iii.) it can be calculated using different
formulations (Gitirana, 2005; Fredlund et al., 2012), such as:
a) Penman (1948) equation,

b) Thornwaite (1948) equation, or
c) Priestley-Taylor (1972) equation.
7.3.1.1 Penman (1948) Equation
The original Penman method is used for the calculation of potential evaporation (PE). The Penman equation uses the routine
weather data as input, namely, relative humidity, air temperature, wind speed, and net radiation.
The potential evaporation at a material-atmosphere boundary can be calculated using the following formulation (Penman,
1948):
Qn +  Ea
PE = [54]
 +
where:
PE = potential evaporation, m/day,
Ea = flux associated with “mixing”; f(u) Cf uv0air (1 – hr), m/day,
f(u) = 0.35 (1. + 0.146 Ww),
Ww wind speed, km/hr,
Cf conversion factor, (i.e., 1 kPa = 0.00750 mHg),
hr = relative humidity in the air above the ground (i.e., hr = uvair/uvoair), obtained
from weather station record,
uvair = water vapor pressure in the air above ground surface, kPa,
uvo air
= saturated vapor pressure at the mean air temperature, kPa,
Γ = slope of saturation vapor pressure versus temperature curve, kPa/oC,
Qn = net radiation at the water surface, m/day,
η = psychrometric constant, (kPa/oC), η = 0.06733 kPa/oC.
The uvair , uvoair and Γ can be calculated from temperature as proposed by Lowe (1977).
The net radiation can be calculated (as suggested by Penman) in the following manner:

( )
0.5 
Qn = (1 − r ) Rc −  (273.15 + Ta )4  0.56 − 0.92 pvsat air  ( 0.10 + 0.90n / N ) [55]

 
where:
Qn = net radiation, m/day,
r = reflection coefficient,
Rc = 0.95 Ra(0.18 + 0.55n/N) = shortwave radiation , m/day,
= coefficient suggested by Penman for evaporation from a wet and bare soil
0.95
as compared to evaporation from an open water surface,
Ra = solar radiation, MJ/m2/day,
 = Stefan Bolzman’s constant, W/m2/K4,
Ta = air temperature, oC,
Pvsatair = vapor pressure of the air above the surface, mmHg, and
= sunshine ratio (actual/possible hours of bright sunshine).
n/N
Note: The vapor pressure, Pvsatair in Equation [55] must be in mmHg.
In addition to the above equation suggested by Penman (1948), the net radiation data can be measured in the field, and is a
typical weather station reading. The net radiation is usually measured in units of MJ/day-m2. Entry of the measured net radiation
is also in units of MJ/day-m2. The software will automatically convert the net radiation from units of MJ/day-m2 into the units
of m/day used in the evaporation calculations in the following expression:
Qn = 106 Qmn ( Lv ) [56]

where:
Qmn = measured net radiation , MJ/day-m2,
Lv = latent heat of evaporation, Lv = 4.187×106 (591 – 0.51Ta), J/m3.
7.3.1.2 Thornthwaite (1948) Equation
Thornthwaite incorporated the variables of length of daylight hours, mean monthly temperature, and an empirical constant into
the calculation of potential evaporation (Fredlund et al., 2012). The daily potential evaporation can be calculated as follows:
at
 L   N   10T 
PE = 0.0005333  d     a  [57]
 12   30   I 
where:
Ld = length of daylight, hr,
N = number of days in the month,
Ta = mean monthly air tempterature, oC,
I = summation for 12 months of the function (Ta/5)1.514 ; i.e.,
1.514
I=   a  ,
12
T
month =1 
 5 
at = complex function of the variable I (i.e., at = (6.75  10-7) I3 – (7.71  10-5)
I2 + (1.79  10-2) I + 0.492.
The Thornthwaite (1948) equation is commonly used to assess the climatic conditions of aridity and humidity based on the
calculation of a moisture index Im. The moisture index Im can be calculated following the modification of the moisture index
equation by Thornthwaite and Mather (1955), and Thornthwaite and Hare (1955):
 P 
I m = 100  − 1 [58]
 PE 
where:
Im = 1955 Thornthwaite moisture index,
P = Total annual precipitation,
PE = total annual potential evapotranspiration calculated as the summation of
the Thornthwaite (1948) monthly potential evaporations.
The climate classification criteria are as shown in Table 2. The moisture index Im is as defined in Equation [58].
Table 2. Criteria for climate classification

Im Category of Climate
> 100 Perhumid

20 – 100 Humid
0 – 20 Moist humid
-33 – 0 Dry subhumid
-67 to -33 Semiarid
-100 to -67 Arid
7.3.2 Actual Evaporation

Actual Evaporation, AE, is difficult to measure directly but it can be calculated from fundamental thermodynamic considerations.
There are several procedures that have emerged for the calculation of the actual evaporative flux. However, there have been
few studies that assess the relationship between each of the suggested solution procedures. As well, the reliability of the
proposed procedures seems to not have been completely verified.
In 1994 Wilson proposed a modification to the well-known Penman (1948) equation for the calculation of Potential Evaporation,
PE. The modified equation has become known as the Wilson-Penman (1994) equation to calculate the actual evaporation. The
Wilson-Penman equation took into consideration the difference in temperature and relative humidity (and therefore vapor
pressure) between the soil surface and the overlying air. The difference in conditions between the air and the water at ground
surface has formed the basis for the Soil-Atmospheric Model which was subsequently implemented into the SoilCover, version
1.0, computer code (1994).
In 1997 a “Limiting Function” type relationship was proposed by Wilson, Fredlund and Barbour. The “Limiting Function” related
Actual Evaporation and Potential Evaporation by scaling the vapor pressures associated with the relative humidity at ground
surface and the relative humidity in the air above ground surface. Inherent in the derivation was the assumption that the air
and soil temperatures were the same.
Wilson (1994) also presented experimental results that showed a unique relationship between total suction at any soil surface
and the ratio of Actual Evaporation to Potential Evaporation, AE/PE. Wilson et al., (1997) presented a unique equation that
passed through the experimental data. As a consequence, there was now another way to empirically relate Actual Evaporation
and Potential Evaporation fluxes.
The above-mentioned relationships give rise to different possibilities for the calculation of Actual Evaporation from the ground
surface. The major difference in the methodologies is related to the assumption regarding the air and soil temperatures. For
example, the soil temperature can be assumed to be equal to the air temperature. This is known to not be the case but there
does not appear to have been a thorough study performed that quantifies the magnitude of the error in calculating AE based
on this assumption.
The Case 1 Solution presented below for Actual Evaporation considers the isothermal case (i.e., no ground thermal flux, Qg =
0). This assumption greatly simplifies the solution for AE since all temperature values are taken as equal to the air temperature
recorded above ground surface. The “Limiting Function” proposed by Wilson et al., (1997) was used to relate Actual Evaporation
and Potential Evaporation. The procedure for solving this case is referred to as the Case 1 Solution.
The second solution procedure considered for calculating Actual Evaporative, AE, utilizes the Wilson-Penman (1994) equation
for the case where ground thermal flux, Qg, is equal to zero beneath the soil surface. However, the soil temperature at ground
surface can be different from the air temperature above ground surface. The solution procedure is referred to as the Case 2
Solution.
The third solution procedure is quite similar to the Case 1 Solution, except that the Actual Evaporative, AE, is approximated
using an empirical expression that is best-fit with Wilson’s 1994 experimental results. This solution procedure is referred to as
the Case 3 Solution.
7.3.3 Case 1 Solution: “Limiting-Function” for AE/PE Model

The following assumptions are made in the isothermal, “Limiting-Function” model.
• There can be liquid and vapor flow through the soil, (i.e., liquid and vapor flow is in response to a hydraulic head
gradient and a vapor pressure gradient, respectively).
• The soil temperature in the entire domain is assumed to be the same, and equal to the air temperature above the
soil surface. (i.e., ground thermal flux is neglected, Qg = 0).
• The ground surface temperature is assumed to be equal to the air temperature.
• Actual evaporation is calculated using the “Limiting Function” proposed by Wilson et al., (1997).
7.3.3.1 Mass Flow Governing Equation

Moisture flow in one-dimension beneath the soil surface is described using the following partial differential equation,
  h  h
y 
(
 kwy + kvh
y
)
− kvh  = − wm2w
 t
[59]
a gDvv  sv 0 hr
kvh = [60]
 w R(273.15 + T )
1.75
 T 
Dv = 2.29  10−5 1 +  [61]
 273.15 
 = ( a )
2/3
= ( n −  L ) 2/3 [62]
− gv 
 w R (273.15 + T ) [63]
hr = e
 = (ua − uw ) +  [64]
where:
h = water head, m,
w = unit weight of water, kN/m3,

kw = hydraulic conductivity, m/s,
kvh = water vapor conductivity by diffusion within the air phase, m/s,
 = soil tortuosity of a dimensionless factor,
a = volumetric air content, a = n - , m3/ m3,
 = volumetric water content, m3/ m3,
n = porosity,
Dv = molecular diffusivity of vapor through soil, m2/s,
ωv = molecular weight of water vapor, 0.018016 kg/mol,
sv0 = saturation vapor density that is dependent of temperature, kg/m3,
hr = relative humidity,
w = water density, kg/m3,
R = universal gas constant, 8.314 J/(mol-K),
T = temperature, oC, and
 = total suction, kPa,
ua = pore-air pressure, kPa,
 = osmotic suction, kPa, and
y = elevation, m.
The osmotic suction in a soil is related to the salt content in the soil. For typical field water content conditions, osmotic suction
may range from 100 kPa to 1000 kPa or more. As the soil dries, the salt contents increase, and the osmotic component increases
(Fredlund, 1991).
It should be noted that PLAXIS LE - Groundwater supports one-dimensional, two-dimensitonal, and three-dimensional
formulations for moisture flow; however, only the one-dimensional partial differential equation is shown (i.e., Equation [59]).
7.3.3.2 Initial soil water contents

The soil moisture can be initialized using an initial set of values for pore-water pressure, water head. It is also possible to
designate the location of the water table (i.e., phreatic line) and assume that hydrostatic conditions exist above and below the
phreatic line.
7.3.3.3 Boundary condition for moisture flow

Using Equation [49] infiltration to the ground surface (i.e., the boundary condition for moisture flow) can be defined as.
qy = P − AE − Roff [65]
surface
where:
qy = moisture flow rate at soil surface, m/day,
P = precipitation flux, m/day,
Roff = water runoff, m/day,
AE = actual evaporation, m/day
7.3.3.4 Actual Evaporation

The “limiting function” proposed by Wilson et al., (1997) is written as follows.
 u soil − u air 
AE = PE  v v 
[66]
 uvo
soil
− uvair 
where:
AE = actual evaporation, m/day,
uv soil
= actual vapor pressure at the soil surface, kPa,
uvosoil = saturated vapor pressure in the soil at the ground surface, kPa, and
uv air
= vapor pressure in the air above the soil surface, kPa.
For the Case 1 Solution, it is assumed that Ts = Ta, and this leads to the vapor pressure in the air being equal to the saturated
vapor pressure in the soil (i.e., uvair
0 = uv 0 ). Equation [66] can also be written in term of relative humidity as follows (Wilson
soil
et al., 1997):
 h − hr 
AE = PE  s  [67]
 1 − hr 
where:
hr = relative humidity of the air above the ground surface, and
hs = relative humidity at the soil surface.
Equations [59], [65], and [66] or [67] together with the initial water content conditions are the equations required to solve
for Actual Evaporation, AE.
7.3.3.5 Determining vapor pressure values throughout each day

The daily air temperature, Ta, and daily relative humidity, hr, of the air above the ground surface are usually measured at a
weather station. Minimum and maximum daily temperature and relative humidity values are generally recorded for each day.
However, usually the hourly values of these variables are also recorded. The amount of data becomes quite excessive when
hourly values are obtained. When minimum and maximum values are given for temperature and relative humidity, then an
assumption can be made regarding their distributions throughout each day.
Figure 11 and Figure 12 shows two options of several possible methods implemented in PLAXIS LE - Groundwater to describe
the daily changing pattern of air temperature and relative humidity based on the daily minimum and maximum value. Figure
11 indicates that the relative humidity has a maximum value at about 6:00 a.m. and a minimum value at about 1:00 p.m. The
time at minimum and maximum value can be specified with other value. In Figure 12, the maximum value of relative humidity
is assumed at the midnight (24:00 a.m.), and minimum value is at noon (12:00 p.m.). The daily changing pattern of air
temperature is in general, opposite to the pattern for relative humidity.
80 40
Relative humidity
Air Temperature
75
35
Relative Humidity (%)

70
30
65
25
60
20
55
6:00 AM 13:00 PM
50 15
0 0.5 1 1.5 2 2.5 3
Time (Day)
Figure 11 Daily changing pattern of air temperature and relative humidity of overlying air
85 32
Relative Humidity
Air Temperature
Relative Humidity and Air Temperature (C)
30
80
28
75
26
70
24
65
22
12:00 PM
60 20
0 0.5 1 1.5 2 2.5 3
Tim e (day)
Figure 12 Symmetric distributions of daily changing of air temperature and relative humidity
The calculated relative humidity at the unsaturated soil surface based on Edlefsen and Anderson (1943) equation may be larger
than the actual value particularly for an unsaturated sand. It can be seen from Figure 13 that when the total suction in soil is
less than about 3000 kPa, the relative humidity at the soil surface approaches 100% (i.e., hs → 1 ). Consequently, AE = PE,
which is not valid for sand soil since the sand may have desaturated at a suction considerably below 100 kPa.
Figure 13 Relationship between relative humidity and total suction using Edlefsen and Anderson (1943) equation
PLAXIS LE - Groundwater provides options to improve the accuracy and stability of evaporative and atmospheric modeling.
Apply Surface Suction Correction
When the Surface Suction Correction option is selected, the total suction that is used to calculate the relative humidity at the
soil surface is adjusted based on an empirical correction factor (Alvenas and Jansson, 1997). In other words, the relative
humidity at soil surface is modified in accordance with the following expression:
u soil  − gv corr 

hs = v = exp   [68]
  w R(273.15 + Ts ) 
soil
uvo
 = (ua − uw ) +  [69]
corr = 10 fcorr [70]
where:
 = total suction, kPa,
ωv = molecular weight of water, 0.018 kg/mol,
w = unit weight of water, 9.807 kN/m3,
g = gravity acceleration, m/s2,
Ts = soil surface temperature, oC,
Ta = air temperature obtained from weather station, oC
ua = pore-air pressure, kPa,
 = osmotic suction, kPa.
corr = correction factor of soil surface suction ranging between 1 to 103.48, and
= empirical number changing from 0 to 3.48. By default, fcorr = 1.2, which
fcorr
corresponds to corr = 15.8.
Figure 13 illustrates the correction factor effect on the relative humidity.
There are two ways the surface suction adjustment factor can be applied: a) Manual input of fcorr value at Evaporation Properties
dialog, or b) estimation of the fcorr value based on soil suction relationship proposed by Fredlund et al. (2016). In the case of
manual input of fcorr value, the same suction correction factor will be applied to all the surface materials whereas the latter
option will calculate different fcorr values for each surface material based on its residual suction.
Fredlund et al., (2016) used the residual suction of the drying SWCC as a reference point to calculate the correction factor, corr.
This correction factor translates the SWCC to a suction value of 3000 kPa on the Lord Kelvin’s curve (Fredlund et al., 2016).
Equation [71] shows the relationship between the residual suction and the empirical number, fcorr and the equation is plotted
in Figure 14 The fcorr adjustment factors for various residual suction values. For coarse-grained soils, the maximum correction
factor, fcorr, is 3.48 and there is no correction factor required for fine-grained soils with a high air-entry value (Fredlund et al.,
2016).
𝑓𝑐𝑜𝑟𝑟 = 3.48 − 𝑙𝑜𝑔 𝜓𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙 [71]

where:
ψresidual = residual soil suction on the drying SWCC, kPa.
4.0
3.5
3.0
2.5
fcorr
2.0
1.5
1.0
0.5
0.0
1 10 100 1000 10000
Residual soil suction (kPa)
Figure 14 The fcorr adjustment factors for various residual suction values
Apply Maximum Gradient Limit at Evaporation Boundaries
When the Apply Maximum Gradient Limit option is selected, it is possible for extremely high gradients to develop at the upper
boundary during evaporative conditions when using the Wilson-Penman climate boundary condition. Extremely high surface
gradients can lead to unreasonable numerical instability. Limiting the gradient to a reasonable maximum value can improve
convergence of climate-based numerical models. A reasonable Gradient Limit might be between 50 to 1000; however, it could
go as high as 10,000.
The vapor pressure of uvsoil used in diffusive flow Equation [66] can be calculated using the following equation.
 gv corr 
uvsoil = uvo
soil soil
hs = uvo exp   [72]
  w R(273.15 + Ts ) 
The saturated vapor pressure is a function of soil surface temperature as given by the following expression (Lowe, 1977;
Gitirana, 2004):
uvsoil 2 3 4 5
0 = a0 + a1Ts + a2Ts + a3Ts + a4Ts + a5Ts [73]
where:
a0 = 0.6183580754
a1 = 0.0411427320
a2 = 0.0017217473
a3 = 0.0000174108
a4 = 0.0000003985
a5 = 0.0000000022
Note : Ts = Ta in this case.
7.3.4 Case 2 Solution: Wilson-Penman (1994) Model

In Case2, the PLAXIS LE - Groundwater model assumes that
• Moisture and vapor flow occurs through the soil.

• Soil temperatures in the entire domain are the same. In other words, the ground thermal flux is neglected (i.e.,
Qg = 0 ).
• The soil temperature at soil surface can be different from the air temperature. The heat exchanged between air
and soil surface follows the convection law as given the closed-form Equations [50] or [51]. However, the surface
temperature is not imposed as a boundary condition for distribution through the underlying soil.
• Actual Evaporation, AE, is calculated using Wilson-Penman (1994) equation.
7.3.4.1 Governing Equations

Moisture flow in one-dimension is defined by equation [74], (note: same as Equation [59]):
  
(
 k wy + kvh
y 
)h
y
− kvh  = − wm2w

h
t
[74]
Since the soil temperature is different from the air temperature, it is necessary to use another equation to determine the soil
temperature. Because this is an isothermal model (i.e., Qg = 0 ), controlled by Equations [50] and [51], a closed-form for the
soil temperature can be written (Wilson, 1994):
1
Ts = Ta + (Qn − AE ) [75]
C f  f (u)
where:
Ts = soil temperature at soil surface, oC,
Ta = air temperature, oC,
Cf = conversion factor (i.e., 1 kPa = 0.00750 mHg),
η = psychometric constant , 0.06733 kPa/oC,
f(u) = function depending speed, f(u) = 0.35 (1. + 0.146 Ww),
Ww = wind speed, km/hr,
Qn = net radiation, m/day, and
AE = actual evaporation, m/day.
7.3.4.2 Initial soil water content profile

The water content in the soil can be initialized using an initial set of values for pore-water pressure or water pressure head. It
is also possible to designate the location of the water table (i.e., the phreatic line) and assume that hydrostatic conditions exist
below and above the water table.
7.3.4.3 Initial soil surface temperature

Temperature can be initialized to the air temperature. The initialized air temperature is:
Ts = Ta [76]

The boundary condition for moisture flow at the ground surface is defined using equation [77].
qy = P − AE − Roff [77]
surface
where:
P = precipitation flux. m/day,
Roff = water runoff, m/day, and
AE = actual evaporation, m/day
7.3.4.5 Actual Evaporation

The Wilson-Penman (1994) equation for actual evaporation, AE, can be re-written as follows:
Qn +  Ea
AE =
 +  / hs [78]
where:
AE = actual evaporation, m/day,
PE = potential evaporation, m/day
Ea = flux associated with “mixing”; f(u) Cf uvair (1/hr – 1/hs), m/day,
f(u) = 0.35 (1. + 0.146 Ww),
Ww wind speed, km/hr,
Cf conversion factor, (i.e., 1 kPa = 0.00750 mHg),
hr = relative humidity in the air above the ground (i.e., hr = uvair/uvoair),
hs = relative humidity at the soil surface (i.e., hs = uvsoil/uvosoil),
uvair = water vapor pressure in the air above ground surface, kPa,
uvo air
= saturated vapor pressure at the mean air temperature, kPa,
uvsoil = vapor pressure in the soil at ground surface, kPa,
uvo soil
= saturated vapor pressure in the soil at ground surface, kPa,
Γ = slope of saturation vapor pressure vs. soil temperature curve (kPa/oC),
Qn = net radiation at the water surface, m/day, and
η = psychrometric constant, kPa/oC, η = 0.06733 kPa/oC.
The governing Equations [74] and [75], initial condition Equation [76], boundary condition Equations [77] and [78] are
essential for solving the isothermal model with atmospheric coupling.
Note: To calculate Actual Evaporation, AE, using Equation [78], the soil temperature at the ground surface must be known.
But when using Equation [75] to calculate the soil temperature, the Actual Evaporation, AE, must be known. In other word,
the Equations [75] and [78] are coupled to each other (i.e., must be solved simultaneously). The coupling Equation [75] is
one of governing equations. Therefore, it must be initialized to a specific value. The Actual Evaporation can be calculated initially
using Equation [78].
7.3.4.6 Determining Vapor Pressure Values throughout each Day

The saturated vapor pressure, u va0ir in Equation [78], is calculated based on air mean temperature.
uvair 2 3 4 5
0 = a0 + a1Ta + a2Ta + a3Ta + a4Ta + a5Ta [79]
where:
Ta = air mean temperature, measured at a weather station.
The parameters of a0, a1, a2, a3, a4 and a5 are previously given in Equation [73]. The water vapor pressure in the air, uvair, is
defined as follows:
uvair = uvair
0 hr [80]
where:
hr = relative humidity of the air above the soil surface.
The relative humidity, hr, of the air above the soil surface is measured at a weather station. The relative humidity, hs, and water
vapor pressure, uvsoil, at the soil surface are calculated using Equations [68] and [72]. It should be noted that soil surface
temperature, Ts, in Equations [68], [72], or [73] are calculated using Equation [75].
7.3.5 Case 3 Solution: “Experimental-Based” for AE/PE Model

For case 3, it is assumed that
• Moisture flow and vapor flow beneath the soil surface are driven by the hydraulic head gradient and the vapor
pressure gradient, respectively.
• Soil temperatures in the model domain are the same, and assumed to be equal to the air temperature above the
soil surface. In other words, the ground thermal flux can be neglected (i.e., Qg = 0 ).
• The soil surface temperature is assumed to be equal to the air temperature.

• The actual evaporation is calculated using an empirical “experimental-based” equation proposed by Wilson et al.,
(1997).
7.3.5.1 Governing Equation

The partial differential equation for the one-dimensional flow of water through a saturated –unsaturated soil is defined as
follows:
  h  h
y 
(
 kwy + kvh
y
)
− kvh  = − w m2w
 t
[81]
7.3.5.2 Initial water content conditions

The water contents of the soil can be initialized using an initial set of values for pore-water pressure or water pressure head. It
is also possible to designate the location of the water table (i.e., the phreatic line) and assume that hydrostatic conditions exist
below and above the phreatic line.

Using the water balance Equation [49], the boundary condition for moisture flow at the ground surface is defined as:
qy = P − AE − Roff [82]
surface
where
P = precipitation flux. m/day,
Roff = water runoff, m/day, and
AE = actual evaporation, m/day.
7.3.5.4 Determine an empirical expression for the ratio of actual evaporation, AE, to potential evaporation, PE
If the soil suction is known at the ground surface, then the rate of evaporation from the ground surface can be estimated from
the empirical “experimental-based” relationship shown in Figure 15.
Figure 15 Relationship between radio of AE/PE and total suction
The ratio of actual evaporation to potential evaporation, AE/PE, can be approximated using the form of the thermodynamic
equilibrium relationship between relative humidity and total suction (Edlefsen and Anderson, 1943). However, a correction
factor, δcorr, must be applied to the calculation of AE and the magnitude of the correction is dependent upon the type of soil
near the ground surface. The ratio of AE/PE has a format similar to that used for Equation [68].
The ratio of AE/PE has a format similar to that used for Equation [68].
 − gv corr 
AE / PE = exp   [83]

 w R (273.15 + Ts 
)
corr = 10 fcorr [84]
where:
AE = Actual evaporation, m/day,
PE = Potential evaporation, m/day,
ψ = total suction (i.e., matric suction plus osmotic suction), kPa,
ωv = molecular weight of water, 0.018 kg/mol,
w = unit weight of water, 9.807 kN/m3,
g = gravity acceleration, m/s2,
Ts = soil surface temperature, oC, and
fcorr = correction variable, and
δcorr = correction factor by which total suction must be multiplied.
The correction factor, δcorr, is computed based on the difference between the residual suction of the soil and a total suction of
3000 kPa. The variable fcorr, is determined based on the shift of the SWCC of the soil and Lord Kelvin’s thermodynamic equilibrium
equation. The fcorr, variable is typically about 1.8 for a coarse sand soil.
To include the relative humidity of the overlying air in equation [83], the equation can be modified to the following format:
 − gv corr 
AE / PE = exp   [85]
  (1 − ha ) w R (Ts + 273.15) 
where:
 = a dimensionless empirical parameter with a suggested value of 0.7,
and
ha = relative humidity overlying air.
Figure 16 shows the predicted values for the ratio of AE to PE when using Equation [67], equation [83], and Equation [85].
In Figure 16, the data for the “Limiting Function (1997)” is calculated using Equation [67]. The data for the Wilson-Penman
(1994) ratio of AE to PE is calculated using Equation [83]. The data for the empirical “experimental-based” ratio of AE to PE is
calculated using Equation [85].
1.2
0.8
Radio of AE/PE
0.6 AE/PE (Wilson Data)

Relative Humidity
Limiting Function (1997)
0.4 Modified empirical expression
0.2
0
1 10 100 1000 10000 100000 1000000
Total Suction (kPa )
Figure 16 Comparison of predicted values of AE/PE using different suggested equations
Note: Equation [85] is currently utilized in PLAXIS LE - Groundwater for obtaining the Case 3 Solution.
BENTLEY SYSTEMS Stream Traces 47 of 53
8 STREAM TRACES
For seepage models, the underground water flow velocity can be expressed using the following formulas. With vx and vy, the
direction of the stream trace can be determined at a specific point.
h
Vx = k x = flux _ x
x
[86]
h
Vy = k y = flux _ y
y
where:
kx = hydraulic conductivity in the x-direction,
ky = hydraulic conductivity in the y-direction,
h = total head,
flux_x = the flow in the x direction,
flux_y = the flow in the y direction.
BENTLEY SYSTEMS Particle Tracking 48 of 53
9 PARTICLE TRACKING
The Particle Tracking feature allows the user to track the advective movement of contaminant particles in a solution. Clicking a
location in the mesh will initiate the plotting of a path of a particle moving from the clicked location. The calculations assume
zero diffusion.
The particle tracking feature is only for transient models. The theory behind the particle tracking is defined as follows:
h
x1 = x0 +  wk x t
x
[87]
h
y1 = y0 +  wk y t
y
BENTLEY SYSTEMS Numerical Implementation 49 of 53
10 NUMERICAL IMPLEMENTATION
The PLAXIS LE - Groundwater module makes use of a finite element solver. The following sections provide a brief overview of
some of the more significant numerical methods used in the software.
Primary advantages of this finite element solver are as follows:
• Fully implicit approach in the solver, which provides for a robust solution of difficult models with convergence
issues,
• 2-noded line elements for 1D analysis, 3-noded triangle and 4-noded quadrilateral as elements for 2D analysis
and 4-noded tetrahedron elements for 3D analysis,
• Automatic generation and control of time steps,
• The implicit integration scheme using the Weighted Residual Method,
• Newton-Raphson convergence iteration schemes, and
• Uses a high performance linear solver engine allowing solution of large models with more than million nodes.
10.1 INTERPOLATION ORDER

Finite element basis is linear. The selected interpolation equation will have an impact on the number of nodes on the resulting
elements.
10.2 TIME STEPPING

Implicit method has been presented in literature as viable methods of solving for nodes while moving forward in time. Explicit
methods use known data to "explicitly" define each nodal value at an advanced time. Such explicit methods usually suffer in
accuracy and stability unless the time step is small. Implicit methods solve "implicitly" for self-consistent values at the advanced
time (i.e., the simultaneous finite element spatial equations are solved at the advanced time). This improves the stability and
accuracy for large time steps, at the cost of a simultaneous solution for all nodes in the mesh.
The implicit method has been selected in the software as the most robust method for general application. The research literature
supports this decision.
10.2.1 Time Integration

The transient equations mentioned above can be spatially discretized and the Galerkin’s principle of weighted residual method
can be applied for the finite element formualations. The numerical integration of the partial differential equations results in a
simpler expression of the equation:
Mh + Kh = q [88]
where：
M = mass matrix,
K = stiffness matrix,
q = flux vector reflecting the boundary conditions,
h = unknown hydraulic head vector, and
h = time derivative of hydraulic heads at the nodal points = h t
The time derivative of Equation [88] can be approximated using a finite difference technique. The relationship between the
nodal heads of an element at two successive time steps can be expressed using the backward difference approximation (implicit
method):
 M M
K +  ht +t = ht + q [89]
 t  t
The finite element flow equation (i.e., Equation [88]) can be written for each element and assembled to form a set of global
flow equations. The set of global flow equations for the whole system is then solved using the sparse linear solver for the
hydraulic heads at the nodal points, h. However, Equation [88] is nonlinear because the coefficients of permeability are a
function of matric suction in unsaturated seepage analysis, which is related to the hydraulic head at the nodal points.
The hydraulic heads are the unknown variables in Equation [88] or [89]. Therefore, Equation [89] must be solved using an
iterative procedure that involves a series of successive approximations. In the first iteration, the coefficients of permeability are
BENTLEY SYSTEMS Numerical Implementation 50 of 53
estimated based on the initial conditions of the hydraulic head or water table in order to calculate the fist set of hydraulic heads
at the nodal points. The computed hydraulic heads are used to calculate the average matric suction in the element. The
coefficient of permeability is adjusted in the subsequent iteration with the new set of matrix suctions. The adjusted set of
coefficients of permeability is then used in the finite element formulations to calculate the new set of hydraulic heads. The above
procedure is repeated until the difference in hydraulic heads at two successive iterations are smaller than a user specified
tolerance.
10.2.2 Convergence Criteria
The convergence criterion used is known as relative residual norm criterion. The solution is deemed to have converged when
the following criterion is satisfied:
h k +1 − h k
 tol [90]
hk
where h k +1 is the solution of Equation [89] at (k+1)th iteration for a particular time step and hk is the solution at kth iteration.
The tolerance (tol) is a user specified parameter. Once a solution for a particular time step satisfies the criterion [90], the
calculation for the time step concludes and moves forward for the next time step.
The convergence rate is highly dependent on the degree on nonlinearity of the permeability function of the materials and spatial
discretization of the problem. A steep permeability function requires more iterations and a larger convergency tolerance. A finer
discretization in both element size and time step will assist in obtaining convergence faster with a smaller tolerance (Fredlund
and Rahardjo, 1993).
10.2.3 Time-Step Control
Staged analysis is obtained by discretizing a stage into several small time-steps and using an incremental analysis. The time-
step for each sub-step is controlled by the following:
• The time discretization set by the user
• The time discretization set by the user for the outputs such as Contour plots and Graph plots or for any other printed
output of the simulations results
• The time discretization of the applied boundary conditions
• The time discretization based on automatic adjustment of the time-step based on the global error in solution in each
subsequent iteration
The adaptive time steps algorithm implemented is based on a slightly modified form of the error control algorithm proposed by
Söderlind (2002). The solver adopts the next time-step considering above controls. If the solution failed to converge within the
maximum iteration set by the user, the time-step is reduced based on the global error in the solution and re-try using a smaller
time step. The above procedure repeats until the solution converges. If the time increment reaches to the “minimum” time
cincrement set by the user in the Stage Settings dialog, the solution is deemed to have failed to converge. If such situation
occurs, the user is advised to review the material properties if there are any materials with a sharp changes in coefficients of
permeability. Similarly, the initial conditions, boundary conditions, and the mesh density can also be reviewed for problems
with convergence issues.
BENTLEY SYSTEMS References 51 of 53
11 REFERENCES
Alvenas, G. and Jansson, P. E. (1997). Model for evaporation, moisture and temperature of bare soil: calibration and sensitivity
analysis. Agricultural and Forest Meteorology, 88:47-56.
Anderson, E.A. (1976). A point energy and mass balance model of a snow cover. Ph.D. Thesis, Department of Civil Engineering,
Stanford University, Palo Alto, CA, USA.
Barbour, S.L., Fredlund, D. G., Gan, J. K-M., and Wilson, G.W. (1991). Prediction of moisture movement in highway subgrade
soils, 45th Canadian Geotechnical Conference, Toronto, ON., Canada.
Bathe, K.J. (1982). Finite element procedures in engineering analysis. Prentice-Hall.
Bear, J. (1972). Dynamics of fluids in porous Media. American Elsevier Publishing Company, New York.
Benson, C. H., Gunter, J. A., Boutwell, G. P., Trautwein, S. J., and Berzanskis, P. H. (1997). Comparison of four methods to
assess hydraulic conductivity. Journal of Geotechnical and Geoenvironmental Engineering, 123(10), 929-937.
Bowles, J.E., (1984). Physical and geotechnical properties of soils. 2nd ed. McGraw-Hill, New York.
Brooks, R., and Corey, T. (1964). Hydraulic properties of porous media. Hydrology Papers, Colorado State University.
Burdine, N. T. (1953). Relative permeability calculations from pore size distribution data. Trans. Am. Inst. Min. Metall. Pet.
Eng., 198, 71–78.
Campbell, J. D. (1973). Pore pressures and volume changes in unsaturated soils. Ph.D. thesis, University of Illinois at Urbana-
Champaign.
Celia, M.A. and Bouloutas, E.T. (1990). A general mass-conservative numerical solution for the unsaturated flow equation.
Water Resources Research, Vol. 26, No. 7, pp. 1483-1496, July.
Chandrakant, S. Desai and John F. Abel, (1972). Introduction to the Finite Element Method. Van Nostrand Reinhold Company,
New York, N.Y.
Chapuis, Robert P., Chenaf, D., Bussiere, B., Aubertin, M., Crespo, R. (2001). A user’s approach to assess numerical codes for
saturated and unsaturated seepage conditions. Canadian Geotechnical Journal, 38: 1113-1126.
Crespo, R., (1993). Modelisation par elements finis des ecoulements a travers les ouvrages de retenue et de confinement des
residus miniers. M.Sc.A. thesis, Ecole Polytechnique de Montreal, Montreal.
Dakshanamurthy, V., and Fredlund, D. G. (1981). A mathematical model for predicting moisture flow in an unsaturated soil
under hydraulics and temperature gradient. Water Resources, vol. 17, No. 3, pp. 714-722.
Dupuit, J., (1863). Etudes theoriques et pratiques sur le mouvement des eaux dans les canaux decouverts et a travers les
terrains permeables. 2nd ed. Dunod, Paris.
Ebrahimi-Biring, N.,Gitirana, Jr., G. F. N., Fredlund., D. G., Fredlund, M. D., and Samarasekera, L. (2004). A lower limit for the
water permeability coefficient, Proceedings of the fifty-seventh Canadian Geotechnical Conference, Quebec, QC, Vol. 1,
pp. 12-19.
Edlefsen, N., and Anderson, A. (1943). Thermodynamics of soil moisture. California Agriculture, 15(2), 31-298.
Fenton, G. A., and Vanmarcke, E. H. (1990). Simulation of random fields via local average subdivision. Journal of Engineering
Mechanics, 116(8), 1733-1749.
Fredlund, D. G. (1991). How negative can pore-water pressures get? Geotechnical news, 9(3), 44-46.
Fredlund, D. G., and Rahardjo, H. (1993). Soil Mechanics for unsaturated soils. John Wiley & Sons, Inc.
Fredlund, D. G., and Xing, A. (1994). Equations for the soil-water characteristic curve. Canadian Geotechnical Journal, 31(4),
521-532.
Fredlund, D. G., Xing, A., and Huang, S. (1994). Predicting the permeability function for unsaturated soils using the soil-water
characteristic curve. Canadian Geotechnical Journal, 31(4), 533-546.
Fredlund M.D. (1996). Design of a Knowledge-Based System for Unsaturated Soil Properties. M.Sc. Thesis, University of
Saskatchewan, Saskatoon, Saskatchewan, Canada.
Fredlund, M. D., Fredlund, D. G., and Wilson, G. W. (2000). An equation to represent grain-size distribution. Canadian
Geotechnical Journal, 37(4), 817-827.
Fredlund, D. G., and Gitirana, Jr., de F. N. G. (2005). Unsaturated soil mechanics as a series of partial differential equations.
In Proceedings of the International Conference on Problematic Soils, Famagusta, Northern Cyprus, pp. 3-30.
Fredlund, D. G., Rahardjo, H., and Fredlund, M. D. (2012). Unsaturated Soil Mechanics in Engineering Practice. John Wiley &
Sons, Inc.
Fredlund, M. D., Tran, D., and Fredlund, D. G. (2016). Methodologies for the Calculation of Actual Evaporation in Geotechnical
Engineering. International Journal of Geomechanics, 16(6), D4016014.
Freeze, R. Allan and John A. Cherry, (1979). Groundwater. Prentice–Hall, Inc., Englewood Cliffs, New Jersey.
Gardner, W. R. (1958). Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation
from a water table. Soil science, 85(4), 228-232.
Gear C.W. (1971). Numerical initial value problems in ordinary differential equations. Prentice Hall.
Granger, R. J. (1989). A complementary relationship approach for evaporation from nonsaturated surfaces. Journal of
Hydrology, 111(1-4), 31-38.
Gray, D. M., (1970). Handbook on the principles of hydrology. Canadian National Committee for the International Hydrological
Decade.
Gitirana, Gilson Jr. (2004). Weather related geo-hazard assessment model for railway embankment stability, PhD Thesis,
University of Saskatchewan, Saskatoon, SK, Canada.
Gustafsson, D., Stähli, M., and Jansson, P. E. (2001). The surface energy balance of a snow cover: comparing measurements
to two differentsimulation models. Theoretical and Applied Climatology, 70(1), 81-96.
Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P. J., and Vachaud, G. (1977). A comparison of numerical simulation models
for one-dimensional infiltration. Soil Science Society of America Journal, Vol. 41, No. 2.
Huang, S., Barbour, S. L., & Fredlund, D. G. (1994, June). A history of the coefficient of permeability function. In Sino-Canadian
Symposium on Expansive Soils/Unsaturated Soils, Chinese Institution for SoilMechanics and Foundation Engineering,
CEES (pp. 57-80).
Hwang, C. T. (1976). Predictions and observations on the behaviour of a warm gas pipeline on permafrost. Canadian
Jame, Y.W. (1977). Heat and mass transfer in freezing soil. Ph.D Thesis, University of Saskatchewan.
Kimbal, B. A., Jackson, R. D., Reginato. R. J., Nakayama, F. S., and Idso, S. B. (1976). Comparison of field-measured and
calculated soil-heat fluxes. Soil Science Society of America Proceedings, Vlo. 40, No. 1, pp. 18-25.
Kozeny, J. (1927). Uber Kapillare Leitung Des Wassers in Boden. Sitzungsber Akad. Wiss., Wien Math.Naturwiss.Kl., Abt.2a,
136:271–306.
Leong, E. C., and Rahardjo, H. (1997). Permeability functions for unsaturated soils. Journal of Geotechnical and
Geoenvironmental Engineering, 123(12), 1118-1126.
Lowe, P. R. (1977). An approximating polynomial for the computation of saturation vapor pressure. Journal of Applied
Meteorology, 16(1), 100-103.
Kunze, R. J., Uehara G. and Graham K. (1968): Factors important in the calculation of hydraulic conductivity. In Proc. Soil Sci.
Soc. Amer., 32, 760-765.
Kuusisto, E. (1980). On the values and variability of degree-day melting factor in Finland. Hydrology Research, 11(5), 235-242.
MEND (1993). SoilCover user’s manual for evaporative flux model. University of Saskatchewan, Saskatoon, SK, Canada.
Mualem, Y. (1976). A catalogue of the hydraulic properties of unsaturated soils. Technion Research and Development
Foundation.
Pentland, J. S., (2000). Use of a general partial differential equation solver for solution of heat and mass transfer problems in
soils. University of Saskatchewan, Saskatoon, SK, Canada.
Penman, H.L., (1948)., Natural evapotranspiration from open water, bare soil and grass. Proc. R. Soc. London Ser. A. 193:
120-145.
Philip, J. R., and de Vries, D. A. (1957). Moisture movement in porous materials under temperature gradients. Transactions,
American Geophysical Union, Vol. 38, No. 2, pp. 232.
Priestley, C. H. B., and Taylor, R. J. (1972). On the assessment of surface heat flux and evaporation using large-scale
parameters. Monthly weather review, 100(2), 81-92.
Ritchie, J. T. (1972). Model for predicting evaporation from a row crop with incomplete cover. Water resources research, 8(5),
1204-1213.
Saxton, K. E. (1982). Mathematical modelling of evapotranspiration on agricultural watersheds. In Proceeding of Modeling

Components of the Hydrologic Cycle. 183-203.
Söderlind, G. (2002). Automatic control and adaptive time-stepping. Numerical Algorithms, 31(1–4), 281–310.
Tetens, V.O. (1930). Uber einige meteorologische. Begriffe, Zeitschrift fur Geophysik. 6:297-309.
Todd, D. K. (1980). Groundwater hydrology. 2nd ed. John Wiley & Sons, New York.
Thornthwaite, C. W. (1948). An approach toward a rational classification of climate. Geographical review, 38(1), 55-94.
Tratch, D. J., (1995). A geotechnical engineering approach to plant transpiration and root water Uptake, University of
Saskatchewan, Saskatoon, SK, Canada.
Tratch, D. J., Wilson, G.W., and Fredlund, D.G. (1995). An introduction to analytical modeling of plant transpiration for
geotechnical engineers. In Proceedings of 48th Canadian Geotechnical Conference, Vancouver, BC, Canada, Vol. 2, pp.
771-780.
Van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil science
society of America journal, 44(5), 892-898.
Vereecken H., J. Maes, J. Feyen, and P. Darius, (1989). Estimating the soil moisture retention characteristic from texture, bulk
density, and carbon content. Soil Science, Volume 148, No. 6
Wilson, W. G. (1990). Soil evaporative fluxes for geotechnical engineering problems. PhD Thesis, Department of Civil
Engineering, University of Saskatchewan, Saskatoon, SK, Canada.
Wilson, W. G., Fredlund, D. G., and Barbour, S. L. (1994). Coupled soil-atmosphere modelling for soil evaporation. Canadian
Wilson, W. G., Fredlund, D. G. and Barbour, S. L. (1997). The effect of soil suction on evaporative fluxes from soil surfaces.
Canadian Geotechnical Journal, 34(4): pp. 145-155.
Wyllie, M. R. J. and Gardner, G. H. F. (1958). The generalized Kozeny-Carman equation II. A novel approach to problems of
fluid flow. World Oil Prod. Sect. 210-228
Zapata, C. E., Houston, W. N., Houston, S. L., and Walsh, K. D. (2000). Soil-water characteristic curve variability. Advances in
Unsaturated Geotechnics, Denver, Colorado, August 5-8, pp. 84-123.

PLAXIS LE - Groundwater - Theory - Manual

Uploaded by

Copyright:

Available Formats

PLAXIS LE - Groundwater - Theory - Manual

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

PLAXIS LE - Groundwater - Theory - Manual

Uploaded by

Copyright:

Available Formats

Last Updated: Thursday, December 16, 2021

Bentley Systems Incorporated

3 PDE’S FOR SEEPAGE ANALYSIS...................................................................................................................................... 11

4 BOUNDARY CONDITIONS IN SEEPAGE .......................................................................................................................... 14

7 SOIL ATMOSPHERE MODELING .................................................................................................................................... 32

7.3.3.3 Boundary condition for moisture flow ................................................................................................................................37

8 STREAM TRACES ........................................................................................................................................................... 47

The primary attributes of the solution process are as follows:

• Fully automatic mesh generation,

2 FUNDAMENTALS OF SEEPAGE THEORY

o Apply conservation laws,

• Soil phases can be described using a continuum mechanics approach,

2.1 CONSERVATION OF MASS

w = density of water, 1000 kg/m3,

2.1.1 Changes in Volume of Stored Water

Vw /Vo = volumetric water content,

S = degree of saturation, and

2.2 FLOW LAWS

2.2.1 Flow of Liquid Water

h = hydraulic head, m, which is equal to uw + y

uw = pore-water pressure, kPa,

2.2.2 Anisotropic Flow of Liquid Water

kwxx = kw1 + kw2 + kw1 − kw2 cos 2

kwxy = kw1 − kw2 sin 2

Figure 3 Principal directions of anisotropy in two-dimensions

3 PDE’S FOR SEEPAGE ANALYSIS

3.1 ONE-DIMENSIONAL SEEPAGE

• hydraulic conductivity function, kw,

3.2 TWO-DIMENSIONAL SEEPAGE

3.3 THREE-DIMENSIONAL SEEPAGE

3.4 UNSATURATED FLOW EQUATION

3.4.1 H-Based Formulation

4 BOUNDARY CONDITIONS IN SEEPAGE

• Natural (or Neumann) boundary condition: flux,

4.1 NATURAL BOUNDARY CONDITION

Natural boundary condition =  ( v  n ) dS

4.1.1 Differences Between Actual and Applied Boundary Flux

4.2 ESSENTIAL BOUNDARY CONDITION

Essential Boundary Condition = h [15]

Forced negative Phreatic

4.2.1 Review Boundary Condition (Drain)

Figure 5 “Review boundary” condition for a free flow surface

A review boundary condition is applied as follows:

4.3 GRADIENT OR UNIT GRADIENT BOUNDARY CONDITION

4.4 INITIAL WATER TABLE

when h>y uw > 0; represents the saturated region

5.1 FLUX SECTION THEORY CALCULATION

The D matrix for an isotropic material is defined as:

The B matrix is defined as:

The H is the total head vector of the element nodes:

The B matrix can be rewritten as:

The normal flux across the flux section is determined as:

Figure 6. Flux section calculation for a 3-noded element

Figure 7. 3-noded element in local (natural) coordinates

6.1 SOIL-WATER CHARACTERISTIC CURVE

6.1.1 Fredlund and Xing (1994) Equation

mf = material parameter which is primarily a function of the residual water

Fitting method: Least squares nonlinear regression