The Slope Cube Modulefor HYDRUS (2D/3D)
Simulating Slope Stress and Stability
in Variably-Saturated Hillslopes

ManualVersion 1.0

Ning Lu, Alexandra Wayllace, and Pan Chen
August 25, 2019
Soil Water Retention, LLC
Madison, Wisconsin, USA© 2019 Ning Lu. All rights reserved.

DISCLAIMER
The Slope Cube module was developed by the SWR, LLC as a supplemental module of the
HYDRUS (2D/3D) software package to simulate transient fields of soil moisture, soil suction,
suction stress, total and effective stresses, and local factor of safety. The computational codeof the Slope Cube module (written in Fortran) was obtained by coupling the standard
computational 2D module of HYDRUS (flow and transport) with an adopted code from a twodimensional hydro-mechanical model (stress and slope stability) developed by Lu et al. (2010).
While this source code is not available to the public (and is proprietary/copyrighted), theoriginal source code in Lu et al. (2010) is not copyrighted, and can be obtained from the seniorauthor. The software has been verified against selected test cases. However, no warranty isgiven that the program is error-free. If you do encounter problems with the code, find errors,or have suggestions for improvement, please contact the senior author at:
Ning Lu
Tel: (303)-995-3980
Soil Water Retention, LLC
www.soilwaterretention.com
Email: ninglu@soilwaterretention.com

iAbstract

Lu, N., Wayllace, A., and Chen, P. The Slope Stress and Stability (Slope Cube) Module for
HYDRUS (3D): Simulating Slope Stress and Stability in Variably-Saturated Soils, Hydrus
Software Series 3.02.0290, PC-Progress.
Rainfall induced shallow landslides are pervasive geologic hazards worldwide, causing
significant damages involving loss of life and properties significant to the GDP of somecountries. Currently, several common research approaches are used to better understand andpredict such geologic hazards, namely statistically or/and physically based analyses. Thismanual documents Slope Stress and Stability, or Slope Cube (2D/3D), a physically basedhydro-mechanical framework focusing on assessment and prediction of the initiation andlocations of landslide occurrences under rainfall conditions. Slope Cube couples variably
saturated flow and stress fields to simulate transient hillslope hydrologic and mechanicalresponses to rainfall events. Three recent advancements are implemented in Slope Cube thatare beyond the traditional physically based landslide models. The first is the utilization of thesuction stress-based effective stress to unify the description of effective stress distribution in
hillslopes under variably-saturated conditions. The second is calculation of deformations due
to changes in suction stress-based effective stress. The third is to employ a recently establishedconcept of the local factor of safety to capture evolution of stress paths toward the failure statein hillslopes. This document aims to describe Slope Cube module for the HYDRUS (2D/3D)
software package. Slope Cube, together with HYDRUS (2D/3D), simulates two or threedimensional transient fields of soil moisture, soil suction, suction stress, total and effective
stresses, deformation, and local factor of safety.
The document serves as a Technical Manual and a User Manual as well as a reference of theGraphical User Interface of Slope Cube related parts of the HYDRUS software package.

ii

Table of Contents
DISCLAIMER.......................................................................................................................... i
Abstract................................................................................................................................... ii
Table of Contents .................................................................................................................. iii
List of Figures .......................................................................................................................... v
List of Tables .......................................................................................................................... ix
1. Introduction ...................................................................................................................... 1
2. Governing Equations for Slope Cube ............................................................................. 2
2.1. Hydrological Model ................................................................................................ 3
2.2. Suction Stress-based Effective Stress Model .............................................................4
2.3. Mechanical Model ................................................................................................... 5
2.4. Equation for deformation moduli ...........................................................................92.5. Local Factor of Safety Model ................................................................................ 10
3. Pre-Processing and Input Data ...................................................................................... 13
3.1. Slope Cube Module Activation............................................................................... 13
3.2. Definition of Material Properties ................................................................... 13
3.3. Definition of Boundary Conditions ................................................................. 14
4. Post-Processing .................................................................................................................. 174.1. The Slope Cube Module Outputs ........................................................................... 17
5. Example Problems for Slope Cube (2D) ......................................................................... 21
5.1. Example 1: Stability Analysis of a Homogeneous Hillslope under a ConstantRainfall Infiltration rate ......................................................................................... 22

.

5.2.

Example 2: Stability Analysis of a Layered Hillslope under a Constant
Infiltration Rate .............................................................................................. 29

5.3.

Example3: Stability Analysis of a Homogeneous Hillslope under a Time-Varying
Heavy Rainfall Event with a Concentrated Mechanical Load ............................. 37

6. Example Problems for Slope Cube (3D) ......................................................................... 44
6.1. Example 1: Stability Analysis of a Homogeneous Hillslope under a ConstantRainfall Infiltration ................................................................................................. 45
6.2.

Example 2: Stability Analysis of a Hillslope under a Time-Varying Heavy
Rainfall Event with a Concentrated Mechanical Load.................................... 55

6.3.

Example 3: Deformation of a Soil Cake due to Change in Water Content,
Comparison of Measured and Predicted Volumetric Change of Soil Cake .......... 63

APPENDIX-A: Finite Element Discretization for 3D Case ..............................................69
References ............................................................................................................................. 79iii

List of Figures
Figure 1. Illustration of the physical processes in Slope Cube module .......................... 2
Figure 2. Illustration of soil water retention curves for typical sand, silt, and clay
(after Lu et al., 2014). ................................................................................... 4
Figure 3. Illustration of hydraulic conductivity functions for typical sand, silt, andclay (after Lu et al. 2014). ............................................................................. 4
Figure 4. Illustration of suction stress characteristic curves for typical sand, silt, andclay (after Lu et al., 2014). ............................................................................ 5
Figure 5. Comparison of the measured data and predicted curve using Eq. (18) for
Young’s modulus of the clayey soil under the variation of water content. .... 10
Figure 6. Conceptual illustration of the local factor of safety (after Lu et al. 2012). .... 11
Figure 7. The Main Processes and add-on Modules dialog window: check the box“Slope Cube” for activation of Slope Cube module ...................................... 13
Figure 8. Slope Cube module: defining materials properties ....................................... 13
Figure 9. Slope Cube module: defining boundary conditions from the “Boundary
Condition” panel in the “Project Navigator – Data” window ........................ 14
Figure 10. Slope Cube module, boundary condition types for: (a) Slope Cube (2D),
and (b) Slope Cube (3D). ........................................................................... 15
Figure 11. Selection of output data for the (a) Slope Cube (2D), and (b) Slope
Cube (3D). ................................................................................................. 17
Figure 12. (a) Customization of the color scales of Slope Cube module, (b) Selectionof contour values in the display of Slope Cube module output by clicking
on the map with the cursor .......................................................................... 19
Figure 13. Example of graphical display of the suction stress (a), and local factorof safety (b) for the observation nodes....................................................... 20
Figure 14. Illustration of hillslope geometry and boundary conditions forhydrological (in black) and mechanical (in gray) models .......................... 22
Figure 15. Contour plots of pressure head [m] and water content evolutions with timefor Example 1a. Results are shown at 0, 3, 6, and 15 hours........................ 23
Figure 16. Contour plots of maximum and minimum effective stress [kPa] evolutions
with time for Example 1a. Results are shown at 0, 3, 6, and 15 hours ........ 24
Figure 17. Contour plots of suction stress [kPa] and field of the local factor of safety
evolutions with time for Example 1a. Results are shown at 0, 3, 6,
iv

and 15 hours ............................................................................................. 25
Figure 18. Contour plots of pressure head [m] and water content evolutions with timefor Example 1b. Results are shown at 0, 30, 60, and 150 hours.................. 26
Figure 19. Contour plots of maximum and minimum effective stress [kPa] evolutions
with time for Example 1b. Results are shown at 0, 30, 60,
and 150 hours .......................................................................................... 27
Figure 20. Contour plots of suction stress [kPa] and field of the local factor of safety
evolutions with time for Example 1b. Results are shown at 0, 30, 60,
and 150 hours ........................................................................................... 28
Figure 21. Comparisons between the failure surface evolution for Example 1aand Example 1b. Local factor of safety <1 is shown in red.
Results are shown at 0, 3, 6, and 15 hours for Example 1a
and at 0, 30, 60, and 150 hours for Example 1b .......................................... 29
Figure 22. Illustration of hillslope geometry and boundary conditions for hydrological(in black) and mechanical (in gray) models for Example 2. ....................... 30
Figure 23. Contour plots of pressure head [m] and water content evolutions with timefor Example 2a. Results are shown at 0, 3, 6, and 15 hours........................ 31
Figure 24. Contour plots of maximum and minimum effective stress [kPa] evolutions
with time for Example 2a. Results are shown at 0, 3, 6,
and 15 hours ............................................................................................ 32
Figure 25. Contour plots of suction stress [kPa] and field of the local factor of safety
evolutions with time for Example 2a. Results are shown at 0, 3, 6
and 15 hours ............................................................................................. 33
Figure 26. Contour plots of pressure head [m] and water content evolutions withtime for Example 2b. Results are shown at 0, 30, 60,
and 150 hours .......................................................................................... 34
Figure 27. Contour plots of maximum and minimum effective stress [kPa]
evolutions with time for Example 2b. Results are shown at 0, 30,
60, and 150 hours…………………………………………………………. 35
Figure 28. Contour plots of suction stress [kPa] and field of local factor of safety
evolutions with time for Example 2b. Results are shown at 0, 30,
60, and 150 hours ...................................................................................... 36
Figure 29. Comparisons between the failure surface evolutions for Example 2aand Example 2b; Local factor of safety <1 is shown in red. Results areshown at0, 3, 6, and 15 hours for Example 2a and at 0, 30, 60, and 150
hours for Example 2b. .............................................................................. 37
Figure 30. Illustration of hillslope geometry and boundary conditions for hydrologicalv

(in black) and mechanical (in gray) models for Example 3. The top layer(in white) represents the sandy soil and the lower layer in grey represents
the silty soil. The concentrated mechanical load is in red ........................... 38
Figure 31. Time-varying rainfall boundary condition used in Example 3 .......................... 38
Figure 32. Contour plots of pressure head [m] and water content evolutions with timefor Example 3. Results are shown at 0, 34, 50, and 70 hours ............................ 39
Figure 33. Contour plots of maximum and minimum effective stress [kPa]
evolutions with time for Example 3. Results are shown at 0, 34, 50, and70 hours……………………………………………………………………..40
Figure34. Contour plots of suction stress [kPa] and field of local factor of safety evolutions with time for Example 3. Results are shown at 0, 34, 50,
and 70 hours………………………………………………………………...41
Figure 35. Contour plots of total displacement [m] evolutions with time for Example 3.
Results are shown at 0, 34, 50, and 70 hours…
.. 42
Figure 36. Contour plots of the local factor of safety < 1 for failure surface evolutions
with time for Example 3. Results are shown at 0, 34, 50, and 70 hours….. 43
Figure 37. Hillslope geometry and boundary conditions for hydrological (in black)and mechanical (in gray) models…………………………………………. 45
Figure 38. Contour plots of pressure head [m] and water content evolutions withtime for Example 1a. Results are shown at 0, 4, 12, and 20 hours……….47
Figure 39. Contour plots of maximum and minimum effective stress [kPa] evolutions
with time for Example 1a. Results are shown at 0, 4, 12, and 20 hours… 48
Figure 40. Contour plots of suction stress [kPa] and field of local factor of safety
evolutions with time for Example 1a. Results are shown at 0, 4, 12, and20 hours………………………………………………………………….. 49
Figure 41. Graphical display of observation nodes. (a) pressure head as a functionof time, (b) water content as a function of time, (c) suction stress as afunction of time, and (d) local factor of safety as a function of time forthe three observation nodes. ………………………………………………50
Figure 42. Contour plots of pressure head [m] and water content evolutions withtime for Example 1b. Results are shown at 0, 40, 120, and 200 hours…… 51
Figure 43. Contour plots of maximum and minimum effective stress [kPa] evolutions
with time for Example 1b. Results are shown at 0, 40, 120,and 200 hours ........................................................................................... 52
Figure 44. Contour plots of suction stress [kPa] and field of local factor of safety
evolutions with time for Example 1b. Results are shown at 0, 40, 120,and 200 hours ............................................................................................ 53vi

Figure 45. Comparisons between the failure surface evolution for Example 1a andExample 1b; Local factor of safety <1 is shown in red. Results are shownat 0, 4, 12, and 20 hours for Example 1a and at 0, 40, 120, and 200 hours
for Example 1b ........................................................................................ 54
Figure 46. Illustration of hillslope geometry and boundary conditions for hydrological(in black) and mechanical (in gray) models for Example 2. ....................... 55
Figure 47. Time-varying rainfall boundary condition used in Example 2 ..................... 55
Figure 48. Contour plots of pressure head [m] and water content evolutions with timefor Example 1a. Results are shown at 0, 34, 50, and 70 hours ................... 56
Figure 49. Contour plots of maximum and minimum effective stress [kPa] evolutions
with time for Example 1a. Results are shown at 0, 34, 50, and 70 hours .... 57
Figure 50. Contour plots of suction stress [kPa] and field of local factor of safety
evolutions with time for Example 1a. Results are shown at 0, 34, 50, and70 hours ..................................................................................................... 58
Figure 51. Graphical display of observation nodes. (a) pressure head as a function oftime, (b) water content as a function of time, (c) suction stress as afunction of time, and (d) local factor of safety as a function of time for thethree observation nodes ............................................................................. 59
Figure 52. Illustration of hillslope geometry and boundary conditions for hydrological(in black) and mechanical (in gray) models for Example 2b. ..................... 60
Figure 53. Contour plots of maximum effective stress [kPa] and total displacement [m]evolutions with time for Example 2b. Results are shown at 0, 32, 56, and72 hours ..................................................................................................... 61
Figure 54. Contour plots of suction stress [kPa] and field of local factor of safety
evolutions with time for Example 1a. Results are shown at 0, 34, 50, and70 hours ..................................................................................................... 62
Figure 55. Soil cake geometry and boundary conditions for hydrological (in black)and mechanical (in gray) models ................................................................ 64
Figure 56. Comparison of the measured and fitted Young’s modulus of Denver
Claystone .................................................................................................. 65
Figure 57. Comparison of the measured and fitted SSCC of Denver claystone............. 66
Figure 58. Comparison between the simulated and measured porosity ........................ 66

Figure A-1. Tetrahedral element for three-dimensional problems ............................... 70vii

List of Tables
Table 1.

List of mechanical parameters for Slope Cube module ................................. 14

Table 2.

List of types of mechanical boundary conditions for Slope Cube (2D) ......... 15

Table 3.

List of types of mechanical boundary conditions for Slope Cube (3D)……... 16

Table 4.

List of output data files for Slope Cube (2D)................................................... 18

Table 5.

List of output data files for Slope Cube (3D)................................................ 18

Table 6.

List of examples for illustration of use of Slope Cube module.. ........................ 21

Table 7.

List of hydro-mechanical parameters defined in Example 1. ............................ 23

Table 8.

List of hydro-mechanical parameters defined in Example 2a and Example 2b…. 30

Table 9.
Table 10.

List of examples for illustration of use of Slope Cube module.. ................ 44

Table 11.

List of hydro-mechanical parameters defined in Example 2b.. .................... 60

Table 12.

Comparison of running time for different sizes of element mesh.. .............. 63

Table 13.

List of hydro-mechanical parameters used in Example 3.. .......................... 64

Table 14.

Steps for the water content change used in Example 3.. .............................. 65

List of hydro-mechanical parameters defined in Example 1........................ 46viii

1. Introduction
Precipitation-induced landslides are one of the most pervasive environmental hazards and poseserious menace for the public safety. Two common research methodologies are used to better
understand and predict shallow landslides: statistically or physically based analyses. While thestatistical methods seek empirical correlations between the occurrence of landslides andenvironmental and geologic factors such as rainfall and geomorphologic characteristics andphysical methods emphasize the quantitative linkages between the occurrences of landslides
and hydrologic and mechanical processes in hillslopes, both methodologies can be applied tohillslope and regional scales with some success. To date, accurate predictions of landslide
initiation and location remain challenging, elusive, and unsettled.
Physically based models (e.g. Montgomery and Dietrich 1994; Baum et al. 2008; Lu and Godt,
2013, Chen et al. 2017) synthetize the interaction between hydrology, geomorphology, and soil
mechanics to improve the predictability of shallow landslides. In general, they include ahydrological component to simulate infiltration and groundwater flow processes, and a soil
slope-stability component to simulate the safety factor of hillslopes. Several models areavailable in the literature with different degrees of complexity. Typical hydrologicalcomponents range from steady state subsurface flow (e.g., Montgomery and Dietrich, 1994;
Lu and Godt, 2008) to transient subsurface flow (e.g., Simoni et al., 2008; Godt et al., 2009).
Common slope stability components are based on limit-equilibrium methods developed
extensively in geotechnical engineering (e.g., Fellenius, 1936; Janbu, 1973; Duncan andWright, 2003), or shear strength reduction analysis (e.g., Smith and Griffiths, 2004). Recently,
a stress field based framework developed by Lu et al. (2012) allows us to compute a scalar
field of factor of safety within a hillslope. Here the factor of safety is defined at each point of
the hillslope, instead of a single factor of safety for the whole slope traditionally employed in
all limit-equilibrium based methods (e.g., Chen et al., 2018b).
We have implemented the hydro-mechanical Slope Stress and Stability (Slope Cube)
framework systematically described in Lu and Godt (2013) as an add-on module of HYDRUS(2D/3D) [Šimunek et al., 2008]. Slope Cube can be used to quantitatively simulate, assess, andpredict the timing and locations of landslide initiation in variably saturated hillslopes.
Slope Cube couples variably-saturated flow and stress fields to simulate transient hillslope
hydrology and stability in response to transient rainfall or precipitation. Three advancements
are implemented in Slope Cube that are beyond the traditional physically based landslideprediction models. The first is the utilization of the suction stress-based effective stress to unify
the description of effective stress distribution and its resulting deformation in hillslopes undervariably-saturated conditions. The second is to calculate deformation due to changes of suction
stress-based effective stress (Chen and Lu, 2018a). The third is to employ the aforementioned
concept of local factor of safety to quantitatively and closely capture the evolution of stresspaths toward the failure state due to rainfall in hillslope environments.
Besides this section for Introduction, this manual includes four other sections: Section 2
presents the fundamental hypotheses and flow and stress governing equation of Slope Cube
module, Sections 3 and 4 illustrate the input data preparation and the modeling visualization,
Section 5 presents 5 model applications using Slope Cube (2D), and Section 6 presents threemodel applications using Slope Cube (3D). An Appendix providing details of the finite element
discretization for the 3D analysis is also included.
12. Governing Equations of Slope Cube
The Slope Cube module, together with HYDRUS, simulates hillslope hydrology and stability
in order to predict the timing and location of rainfall-induced landslides. Gravitational-induced
stresses, rainfall-induced soil moisture and soil suction variations are concurrently simulated
by a coupled variably-saturated flow and stress field finite element method. The componentsof the hydro-mechanical framework are fully defined and validated in the textbook Hillslope
Hydrology and Stability (Lu and Godt, 2013), and are illustrated in Figure 1. The spatial andtemporal distributions of suction-stress based effective stress and the factor of safety at eachpoint throughout a hillslope are computed by taking into account of the hillslope’s
geomorphology, hydrology, and stress, strain and deformation boundary and initial conditions.

Figure 1. Illustration of the physical processes in Slope Cube module.
The hydrological and mechanical responses of hillslope under rainfall conditions are simulated
by combining Richards equation with a set of equations of a classical linear elasticity thatincorporate effective stresses (Lu et al., 2012; Lu and Godt, 2013). Knowledge of the soil water
retention curve (SWRC), hydraulic conductivity function (HCF), and mechanical properties ofhillslope materials is needed and is defined as the physical properties in the afore-mentioned
governing equations. The SWRC describes the relation between soil (volumetric) moisture
content (θ) and pressure head (h [m]). The HCF gives the relation between unsaturated
hydraulic conductivity (k [L/T]) and moisture content (θ) or pressure head (h). Lastly, acomponent including suction stress and local field safety factor is integrated into a hillslope
stability model (Lu et al., 2010) in order to analyze and predict effective stress, deformations,and factor of safety throughout the entire hillslope.
The module is able to account for different soil layers with their own mechanical andhydrological properties such as SWRC parameters, mechanical properties of internal friction
angle, cohesion, elastic modulus, Poisson’s ratio, and specific gravity.

22.1. Hydrological Model
The numerical solution of transient unsaturated flow is described in the HYDRUS (2D/3D)
model [Šimunek et al., 2008]. It solves the 2D and/or 3D Richards equation (Equation 1):∂θ ( h )= ∇ ⋅ k ( h ) ∇H + Ω(1)∂twhere θ [-] is the volumetric water content, t [T] is the time, h [L] is the water pressure or
suction head, hereafter the suction head can be expressed by soil suction, h = −ψ ( ρ w g ) as
the air pressure is neglected. H is the total head equal to the suction head h plus the elevation
head z [L], and Ω is the flux from a source or to a sink. k ( h ) is the hydraulic conductivity
function (HCF); and θ (h) is the SWRC. The vertical dimension is z, while x and y encompass
the horizontal plane. In 3D, Eq. (1) can be written as follows:
∂θ
∂
∂h
=
[k ( KijA+ KizA )] − Ω ( i, j = x, y, z )∂t ∂xi∂x j(2)

1 0 0 
1 0  , ( i, j = x, y, z ).
where0 0 1 
The HYDRUS model implements an array of SWRC and HCF models [Šimunek et al., 2008].
The van Genuchten model [van Genuchten, 1980] is reported as an example of SWRC in
Equation (3) and the Mualem model [Mualem, 1976] as an example of HCF in Equation (4):θs − θ rθ ( h=(3)) θr +
1−1/ n1 + (α h )n 
KijA = 0l 

1/ m m  θ − θr     θ − θr   k (θ ) k sat = 1 − 1 −    θ s − θ r     θ s − θ r   



2

(4)

where θs is the saturated volumetric water content [-]; θ r is the residual volumetric watercontent; α [1/L] is related to the inverse of the air entry value suction, and n [-] is a pore-size
distribution index. k sat is the saturated hydraulic conductivity coefficient, m and l model
parameters, m = 1 − 1 n and l a pore-connectivity parameter commonly assumed to be 0.5 byMualem (1976).
SWRC and HFC show different patterns depending on the soil types. Figures 2 and 3 show
SWRC and HFC for clay, silt or sand based on Equations (3) and (4).

3

Figure 2. Illustration of soil water retention curves for typical sand, silt, and clay (after Lu etal., 2014).

Figure 3. Illustration of hydraulic conductivity functions for typical sand, silt, and clay (after
Lu et al. 2014).
All processes related to variably-saturated water flow are described in detail in the HYDRUS(2D/3D) documentation [Šimunek et al., 2008], and will not be repeated here.
2.2. Suction Stress-based Effective Stress Model
The state of effective stress at each point in variably-saturated hillslopes is calculated by using
the unified effective stress theory proposed and validated by Lu and Likos (2004; 2006) andLu et al. (2010). It has the identical form to Terzaghi’s effective stress for saturated soil:
σ ′= σ − σ s(5)4where σ ′ is effective stress, σ is the total stress, and σ s is the suction stress.
Suction stress is a characteristic function of hillslope materials, and is called suction stresscharacteristic curve (SSCC). The SSCC of a material can be expressed in a sole function of
matric suction ( ua - uw ) through a common set of hydrologic parameters defined in the van
Genuchten’s (1980) soil water retention model (VG model):(ua − uw )(6a)
σs = −(1 + [α (ua − uw )]n )( n −1)/ n
or a sole function of the effective degree of saturation Se :Sσs =− e ( Se1/(1− n ) − 1)1/ n
α(6b)

The SSCC exhibits different patterns and magnitudes for different types of soil. Figure 4 shows
the patterns of the SSCC for three representative soils, namely, sand, silt, and clay. The
hydrologic parameters correspond to these three soils are also shown in Figure 4. The SSCC insandy and silty soils shows a non-monotonic variation pattern with saturation and reaches is its
maximum value of zero at both residual and full saturation and minimum value at somewhere
between. In contrast, the SSCC in clay soil shows a monotonic behavior and reaches its
minimum value at the residual saturation. For these three representative soils, the sand has its
minimum suction stress -1.5 kPa at the effective degree of saturation Se = 0.70, the silt has its
minimum suction stress -10 kPa at the effective degree of saturation Se = 0.50, and the clay
has its minimum suction stress -100 kPa at the effective degree of saturation Se = 0.0.

Figure 4. Illustration of suction stress characteristic curves for typical sand, silt, and clay
(after Lu et al., 2014).
2.3. Mechanical Model

5

The mechanical model computes the soil stresses and displacements by using the field of moist
unit weight as the gravity term and the suction stress as a component of the effective stress.
For the two-dimensional analysis, the model adopts the finite element code FEM2D (Reddy,
1985) for plane stress linear elasticity and extends the capabilities to simulate gravity, its timedependency due to transient unit weight, and implementation of suction stress in the effective
stress and displacement calculations. The model was then further extended to perform a similar
analysis in a three-dimensional field in SLOPE CUBE(3D); a detailed description of thisextension is provided in the Appendix A.
Stress at a point in hillslopes is a tensor. In two-dimensional space, it has three independent
components and is often expressed as:𝜎𝜎𝑥𝑥𝑥𝑥𝜎𝜎𝑖𝑖𝑖𝑖 = � 𝜎𝜎𝑧𝑧𝑧𝑧

𝜎𝜎𝑥𝑥𝑥𝑥𝜎𝜎𝑧𝑧𝑧𝑧 �(7a)where subscripts i and j run from x to z with x in the horizontal direction and z in the verticalor gravitational direction, σxx is the normal total stress in the horizontal direction, σzz is thenormal total stress in the vertical direction, and σxz is the shear stress in the x-z plane.
In a three-dimensional space, the tensor has 6 independent components and is expressed as:

σ xx σ xy σ xz 
σ ij = σ yx σ yy σ yz σ zx σ zy σ zz 

(7b)where subscripts i and j run from x to z with x, y in the horizontal direction and z in the verticalor gravitational direction, σ xx , σ yy are the normal total stresses in the horizontal plane, σ zz
is the normal total stress in the vertical direction, σ xy , σ yx is the shear stress in the x-y plane,

σ xz , σ zx is the shear stress in the x-z plane, σ yz , σ zy is the shear stress in the y-z plane.
Equivalent tensors for effective stresses in 2D and 3D respectively, are:𝜎𝜎 ′𝜎𝜎𝑖𝑖𝑖𝑖′ = � 𝑥𝑥𝑥𝑥′
𝜎𝜎𝑧𝑧𝑧𝑧′𝜎𝜎𝑥𝑥𝑥𝑥′ �
𝜎𝜎𝑧𝑧𝑧𝑧

σ xx′ σ xy′ σ xz′ σ ij′ = σ ′yx σ ′yy σ ′yz  ′′ σ zz′ σ zx σ zy

(8a)

(8b)

Strain at a point in hillslopes is also a tensor. In two-dimensional space it has 3 independent
components:

ε xz 
εε ij =  xxε zx ε zz 6(9a)

where εxx is the normal strain in the horizontal direction, εzz is the normal strain in the verticaldirection, and εxz is the shear strain in the x-z plane.
In a three-dimensional space, the tensor has 6 independent components and is expressed as:ε xx ε xy ε xz (9b)ε ij = ε yx ε yy ε yz ε zx ε zy ε zz 
where ε xx , ε yy are the normal strain in the horizontal plane, ε zz is the normal strain in thevertical direction, ε xy , ε yx are the shear stain in the x-y plane, ε xz , ε zx are the shear strain in
the x-z plane, and ε yz , ε zy are the shear strain in the y-z plane.
In linear elasticity theory, strain can be defined by the first derivatives of displacements asfollows:
1 𝜕𝜕𝑢𝑢In two-dimensions:

𝜕𝜕𝑢𝑢

𝜀𝜀𝑖𝑖𝑖𝑖 = � 𝜕𝜕 𝑖𝑖 + 𝜕𝜕 𝑗𝑗 �
2

𝑗𝑗


∂u x

u
∂
1  ∂u∂xj  = 
ε ij =  i +2  ∂j∂i   1  ∂u x ∂u z  2  ∂z + ∂x 


𝑖𝑖

1  ∂u x ∂u z +
2  ∂z∂x 
∂u z


∂z


(10)

(11a)

In three-dimensional space:∂u x
1 ∂u x ∂u y
1 ∂u x ∂u z 
+
+
(
)
(
)

∂x
∂x
∂x 
2 ∂y
2 ∂z

 1 ∂u ∂u∂u y
1 ∂u ∂u
1 ∂u z ∂u y +()ε ij = ( i + j ) =  ( x + y )
∂i
∂x
∂y
∂z 
2 ∂j
2 ∂y
 2 ∂y 1 ∂u ∂u∂u z
1 ∂u z ∂u y
 ( x + z)

+
(
)
∂x
∂z
∂z
2 ∂y
 2 ∂zFor example, the Cartesian coordinates for the three-dimensional case are:∂u
ε xx = x∂x∂u yε yy =∂y∂u
ε zz = z
∂z

7

(11b)

(11c)
(11d)
(11e)

1 ∂u x ∂u y
(
)
+
(11f)
2 ∂y
∂x
1 ∂u y ∂u z
ε=
ε=
(
+
)
(11g)
yz
zy
2 ∂z
∂y
1 ∂u x ∂u z
(11h)
+
(
)
ε=
ε=
xz
zx
∂x
2 ∂z
where u x , u y are the displacements in the horizontal plane, u z is the displacement in thevertical direction.

ε=
ε=
xy
yx

In linear elasticity, stress and strain can be related by Hooke’s law through the introduction ofmaterial’s elastic properties:(12)
For example, the Cartesian coordinates in the three-dimensional case are:′𝜎𝜎𝑥𝑥𝑥𝑥=′𝜎𝜎𝑦𝑦𝑦𝑦=′𝜎𝜎𝑧𝑧𝑧𝑧
=

2𝐺𝐺𝐺𝐺

1−2𝜇𝜇
2𝐺𝐺𝐺𝐺

1−2𝜇𝜇
2𝐺𝐺𝐺𝐺

1−2𝜇𝜇

(𝜀𝜀𝑥𝑥𝑥𝑥 + 𝜀𝜀𝑦𝑦𝑦𝑦 + 𝜀𝜀𝑧𝑧𝑧𝑧 ) + 2𝐺𝐺𝐺𝐺𝑥𝑥𝑥𝑥(𝜀𝜀𝑥𝑥𝑥𝑥 + 𝜀𝜀𝑦𝑦𝑦𝑦 + 𝜀𝜀𝑧𝑧𝑧𝑧 ) + 2𝐺𝐺𝐺𝐺𝑦𝑦𝑦𝑦(𝜀𝜀𝑥𝑥𝑥𝑥 + 𝜀𝜀𝑦𝑦𝑦𝑦 + 𝜀𝜀𝑧𝑧𝑧𝑧 ) + 2𝐺𝐺𝐺𝐺𝑧𝑧𝑧𝑧

(13a)
(13b)
(13c)

′ σ=
′yx 2Gε xy
σ=
xy

(13d)

′yz σ=′ 2Gε yz
σ=
zy

(13e)

′ σ=
′ 2Gε zx
σ=
zx
xz

(13f)

where δ ij is the identity tensor, µ is the Poisson’s ratio reflecting the change in the normal
strain in a specific direction due to the change in a unit strain orthogonal to it, and G is theshear modulus reflecting the change in shear stress due to the change in a unit shear strain.
Under the static force equilibrium conditions, i.e., no inertial force involved, at any point within
hillslopes, the momentum balance in both the horizontal and vertical directions should besatisfied, leading to the following force equilibrium equations in index notation:
To solve the elastic problem in three-dimensional conditions, the equilibrium equations need
to be constructed. Under the equilibrium conditions, at any point within hillslopes, themomentum balance in both the horizontal and vertical directions should be satisfied, leading
to the following force equilibrium equations in index notation:′
𝜕𝜕𝜎𝜎𝑖𝑖𝑖𝑖

𝜕𝜕𝜕𝜕+ 𝑓𝑓𝑖𝑖 = 0
8(14)where f is the body force vector that only has the non-zero component ρg in the z direction,
ρ is the density of the variably-saturated material and g is the gravitational constant. The
quantity ρg is also called unit weight of hillslope material.
For example, the Cartesian coordinates in the three-dimensional case are:∂σ x′ ∂σ xy ∂σ zx+++ fx =
0
∂x
∂y
∂z
∂σ xy

∂σ ′y

(15b)

∂σ zx ∂σ yz ∂σ z′+++ fz =
0
∂x
∂y
∂z

(15c)

∂y

+

∂σ yz

+ fy =
0

∂x

+

(15a)

∂z

Implementing the unified effective stress Eq. (5) into Eq. (14) yields:
𝜕𝜕𝜕𝜕𝑖𝑖𝑖𝑖

𝜕𝜕𝜎𝜎 𝑠𝑠= −𝑓𝑓𝑖𝑖 + 𝜕𝜕𝜕𝜕 = 0
𝜕𝜕𝜕𝜕

(16)

For example, the Cartesian coordinates for the three-dimensional case are:∂σ x ∂σ xy ∂σ zx
∂σ s
(17a)
+
+
=
− fx +∂x∂y∂z∂x∂σ xy ∂σ y ∂σ yz
∂σ s
(17b)
+
+
=
− fy +∂x∂y∂z∂y∂σ zx ∂σ yz ∂σ z
∂σ s
(17c)
+
+
=
− fz +∂x∂y∂z∂zEquations 10, 12, and 16 completely define the variables in a point in a hillslope. In the twodimensional case this involves 8 equations with 8 variables (3 total stresses, 3 strains, and 2
displacements). Conversely, in the three-dimensional case there are 15 equations with 15
variables (6 total stresses, 6 strains, and 3 displacements). The Cartesian equations for this caseare provided in equations (11c-11h), (13a-13f), and (17a-17c). The numerical solution of thesevariables is implemented in the coupled hydro-mechanical FE code with proper boundary
conditions of displacements and stresses.
2.4. Equation for deformation modulus
The Young’s modulus of a soil is a function of its volumetric water content or its matric suction(Lu and Kaya 2013). Thus, it is important to consider it when quantifying the deformation of
unsaturated soils under drying and wetting conditions. Based on experimental data, Lu andKaya (2014) proposed an empirical power law equation for Young’s moduli as a function of
volumetric water content:θ − θd m
E (θ ) =Ed + ( Ew − Ed )()(18)
θ w − θd9

where E is the Young’s modulus; for the drying state E = Ed when θ = θ d usually measured
at the residual saturation state; for the wetting state, E = Ew , when θ = θ w , usually measured
at the fully saturated state; m is an empirical parameter. The evolution of Young’s modulus ofa soil predicted by Eq. (1) can be seen in Figure 5.

Figure 5 Comparison of measured data and predicted curve using Eq. (18) for Young’s
modulus of a clayey soil under drying conditions.
Eq. (18), can be used when describing the stress-strain relationship of an element.
Combining the hydrological model, the mechanical model that implements variable Elastic
modulus and suction stress-based effective stresses, the deformation of a soils can be
determined for drying and wetting conditions.
2.5. Local Factor of Safety Model
Slope stability analysis has been traditionally performed with limit equilibrium methods (e.g.,
Fellenius, 1936; Bishop, 1955; Lam and Fredlund, 1993). The limit equilibrium methods
greatly simplify slope stability analysis with two propounding assumptions: presuming failuresurface geometry, and presuming failure occurrence simultaneously at every point along the
failure surface. While these two assumptions allow straightforward definition and computationof the factor of safety for a given slope; driving force (or moment) vs. resistance force (ormoment) along the failure surface, the limit equilibrium methodologies could not identify thelocation of failure initiation and predict real failure surface.
With the advent of modern computational methods for the numerical solutions of stress fieldby methods such as finite elements and finite differences, numerical algorithms for shear
strength reduction analysis in the context of finite element methods (e.g. Matsui and San, 1992;
Smith and Griffiths, 2004) have been developed to better calculate the driving forces andresistance forces along the failure surface. Because the afore-mentioned two assumptions arestill embedded in the strength reduction methodology, the use of a single FS as slope stability
indicator for the whole hillslope does not allow the description of the spatial and temporaldevelopment and evolution the failure surface.

10

Another major limitation to the traditional slope stability analysis methodologies is the use of
Terzaghi’s effective stress equation, which could not account for the development andevolution of effective stress under unsaturated conditions. This problem is overcome by theuse of the suction stress-based effective stress described in Section 2.2.
The concept of local factor of safety (LFS) presented in Lu et al., (2012) avoids theaforementioned two assumptions. Instead of presumed failure surface, the state of stress at eachpoint within a hillslope is assessed by the local factor of safety. Assuming the validity of MohrCoulomb failure criterion, the LFS is defined as the ratio of the potential Coulomb stress to thecurrent Coulomb stress, as illustrated in Figure 6:
τ*
(19)
LFS =
τ
where τ is the shear stress or the current Coulomb stress, and τ * is the shear strength or the
potential Coulomb stress.

Figure 6. Conceptual illustration of the local factor of safety (after Lu et al. 2012).
Assuming that the current state of stress at a point is represented by the solid Mohr circle(Figure 6), the potential Coulomb stress τ * can be uniquely identified as the intercept of theextension of the Coulomb stress with the Mohr circle envelope, shown as point B in Figure 6.
By the similarity of triangles ACD and ABE, it can be shown that:
AE BE BG τ *(20)=== = LFSAD CD CH τ
σ ′ + σ 3′c′(21)
BE AB=
sin φ ′ () sin φ ′=+ 1
tan φ ′
2
which leads to:cos φ ′
[2c′ + (σ 1′ + σ 3′ ) tan φ ′]σ 1′ − σ 3′
Inserting the unified effective stress Eq. (5) into Eq. (21) yields:
LFS F ( x, =y, z , t )
=

11

(22)

cos φ ′(23)[2c′ + (σ 1 + σ 3 − 2σ s ) tan φ ′]σ1 − σ 3
Note that for both 2D and 3D analyses, σ1 ’ is the major principal effective stress and σ3 ’ is theminor principal effective stress.
The stress paths under infiltration conditions in hillslopes follow greatly the leftward shifting
pattern. This is because the size of the Mohr circle, which is determined by the difference inmajor and minor principal total stresses, varies much less than the suction stress, which controls
the center of the Mohr circle during infiltration and seepage processes.
=
LFS F ( x,=y, z , t )Based on the matric suction and the soil moisture fields obtained by solving the Richards
equation, and the total stress field obtained by solving the elasticity problem, the effective stresscan be determined. Hence, the LFS field for a given hillslope can be obtained quantitatively.
Note that since the major and minor principal stresses are used, the LFS obtained in the threedimensional analysis is the smallest number regardless of the plane analyzed.

12

3. Preprocessing and Input Data
3.1 Slope Cube Module Activation
The Slope Cube module is activated when the “Slope Cube” check box is selected in theMain Processes and add-on Modules window (Figure 7).

Figure 7. The Main Processes and add-on Modules dialog window: check the box “Slope
Cube” for activation of Slope Cube module.
3.2 Definition of Materials Properties
The soil mechanical parameters for the Slope Cube module are specified in the Slope Cube
Parameters window (Figure 8). The definitions of these parameters are given in Table 1. All
parameters could be different for different soil layers.

Figure 8. Slope Cube module: defining materials properties.
13

Table 1. List of mechanical parameters for Slope Cube module.
Soil Parameter

Notation

Units

Specific gravity

Gs

dimensionless

Effective cohesion

c′

[kPa]

Effective internal friction angle

φ′[degree]Poisson ratiov

dimensionless

Ez/Ex ratio-

dimensionless

Young modulus at drying state

Ed

[kPa]

Young modulus at wetting state

Ew

[kPa]

Water content at drying state

cm3/cm3

Water content at wetting stateθdθwPower law parametermdimensionless

cm3/cm3

3.3 Definition of Boundary Conditions
The Slope Cube boundary conditions are specified by selecting “Slope Cube” in the “Boundary
Conditions” panel of the “Project Navigator-Data” window and defining a new boundary
condition (Figure 9).
In order to assign a new boundary condition, the user selects “New Boundary Condition” inthe “Edit Bar - Boundary Conditions” located on the right-hand side of the graphical userinterface (Figure 9). Different types of boundary conditions will appear in the “Edit Bar Boundary Conditions” window (Figure 10a, 10b). The detailed description of boundary
conditions for analysis in 2D and 3D is given in Tables 2 and 3, respectively.

Figure 9. Slope Cube module: defining boundary conditions from the “Boundary Condition”
panel in the “Project Navigator – Data” window.
14

(a)

(b)
Figure 10. Slope Cube module, boundary condition types for: (a) Slope Cube (2D), and (b)Slope Cube (3D).
Table 2. List of types of mechanical boundary conditions for Slope Cube (2D).
Parameter
Free displacement
X-displacement
Z-displacement
XZ-displacement
X-forceZ-force
XZ-force

Description
Free displacement condition
A displacement in the x-direction is specified
A displacement in the z-direction is specified
A displacement in the x- and z-directions is specified
A boundary force in the x-direction is specified
A boundary force in the z-direction is specified
Boundary forces in the x- and z-directions are specified

15

Table 3. List of types of mechanical boundary conditions for Slope Cube (3D).
Parameter
Free displacement
X-displacement
Y-displacement
Z-displacement
XY-displacement
YZ-displacement
XZ-displacement
XYZ-displacement
X-forceY-forceZ-force
XY-force
YZ-force
XZ-force
XYZ-force

Description
Free displacement condition
A displacement in the x-direction is specified
A displacement in the y-direction is specified
A displacement in the z-direction is specified
A displacement in the xy-direction is specified
A displacement in the yz-direction is specified
A displacement in the xz-direction is specified
A displacement in the xyz-direction is specified
A boundary force in the x-direction is specified
A boundary force in the y-direction is specified
A boundary force in the z-direction is specified
A boundary force in the xy-direction is specified
A boundary force in the yz-direction is specified
A boundary force in the xz-direction is specified
A boundary force in the xyz-direction is specified

16

4 Post-processing
4.1 The Slope Cube Module Outputs
The output files from the Slope Cube are described in detail in Tables 4 and 5, and shown in
Figure 11. Besides the hydrological variables computed by HYDRUS (3D), the Slope Cube
provides the fields of the minimum and maximum total and effective stresses, suction stress,
total deformations, and local factor of safety. The model outputs are similar to the HYDRUS
standard variables such as pressure head, water content. These variables can be visualized in
the graphical user interface, edited (Figure 12a) or queried using the cursor (Figure 12b). If
observation points were added in the model, graphical output of water content, pressure head,suction stress, maximum and minimum total and effective stresses, deformations, and local
factor of safety as a function of time are provided for the observed point (Figure 13).

(a)

(b)

Figure 11. Selection of output data for the (a) Slope Cube (2D), and (b) Slope Cube (3D).

17

Table 4. List of output data files for Slope Cube (2D).
File name
S3_SIG1.out
S3_SIG3.out
S3_SUCTSTRESS.out
S3_ESIG1.out
S3_ESIG3.out
S3_LFSF.out
S3_SXZ
S3_DXZ.out

Description
Nodal values of the maximum total stress in the domain
Nodal values of the minimum total stress in the domain
Nodal values of the suction stress in the domain
Nodal values of the maximum effective stress in the domain
Nodal values of the minimum effective stress in the domain
Nodal values of the local factor of safety in the domain
Nodal values of the shear stress in the XZ-plane in the domain
Nodal values of the displacement in the domain

Table 5. List of output data files for Slope Cube (3D).
Parameter
S3_SIG1.out
S3_SIG2.out
S3_SIG3.out
S3_ESIG1.out
S3_ESIG2.out
S3_ESIG3.out
S3_SUCTSTRESS.out
S3_LFSF.out
S3_SXY
S3_SYZ
S3_SZX
S3_DXYZ.out
S3_VOL.out

Description
Nodal values of the maximum total stress in the domain
Nodal values of the intermediate total stress in the domain
Nodal values of the minimum total stress in the domain
Nodal values of the maximum effective stress in the domain
Nodal values of the intermediate effective stress in the domain
Nodal values of the minimum effective stress in the domain
Nodal values of suction stress in the domain
Nodal values of the local factor of safety in the domain
Nodal values of the shear stress in the XY-plane in the domain
Nodal values of the shear stress in the YZ-plane in the domain
Nodal values of the shear stress in the XZ-plane in the domain
Nodal values of the displacement in the domain
Volume changes of the geometry

18

(a)

(b)
Figure 12. (a) Customization of the color scales of Slope Cube module, (b) Selection of contour
values in the display of Slope Cube module output by clicking on the map with the cursor.

19

(a)

(b)
Figure 13. Example of graphical display of the suction stress (a), and local factor of safety (b)for the observation nodes.

205. Example Problems for Slope Cube (2D)The Slope Cube module supplements the standard HYDRUS program, which includes thevariably-saturated water flow model that has been widely tested earlier. Therefore, the example
problems presented in this section concentrate on the demonstration of the new the Slope Cube
module.
The module and most of its features have been tested on five problems with increasing
complexity (Table 6). All the problems use the identical 23 degree-slope hillslope, which issimilar to the previous studied slope by others (e.g., Bogaard, 2002; Shao et al., 2015). The
slope stability analyses of such slope are conducted under different rainfall or load conditions.
The first example, described in Section 5.1 (Example 1; see Table 6) demonstrates the stability
analysis of a hillslope with homogeneous soil subjected to a constant infiltration rate. Two
combinations of rainfall intensity and duration were simulated: a heavy rainfall with anintensity of 0.02 m/h for a short duration of 15 h (S1), and a low rainfall with an intensity of
0.02 mm/h for a long duration of 150 h (S2). The scenario S1 is reported as Example 1a andscenario S2 as Example 1b.
The second example presents the stability analysis of a two-layer hillslope. Each layer has its
own hydrological and mechanical parameter set. Example 2a considers the rainfall scenario S1and Example 2b the rainfall scenario S2.
Finally, in Example 3, a stability analysis was carried out by applying a time-varying heavy
rainfall with a concentrated mechanical load at the top of the hillslope.
Table 6. List of examples for illustration of use of Slope Cube (2D) module.
Example nameExample 1a
Example 1b
Example 2a
Example 2b
Example 3

Description
Stability analysis of a hillslope with homogeneous sandy soil
subjected to a constant infiltration rate of 0.02 m/h for 15 hours.
Stability analysis of a hillslope with homogeneous sandy soil
subjected to a constant infiltration rate of 0.002 m/h for 150 hours.
Stability analysis of a layered hillslope subjected to a constantinfiltration rate of 0.02 m/h for 15 hours
Stability analysis of a layered hillslope subjected to a constantinfiltration rate of 0.002 m/h for 150 hours.
Stability analysis of a layered hillslope subjected to a time-varying
heavy rainfall event with a concentrated mechanical load

21

5.1 Example 1: Stability Analysis of a Hillslope with Homogeneous Sandy Soil under aConstant Infiltration Rate
In this example, stability analysis was carried out for a homogeneous sandy hillslope subjected
to two different rainfall scenarios S1 and S2 presented as Example 1a and Example 1b,
respectively. The variably-saturated water flow is described by Equation (1) with the van
Genuchten model for the SWRC, and the Mualem model as the HCF. Time-varying fields oftotal stress, effective stress, and local factor of safety are described respectively by Equations(10)-(15), and (19).
The hillslope geometry and boundary conditions for the Water Flow and the Slope Cubemodels are presented in Figure 14. We used hydrostatic conditions as the initial conditions withthe water table set at 2 m depth. The hillslope hydro-mechanical parameters are specified inTable 7. Model results are illustrated as the transient fields of pressure head, water content,
minimum and maximum effective stress, suction stress and local factor of safety for the
Example 1a (Figures 15-17), and for Example 1b (Figures 18-20). Finally, Figure 21 shows thecomparison between the failure surface evolutions for Example 1a and Example 1b. The stepby-step procedure to implement Example 1a and Example 1b is presented in the Slope Cube
Tutorial Section 1.

No flow
Z-displacement=0

Figure 14. Illustration of hillslope geometry and boundary conditions for hydrological (in black)and mechanical (in gray) models.

22

Table 7. List of hydro-mechanical parameters defined in Example 1.
Soil

Sand

Parameter

Value

[1/m]

4.0

n [-]

1.80

ks [m/h]
[-]
[-]

0.21
0.045
0.43

c [kPa][ ]140

Gs [-]

2.6

Ex=Ey [GPa]
[-]

10
0.33

Figure 15. Contour plots of pressure head [m] and water content evolutions with time for Example
1a. Results are shown at 0, 3, 6, and 15 hours.

23

Figure 16. Contour plots of maximum and minimum effective stress [kPa] evolutions with timefor Example 1a. Results are shown at 0, 3, 6, and 15 hours.

24Figure 17. Contour plots of suction stress [kPa] and field of local factor of safety evolutions withtime for Example 1a. Results are shown at 0, 3, 6, and 15 hours.

25

Figure 18. Contour plots of pressure head [m] and water content evolutions with time for
Example 1b. Results are shown at 0, 30, 60, and 150 hours.

26

Figure 19. Contour plots of maximum and minimum effective stress [kPa] evolutions with timefor Example 1b. Results are shown at 0, 30, 60, and 150 hours.

27

Figure 20. Contour plots of suction stress [kPa] and field of local factor of safety evolutions withtime for Example 1b. Results are shown at 0, 30, 60, and 150 hours.

28

Failure surface evolution for Example 1a

Failure surface evolution for Example 1b

Figure 21. Comparisons between the failure surface evolution for Example 1a and Example 1b;
Local factor of safety <1 is shown in red. Results are shown at 0, 3, 6, and 15 hours for Example
1a and at 0, 30, 60, and 150 hours for Example 1b.

5.2 Example 2: Stability Analysis of a Layered Hillslope under a Constant Infiltration Rate.
In this example, a stability analysis was carried out for a layered hillslope subjected to two
different rainfall scenarios S1 and S2 presented as Example 2a and Example 2b, respectively.
The hillslope consists of a 2 m-thickness top sandy layer and a deeper silty layer. The hillslope
geometry and boundary conditions for the water flow and the Slope Cube models are presented
in Figure 22. The hydrostatic conditions were used as the initial conditions with the water table
set at 2 m depth. The hillslope hydro-mechanical parameters are specified in Table 8 for bothsoil layers. Model results are illustrated as the transient fields of pressure head, water content,
minimum and maximum effective stress, suction stress and local safety factor for the Example
2a (Figures 23-25) and for Example 2b (Figures 26-28). Finally, Figure 29 shows thecomparison between the failure surface evolution for Example 2a and Example 2b. The step29

by-step procedure to implement Example 2a and Example 2b are presented in the Slope Cube
Tutorial Section 2.

sand
silt

No flow
Z-displacement=0

Figure 22. Illustration of hillslope geometry and boundary conditions for hydrological (in black)and mechanical (in gray) models for Example 2.
Table 8. List of hydro-mechanical parameters defined in Example 2a and Example 2b.
Soil Parameter

Value

Soil

Sand

Silt

[1/m]

4.0

1.6

n [-]

1.80

1.37

ks [m/h]
[-]
[-]

0.21

0.002

0.045

0.034

0.430

0.46

c [kPa]
[ ]

1.0

5.0

40

35

Gs [-]

2.6

2.7

Ex=Ey [GPa]
[-]

10

10

0.33

0.33

30

Figure 23. Contour plots of pressure head [m] and water content evolutions with time for Example
2a. Results are shown at 0, 3, 6, and 15 hours.

31

Figure 24. Contour plots of maximum and minimum effective stress [kPa] evolutions with timefor Example 2a. Results are shown at 0, 3, 6, and 15 hours.

32Figure 25. Contour plots of suction stress [kPa] and field of local factor of safety evolutions withtime for Example 2a. Results are shown at 0, 3, 6, and 15 hours.

33

Figure 26. Contour plots of pressure head [m] and water content evolutions with time for Example
2b. Results are shown at 0, 30, 60, and 150 hours.

34

Figure 27. Contour plots of maximum and minimum effective stress [kPa] evolutions with timefor Example 2b. Results are shown at 0, 30, 60, and 150 hours.

35

Figure 28. Contour plots of suction stress [kPa] and field of local factor of safety evolutions withtime for Example 2b. Results are shown at 0, 30, 60, and 150 hours.

36Failure surface evolution for Example 2a

Failure surface evolution for Example 2b

Figure 29. Comparisons between the failure surface evolutions for Example 2a and Example 2b;
Local factor of safety <1 is shown in red. Results are shown at 0, 3, 6, and 15 hours for Example
2a and at 0, 30, 60, and 150 hours for Example 2b.

5.3 Example 3: Stability Analysis of a Layered Hillslope under a Time-varying Heavy Rainfall
Event with a Concentrated Mechanical Load.
Example 3 demonstrates the stability analysis of a layered hillslope subject to the hourly
averaged rainfall measured during the Hurricane Frances that hit North Carolina in September
2004 (Wooten et al., 2008). Furthermore, the example includes a concentrated mechanical loadof 20 kN applied on the crest of the hillslope as shown in Figure 30. The hillslope geometry
and boundary conditions for the water flow and the Slope Cube models are presented in Figure
30. The hydrostatic conditions are used as the initial hydrologic conditions with the water table
set at 2 m depth. The hillslope hydro-mechanical parameters are specified in Table 8. The
rainfall time series is used as the top hydrologic boundary condition and is shown in Figure 31.

37

The results are presented in the transient fields of pressure head and water content (Figure 32),
minimum and maximum effective stress (Figure 33), and suction stress and field of local factorof safety (Figure 34). The stress concentration around the crest area can be seen in Figures 3334. Figure 35 shows the total displacement evolution with time. Finally, Figure 36 shows thefailure surface evolution with time. The mechanical effect to the slope stability can bewitnessed near the crest of the hillslope in Figure 36.

sand
silt

No flow
Z-displacement=0

Figure 30. Illustration of hillslope geometry and boundary conditions for hydrological (in black)and mechanical (in gray) models for Example 3. The top layer (in white) represents the sandy soil
and the lower layer in grey represents the silty soil. The concentrated mechanical load is in red.

Figure 31. Time-varying rainfall boundary condition used in Example 3.

38

Figure 32. Contour plots of pressure head [m] and water content evolutions with time for Example
3. Results are shown at 0, 34, 50, and 70 hours.

39Figure 33. Contour plots of maximum and minimum effective stress [kPa] evolutions with timefor Example 3. Results are shown at 0, 34, 50, and 70 hours.

40Figure 34. Contour plots of suction stress [kPa] and field of local factor of safety evolutions withtime for Example 3. Results are shown at 0, 34, 50, and 70 hours.

41Total displacement evolution [m] for Example 3

Time: 0 [h]Time: 34 [h]Time: 50 [h]Time: 70 [h]Figure 35. Contour plots of total displacement [m] evolutions with time for Example 3.
Results are shown at 0, 34, 50, and 70 hours.

42Failure surface evolution for Example 3

Figure 36. Contour plots of the local factor of safety < 1 for failure surface evolutions with timefor Example 3. Results are shown at 0, 34, 50, and 70 hours.

43

6. Example Problems for Slope Cube (3D)The Slope Cube module supplements the standard HYDRUS program, which includes thevariably-saturated water flow model that has been widely tested earlier. Therefore, the example
problems presented in this section concentrate on the demonstration of the new the Slope Cube
module.
The module and most of its features have been tested on three problems with increasing
complexity (Table 9). All the problems use the identical a 30-60-degree homogeneous
hillslope. The slope stability analyses of such slope are conducted under different rainfall or
load conditions.
The first example, described in Section 6.1 (Example 1; see Table 9) demonstrates the stability
analysis of a hillslope with homogeneous soil subjected to a constant infiltration rate. Two
combinations of rainfall intensity and duration were simulated: a heavy rainfall with anintensity of 0.05 m/h for a short duration of 20 h (S1), and a low rainfall with an intensity of
0.005 mm/h for a long duration of 200 h (S2). The scenario S1 is reported as Example 1a andscenario S2 as Example 1b.
The second example presents a stability analysis that was carried out by applying a timevarying heavy rainfall with a concentrated mechanical load at the top of the hillslope. Scenario
1 models a homogenous slope, while scenario 2 models a layered slope and the mesh has alarger number of elements.
Finally, the fourth example models a drying cake test, which consists on allowing a cylindricalsoil sample to deform due to changes in water content and therefore, changes in suction stress.
Table 9. List of examples for illustration of use of Slope Cube module.
Example nameExample 1a
Example 1b
Example 2a
Example 2b
Example 3

Description

Stability analysis of a hillslope with homogeneous sandy soil
subjected to a constant infiltration rate of 0.05 m/h for 20 hours.
Stability analysis of a hillslope with homogeneous sandy soil
subjected to a constant infiltration rate of 0.005 m/h for 200 hours.
Stability analysis of a homogeneous hillslope subjected to heavy
rainfall and a mechanical load.
Stability analysis of a homogeneous hillslope subjected to heavy
rainfall and a mechanical load.
Deformation of a Soil Cake due to Change in Water Content,
Comparison of Measured and Predicted Volumetric Change ofSoil Cake

44

6.1. Example 1: Stability Analysis of a Homogeneous Hillslope under a Constant Infiltration
Rate
In this example, stability analysis was carried out for a homogeneous sandy hillslope subjected
to a high infiltration rate presented as Example1. The variably-saturated water flow is describedby Equation (1) with the van Genuchten model for the SWRC, and the Mualem model as the
HCF. Time varying fields of effective stress, suction stress, and local factor of safety aredescribed respectively by Equations (10) - (11), (16) and (20).
The hillslope geometry and boundary conditions for the Water Flow and the Slope Cubemodels are presented in Figure 37. Initial conditions were hydrostatic with the water table set
at bottom. Two different rainfall scenarios are presented as Example 1a and Example 1b
respectively. The hillslope hydro-mechanical parameters are specified in Table 10. Model
results are illustrated as the transient fields of pressure head, water content, minimum andmaximum effective stress, suction stress, and local factor of safety for the Example 1a (Figures
38-41) and for Example 1b (Figures 42-45). Finally, Figure 23 shows the comparison between
the failure surface evolutions for Example 1a and Example 1b. The step-by-step procedure to
implement Example 1a and Example 1b is presented in the Slope Cube Tutorial Section 1.

Figure 37. Hillslope geometry and boundary conditions for hydrological (in black) andmechanical (in gray) models.

45

Table 10. List of hydro-mechanical parameters defined in Example 1.
Soil Parameter

Value

Soil type

Silt

[1/m]

3.6

n [-]

1.56

ks [m/h]
[-]

0.10
0.078

[-]

0.430

Gs [-]

2.6

c' [kPa]

12.0

φ' [ ]

45.0

µ [-]

0.30

Ez/Ex

1.0

Ed [kPa]

10000

Ew [kPa]

10000

θd

0.1

θw

0.35

m

0.25

46

Pressure head [m]Water content [-]Time: 0 [h]Time: 0 [h]Time: 4 [h]Time: 4 [h]Time: 12 [h]Time: 12 [h]Time: 20 [h]Time: 20 [h]Figure 38. Contour plots of pressure head [m] and water content evolutions with time for
Example 1a. Results are shown at 0, 4, 12, and 20 hours.
47

Maximum effective stress [kPa]

Minimum effective stress [kPa]

Time: 0 [h]Time: 0 [h]Time: 4 [h]Time: 4 [h]Time: 12 [h]Time: 12 [h]Time: 20 [h]Time: 20 [h]Figure 39. Contour plots of maximum and minimum effective stress [kPa] evolutions withtime for Example 1a. Results are shown at 0, 4, 12, and 20 hours.

48Suction stress [kPa]

Field of local factor of safety [-]Time: 0 [h]Time:Time: 00 [h][h]Time: 4 [h]Time: 4 [h]Time:Time:4 4[h][h]Time: 12 [h]Time: 12 [h]Time: 20 [h]Time: 20 [h]Figure 40. Contour plots of suction stress [kPa] and field of local factor of safety evolutions
with time for Example 1a. Results are shown at 0, 4, 12, and 20 hours.

49

(a)

(b)

(c)

(d)

Figure 41. Graphical display of observation nodes. (a) pressure head as a function of time, (b)water content as a function of time, (c) suction stress as a function of time, and (d) local factorof safety as a function of time for the three observation nodes.

50

Pressure head [m]Water content [-]Time: 0 [h]Time: 0 [h]Time: 40 [h]Time: 40 [h]Time: 120 [h]Time: 120 [h]Time: 200 [h]Time: 200 [h]Figure 42. Contour plots of pressure head [m] and water content evolutions with time for
Example 1b. Results are shown at 0, 40, 120, and 200 hours.

51Maximum effective stress [kPa]

Minimum effective stress [kPa]

Time: 0 [h]Time: 0 [h]Time: 40 [h]Time: 40 [h]Time: 120 [h]Time: 120 [h]Time: 200 [h]Time: 200 [h]Figure 43. Contour plots of maximum and minimum effective stress [kPa] evolutions withtime for Example 1b. Results are shown at 0, 40, 120, and 200 hours.

52Suction stress [kPa]

Field of local factor of safety [-]Time: 0 [h]Time: 0 [h]Time: 40 [h]Time: 40 [h]Time: 120 [h]Time: 120 [h]Time: 200 [h]Time: 200 [h]Figure 44. Contour plots of suction stress [kPa] and field of local factor of safety evolutions
with time for Example 1b. Results are shown at 0, 40, 120, and 200 hours.

53Field of local factor of safety [-]Field of local factor of safety [-]Time: 0 [h]Time: 0 [h]Time: 4 [h]Time: 40 [h]Time: 12 [h]Time: 120 [h]Time: 20 [h]Time: 200 [h]Figure 45. Comparisons between the failure surface evolution for Example 1a and Example
1b; Local factor of safety <1 is shown in red. Results are shown at 0, 4, 12, and 20 hours for
Example 1a and at 0, 40, 120, and 200 hours for Example 1b.
54

6.2. Example 2: Stability Analysis of a Homogeneous Hillslope under a Time-Varying
Heavy Rainfall Event with a Concentrated Mechanical Load
Example 2 demonstrates the stability analysis of a homogeneous hillslope subject to the hourly
averaged rainfall measured during the Hurricane Frances that hit North Carolina in September
2004 (Wooten et al., 2008). Furthermore, the example includes a concentrated mechanical loadof 20 KN applied on the crest of the hillslope as shown in Figure 46. The hillslope geometry
and boundary conditions for the water flow and the Slope Cube models are presented in Figure
46. The hydrostatic conditions are used as the initial hydrologic conditions with the water table
set at the bottom. The hillslope hydro-mechanical parameters are specified in Table 10. The
rainfall time series is used as the top hydrologic boundary condition and is shown in Figure 47.
The results are presented in the transient fields of pressure head and water content (Figure 48),
minimum and maximum effective stress (Figure 49), and suction stress and field of local factorof safety (Figure 50). The information of observed nodes can be seen in Figure 51.
20 KN

sand
silt

Figure 46. Illustration of hillslope geometry and boundary conditions for hydrological (inblack) and mechanical (in gray) models for Example 2.

Figure 47. Time-varying rainfall boundary condition used in Example 2.

55

Pressure head [m]Water content [-]Time: 0 [h]Time: 0 [h]Time: 34 [h]Time: 34 [h]Time: 50 [h]Time: 50 [h]Time: 70 [h]Time: 70 [h]Figure 48. Contour plots of pressure head [m] and water content evolutions with time for
Example 1a. Results are shown at 0, 34, 50, and 70 hours.

56Maximum effective stress [kPa]

Minimum effective stress [kPa]

Time: 0 [h]Time: 0 [h]Time: 34 [h]Time: 34 [h]Time: 50 [h]Time: 50 [h]Time: 70 [h]Time: 70 [h]Figure 49. Contour plots of maximum and minimum effective stress [kPa] evolutions withtime for Example 1a. Results are shown at 0, 34, 50, and 70 hours.

57Suction stress [kPa]

Field of local factor of safety [-]Time: 0 [h]Time: 0 [h]Time: 34 [h]Time: 34 [h]Time: 50 [h]Time: 50 [h]Time: 70 [h]Time: 70 [h]Figure 50. Contour plots of suction stress [kPa] and field of local factor of safety evolutions
with time for Example 1a. Results are shown at 0, 34, 50, and 70 hours.

58

(a)

(b)

(c)

(d)

Figure 51. Graphical display of observation nodes. (a) pressure head as a function of time, (b)water content as a function of time, (c) suction stress as a function of time, and (d) local
factor of safety as a function of time for the three observation nodes.
The second scenario of this example demonstrates the stability analysis of a layered hillslope:silt on top and clay underneath. The hillslope hydro-mechanical parameters are specified inTable 13.

59

20 KN
25.78

Back surface, Y=36 m:
No flow
Top surface:u y =0AtmosphericFree displacement

23.78

Left surface, X=0 m:
No flow
ux=0

55°

Silt

27.73

Clay
Front surface, Y=0 m:
No flow
uy=0

Z
Y
O

X

Unit: m

36

15
13

26°

0

Bottom surface, Z=0 m:
pw =0
uz=0

36.52

Figure 52. Illustration of hillslope geometry and boundary conditions for hydrological (inblack) and mechanical (in gray) models for Example 2b.
Table 11. List of hydro-mechanical parameters used in Example 2b
Soil Parameter

Value

Soil type

Silt

Clay

[1/m]

3.6

1.9

n [-]

1.56

1.31

ks [m/h]
[-]

0.10

0.01

0.078

0.095

[-]

0.430

0.410

Gs [-]

2.6

2.7

c' [kPa]

22.0

30.0

φ' [ ]

50.0

35.0

v

0.20

0.25

Ez/Ex

1.0

1.0

Ed [kPa]

15000

10000

Ew [kPa]

15000

10000

θd

0.10

0.10

θw

0.35

0.35

m

0.25

0.25

60

Maximum effective stress [kPa]XYZ displacement [m]Time: 0 [h]Time: 0 [h]Time: 32 [h]Time: 32 [h]Time: 56 [h]Time: 56 [h]Time: 72 [h]Time: 72 [h]Figure 53. Contour plots of maximum effective stress [kPa] and displacement [m] evolutions
with time for Example 4. Results are shown at 0, 32, 56, and 72 hours.

61Suction stress [kPa]

Field of local factor of safety [-]Time: 0 [h]Time: 0 [h]Time: 32 [h]Time: 32 [h]Time: 56 [h]Time: 56 [h]Time: 72 [h]Time: 72 [h]Figure 54. Contour plots of suction stress [kPa] and field of local factor of safety [-] evolutions
with time for Example 4. Results are shown at 0, 32, 56, and 72 hours.

62

The comparison of running time for the example is also done under different sizes of element
mesh, seen in Table 14.
Table 12. Comparison of running time for different sizes of element mesh.
Geometry

Max: 36.52×36×25.78 meters

Modeling time

80 hours

Print time count10
Element number

Running time (s)171,499
Size=0.96 m（Automatic）

347,912
Size=0.75 m

708,491
Size=0.60 m

218.7

496.5

1108.66.3. Example 3: Deformation of a Soil Cake due to Change in Water Content, Comparisonof Measured and Predicted Volumetric Change of Soil Cake
The third example simulates an experiment using the drying cake method (Lu and Kaya, 2013).
Results from the numerical simulation are compared to experimental results obtained by Dongand Lu in 2017. The set up consists of a soil cake made out of claystone of 2.00 cm in heightand 7.62 cm in diameter; the sample was compacted to a porosity of 0.47, initially fully
saturated, and subjected to drying conditions. Volumetric water content decreased from 0.47
to 0.21; changes in diameter and height were measured at ten different times and are presented
in Table 13. The geometry and boundary conditions are presented in Figure 55; hydrologicalconditions are labeled in black, while mechanical boundary conditions are labeled in gray. Soilparameters are provided in Table 14. Young modulus was obtained experimentally at each stepby applying a constant strain and measuring the stress in the soil sample. For the numericalsimulation, the Young modulus values at different water contents were obtained using thepower law equation developed by Lu and Kaya in 2013. A comparison between the measured
and predicted Young modulus values is presented in Figure 56.

63

Figure 55. Soil cake geometry and boundary conditions for hydrological (in black) andmechanical (in gray) models.
Table 13. List of hydro-mechanical parameters used in Example 3.
Soil Parameter

Value

Soil type
α [1/m]n [-]Denver claystone
0.29431.29

ks [m/h]

0.0003

θ r [-]
θs [-]c′ [kPa]φ′ [ ]

0.04
0.47
10.0
30.0

Gs [-]

2.7

Ed [KPa]2391

θd [-]

0.04

Ew [KPa]196

θ w [-]
m [-]
µ [-]

0.47
1.59
0.25

64

Table 14. Steps for the water content change used in Example 3.
Volumetric watercontent
0.470
0.457
0.442
0.421
0.401
0.366
0.344
0.314
0.282
0.247
0.206

Step
0
1
2
3
4
5
6
7
8
9
10

Figure 56. Comparison of the measured and fitted Young’s modulus of Denver claystone.
The measured and predicted suction stress characteristic curves (SSCC) are provided in Figure
57. Measured values were obtained using the Drying Cake method (Dong and Lu in 2017), andpredicted values were obtained with Slope Cube 3D.

65

104
Suction stress (-kPa)

Denver claystone

103

SSCC

102

101

100
0.0

0.1

0.2

0.3

0.4

0.5

Volumetric water content
Figure 57. Comparison of the measured and predicted SSCC of Denver claystone.

Denver claystone

0.5

Porosity

0.4

0.3

0.2

Measured
Predicted

0.1
0.0

0.1

0.2
0.3
0.4
Volumetric water content

0.5

Figure 58. Comparison between the simulated and measured porosity.
Finally, Figure 58 displays a comparison of porosity changes measured with a particle image
velocimetry (PIV) method and predicted using SlopeCube3D. Measured values were obtainedusing the Drying Cake method (Dong and Lu in 2017), and predicted values were obtainedwith Slope Cube 3D.

66

APPENDIX-A

67

APPENDIX-A: Finite Element Discretization for 3D Case
The finite element method was used to solve the elastic problem under three-dimensional (3D)
conditions. The variables on the continuous medium can be described by the variables on thenodes of elements. Hence, the finite element discretization technique needs to be used totransfer the variables to the nodes. Then all the variables can be solved using the 3D finite
element method (e.g., Zienkiewicz et al. 2005; Tirupathi and Ashok 2002). The detailed
discretization method and implemented process is described in this section.
The three invariants of stress tensor can be determined by the effective stress components inthree directions,′ + σ ′yy + σ zz′
I1 = σ xx(A-1-a)′ σ ′yy + σ ′yyσ zz′ + σ xx
′ σ zz′ − (σ ′yz ) 2 − (σ xz′ ) 2 − (σ xy′ )2I 2 = σ xx
′ σ ′yyσ zz′ + 2σ ′yzσ xz′ σ xy′ − σ xx′ (σ ′yz ) 2 − σ ′yy (σ xz′ ) 2 − σ zz
′ (σ xy′ )2I 3 = σ xx

(A-1-b)(A-1-c)In addition:
I12
− I2
3
I
II
b=
−2( 1 )3 + 1 2 − I 3
3
3
a
c=2
3
1
3b
=
θ
cos −1 (− )
3
ac
=
a

Thus, the three principal effective stresses are:I1
+ c cos θ
σ=
1′
3
I1
2π
+ c cos(θ +σ 2′ =
)
3
3
I1
4π
)
σ 3′ =+ c cos(θ +33(A-2-a)(A-2-b)(A-2-c)(A-2-d)(A-3-a)(A-3-b)(A-3-c)A.1. Displacement Function

For the tetrahedral element (Figure A-1), the displacement field on any node in the element canbe described by using a linear displacement function, given as follows:
u =α1 + α 2 x + α 3 y + α 4 z(A-4-a)v =α 5 + α 6 x + α 7 y + α8 z
w =α 9 + α10 x + α11 y + α12 z(A-4-b)(A-4-c)The displacement function should also be satisfied by the displacements in the four nodes of
the tetrahedral element.
For node p,up =α1 + α 2 x p + α 3 y p + α 4 z p(A-5-a)68

For node i,vp =α 5 + α 6 x p + α 7 y p + α8 z p(A-5-b)wp =α 9 + α10 x p + α11 y p + α12 z p(A-5-c)ui =α1 + α 2 xi + α 3 yi + α 4 zi

(A-6-a)vi =α 5 + α 6 xi + α 7 yi + α8 zi
wi =α 9 + α10 xi + α11 yi + α12 zi
For node j,

For node m,(A-6-b)(A-6-c)uj =α1 + α 2 x j + α 3 y j + α 4 z j(A-7-a)vj =α 5 + α 6 x j + α 7 y j + α8 z j(A-7-b)wj =α 9 + α10 x j + α11 y j + α12 z j(A-7-c)um =α1 + α 2 xm + α 3 ym + α 4 zm

(A-8-a)vm =α 5 + α 6 xm + α 7 ym + α8 zm
wm =α 9 + α10 xm + α11 ym + α12 zm(A-8-b)(A-8-c)Figure A-1. Tetrahedral element for three-dimensional problems
The parameters α i ( i = 1,12 ) can be solved by the twelve equations in Eqs. (A-5 to A-8) given
the node displacements and coordinates of the four nodes in three directions. After that, thedisplacement on any node within an element can be expressed:
1=[(ai + bi x + ci y + di z )uiu6V+ (a j + b j x + c j y + d j z )u j
(A-9-a)+ (am + bm x + cm y + d m z )um+ (a p + b p x + c p y + d p z )u p ]Similarly,

69

=
v

1
[(ai + bi x + ci y + di z )vi
6V+ ( a j + b j x + c j y + d j z )v j

(A-9-b)+(am + bm x + cm y + d m z )vm+ ( a p + b p x + c p y + d p z )v p ]w=1[(ai + bi x + ci y + di z ) wi
6V
+ (a j + b j x + c j y + d j z ) w j(A-9-c)+(am + bm x + cm y + d m z ) wm+(a p + bp x + c p y + d p z ) w p ]where V is the volume of a Tetrahedral element, 6V = det

determinant of the matrix.
The parameters in Eq. (A-9) can be determined as follows:
xj yj zjai = det xm ym zm

xp yp z p
1 yj zjbi = − det 1 ym zm

1 yp zp
xj 1 zjc = − det x 1 z
m
m
 i

xp 1 z p
xj yj 1di = − det xm ym 1
xp yp 1


70

1 xi
1 xj

yi
yj

zi
zj

1 xm
1 xp

ym
yp

zm
zp

. Here det is the(A-10-a)
xi yi zi
a j = − det xm ym zm

xp yp z p
1 yi zib j = det 1 ym zm

1 yp zp
xi 1 zic = − det x 1 zmm
 j
xp 1 z p
xi yi 1d j = − det xm ym 1
xp yp 1


(A-10-b)


xi yi ziam = − det x j y j z j
xp yp z p
1 yi zib = − det 1 y
zj
j
 m

1 yp zp
xi 1 zicm = det x j 1 z j
xp 1 z p
xi yi 1d m = − det x j y j 1

xp yp 1


(A-10-c)71
xi yi zia p = − det x j y j z j
xm ym zm

1 yi zib p = − det 1 y j z j
1 ym z m
xi 1 zic = − det x 1 z
j
j
 p

xm 1 zm
xi yi 1d p = det x j y j 1xm ym 1


(A-10-d)

If shape functions,1(a + b x + c y + d z ) ,N=
i 6V i iii
1(a j + b j x + c j y + d j z ) ,
N=
j
6V
1
(am + bm x + cm y + d m z ) ,
N=
m
6V
1
(a p + bp x + c p y + d p z )N=p6Vare assumed, equations A-9 can be expressed as:
u =Ni ui + N j u j + N mum + N p u p(A-11-a)v =Ni vi + N j v j + N m vm + N p v p(A-11-b)w =Ni wi + N j w j + N m wm + N p w p(A-11-c)Here Ni , N j , N m , and N p can be called shape function.

A.2. Strain of Element
The strain of the element can be determined by the displacement of the element as follows:72ε x  
ε y ε z ={ε } =γ xy γ  yz 
γ zx 

 ∂Ni

 ∂x

 0


 0

 ∂Ni

 ∂y

 0


 ∂N i
 ∂z

 ∂u

 ∂x



 ∂v

 ∂y



 ∂w

 ∂z

=


 ∂u + ∂v  ∂y ∂x  ∂v + ∂w  ∂z ∂y  ∂w + ∂u  ∂x ∂z 

∂
 ∂x

0


0


∂
 ∂y

0


∂
 ∂z

∂N j

∂
∂y
0
∂
∂x
∂
∂z
0

0

0

∂Ni
∂y

0

0

0

∂Ni
∂z

0

0

∂Ni
∂x

0

∂N j

∂N j

∂y

∂Ni
∂z

∂Ni
∂y

∂x
∂N j

0

∂N i
∂x

 bi
ε x 

 
0
ε

y
 
ε  1  0z ={ε } =

 ci
γ
6
V
 xy 

γ 
0
yz
 

γ zx 
 di

0

∂x

0

∂N j
∂y

∂z

0

0

0

0

0

∂N m
∂y

0

0

0

0

∂N m
∂z

0

0

∂N m
∂y

∂N m
∂x

0

∂N p

∂N p

∂y

0

∂N m
∂z

∂N m
∂y

∂x
∂N p

∂N m
∂z

0

∂N m
∂x

∂z
0
∂N j
∂y
∂N j

0

∂N p

∂N m
∂x

∂N j

∂z

∂N j


0

0


∂
u 
∂z    v  =0   w


∂
∂y 

∂
∂x 

0

∂x

0

∂x

∂N p
∂y

0

∂z

∂N p

0

∂z

0

0

bj

0

0

bm

0

0

bp

0

ci

0

0

cj

0

0

cm

0

0

cp

0

di

0

0

dj

0

0

dm

0

0

bi

0

cj

bj

0

cm

bm

0

cp

bp

di

ci

0

dj

cj

0

dm

cm

0

dp

0

bi

dj

0

bj

dm

0

bm

dp

0

= [ B ]{u}

73

ui 
 
0  vi  w 
 i 
0  u j  ∂N p  v j  ∂z   w j  u 
0  m  vm  
∂N p   wm ∂y  u p  ∂N p  v  p∂x   w 
 p

ui 
 
vi 
 wi  
0  u j 0  v j  
d p   w j  
0  um c p  vm  b p   wm 
u 
 p
v p  
 w p (A-12-a)(A-12-b)where [B] is the strain matrix
 bi 0 0 ci
1 0 0[ B] =  c b
6V ii 0 di di 0

0

bj

0

0

bm

0

0

bp

0

0

0

cj

0

0

cm

0

0

cp

di

0

0

dj

0

0

dm

0

0

0

cj

bj

0

cm

bm

0

cp

bp

ci

0

dj

cj

0

dm

cm

0

dp

bi

dj

0

bj

dm

0

bm

dp

0

v
v
0
1 − vv
0
σ x 
 v 1− v  v
v 1− v0
σ y σ 
1 − 2v
 0
E
 z
0
0
=
=
σ
{ } τ 2 xy  (1 + v)(1 − 2v) 
τ 
 0
0
0
0
 yz 
τ zx 

0
0
0
 0


0
0
0

0

0
dp 

0

cp 

b p 

(A-13)

A.3. Stress and Strain Relation0
1 − 2v20
 ε x  ε y 
 
ε 
z 
0  =
 γ  [ D]{ε }  xy 
0  γ yz  
1 − 2v  γ zx 2 (A-14)
000where [D] is the stiffness matrix,vv
0
1 − v
 v 1− vv
0 v
v 1− v0
1 − 2v
 0E
00[ D] =
2(1 + v)(1 − 2v) 
 0
0
0
0


0
0
0
 0


0
0
0
0
1 − 2v
2
0






0 


0 

1 − 2v 

2 
0
0
0

(A-15)

A.4. Element Stiffness Matrix[ K e ] = ∫ [ B e ]T [ D e ][ B e ]dxdydz(A-16-a)[ K e ] = [ B e ]T [ D e ][ B e ]V

(A-16-b)

74

 bi
0

0

bj
0

T
1 0
where  B e  = 
6V bm
0

0

bp
0

 0
 bi

0

1 0
 Be  =  6V  c
i

0

 di

cp

dm
cm
0

0

bp

dp

0

dp

0

cp

di 
0 
bi 

dj 
0

bj 
,
dm 
0

bm 

dp 
0

b p 

0

0
0

ci
bi
0

0

ci
0
0

di
0

cj

ci
0

cj

0

bj

dj

0

dj

0

cj

0

0
0

cm
bm
0

0

cm
0
0

dm
0

cp

di

0

0

bj

0

0

bm

0

0

bp

0

ci

0

0

cj

0

0

cm

0

0

cp

0

di

0

0

dj

0

0

dm

0

0

bi

0

cj

bj

0

cm

bm

0

cp

bp

di

ci

0

dj

cj

0

dm

cm

0

dp

0

bi

dj

0

bj

dm

0

bm

dp

0

v
v
0
1 − v
 v 1− vv
0 v
v 1− v0
1 − 2v
 0
E00
 De  =2  (1 + v)(1 − 2v) 

 0
0
0
0


0
0
0
 0


0
0

A.5. Body Force75001 − 2v
2
0

0 
0 
0 

0 


0 

1 − 2v 

2 

0

0
dp 
,
0

cp 

b p 

bxe0  e 
by 0  e 
bz 0  e 
bx 0 V 
[ f e ] = bye 0 6 ebz 0  e 
bx 0 be 
 y0 
be 
 z0 

(A-17)

A.6. Boundary Conditions
t x = σ x nx + τ xy n y + τ zx nz
t y = τ xy nx + σ y n y + τ yz nz
t z = τ zx nx + τ yz n y + σ z nz
u = u , v = v , and w = w

on Γ 2(A-18-a)on Γ1(A-18-b)

A.7. Partial Differential Equation Discretization
Eqs. (17) can be expressed by the matrix form,
σ
∂
∂
∂
∂  x  00 ∂x 0 ∂xσ∂y∂z   y 




f
x
σ 


  
∂
∂
∂
z 
− fy  +  0
00  =
0
σ
∂
∂
∂
y
x
z
xy
  

 
 fz  ∂∂∂  σ yz 
0
0
0
 
0
∂z
∂y ∂x  σ  zx 

Inserting Eq. (A-14) into Eq. (A-19) yields:

76

0

0

∂
∂y

0

0

∂
∂z

∂
∂y
∂
∂x
0

0
∂
∂z
∂
∂y

s
∂  σ  ∂z  σ s 

 s

0  σ  (A-19)
 0 ∂  0  ∂x  0 
 

∂

 ∂x

0


0


0

0

∂
∂y

0

0

∂
∂z

∂
∂y
∂
∂x
0

0
∂
∂z
∂
∂y

 ∂u

 ∂x



 ∂v



∂
∂
∂y



 ∂x
∂z 
 ∂w

 fx  
 ∂z


  


=
− fy  +  0
0  [ D]  ∂u + ∂v    fz   ∂y ∂x 
∂



0
∂x 

 ∂v + ∂w  ∂z ∂y  ∂w + ∂u  ∂x ∂z 



0

0

∂
∂y

0

0

∂
∂z

∂
∂y
∂
∂x
0

0
∂
∂z
∂
∂y

s
∂  σ ∂z  σ s 

 s

0  σ  0 ∂  0  ∂x  0  (A-20)
Expressed in matric form
σ s 
 s
σ 
 fx u   T  s
−  f y  + [ B ] σ 
[ B ]T [ D ][ B ] v  =
 w
 
0 
 
 fz 
0 
 
0 



(A-21)

where [ K e ] = [ B ] [ D ][ B ] is element stiffness matrix.
T

A.8. Finite Element Form
The global stiffness matrix is determined by assembling the element stiffness matrix of all theelements. In a similar manner, the global force vector is determined by assembling the element
force vector of all the elements. Thus, the global equations can be expressed as:(A-22)[ K ]{δ } = {F }

u 
 
Addressing the initial conditions and solving Eq. (A-22), the node displacement {δ } = v  w
 
can be obtained by the iterative numerical method.

77

References
Baum, R. L., Godt, J. W., and Savage, W. Z. (2010). Estimating the timing and location ofshallow rainfall induced landslides using a model for transient, unsaturated infiltration.
Journal of Geophysical Research: Earth Surface, 115, F03013, 1-26.
Bishop, A. W. (1955), The use of slip circle in the stability analysis of slopes. Geotechnique,5(1), 7–17.
Bogaard, T. (2002). Analysis of Hydrological Processes in Unstable Clayey Slopes. Ph.D.
dissertation, Universiteit Utrecht, Utrecht.
Chen, P., Mirus, B., Lu, N., and Godt. J. W. (2017). Effect of hydraulic hysteresis on thestability of infinite slopes under steady infiltration. J. Geotech. Geoenviron. Eng., 143(9),04017041, 1-10.
Chen, P., and Lu, N. (2018a). Generalized equation for soil shrinkage curve. J. Geotech.
Geoenviron. Eng., 144(8), 04018046, 1-10.
Chen, P., Lu, N., Giuseppe Formetta, Godt, J. W. and Wayllace, A. (2018b). Tropical storminduced landslide potential using combined field monitoring and numerical modeling. J.
Geotech. Geoenviron. Eng., 144(11): 05018002, 1-12.
Dong, Y., and Lu, N. (2017). Measurement of suction-stress characteristic curve under drying
and wetting conditions. Geotech. Test. J., 40(1), 1–15, doi:10.1520/GTJ20160058.
Fellenius, W. (1936). Calculation of the stability of earth dams, in Transactions of the 2nd
Congress on Large Dams, vol. 4, pp. 445–463, Int. Comm. on Large Dams, Washington,
D.C.
Janbu, N. (1973). Slope stability computations. Embankment-dam Engineering, in Casagrande,
vol. R, edited by C. Hirschfeld and S. J. Poulos, pp. 49–86, John Wiley, New York.
Lam, L., and Fredlund, D. G. (1993). A general limit equilibrium model for three-dimensionalslope stability analysis. Canadian Geotechnical Journal, 30(6), 905-919.
Lu, N., and Godt, J. W. (2013). Hillslope Hydrology and Stability. Cambridge University Press.
Lu, N., SenerıKaya, B., Wayllace, A., and Godt, J. W. (2012). Analysis of rainfall induced
slope instability using a field of local factor of safety. Water Resources Research, 48,
W09524, 1-14.
Lu, N., Wayllace, A., and Godt, J. W. (2010). A hydro-mechanical model for predicting
infiltration-induced landslides, The 3rd USGS Modeling Conference, June 7-11, 2010,
Geological Society of America No. 3-3.
Lu, N., and Likos, W. J. (2004). Unsaturated Soil Mechanics. Wiley.
Lu, N., and Likos, W. J. (2006). Suction stress characteristic curve for unsaturated soil. Journalof Geotechnical and Geoenvironmental Engineering, 132(2), 131-142.
Lu, N., and Kaya, M. (2013). “A drying cake method for measuring suction stress characteristic
curve, soil-water retention, and hydraulic conductivity function.” Geotech. Test. J., 36, 1–
19, doi:10.1520/GTJ20120097.
Lu, N., and Kaya, M. (2014). “Power law for elastic moduli of unsaturated soil.” J. Geotech.
Geoenviron. Eng., 140(1), 46–56.
Lu, N., Kaya, M., and Godt, J. W. (2014). Interrelations among the soil-water retention,
hydraulic conductivity, and suction-stress characteristic curves. Journal of Geotechnicaland Geoenvironmental Engineering, 140(5), 04014007, 1-10.
Matsui, T., and San, K. C. (1992). Finite element slope stability analysis by shear strength
reduction technique, Soils and Foundations, 32(1), 59–70.
Montgomery, D. R., Dietrich, W. E. (1994). A physically based model for the topographic
control on shallow landslide. Water Resources Research, 30, 83–92
Mualem, Y. (1976). A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resource Research, 12(3), 513-522.
78

Reddy, J. N. (1993). An Introduction to the Finite Element Method (Vol. 2, No. 2.2). New
York: McGraw-Hill.
Simoni, S., Zanotti, F., Bertoldi, G., and Rigon, R. (2008). Modeling the probability ofoccurrence of shallow landslides and channelized debris flows using GEOtop-FS,
Hydrological Processes, 22, 532-545,
Shao, W., Bogaard, T. A., Bakker, M., and Greco, R. (2015). Quantification of the influenceof preferential flow on slope stability using a numerical modelling approach. Hydrology and
Earth System Sciences, 19(5), 2197-2212.
Šimunek, J., M. Th. van Genuchten, and M. Šejna. (2008). Development and applications of
the HYDRUS and STANMOD software packages and related codes, Vadose Zone Journal,
7(2), 587-600.
Smith, I. M., and D. V. Griffiths (2004), Programming the Finite Element, 4th ed., John Wiley,West Sussex, U. K.
Tirupathi, R. C., and Ashok, D. B. (2002). Introduction to finite elements in engineering (3thedition). P473, Prentice Hall, USA.
van Genuchten, M. T. (1980). A closed-form equation for predicting the hydraulic conductivityof unsaturated soils. Soil Science Society of America Journal, 44(5), 892-898.
Witt, A. C. (2005). A brief history of debris flow occurrence in the French Broad River
Watershed, western North Carolina. NC Geogr, 13, 58-82.
Wooten, R. M., Gillon, K. A., Witt, A. C., Latham R. S., Douglas, T. J., Bauer, J. B.,
Fuemmeler, S. J., Lee, L. G. (2008). Geologic, geomorphic, meteorological aspects of debris
flows triggered by hurricanes Frances and Ivan during September 2004 in the southern
Appalachian Mountains of Macon County, North Carolina (southeastern USA). Landslides
5(1):31–44.
Zienkiewicz, O., Taylor, R., Zhu, J. Z. (2005). The finite element method: Its basis andfundamentals (6th edition). p752, Butterworth-Heinemann, British.

79