Computers and Geotechnics 88 (2017) 222–241

Contents lists available at ScienceDirect

Computers and Geotechnics
journal homepage: www.elsevier.com/locate/compgeo

Research Paper

The role of uncertainty in bedrock depth and hydraulic properties on thestability of a variably-saturated slope
Guilherme J.C. Gomes a,b, Jasper A. Vrugt b,c, Eurípedes A. Vargas Jr. a,⇑, Julia T. Camargo d, Raquel Q. Velloso a,
Martinus Th. van Genuchten e,fa

Department of Civil Engineering, Pontifical Catholic University of Rio de Janeiro, Rio de Janeiro, Brazil
Department of Civil and Environmental Engineering, University of California, Irvine, CA, USA
Department of Earth System Science, University of California, Irvine, CA, USAd
Department of Civil and Environmental Engineering, Stanford University, Stanford, CA, USAe
NIDES Interdisciplinary Center for Social Development, Federal University of Rio de Janeiro, Rio de Janeiro, BrazilfDepartment of Earth Sciences, Utrecht University, Utrecht, The Netherlandsbca r t i c l e

i n f oArticle history:
Received 30 November 2016
Received in revised form 16 February 2017
Accepted 22 March 2017

Keywords:
Bedrock depth
Hydraulic parameters
Slope stability
Unsaturated flow
Limit analysis
Safety factora b s t r a c t
We investigate the uncertainty in bedrock depth and soil hydraulic parameters on the stability of avariably-saturated slope in Rio de Janeiro, Brazil. We couple Monte Carlo simulation of a threedimensional flow model with numerical limit analysis to calculate confidence intervals of the safety factor using a 22-day rainfall record. We evaluate the marginal and joint impact of bedrock depth and soil
hydraulic uncertainty. The mean safety factor and its 95% confidence interval evolve rapidly in responseto the storm events. Explicit recognition of uncertainty in the hydraulic properties and depth to bedrock
increases significantly the probability of failure.
Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction
Hilly landscapes are typically mantled with soil and underlain
by a weathered bedrock zone that may extend tens of metersbeneath the surface before reaching fresh bedrock [1]. The design,
construction, and maintenance of buildings, infrastructure, andother man-made structures in such environments depends critically on the stability of the soil and underlying bedrock. It is
well-known that precipitation exerts a strong control on hillslope
stability via infiltrating water, which elevates pore pressures
within the soil mantle and reduces the shear strength. Thus, themore water is stored in the soil mantle, the more susceptible a hillslope is to landsliding. Numerous contributions to the geotechnicalliterature have investigated the triggering mechanisms and probability of occurrence of rainfall-induced landslides. This body ofwork has shown that the stability of hillslopes is controlled by⇑ Corresponding author.
E-mail addresses: guilhermejcg@yahoo.com.br (G.J.C. Gomes), jasper@uci.edu(J.A. Vrugt), vargas@puc-rio.br (E.A. Vargas Jr.), juliatcamargo@gmail.com(J.T. Camargo), raquelvelloso@puc-rio.br (R.Q. Velloso), rvangenuchten@hotmail.
com (M.Th. van Genuchten).
http://dx.doi.org/10.1016/j.compgeo.2017.03.016
0266-352X/Ó 2017 Elsevier Ltd. All rights reserved.

many different factors including surface topography [2], depth tounweathered bedrock [3–6], the hydraulic properties [7,8] andshear strength [9–11] including their spatial variability [12,13].
The stability of a hillslope is commonly characterized by a single integrated measure of its load carrying capacity. This measure,coined the factor of safety or safety factor (SF), requires detailed
computational analysis since field experiments are often impractical, time-consuming, labor-intensive and expensive. Many different methods have been developed in the geotechnical literature
to compute the SF of a slope, ranging from simplified infinite slope
approaches to more advanced limit equilibrium methods andsophisticated numerical procedures. These methods differ in theirunderlying assumptions and rigor, and consequently, may provide
conflicting results. Slope stability analyses are further complicatedby our inability to characterize adequately the hillslope interior.
Indeed, soil properties and/or related variables deemed necessary
for the different methods cannot be measured in the field with infinite precision, and consequently, the computed SF values are subject to considerable uncertainty [14].
Conventional (deterministic) methods for slope stability analyses calculate the SF without recourse to the underlying uncertainty
in soil and hillslope properties [15,16]. Probabilistic methods allow

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

designers to address issues beyond those that can be dealt withusing deterministic methods. Such methods allow us to calculateconfidence intervals of the SF due to uncertainty in soil properties
and stratigraphic conditions, data errors, and model structural
inadequacies [15,17]. The impact of uncertainty in the soil properties on the SF is relatively easily quantified using statistical methods [15], and has been explored by different authors [7,9,11,18,19].
This includes treatment of uncertainty in the porosity [20], specificweight [11,17], cohesion [11,17], and friction angle [11,17] of thesoil mantle, and the spatial variability of the hydraulic conductivity[7,20]. Data errors arise from improper, incomplete and/or inaccurate measurement methods [15]. For example, the subsurface
stratigraphy can exhibit considerable spatial variability, which, if
poorly characterized, compromises the accuracy of slope stability
studies [3,4]. What is more, spatiotemporal variability of the moisture content and distribution of the hillslope interior controls slope
stability, yet is difficult to characterize adequately with porous
flow simulators due to large uncertainties in the hydraulic properties and boundary conditions of the considered soil-bedrockdomain [15,21].
Geotechnical stability analyses under transient soil-water conditions require knowledge of the soil hydraulic properties, namely
the water retention function (or soil-water characteristic curve,
SWCC) and the unsaturated soil hydraulic conductivity function,
HCF. These two constitutive relationships generally exhibit considerable spatial variability, and are very time consuming, labor intensive, and costly to measure at the scale of a hillslope [22,23].
Whereas several authors have investigated the effect of the saturated hydraulic conductivity on slope (in) stability [12,24–26], relatively few contributions in the geotechnical literature have
explored properly the coordinated impact of SWCC and HCF uncertainty on the SF values derived from slope stability studies[4,6,18,27]. We agree wholeheartedly with Liang and Uchida [23]that a detailed characterization of the temporal and spatial variability of the moisture content of the soil mantle is warranted inslope stability studies. Such studies should also include proper
recognition of uncertainty in the soil hydraulic properties [8].
Liang and Uchida [23] evaluated, in a recent study, the impact ofsoil hydraulic parameter uncertainty on distributed moisture contents simulated with a three-dimensional variably-saturated flow
model. Using a detailed depth to bedrock map of a small catchment
in Japan and different parametrizations of the soil hydraulic functions, these authors found that saturation develops predominantly
at the soil-bedrock interface. The bedrock surface connects sparsely saturated regions [28], and thus determines, along with regolith thickness, the water pressure within the soil mantle’s pores.
Indeed, depth to bedrock is a key variable that controls subsurface
flow [29], and triggers landslides during rainfall events [30–32].
As hillslope interiors are costly, labor-intensive, and difficult toaccess, most slope stability studies use a relatively simple description of the bedrock surface [3,6,7,11]. Resulting topographic maps
often do not do justice to complex field measurements of the bedrock depth, which often demonstrate significant spatial variability[3,33] with a geometry that is difficult to characterize adequately
with some closed-form mathematical expression, while hydraulicand strength parameters can vary abruptly at the soil-bedrock
interface. Whereas some authors have used high-resolution bedrock depth maps to assess hillslope stability [28,32,34,35], existingstudies in the literature do not properly recognize the effect of bedrock depth uncertainty on slope stability assessments.
The purpose of the present study is to evaluate the marginal andjoint impact of bedrock depth and hydraulic uncertainty on thestability and probability of failure of a natural soil-mantled hillslope in Rio de Janeiro, Brazil. The soils in our experimental watershed are unsaturated during most of the year, but this state canrapidly change during periods of heavy rainfall. Therefore, [3]223advocated monitoring variations in pore water pressures (bothnegative and positive) at specific sites during the year in order toassess their impact on slope stability. We use the calibrated geomorhologic model of [33] to generate plausible maps of the bedrock depth. The model was calibrated against a rich data set of
distributed bedrock depth measurements using Bayesian inferencewith the DREAM algorithm [36,37]. Posterior maps of the simulated bedrock depth topography serve as input to a threedimensional finite-element (FE) water flow model of the
bedrock-soil domain, and are used to evaluate hillslope stability
for a 22-day rainfall record via numerical limit analysis [38] using
the ensemble of simulated transient pore water pressures. Our
analysis will also consider soil hydraulic parameter uncertainty
derived from pedotransfer functions using the Rosetta program[39]. We are particularly interested in the coordinated impact ofboth sources of uncertainty (bedrock depth and soil hydraulicproperties) on the mean SF and its associated 95% confidence interval instead of looking at their marginal impact. The framework presented herein allows geotechnical engineers to address slope
stability issues beyond those that can be addressed using deterministic methods, and provides a basis for reliability analyses ofgeotechnical hazards.
The remainder of this paper is organized as follows. Section 2
reviews the different building blocks of our integrated fluid flow
and numerical limit analysis methodology. This framework is usedherein to assess the SF values and corresponding failure probabilities. In Section 3 we are particularly concerned with quantificationof soil hydraulic and bedrock depth uncertainties. This is followed
in Section 4 with a detailed description of three different case studies of the same hillslope. Section 5 presents results of the threecase studies, whereas Section 6 discusses the implications of ourfindings to hydrologic and geotechnical modeling. Finally, Section 7
concludes this manuscript with a summary of our main findings.

2. Modeling framework
We developed an integrated framework for stochastic slope stability analyses. The probabilistic framework couples a Monte Carloalgorithm with a three-dimensional variably-saturated flow modelof the soil mantle and assesses hillslope stability via numericallimit analysis (NLA). The flow model solves numerically in threedimensions Richards’ equation and simulates transient values ofthe pressure head, wðÞ, for the spatial coordinates, x and y, anddepth, z, of the hillslope and time, t, of the supplied rainfall data
record. The values of wðx; y; z; tÞ in the center of each tetrahedral
element are used to characterize the soil-mantle’s water pressures.
These locations are indicated schematically with the red dot in themesh of Fig. 1. The evolving pore water pressures then serve as
input to NLA to solve numerically for the transient stress andvelocity distributions of the soil-mantle, which in turn are used
to calculate the SF during the rainfall record. The Monte Carlo algorithm is used to quantify the impact of uncertainty in the bedrock
depth and the soil hydraulic properties on hillslope stability byrepeated numerical solution of the NLA equations for the ensembleof pore water pressures simulated with the Richards-based flow
model.
Our framework is coded in MATLAB and integrates the twonumerical models so that users do not have to port data fromone model to the next [10], thus simplifying considerably slope
stability analyses. The MATLAB code has a built-in routine which
generates automatically prismatic and tetrahedral FE meshes forthe hillslope water flow model and NLA equations, respectively.
The code furthermore has pre- and postprocessing capabilities tosimplify model setup and visualization of the results. Fig. 1 provides a schematic overview of our coupled framework. Users have

224

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

Fig. 1. Schematic overview of the different building blocks of our stochastic slope stability analysis framework in MATLAB. The middle column details the different modelingsteps of our framework, whereas the left and right columns summarize the input and output data of the Richards-based hillslope flow model and NLA equations, respectively.
This setup of the consecutive modeling steps is in support of the Monte Carlo simulations to evaluate the impact of bedrock depth and soil hydraulic parameter uncertainty on
slope stability.

to supply several inputs to execute our framework, including slope
geometry (this includes a detailed topographic map of the bedrocksurface), initial and boundary conditions used by the hillslope
water flow model and NLA, and the hydraulic and strength parameters of the domain of interest. Users can also change default settings of algorithmic variables used for mesh creation andrefinement, as well as Monte Carlo simulations of the bedrock surface and soil hydraulic parameters. The framework also has built-in
capabilities for distributed computing, which speeds up significantly stability analyses by evaluating the Monte Carlo realizationsof the bedrock depth topography and soil-mantle’s hydraulic properties on different processors.
The next subsections describe the different modeling steps ofthe flowchart presented in Fig. 1, and summarize the spatial discretizations and initial and boundary conditions that are used tosolve numerically the governing equations of water flow and slope
stability. Further details are presented in Appendices A and B. In
Section 3 we then discuss our treatment of the bedrock topography
and soil hydraulic properties and their associated uncertainty insupport of probabilistic slope stability analyses.

2.1. Variably saturated model description
In this subsection we describe the Richards-based fluid flow
model used to simulate the transient pore water pressures of thesoil-mantled hillslope in Rio de Janeiro. The simulated pore waterpressures serve as input to NLA to compute confidence intervals of
the SF.
2.1.1. Governing equation
The hillslope water flow model solves the three-dimensionalpressure head formulation of Richards’ equation for the spatial
coordinates, x; y and z (L) of the soil-mantle given by

CðwÞ





@w
@w
@
@w
@
@w
þ
þ Ss SðwÞ
¼
KðwÞ
KðwÞ
@t
@t @x
@x
@y
@y


@
@w
þ
KðwÞ
þ1 ;
@z
@z

ð1Þ

where w (L) is the pressure head, CðwÞ (L1 ) denotes the soil watercapacity function, SðwÞ (-) signifies the degree of saturation, Ss(L1 ) represents the specific storage, KðwÞ (L/T) is the HCF, and t(T) denotes time. The degree of saturation is defined asSðwÞ ¼

hðwÞ
;
hs

ð2Þ

where hðwÞ (L3 =L3 ) is the moisture content, and hs (L3 =L3 ) denotes
the saturated water content. Note that Ss can be omitted in ourapplication since it does not affect the simulated moisture distribution and pore water pressures.

225

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

2.1.2. Temporal and spatial discretization
Numerical solutions of Eq. (1) require the users to define a suitable spatial and temporal discretization for the domain of interest.
We use herein a mesh consisting of linear prismatic elements composed of six nodes. Fig. 2a presents a three-dimensional plot of the
FE mesh of the experimental hillslope in Rio de Janeiro, which will
be discussed in more detail later. Unfortunately, one cannot build
the mesh without detailed knowledge of the surface and bedrock
topography. The model domain was generated in our program byimporting x; y; z coordinates for the surface and bedrock topography. The MATLAB toolbox we developed herein allows users todefine the density of the surface mesh (number of nodes/elements). The toolbox also allows users to define the number ofnodes in the z (vertical) direction between the hillslope surface
and the bedrock. The structured mesh presented in Fig. 2a is madeup of 9600 elements with 5797 nodal points. This mesh provides agood compromise between the CPU costs of the BFGS solver andthe corresponding accuracy of the simulated distributed pressure
head values. Indeed, as shown in Fig. 3, the use of a larger numberof elements (and thus nodes) hardly changes the numerical resultsof our flow model. Since our measurements of the soil mantle did
not reveal the presence of a well-defined stratigraphy with multiple layers [42], we conveniently assumed the hillslope interior tobe homogeneous. Nonetheless, our MATLAB flow model can solve
for water movement in heterogeneous soil-bedrock domains, if
deemed appropriate.
The BFGS solver implements a variable time step, Dt, to solve
numerically Eq. (1). This time step is reduced when convergence

0-100

ψ (cm)

Unfortunately, no analytic solutions exists of Eq. (1) for sloping
soils with a complicated stratigraphy and realistic boundary conditions. We therefore resort to numerical solutions, and use a fully
implicit finite-difference scheme with variable time steps to simulate spatially distributed values of the pressure head and moisture
content. The quasi-Newton Broyden-Fletcher-Goldfarb-Shanno(BFGS) method [40,41] was used to solve numerically in time Eq.
(1) at relevant x; y; z locations of the hillslope. This method calculates the Jacobian (derivative) matrix only once at each time,
thereby reducing considerably the CPU cost for highly discretized
FE meshes. The cited publications detail the numerical implementation of the BFGS method.

-200
-300
-400

22143 nodes; CPU time = 1786 s
5797 nodes; CPU time = 424 s

-500
0

100

200

300

400

500

600

t (hours)
Fig. 3. Simulated pressure head values, wðx; y; z; tÞ, at some random node in the soilmantle of our hillslope using a mesh with 5797 nodal points (solid line) and a much
more refined mesh composed of 22143 nodes (dashed line). The simulations appear
almost identical, yet the CPU cost of the densest mesh was more than four times
larger than its counterpart with 5797 nodes.

criteria of the simulated wðx; y; z; tÞ values are not satisfied withthe current value of Dt. Our flow model relates the increment of
Dt to the number of iterations used during the previous time step.
The initial, minimum and maximum time steps were set to0:01; 109 and 1 h, respectively. Based on the work of Paniconi
et al. [43] we did set the convergence threshold of w at 0.1 m.
2.1.3. Boundary and initial conditions
Numerical solutions of Richards’ equation were obtained by
imposing a precipitation rate, Q (L/T), across the surface boundary(C) as follows

qn ¼ Q ðtÞ at C;

ð3Þ

where qn represents the incoming flux, normal to the surface of the
FE mesh. Rainfall was recorded at a station near the study site
(Fig. 4). The measured hyetograph displayed three consecutive
storm events (days 2–4, days 7–13 and days 18–21), which were

Fig. 2. Examples of meshes composed of (a) prismatic and (b) tetrahedral elements for our soil-mantled hillslope. (a) The boundary conditions of Eq. (1) are separately
defined at the edges of the mesh of our flow domain. The precipitation data record, Q ðtÞ, is used as upper boundary condition for the surface of our structured mesh, whereas ano-flow boundary condition, qn ¼ 0 is assumed at the soil-bedrock interface, and at the upslope and lateral faces of the hillslope domain. Finally, a seepage face boundary
condition is used at the downslope face. (b) The boundary conditions of Eq. (A.4) are defined at the edges of the tetrahedral mesh used by the NLA. The symbols u; v and w
signify the strain velocities in the directions of the three spatial coordinates, x; y, and z of the hillslope, respectively. Their values follow logically from the assumed boundary
conditions, that is, a zero velocity equates to u ¼ v ¼ w ¼ 0.

226

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

time (days)

precipitation (mm/h)

0

2

4

6

8

10

12

14

16

18

20

22

7
6
5
4
3

soil-moisture state in the flow model to assess slope stability during the following 22-days. Of course, one can debate the choice of
our spin-up period. We note that, despite recent research (e.g. see[45]), no consensus exists on what constitutes a proper spin-up
time and period for particular applications. Without formal guidelines, our length of spin-up period is in agreement with [46]. It is
important to stress here that the sensitivity of our flow model tothe initial soil moisture state diminished rapidly during the 22day period because of the rather large intensity of the rainfall
events of the hyetograph (see Fig. 4).

2
1
0

100

200

300

400

500

time (h)Fig. 4. Atmospheric boundary conditions used in all of our numerical simulations.
The hyetograph envelopes a period of 22 days with hourly measurements of rainfall
rates (mm/hour) from a nearby meteorological station.

2.1.4. Soil hydraulic properties
Numerical solutions of Eq. (1) require knowledge of the hydraulic properties of the soil mantle. We use for this purpose the soil
water characteristic curve (SWCC) and hydraulic conductivity
function (HFC) of van Genuchten-Mualem (VGM) [47,48], which
are given by

Se ðwÞ ¼
separated well in time. The last of these three events exhibited thelargest intensity and volume of rainfall. In this study, we assumed
the infiltration flux at the surface nodes to be equivalent to measured rainfall rates. If, in response to this rainfall, the simulatedpressure heads of the surface nodes became larger than zero, thenthe excess water was removed instantaneously by surface runoff
and the pressure heads were set to zero [44]. If, at the end of sometime step, Dt, the infiltration rate at some surface node was larger
than the prescribed precipitation flux, then the infiltration flux
was set equal to the measured hourly rainfall rate, Q ðtÞ. We refer
to Fig. 2a for definition of the boundary conditions of the hillslope.
In summary, a seepage face boundary condition was assumed along
the downslope face. The remaining boundary nodes of the soilbedrock interface and the lateral and upslope faces assumed a noflow boundary conditions, that is, qn ¼ 0.
In the absence of detailed information about the soil-mantle’spressure head values at the start of the 22-day experimental period, we preceded our simulations with a one-year spin-up period
to derive reasonable initial values of wðx; y; z; 0Þ. We started thespin-up period with wðx; y; zÞ ¼ 500 cm and forced the flow
model with a yearly record of measured daily rainfall rates froma nearby meteorological station in Rio de Janeiro (see Fig. 5). At
the end of the 365-day spin-up period, the highest moisture contents (and thus pressure head values) were found at places of theslope with a relatively thin regolith depth. This simulated distribution of pressure head values was subsequently used as the initial

precipitation (mm/day)

60
50
40
30
20
10
0
0

50

100

150

200

250

300

350

time (days)
Fig. 5. Yearly record of daily precipitation rates used during the spin-up period. The
data originated from a nearby meteorological station.

hðwÞ  hr¼
hs  hr(ð1 þ jawjn Þ1mfor w 6 0
for w > 0

n
h
im o2
KðwÞ ¼ K s Se ðwÞL 1  1  Se ðwÞ1=m

ð4Þ

ð5Þ

where Se (-) is effective saturation, hr (L3 =L3 ) denotes the residual
water content, a (L1 ) is a proxy of the reciprocal of the air-entry
value of the soil, n and m are unitless shape parameters of the SWCC
and HCF, K s (L/T) denotes the saturated hydraulic conductivity and L(-) is a pore-connectivity parameter. We make the common
assumptions that m ¼ 1  1=n and L ¼ 0:5, and use the remaining
parameters, fhs ; hr ; a; n; K s g, to characterize the hydraulic propertiesof the soil mantle.
2.2. Numerical limit analysis
Numerical limit analysis (NLA) was used to evaluate thegeotechnical stability of the hillslope under variably-saturated
conditions. NLA as implemented in our study uses as input the porewater pressures simulated by the Richards-based flow model (seeFig. 1) and returns the corresponding estimates of the collapse
load, that is, the maximum load that can be sustained by the hillslope. We briefly summarize the main elements of NLA, and referreaders to the excellent review by Sloan [49] for further details
on limit analysis approaches.
We use the mixed-formulation, second-order cone programming implementation of [38] to solve for the stress and velocity
distributions in the three-dimensional finite element mesh of ourhillslope. The mixed-formulation assumes stresses to be constantwithin each FE, while velocities in the u; v and w directions are varied linearly between neighboring nodes [50]. The second-ordercone programming [51] exploits the underlying convexity of the
optimization problem and searches for the collapse factor that satisfies the equivalent static and kinematic principles of limit analysis [52].
A summarizes the NLA method described by Camargo et al. [38]and adopted herein for computation of the SF of the hillslope during the 22-day rainfall record.

3. Treatment of bedrock depth and soil hydraulic uncertainty
This section characterizes the bedrock depth topography andthe hydraulic functions of the soil mantle. We are particularly concerned with proper treatment of their uncertainty.

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

3.1. Uncertainty of bedrock topography
In a recent paper, Gomes et al. [33] developed a geomorphologic
model to predict the spatial distribution of depth to bedrock (DTB)
from high-resolution topographic data, numerical modeling andBayesian analysis. We used the DTB model of [33] and refer interested readers to this publication for further details. Bayesian analysis was used to reconcile the DTB model with distributed bedrock
depth measurements. We used Markov chain Monte Carlo simulation (MCMC) with the DREAM algorithm [36,37] for posterior inference in the Bayesian framework. This approach searched the DTB
model parameter space for posterior solutions that ”best” honor
the observed bedrock depth data. The posterior mean bedrock
topography simulated with the DTB model was shown to be ingood quantitative agreement with point observations of the regolith depth. The associated prediction intervals appeared rather
large, yet enveloped the vast majority of bedrock depth observations. Note that [8] used the DREAM algorithm to evaluate slope
stability under unsaturated conditions. We refer to this publication
for a geotechnical application of Bayesian analysis. B summarizes
how we derived the posterior mean bedrock topography and asample set of surfaces that quantify its associated 95% prediction
intervals.
DTB measurements were carried out in two adjacent watersheds in the Papagaio river basin, Rio de Janeiro, Brazil. The two227

watersheds studied extensively by several authors in the literature
after mass movements occurred in 1996 (see for instance [53]). The
geographic location of the field site is depicted with a red cross in
Fig. 6, which also displays (right hand side) the topography of thetwo experimental basins. The white dots indicate the measurement locations of the bedrock depth using a light dynamic cone
penetrometer. The eight black circles represent locations where
additional soil samples were collected via drilling with a mechanical auger. We refer to [42] for a detailed description of the field
campaign and measurements.
We next illustrate the results obtained with the DTB model, anddisplay in Fig. 7 the main morphological features of the experimental hillslope. The top panel presents a three-dimensional view ofthe surface elevation of the hillslope as derived from highresolution topographic data. The elevation of the hillslope ranges
between 280 and 340 m above mean sea level, with slopes thatare steepest near the crest (hilltop) and become increasingly planar
downslope towards the drainage valley. The color map that is projected onto the surface topography depicts the posterior mean
depth to (unweathered) bedrock simulated with the DTB model.
The DTB model produced a rather smooth bedrock surface fromthe hilltop to the drainage channel in the valley. The bottom panelof Fig. 7 presents a two-dimensional cross section of the hillslopeand depicts the topography of the surface (solid black line) andDTB simulated posterior mean bedrock surface (solid red line) ofFig. 6. Geographic overview (left) and shaded relief image (right) of the Papagaio river basin experimental field site (red cross) in Rio de Janeiro, Brazil. The white and black
circles on the relief image indicate, respectively, the measurement locations of the dynamic cone penetrometer and the places were drilling was used with a mechanical auger
to collect soil samples. Transect AA’ was used to summarize the results of the DTB model. (For interpretation of the references to colour in this figure legend, the reader isreferred to the web version of this article.)

228

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241
Table 1Comparison of the mean (top row), maximum (middle row) and minimum (bottomrow) gradients of the surface topography (middle column) and posterior mean soilbedrock interface simulated with the DTB model (right column).
Attribute

Surface

Bedrock

Mean gradient (°)
Max. gradient (°)
Min. gradient (°)

30
44
9

23
43
3

3.2. Uncertainty in the soil hydraulic properties

Fig. 7. Topography and bedrock surface of the hillslope under investigation. (a)Color projection of the posterior mean bedrock depth simulated with the DTB
model. The vertical bars along the axes denote the elevation of the hillslope. (b)Two-dimensional plot of the surface topography (solid black line) and posteriormean bedrock surface (solid red line) of transect AA’. The elevations refer to thebase of the hillslope (drainage valley). The dark and light gray regions denote the95% prediction uncertainty intervals of the DTB model due to parameter uncertainty and total uncertainty. (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)transect AA’. The dark gray region in this graph represents the 95%
prediction intervals of the DTB model due to parameter uncertainty. This uncertainty was derived from the M posterior realizations stored in the three-dimensional array h, as detailed inAppendix B.
Note that in Fig. 7 we also portrayed the 95% ranges of the total
DTB model uncertainty (gray region). These enlarged prediction
intervals of the bedrock depth summarize the combined impactof DTB parameter and model uncertainty, respectively, and arecomputed by adding to the simulated bedrock depths of each posterior realization of 
h a random draw from the model error distribution. This step introduces irregularities in the soil-bedrock interface
and characterizes structural inadequacies (epistemic errors) of the
DTB model. We made herein the common assumption of a zeromean normally distributed model error with standard deviation
equal to the root mean square error of the posterior mean DTB
model fit. The resulting ”perturbed” ensemble of M different bedrock depth maps summarizes the uncertainty in the bedrock topography underneath the hillslope and is now used in our Monte Carloframework for probabilistic slope stability analysis.
Table 1 lists summary statistics of the gradient of the (measured) hillslope topography (middle column) and the DTB simulated posterior mean bedrock surface (right column). Note thatthe surface topography exhibited the largest gradients with mean
and minimum values that are larger than their counterparts ofthe bedrock surface. The maximum gradients of the hillslope surface and soil-bedrock topography were very similar (43–44°).

Now that the soil-bedrock interface has been delineated, we areleft with characterization of the hydraulic properties of the soil
mantle. In the absence of soil water retention and unsaturated
hydraulic conductivity data, we used the Rosetta program of [39]to estimate values of the VGM parameters in the soil hydraulicfunctions given by Eqs. (4) and (5). Rosetta uses a set of hierarchical pedotransfer functions to estimate the VGM parameters frombasic soil data, such as soil textural class, sand, silt, and clay percentages, bulk density, and direct water retention observations.
We followed the methodology developed by [22] and quantified
the uncertainty in the hydraulic properties (SWCC and HCF) of
the soil mantle as follows.
First, 52 measurements of soil texture (sand, silt and clay %) ateight different locations and depths (see black circles in Fig. 6)
were used to compute the mean texture, a, of the soil mantle
and its associated 3  3 covariance matrix, b. We drew randomly
a large sample of P texture samples from N 3 ða; bÞ, a three-variate
normal distribution with mean a and covariance matrix b, andevaluated each of the P samples using Rosetta. This did result ina large ensemble of P realizations of the VGM soil hydraulic parameters, fhs ; hr ; a; n; K s g. Table 2 summarizes the results of our Rosetta
ensemble, including mean values of the soil hydraulic parameters(third column), their standard deviations (fourth column), and linear correlation coefficients (last five columns). The negative correlation coefficient (-0.86) between log10 ðaÞ and log10 ðnÞ is inagreement with field data and supported by statistical analysis ofthese two soil hydraulic parameters [54,18]. Indeed, Phoon et al.
[18] argued that this correlation between a and n should behonored otherwise this might lead to a SWCC that does not characterize accurately the retention function of the soil domain underinvestigation. A similar conclusion was reported independently
by Scharnagl et al. [22] who showed that a proper treatment of
the correlation structure between the soil hydraulic parameters
is a prerequisite to accurate sampling of the water retention functions. The relatively low correlation between the VGM parameter
pairs hr  hs and log10 ðaÞ  log10 ðK s Þ is in agreement with previous
findings (e.g., [22]).
The P different ROSETTA predictions of the VGM parameters
were used next to characterize the uncertainty in the hydraulicproperties of the soil-mantle. Alternatively, one could generate P
realizations of the VGM parameters by drawing randomly from afive-variate normal distribution with mean and covariance matrix
from the P Rosetta realizations. We illustrate the results of this second approach in Fig. 8, which shows the mean (red line) and 95%
confidence intervals (gray region) of the water retention (top)and hydraulic conductivity (bottom) functions. The two functions
suggest an excellent water holding capacity of the soil mantle, withpressure heads that must be below 80 cm for air to enter the
pores. This finding is consistent with our soil textural data, which
demonstrate relatively high percentages of silt and clay [42]. The
confidence interval of the SWCC is almost perfectly symmetrical,
being similarly large at saturation and the residual moisture content, and smallest in the interval between the two inflection points(moisture contents between 0.2 and 0.4 cm3/cm3). Note that the229

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241
Table 2
Mean values, standard deviations and correlation coefficients of the Rosetta [39] predicted soil hydraulic parameters.
Parameter

Unit

hr
hs
log10 ðaÞlog10 ðnÞlog10 ðK s Þ

cm3 cm3
cm3 cm3
cm1
–
cm/h

Mean

Standard deviation

0.05770.4337
2.27190.2146
0.05830.01530.01590.07480.01400.0128θ (cm3 /cm3 )hr

hslog10 ðaÞlog10 ðnÞlog10 ðK s Þ

1.00
0.09
0.15
0.58
0.97

1.00
0.24
0.04
0.12

1.00
0.86
0.10

1.00
0.57

1.00

Table 3
Constant values of the geotechnical parameters used in our numerical simulations.

0.5
0.4

Type

Parameter

Value

Unit

Soil and water
Properties

cw
cnat
csat

10
18
20

kN/m3
kN/m3
kN/m3

Strength
Parameters

c0
/0
/b

5
30
15

kPa

0.3
0.2

mean
95% parameter uncertainty

0.1
0 4
-10

-10

3

-10

2

-10

1

-10

0

ψ (cm)(a)10

K(ψ) (cm/h)Correlation coefficients

10

2




of the stress field, r in Eq. (A.3). Table 3 lists the values of the specific weights and strength parameters used herein. We conveniently
assumed constant values of the strength parameters in all of ournumerical simulations. This approach ignores spatial variabilities
in c0 and /0 [11], and may lead to inadequate characterization of
the stability of hillslopes during rainfall events [8,14]. Nevertheless, the use of constant strength parameters simplifies investigations into the impact of bedrock depth uncertainty and/or soil
hydraulic uncertainty in slope stability studies.

0

4. Case studies

10
10
10

-2

-4

-6

-10

4

-10

3

-10

2

-10

1

-10

0

ψ (cm)(b)
Fig. 8. 95% prediction intervals of the hydraulic functions of the soil mantle derived
with Rosetta using 52 different soil textural samples from eight different
measurement locations on the hillslope in Rio de Janeiro, Brazil: (a) water retentionfunction, and (b) unsaturated hydraulic conductivity function. The hydraulicfunctions of the mean VGM parameters are separately indicated with solid redlines. (For interpretation of the references to colour in this figure legend, the reader
is referred to the web version of this article.)confidence interval at w ¼ 0 cm is determined only by the value ofhs , whereas the ranges at w ¼ 104 cm are determined by the values of hr and n. The HCF exhibits asymmetrical confidence intervals
that become increasingly smaller (even on a logarithmic scale)with decreasing pressure heads as the soil dries out. Uncertainty
in the HCF is determined mostly by K s .
3.3. Strength parameters
Numerical simulations with the flow model and NLA also
require estimates of the strength parameters, c0 ; /0 and /b in Eq.
(A.1). Additionally, values of the specific weights of natural soil,
cnat , saturated soil, csat , and water, cw , are required for computation

The MATLAB framework with its various elements and settings
pertained to the variably-saturated flow model and NLA was used
next to assess the impact of uncertainty in the bedrock depth andhydraulic parameters on the stability of the soil-mantled hillslopein Rio de Janeiro. A topographic map of the hillslope was shown
previously in Fig. 7, including a three-dimensional color projectionof the underlying bedrock surface. We also presented in Fig. 2 two
example FE meshes of the hillslope that will be used by the waterflow model and NLA. We describe below three different case studies to illustrate our findings. The first two case studies examine
only the marginal (individual) impact of bedrock depth uncertainty(case study 1) and soil hydraulic parameter uncertainty (case study2) on slope stability. The third case study evaluates their joint or
combined impact on hillslope stability.
4.1. Case study 1: Uncertainty in bedrock topography
This first case study considered the impact of bedrock topographic uncertainty on slope stability. We necessarily assumed
mean values of the hydraulic parameters of the soil mantle (seeTable 2). A total of 100 different bedrock depth maps were drawn
at random from the ”perturbed” ensemble of posterior DTB-model
simulations. This number of realizations was deemed sufficient tocharacterize the uncertainty in the bedrock topography underneath the hillslope (see Fig. B.1). The MATLAB framework createdautomatically, for each bedrock depth map, two FE meshes: onecomposed of prismatic elements and one comprising tetrahedral
elements.
To provide a better understanding of the variations in the bedrock topography of the Monte Carlo ensemble, consider Fig. 9
which displays graphically the regolith surface of two plausible

230

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

A’A

hilltopsideslopesurface⏐⏐⏐soil-bedrock interface

↑
)
E(

base

(a)

(b)

Fig. 9. Two different realizations of the bedrock surface for transect AA’ of the hillslope in Rio de Janeiro. The posterior mean soil-bedrock interface (EðhÞ) simulated with the
DTB model is shown separately in each panel by the dashed blue line. The right plot (b) clarifies the position of the landscape elements ‘‘hilltop”, ‘‘sideslope” and ‘‘base” usedherein to discuss our findings. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)realizations, coined H1 and H2 . The two samples are characterizedwith the transect AA’ depicted previously in Fig. 7. The posteriormean bedrock surface simulated with the DTB model is indicatedseparately in each plot using the dashed blue line. We used green(a) and red (b) colors to differentiate between the bedrock topographies of H1 and H2 . The shape of the two surfaces in Fig. 9 will
enhance our understanding of the simulated flow patterns and failure mechanisms. Note that the two bedrock surfaces exhibit manyof the common features found in Rio de Janeiro [3]. To simplify
interpretation we will use the landscape elements ”hilltop”,”sideslope” and ”base” to discuss our findings. Their position isindicated in Fig. 9b.
The 100 prismatic FE meshes served as input to the soil waterflow model to create an ensemble of transient pressure head valuesof the hillslope interior during the 22-day experimental period. The
100 transient simulations generated input to the tetrahedral FE
meshes of the NLA to evaluate slope stability.

4.2. Case study 2: Uncertainty in the soil hydraulic parameters
In this study case, the uncertainty in the VGM parameters wasevaluated by defining a prior distribution of the soil hydraulicparameters. This informative prior distribution can significantly
reduce the uncertainty of the soil hydraulic parameters in backanalyses [8] and avoid unrealistic water retention and hydraulicconductivity functions [18]. A total of 100 samples of the soil
hydraulic parameters were generated by drawing from a multivariate normal distribution N ðlp ; Rp Þ with mean lp and covariancematrix, Rp (Section 3.2). We ignore uncertainty in the bedrock
depth and assume the posterior mean bedrock topography, EðhÞ,
as illustrated previously in Fig. 9.

Utilizing the same method as adopted in the previous subsection, two samples of soil hydraulic parameters were drawn
from N ðlp ; Rp Þ to help us better understand the effects of the
SWCC and HCF on the Monte Carlo results. The two samples arereferred to hereafter as p1 and p2 , respectively, and assume different values of the soil hydraulic parameters p ¼ fhs ; hr ; a; n; K s g, aslisted in Table 4. For completeness, the last row of Table 4 lists
the mean values of the soil hydraulic parameters, lp . The saturated
conductivity of p1 (K s ¼ 4:32 cm/h) was found to be higher than
the sample mean (K s ¼ 1:15 cm/h). The opposite is true for p2
which exhibited a much lower saturated conductivity(K s ¼ 0:76 cm/h).
4.3. Case study 3: Combined uncertainty in bedrock depth and soil
hydraulic functions
Case study 3 treats the combined effect of uncertainty in the
bedrock surface and the soil hydraulic parameters. The M realizations of the bedrock topography now have different hydraulicparameter values of the soil mantle by drawing from the multivariate normal prior distribution of the VGM parameters (Section 4.2).
Specifically, we illustrate our findings with the help of two
bedrock-soil hydraulic parameter samples, coined S13 and S23 .
These two samples use the bedrock surfaces, H1 and H2 of Fig. 9
and the hydraulic parameter values, p1 and p2 of Table 4, i.e.,
S13 ¼ fH1 ; p1 g and S23 ¼ fH2 ; p2 g, where the subscripts of S13 referto sample 1 of the case study 3, etc.
Table 5 summarizes the three case studies used to evaluateslope stability in the presence of bedrock depth and soil hydraulicparameter uncertainty. The table also lists several samples used toclarify our findings. Results are discussed in the next section.

Table 4
Values of the VGM soil hydraulic parameters used in case study 1 (lp ), and their values of samples p1 and p2 , used in case studies 2 and 3.
Samplehs (cm3 cm3)hr (cm3 cm3)a (cm1)n (–)

K s (cm/h)

p1
p2

0.453
0.425
0.429

0.036
0.063
0.056

0.00560.00420.0055

1.66
1.68
1.63

4.32
0.76
1.15

lp

Table 5
Summary of the three case studies and their samples used to illustrate our findings.
Case

Source of uncertainty

Samples

Bedrock surface

Hydraulic parameters

1

2


h

lp

S11 ¼ fH1 ; lp g

S21 ¼ fH2 ; lp g

2

Eð
hÞ

N ðlp ; Rp ÞhÞ; p1 g
S12 ¼ fEð

S22 ¼ fEð
hÞ; p2 g

3hN ðlp ; Rp Þ

S13 ¼ fH1 ; p1 g

S23 ¼ fH2 ; p2 g

1

231

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

5. Results
5.1. Uncertainty in the safety factor
Fig. 10 shows the evolution of the mean SF (solid red line) and
its associated 95% confidence intervals (gray region) for case study
1 (a), case study 2 (b) and case study 3 (c). To better interpret theresults, we also plotted separately in the figure the measured hyetograph and listed the precipitation amounts on the left y-axis.
Fig. 10a shows that the first storm between t ¼ 0 and t ¼ 8(days) decreased the mean SF steadily from values of 1.70 to about
1.60, in response to increased storage of water in the soil mantle.
The associated 95% intervals of the SF indicated a small controlof the bedrock topography on slope stability during the first storm
event. The second storm event between days 9 and 13 exerted far
more control on slope stability. The SF plummeted after the peak
rainfall on day 11 from 1.60 to about 1.25 on day 13. The second
storm event also increased substantially the width of the 95%
intervals of the SF. During the next four consecutive dry days (days
14–17), the stability of the hillslope improved somewhat with amean SF on day 17 of about 1.35, while the associated 95% ranges
that exhibited a smaller spread. The third and last storm event(days 18 to 22) was most severe, with its relatively large precipitation amount threatening the stability of the hillslope. Indeed, themean SF of the hillslope deteriorated rapidly from values of about
1.35 on day 18 to 0.95 at the end of the 22-day simulation period.

The hillslope had become unstable and prone to mass movement
and failure during the last few days of this third storm event.
The 95% intervals of the SF became relatively low, yet exhibited amuch larger spread on day 18 (1.2–1.4) than on day 22 (0.9–1.0).
Thus, the bedrock topography exerted far more control on theslope stability on day 18 than on day 22.
We next analyzed how uncertainty in the soil hydraulic parameters affected the simulated SF values (Fig. 10b). The mean SF, represented by the red line, was very similar to the values obtained forcase study 1 (Fig. 10a). However, during the first storm event
between 0 and 5 days, the 95% confidence intervals exhibit very
tight bounds, thus suggesting no control of the soil hydraulicparameters on slope stability. We will revisit this assumption later.
Fig. 10b shows that, after some precipitation threshold is reached,
the 95% intervals of the SF increase. This occurs after about t ¼ 5
days. Between the second and third precipitation events, thebounds reached their maximum spread. By comparison, Fig. 10aand b show that, a priori, the effects of both sources of uncertainty(bedrock depth and soil hydraulic parameters) are quite similar.
Fig. 10c presents plots of the mean SF and its uncertainty
bounds for case study 3. While the mean values are very similaras those of the other two case studies, the 95% confidence intervals
are now slightly greater. The uncertainty during the first rainfall
event (until t ¼ 5 days) seems related more to the variability inbedrock depth, rather than the soil hydraulic parameters. Indeed,
for this simulation period, Fig. 10c and a depict equivalent uncer-

time (days)

time (days)
4

6

8

10 12 14 16 18 20 22

1

0

2.0

4

1.4

5
1.2

6
7

(a)

precipitation (mm/h)

1.8
1.6

4

6

8

10 12 14 16 18 20 22

1

2
3

2

1.8

2
1.6

3
4

1.4

5
1.2

6
7

1.0

8

2.0

(b)

safety factor

2

safety factor

precipitation (mm/h)

0

1.0

8
0

100

200

300

400

500

0

100

200

300

400

500

time (h)time (h)time (days)
2

4

6

8

10 12 14 16 18 20 22

1

2.0
1.8

2
1.6

3
4

1.4

5
1.2

6
7

(c)

safetyfactor

precipitation (mm/h)

0mean SF
sample S13
sample S23
95% confidence intervals
b precipitation

1.0

8
0

100

200

300

400

500

time (h)Fig. 10. Uncertainty in the safety factor (SF): 95% confidence intervals (gray region) and the mean SF (red line) for case study 1 (a), case study 2 (b) and case study 3 (c), areshown on the right axes. The left axes show hourly precipitation amounts. The SF values of samples S13 and S23 are indicated in (c). The vertical dashed blue lines in the plots
represent the time levels that were analyzed. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

232

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

tainty limits. After this time, the multivariate normal prior distribution of the soil hydraulic parameters also affects the confidence
intervals of the SF, which are now wider than those of the first twocase studies. Notice also the relatively large differences between
the SF traces of samples S13 and S23 . These two realizations differ
substantially in their assessment of slope stability after the second
storm event, and will be investigated in more detail later.
5.2. Pressure head and failure zone dynamics
We now investigate how uncertainty in the bedrock depth andsoil hydraulic parameters affects pressure head dynamics within
the hillslope interior and the failure zones at collapse. Fig. 11
depicts graphically for samples S13 (top part) and S23 (bottom part)the simulated pressure values at t ¼ 5; 12; 16 and 19 days (left-toright) at the soil-bedrock interface (top panels) and along crosssection AA’ (middle panels). The two bottom panels in both graphs
display the velocity vectors at collapse. The samples S13 and S23
exhibit contrasting bedrock topographies (previously depicted inFig. 9) and soil hydraulic parameter values (listed in Tables 4 and
5). The snapshots from left to right match exactly the times ofthe vertical dashed blue lines in Fig. 10 and characterize different
saturation conditions of the soil-mantle in response to the successive rainfall events. To simplify graphical interpretation, the pressure head values share a common color bar. We first discuss theresults of sample S13 displayed in the top three panels, followedby those of S23 in the bottom three panels.

The top panel of four plots of Fig. 11a displays simulated pressure heads of sample S13 at the soil-bedrock interface. Initially, at
t ¼ 5 days, the soil-bedrock interface appears almost entirely
unsaturated, with the exception of a small area of saturation at
the bedrock surface in the middle part of the hillslope. This saturation condition is a leftover of the one-year spin-up period and notcaused by the first storm event. In fact, the yellowish-green color isdominant which indicates negative w values for large portions of
the soil-bedrock interface. Results show that the size of the saturated area increased slowly and upward along the drainage channel during the next three time snapshots in response toantecedent rainfall. This finding is not surprising as moisture contents are expected to increase most in shallow regions [29], where
the distance to the impermeable layer is relatively small, thus
causing a rapid build-up of infiltrated rainwater. Note that thepressure head values within the saturated areas of the soilbedrock interface tended to increase between t ¼ 5 and t ¼ 16
days (darker blue).
While the top panel of Fig. 11a highlights important changes tothe extent of the saturated zone at the soil-bedrock interface, themiddle panel visualizes the evolution of the saturated area alongcross section AA’. At t ¼ 5 days the soil is mostly unsaturated,
except for a small area of the transect where the regolith is thinnest. After the second storm ceased, at t ¼ 12 days, the antecedent
rainfall had caused the area of saturation of AA’ to expand from thesideslope to the base. This area remained saturated until at least
the last snapshot at t ¼ 19 days. Note that near-saturated zones

Fig. 11. Pressure head and failure zone dynamics of samples S13 (top part) and S23 (bottom part): Distributed values of the simulated pressure heads, wðx; y; z; tÞ at the soilbedrock interface (top panel) and along cross section AA’ (middle panel) at t ¼ 5; 12; 16 and 19 days (plots from left-to-right). For completeness we also displayed in thebottom panel the velocity vectors at collapse.

233

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241(light blue areas) developed initially only at the soil-bedrock
interface.
The bottom panel of Fig. 11a presents a two-dimensional projection (cross section AA’) of the velocity vectors of sample S13 atcollapse. The four snapshots in time, as derived from NLA, visualize
the possible collapse of the hillslope under different saturation
conditions [49]. At t ¼ 5 days the velocity vectors indicate an
approximately circular failure surface at the hilltop. The soil moisture conditions at this time were clearly unsaturated (see middle
panel). In the collapse region, topographic gradients were relatively high, with failure occurring due to loss of suction (to be discussed in more detail later). At the end of the second storm event,at t ¼ 12 days, the shallow soil of the sideslope had become fully
saturated, thereby inducing failure in the region of minimum regolith thickness. Then, after a few drier days at t ¼ 16 days, the velocity vectors extended uphill to areas with the highest topographic
gradients near the hilltop. Subsequently, at t ¼ 19 days, after two
more days of rainfall (third storm), the failure zone concentrated
again on the thin regolith zone of the sideslope.
The results of sample S23 in the bottom three panels appear very
similar to those presented previously for S13 but with two main
exceptions. First, the pressure head values at the soil-bedrock
interface demonstrate an enlarged area of saturation which covered a much larger portion of the drainage channel and the sideslope. This is most evident in the snapshots of the top and middle
panels at t ¼ 16 and t ¼ 19 days, and can be explained in part by
the much lower saturated conductivity of sample S23(K s ¼ 0:76 cm/h) compared to S13 (K s ¼ 4:32 cm/h). The lower conductivity limits rainfall infiltration, thereby promoting nearsaturated conditions at the hilltop and enhanced storage of wateralong the sideslope and base (see middle panel in Fig. 11b at
t ¼ 12 days). This enlarged area of saturation can be explained further by the bedrock topography of sample S23 , which differs in two
important ways from S13 . The deeper mantle at the base of the hillslope acts as a storage area for water [29], with higher values of thesoil water pressure head, particularly at the soil-bedrock interface.
The thinner soils upslope from the inflection point to the hilltop
expand uphill the area of saturation to cover a larger portion of

the drainage channel. Second, the collapse vectors (red arrows inbottom panel) of S23 at t ¼ 16 days do not extent to the hilltop
as is the case for sample S13 but remain confined to a small areaof the soil mantle closest to the underlying bedrock. This finding
is best understood by considering more closely the pressure head
distributions and flow patterns within the slope interior and failurezone.
The differences in the extent of saturated area of samples S13(squares) and S23 (circles) are illustrated in Fig. 12, which presents
traceplots of the fraction (expressed as percentage, %) of saturated
nodes at the soil-bedrock interface (solid symbols) and at the soil
surface (open symbols) during the 22-day study period. The simulated traces confirm our findings of Fig. 11 in that sample S23
exhibited the largest fraction of saturated nodes at the soilbedrock interface as a consequence of its relatively low saturated
conductivity and thin soil mantle in the upper part of the hillslope.
The low saturated conductivity of the soil mantle of S23 has important consequences since it slows down the vertical infiltration of
rainwater, thereby leading to saturated conditions of a significant
fraction of the surface nodes during intense storm events. In thedays after each storm event, the saturated surface nodes of S23
dry slowly as water moves to deeper nodes within the soil mantle.
In comparison, far fewer (surface) nodes of sample S13 achieved
saturated conditions during the course of our 22-day rainfall
record due to the higher saturated conductivity of the soil mantle.
We next investigated in more detail the pressure head values
within the failure zone of our hillslope. The failure zone may
migrate as a result of the variable moisture state of the soil mantle.
We used the velocity vectors of Fig. 11 (bottom two panels) todelineate the FE nodes of the evolving failure zones of samples
S13 and S23 . Fig. 13 shows a plot of their mean pressure head, wf ,during the 22-day rainfall recorded, marked by open squares(S13 ) and solid circles (S23 ) symbols. The simulated trajectories provide further insight into the impact of the soil hydraulic properties(notably the saturated conductivity) and bedrock surface topography on slope stability. The trajectories of the two samples differ
substantially, yet demonstrate the steadily increasing amounts ofwater being stored in their respective failure zones. At the end of
the simulation, both samples simulate positive mean pressure

time (days)
0

2

4

6

8

10

12

14

16

18

20

time (days)

22

100

0

70

6

8

10

12

14

16

18

20

22

0

-100

60
50

ψ f (cm)saturated nodes (%)

80

4

100

bedrock (S13 )bedrock (S23 )surface (S13 )surface(S2383)
bedrock S11

90

2

40
30

-200

-300

20

sample S13

-400

10

sample s33se
S23
amostra

0
0

100

200

300

400

500

time (h)Fig. 12. Time series plot of the percentage (%) of saturated nodes at the soil-bedrock
interface (solid symbols) and at the soil surface (open symbols) of transect AA’during the 22-day experimental period for samples S13 (squares) and S23 (circles).
The two realizations differ in terms of their bedrock topography and hydraulicproperties of the soil mantle. The four vertical dashed blue lines match the times ofthe snapshots depicted previously in Fig. 11.

-500
0

100

200

300

400

500

time (h)Fig. 13. Evolution of the mean pressure head, wf , in the failure zones of samples S13(open squares) and S23 (solid circles) for transect AA’ of our hillslope. The verticaldashed blue lines signify the previously analyzed simulation times. (For interpretation of the references to colour in this figure legend, the reader is referred to the
web version of this article.)

234

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

heads due to the large rainfall depths of the three consecutive
storm events. In general, the mean pressure head of the failurezone increases in response to rainfall infiltration during and/or
immediately after a storm event, and remains constant ordecreases somewhat during drier days. The first storm between
days 1 and 6 of the 22-day hyetograph may have been sufficient
to saturate the nodes of the failure zone of sample S23 as illustrated
by the simulated positive mean pressure heads at t ¼ 6 days. The
second storm event was needed to saturate the nodes of the failurezone of S13 , with the mean pressure head becoming positive for the
first time at t ¼ 13 days. The next 120 h after this second storm
event barely changed the conditions within the failure zone of
S23 with FE nodes that continuing to be saturated as evidenced
by the simulated positive pressure heads. During this same 5-day
period the mean pressure head of the failure zone of S13 decreased
from positive values of about +40 cm at t ¼ 13, the day after thesecond storm, to negative values of about 100 cm at the end of
day 17. The rainfall from the third and last storm of the hyetographbetween days 18 and 21 had the most impact on sample S13 in thatthe mean pressure head of its failure zone that changed rapidly on
day 18 from negative to positive values. The failure zone continued
to be saturated during the remaining days of the 22-day record.
The mean pressure head of sample S23 was not much affected by
the third storm since the nodes of its failure zone were already saturated at t ¼ 18 prior to the start of this last precipitation event.
Altogether, the results of the two samples demonstrate a more
rapid build-up of water in the failure zone of sample S23 , withnodal pressure heads quickly attaining large positive values. As
discussed before, this finding is explained by the low conductivityof the soil mantle of S23 at saturation, and thinner depths of its
regolith uphill from the inflection point of the hillslope. Nevertheless, both samples simulated a saturated failure zone at the end of
the 22-day hyetograph. As previously shown in Fig. 10c, samples
S13 and S23 produced SF values that fell within its 95% ranges. In
other words, the two trajectories of Fig. 13 are a plausible characterization of the pressure head values within the failure zone of thehillslope.
In summary, in this section we have used two samples with
contrasting permeabilities and bedrock topographies to analyze

the spatiotemporal dynamics of saturated areas within the hillslope interior. The transient pore water pressures of these two
samples led to different locations and mean pressure heads ofthe failure zone within the soil. The next section analyzes moreclosely how regolith depth, bedrock topography and soil hydraulicproperties impact the simulated pressure heads, and thus the porewater pressures, of the soil mantle, and their control on the SF of
the hillslope.
5.3. Impact of regolith depth, bedrock topography and soil hydraulicproperties on hillslope stability
The simulations thus far have demonstrated that a deeper soil
mantle accelerates the rapid buildup of water (unless saturated
conductivity of the soil mantle is low). This finding is not surprising since deeper soils can store more water. Our simulationsindeed showed that moisture content will increase most quickly
in areas with a shallow bedrock surface, where a temporary watertable may developed in response to antecedent rainfall. This finding is in agreement with previous literature [55,29]. Fig. 14 presents scatterplots of the relationship between the depth of thesoil mantle (top panel) and the pressure head at the soil-bedrock
interface at t ¼ 5; 12; 16 and 19 days. The data points pertain to
the case studies 1 (sample S11 , solid squares) and 3 (sample S13 ,blue triangles). As listed in Table 5, these samples have the same
bedrock surface (H1 ) but different values of the hydraulic parameters. We used Pearson’s correlation coefficient, q 2 ½1; 1, to quantify the strength of the relationship between the scatter in eachsnapshot of Fig. 14. The stronger the correlation between the variables, the closer the value of q is to 1 or +1. Note that each scatterplot displays the values of q1 and q3 for the samples of case studies
1 (S11 ) and 3 (S13 ), respectively. The results in the top panel ofFig. 14 confirm the presence of a negative relationship between
the depth of the regolith layer and the corresponding pressure
head at the soil-bedrock interface. Indeed, the larger the depth of
the regolith, the lower the value of the pressure head at the bottomof the soil mantle. This is noticeable in the snapshots at all fourtimes. The correlation between bedrock depth and the pressure
head of the soil-bedrock interface increased with time for bothFig. 14. Relationship of regolith depth (top panel) and bedrock topography (bottom panel) to changes in the pressure head (w) distribution at the soil-bedrock interface:
Pearson’s correlation coefficients of samples S11 (solid squares) and S13 (blue triangles) illustrate different patterns of case studies 1 and 3, respectively.

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

samples during our 22-day study. Nevertheless, the scatter appears
somewhat complicated as reflected by the banana-shaped correlation patterns. This suggests that the pressure head at the bedrocksoil interface may be controlled in part by factors other than only
the depth of the overlying soil mantle.
Small depressions in the bedrock surface (concave or planar
areas) are known to impact the storage, distribution, and movement of water. For example, Liang and Uchida [23] contend that
subsurface saturation is governed by three variables, i.e., slope gradient, contributing area, and topographic wetness index. We tested
this hypothesis by characterizing the roughness of the bedrock surface using the index proposed by [56]. The bottom panel of Fig. 14
summarizes the results of our analysis. The plots display theroughness of the bedrock topography as a function of the pressure
head at the soil-bedrock interface. In general, the larger the valueof the bedrock roughness, the larger the variations in bedrock235

topography of neighboring cells or nodes. Per definition, a flat surface will have a roughness index of zero. The four scatter plots inthe bottom panel demonstrate the presence of a negative relationship between bedrock roughness and the pressure head at the soilbedrock interface. The larger the variations in bedrock topography,
the lower the pressure head values at the soil-bedrock interface.
The strength of this relationship does not change much in time,with correlation coefficients that are very similar for both samples(i.e., S11 and S13 ). We conclude that the pressure head at the soilbedrock interface is moderately affected by the shape of the underlying bedrock surface. The pressure head distribution of the soil
mantle hence is not only controlled by regolith thickness but also
by small-scale variability of the bedrock topography which affects
water flow and redistribution.
We next analyzed how the soil hydraulic parameters
p ¼ fhs ; hr ; a; n; K s g control slope stability. Fig. 15 presents two-Fig. 15. Two-dimensional scatter plots of the relationship between the safety factor (SF) and the hydraulic parameters hs ; hr ; a; n, and K s . Each panel also shows values ofPearson’s correlation coefficient for case studies 2 (red circles) and 3 (blue triangles). (For interpretation of the references to colour in this figure legend, the reader is referredto the web version of this article.)

236

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

dimensional scatter plots of relationships between p and the SF forcase studies 2 (sample S12 , red circles) and 3 (sample S13 , blue triangles). Note that these case studies differ only in terms of the bedrock surface used in the simulations (Table 5). Correlation
coefficients for case studies 2 (q2 , sample S12 ) and 3 (q3 , sample
S13 ) are indicated in each plot. As shown, at t ¼ 5 days, the SF
was not affected by any of the soil hydraulic parameters, certainly
not for case study 2, thus suggesting a lack of initial control of thehydraulic properties on slope stability. This finding is consistent
with the results of Fig. 10b. The third case study (S13 ) also showeda negligible correlation between the soil hydraulic parameters andthe SF. Since both samples use the same hydraulic parameters of
the soil mantle, the scatter of case study 3 is due to the use of a different bedrock surface. The scatter is persistent, and visible in thesnapshots at all four times. As illustrated in Fig. 15, the soil hydraulic parameters, hr and K s are strongly correlated with SF in casestudy 2 and moderately correlated in the case study 3. The negative
feedback between hr and K s was previously listed in Table 2. It
seems that the SF increases for larger values of the saturated
hydraulic conductivity, K s , of the soil mantle. The results in Sec-10case 1CVSF (%)

8case 2case 3

6
4
2

0
0

100

200

300

400

500

time (h)Fig. 16. The coefficient of variation of the safety factor (CVSF ) for the three different
case studies.

tion 5.2 indicated that the less permeable soil mantle induces
zones of positive pore pressures at the more steeper zones of thehillslope, thereby reducing the SF. On the other hand, the more
permeable soil mantle could increase drainage to the soilbedrock interface, and thus cause water to eventually move laterally to the drainage channel. This observation is in agreement with
Lanni et al. [55]. Hence, the planar (saturated) regions of the drainage channel are less susceptible to slope failure, thus leading to ahigher SF. This result also agrees with the analysis presented inFig. 11.
Based on the findings of Figs. 14 and 15, we suggest that theuncertainty in the SF during storm events is determined by complex nonlinear interactions between regolith depth, bedrock surface topography and the permeability of the soil mantle. We next
address the combined impact of bedrock depth and soil hydraulicparameter uncertainty on the probability of failure of the hillslope.
5.4. Probability of failure
Fig. 16 provides the coefficient of variation of the safety factor(CVSF ) with time for the three case studies. The coefficient of variation is a common and convenient statistical metric because it isdimensionless. The metric is given by the quotient of the standard
deviation and the corresponding mean of the ensemble of safety
factors, lSF . As expected, our results agree well with the confidence
intervals for the SF (Fig. 10). During the initial stages of the simulation, CVSF for case study 1 (black line in Fig. 16) was dominant,although CVSF values appeared rather small (around 1%). CVSF values for case study 3 (blue line) were in excellent agreement with
their values for case study 1, whereas the second case study (redline) produced negligible values. At t ¼ 120 h, CVSF for case study
3 quickly exceeded the other cases, and remained always higher
until the end of simulation. These results confirm the importanceof the combined effect of bedrock depth and soil hydraulic parameter uncertainty on probabilistic slope stability analyses. If the various sources of uncertainty are not considered simultaneously,
non-conservative CVSF may results.
We conclude this section with Fig. 17, which displays histograms of the SF for case studies (a) 1, (b) 2, and (c) 3 derived from
our Monte Carlo realizations at t ¼ 19 days, a few days before theend of the 22-day simulation period. The mean SF is separately

Fig. 17. Histograms of the marginal distribution of the SF at t ¼ 19 days for case studies (a) 1, (b) 2, and (c) 3. The mean value of the SF is indicated separately in eachhistogram by the blue cross symbol. The dashed gray lines delineate the boundaries between stability and instability (SF of one). The bottom panels present the corresponding
quantile-quantile plots for each case study. These plots demonstrate that the SF follows an approximately log-normal distribution. The reference one-to-one lines are shownin red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

237

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

indicated in each graph by the blue cross symbol. The verticaldashed lines in each histogram delineate stability and instability(SF values of unity). As expected, the histogram of the third casestudy exhibited the most dispersion in response to the combined
treatment of bedrock depth and soil hydraulic parameter uncertainty (Figs. 10 and 16).
The shape of the SF distributions (Fig. 17) and the corresponding coefficient of variation (Fig. 16) provide a means for evaluating
an alternative measure of safety: the probability of failure, P f . For agiven SF distribution, the probability, P, of hillslope failure isreported as Pf ¼ PðSF 6 1Þ, i.e., P f 2 ½0; 1 quantifies the probability
that the soil mass will fail (colapse). By assuming the distributions
shown in Fig. 17 for the SF, Pf can be estimated from the lognormal reliability index, bLN , so that Pf ¼ 1  UðbLN Þ, where U
denotes the cumulative distribution function of the standard normal variate [57]. bLN can be calculated as follows [14]:



ln
bLN ¼

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

lSF = 1 þ CV2SF

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
lnð1 þ CV2SF Þ:ð6ÞFrom Eq. (6), one can note that for a fixed lSF , the larger the value of
CVSF , the lower the value of bLN , and the larger the probability of
failure, P f .
Table 6 reports P f -values for the three case studies at threepoints in time (i.e., 12; 16, and 19 days). While the mean SF (thirdcolumn) was the same for all case studies, the values of CVSF
(fourth column) differed noticeably, which caused the values ofP f in the last column to vary between the three case studies. At
t ¼ 12 days, P f increased an order of magnitude from 0:2  105
to 0:4  104 for case study 1, to 0:2  104 for case study 2, andthen to 0:2  103 for case study 3. At t ¼ 16 days, the lSF values
remained rather constant and the already negligible probabilityof failure of the hillslope has further decreased for each case study,
due to smaller values of CVSF . Finally, at t ¼ 19 days, the listed Pf ’sdemonstrate an increased instability of the hillslope, primarily inresponse to smaller lSF values. Note also that Table 6 indicatedthe largest probability of failure for case study 3. Indeed, the valuesof Pf of this study always exceeded their counterparts of the othertwo studies. This discrepancy is orders of magnitude at t ¼ 12 and16 days, but smaller at 19 days. We thus conclude that the joint
impact of bedrock depth and soil hydraulic uncertainty impacts
significantly the probability of failure.
6. Discussion
The integrated hillslope water flow and numerical limit analysis
methodology used herein enables explicit treatment of uncertainty
in key inputs of slope stability studies of unsaturated soils. We

used this Monte Carlo approach and evaluated in this paper the
impact of bedrock depth uncertainty (case study 1), soil hydraulicuncertainty (case study 2), and their combined uncertainties (casestudy 3) on the stability of a variably-saturated hillslope in Rio deJaneiro, Brazil. We investigated in detail the marginal and joint
impact of these uncertainties on distributed pressure head dynamics of the soil mantle, the mean and confidence intervals of the SF,
failure mechanisms, and the probability of slope failure. Joint treatment of bedrock depth and hydraulic parameter uncertainty of thesoil mantle allowed us to investigate how the balance between
regolith depth, bedrock topography and parameters of thevariably-saturated flow model can explain complex interactions
that drive pore pressures and failure zones. One of the major
strengths of our proposed method is that, instead of merely assuming a known and deterministic bedrock surface, one can propagate
its uncertainty in geotechnical analysis.
The uncertainty of the bedrock surface was characterized using
Bayesian analyses coupled with MCMC simulations. The results
allowed us to determine more appropriate confidence limits for
the SF. Moreover, having a more accurate bedrock depth map also
makes it easier to characterize the impact of soil hydraulic parameters on rainfall-induced landslides. The SF confidence bounds forcase study 1 (Fig. 10a) indicate a moderate sensitivity of different
stratigraphic profiles when the soil-bedrock interface is predominantly unsaturated. This result is important since past studies indicated that slope failure may occur while some suction is stillpresent in the soil mass [9]. Our simulations revealed that shallow
regolith zones may increase pressure heads quickly because thetravel distance for the infiltration process is short. However, highpositive pore water pressures will only reduce the SF if they
develop at bedrock depressions or planar surfaces. Therefore, convex bedrock shapes at the sideslope were effective to drainage andpore pressure dissipation.
The use of a multivariate normal distribution for the van
Genuchten-Mualem (VGM) parameters restricts the soil mantle’s
hydraulic characterization to realistic retention and unsaturated
permeability functions as the Monte Carlo sampling honors, via
the covariance matrix, the multidimensional correlation structureof the soil hydraulic parameters. This finding is not new and has
been previously reported [18,22,54]. In practice, however, thisapproach has not been explored adequately in stability analyses.
Our findings suggest that the multivariate normal distributions of
the soil hydraulic parameters sampled in a Monte Carlo approach
can lead to significant uncertainty in the SF when positive pressure
heads start to develop in the soil mantle. Samples from our priordistribution with higher K s values generally led to more drainage
and contributed to the upper bound of the SF uncertainty. Samples
with a smaller K s produced patches of pressure heads beyond thedrainage channel, at high-gradient bedrock surface zones, thus
reducing the SF.

Table 6
Probabilities of hillslope failure, P f , for three different measurement times (first column) and three case studies (second column) analyzed in this paper. We separately report thecorresponding values of lSF and CVSF . The listed values of Pf assume a log-normal distribution of the SF at each time.
t (days)

Case

lSF

CVSF (%)

Pf

12

1

1.29

5.50

0:2  105

2

1.28

6.14

0:4  104

3

1.29

7.28

0:2  103

1
2
3

1.33
1.31
1.32

4.12
5.28
5.90

0
0

1

1.19

4.93

0:2  103

2

1.18

4.65

0:2  103

3

1.19

7.23

8:8  103

16

19

0:1  105

238

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

In case study 3, the combined effect of uncertainties showed an
overall increase in the variability of the SF as compared to casestudies 1 and 2 (Fig. 16). The large coefficient of variation at theend of the simulations indicated that the joint uncertainty of bedrock depth and the soil hydraulic parameters produced more conservative values of CVSF . While the mean SF was about the same forthe three case studies (Table 6), the probability of failure for casestudy 3 was always higher than that for the other two studies. Thiscan be attributed to the interactions between regolith depth, bedrock topography and soil hydraulic parameters as detailed in Section 5.3. These three factors determine patches of saturation that
may increase pore pressure and, in turn, influence the SF values
derived from NLA. Ignoring uncertainty in bedrock depth and thehydraulic properties of the soil mantle would hence result in overly
optimistic, or non-conservative, estimates of the probability of
failure.
Despite our efforts to address two important sources of uncertainty on probabilistic analyses of rainfall-induced landslides,
some other sources of uncertainty have been conveniently ignored.
For instance, some authors have pointed out that landslides arecontrolled by critical rainfall intensity and duration [7], anisotropic
hydraulic conductivity [7,58], slope geometry [2–6], initial saturation [3,4,7], recharge from bedrock [3,58], and strength parameters[11] on landslides. These studies, however, have been restricted tosimplified geometries. A large computational burden would be
required to analyze all of these variables within a transient
three-dimensional framework as proposed herein. In this paper
we assumed that soil properties are homogeneous and isotropic
and that bedrock is impermeable. This simplification is particularly
convenient for large hillslopes for which characterization of subsurface heterogeneity and anisotropy is quite difficult, if not
impossible. However, flow through the soil-bedrock interface[3,23,27], rainfall intensity, duration and frequency [7], as well soil
macro-permeability [9], if available, should be considered. As
shown in this study, the treatment of uncertainty in the bedrocksurface is paramount in stochastic modeling of landslides, butother sources of error could be equally important.

7. Summary and conclusions
In this paper we presented a stochastic framework to investigate the role of uncertainty in the bedrock surface and hydraulicproperties on the geotechnical stability of a variably-saturated,soil-mantled hillslope. This framework couples Monte Carlo simulation of a three-dimensional finite element water flow model ofthe bedrock-soil domain with numerical limit analysis to computeconfidence intervals of the hillslope’s safety factor (SF) duringsome historical period of rainfall events.
Three case studies were used to illustrate our methodology witha regolith depth model of a natural unsaturated hillslope in Rio deJaneiro, Brazil. Our depth-to-bedrock model predicted accurately
distributed measurements of the bedrock depth, thus providing ameans to quantify, via Bayesian analysis, bedrock topographic
uncertainty. The first case study evaluates this uncertainty of the
bedrock depth on slope stability using a fixed hydraulic characterization of our soil-hillslope domain. Our results demonstrate thatthe confidence intervals of the SF are affected by bedrock depth
uncertainty, including a moderate variability in the SF during periods in which the soil-bedrock interface is unsaturated. The secondcase study focused on uncertainty in the soil hydraulic parameters,
while keeping the posterior mean bedrock surface derived with
Bayesian analysis as a fixed slope geometry. The SF showed negligible variability when the soil mantle was mostly unsaturated, but
the confidence intervals increased dramatically with saturation. In
the third case study, the combined effects of both sources of uncer-

tainties produced broader confidence bounds for the SF. This finding suggests that our Monte Carlo approach improves
characterization of hillslope stability and the reliability of the output of hydro-mechanical models, such as safety factors. The larger
probability of failure obtained in the third case study further supports this conclusion.
Altogether, the results of our case studies advocate the need foran explicit treatment of bedrock depth uncertainty in hillslope stability analysis studies. Without this treatment, reliability analysis
can lead to overly optimistic results by underestimating considerably the actual probability of hillslope failure. Our paper showed
how to incorporate bedrock depth uncertainty in geotechnicalanalyses using a single integrated framework of Monte Carlo simulation, finite element soil water flow modeling and numericallimit analysis.
Conflict of interest
There is no conflict of interest in this study.
Acknowledgements
We appreciated the useful comments of the two referees, which
helped us to further enhance this manuscript. The authors
acknowledge financial support from the Brazilian National Council
for Scientific and Technological Development, CNPq. The firstauthor is also grateful to the Highway Department of Espírito Santo
State in Brazil (DER-ES) for their help during this research. We
would like to thank Alonso Carbono for his help with the compilation of the hillslope water flow model. The MATLAB code of thestochastic slope stability framework presented herein is available
upon request from the first author. This includes access to the different codes and scripts used in the three case studies of this paper.
Appendix A. Numerical limit analysis
We use numerical limit analysis (NLA) to calculate the SF of thevariably-saturated hillslope. This Appendix A provides a brief summary of the NLA method as described by Camargo et al. [38]. We
refer interested readers to this publication for further details.
A.1. Yield criterion
In this study we used the Drucker-Prager failure criterion[59,60] as the yield condition. The Drucker-Prager model uses
two stress tensor invariants and two material properties. Because
we assume plane strain conditions, the material properties canbe directly related to the strength parameters of Mohr-Coulomb[38], which are the effective cohesion, c0 , and the friction angle,/0 . The shear strength in variably-saturated soils also depends on
the soil suction, which affects c0 . We therefore adopt the extended
Mohr-Coulomb criterion [61], and make use of the effective stressprinciple to delineate two independent stress state variables,
namely the net normal stress, r  ua , and the suction, ua  uw asfollows

s ¼ c0 þ ðua  uw Þ tan /b þ ðr  ua Þ tan /0¼ cap þ ðr  ua Þ tan /0 ;ðA:1Þ

where ua (MT2L1) and uw (MT2L1) signify the air pressure in thevoids and the pore water pressure, respectively, /b ( ) is a proxy(called angle) for the rate of change in matric suction, while r(MT2L1) represents the total stress. The first two terms at theright-hand-side of Eq. (A.1) measure the combined effects of soil
cohesion and soil suction on the shear strength. The sum of these

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

two terms is equivalent to the apparent cohesion, cap . The parameters cap and /0 are thus related to the Drucker-Prager soil material.
A.2. Optimization problem
The collapse factor, k, was evaluated by imposing on the statiTcally admissible stress field, ri ¼ rxx ; ryy ; rzz ; sxy ; syz ; sxz
ofdimension 1  6 (where T denotes transpose), constraints of equilibrium throughout the hillslope domain and the yield criterion of
Drucker-Prager. This is equivalent to the following deterministic
and convex optimization problem

maximize ksubjet toequilibrium conditioni ¼ 1; . . . ; Ne

yield criterion;
ðA:2Þ

where G denotes the 3Nn  6Ne equilibrium matrix of the N n nodes
and Ne elements of the finite-element mesh of the hillslope, f is the
3N n  1 vector of body forces associated with the weight of the soilmaterial at each element of the mesh, and gðri Þ signifies the yield
condition.
The convexity of the collapse load optimization problem of Eq.
(A.2) can be exploited using second-order cone programming[49,51]. In this approach, the 6N e  1 vector r is transformed linearly into an auxiliary vector, q, composed of stresses in a conic
quadratic space [51]q ¼ D  r þ d;

ðA:3Þ

where the matrix D and the vector d store the coefficients of the
Drucker-Prager yield criterion. For each element of the mesh, Di isi

a 6  6 matrix, while vector d has dimensions 6  1. In a conic
quadratic form, the collapse factor optimization problem now reads[38]

maximize ksubjet to

A.4. Safety factor (SF)A shear strength reduction technique was used to compute theSF of the hillslope during our 22-day rainfall record. This method
decreases iteratively the values of the strength parameters, cap
and /0 , via a reduction factor until the collapse factor, k, approaches
zero. We invoked Newton’s method to solve for the root of the
reduction factor at many different simulation times [38]. This
resulted in a 22-day SF record of the variably-saturated hillslope.
To maximize search efficiency of Newton’s method, we initialized
the reduction factor with the root (its optimized solution) of theprevious simulation time.

collapse factorG  r ¼ kf
gðri Þ 6 0;239

G  ðD1  q  D1  dÞ ¼ f þ kf h ;qi 2 K; i ¼ 1; . . . ; Ne ;ðA:4Þ

where the 3Nn  1 vector f h stores the external nodal forces acting
on the soil mass of the hillslope, (net effect of horizontal and fictitious body forces), and where each auxiliary (elementary) vector,
qi , is restricted to the conic quadratic set, K. We solved numericallyfor the collapse factor in Eq. (A.4) using the MOSEK optimization
toolbox in MATLAB [62]. This toolbox gives users access to a large
array of powerful optimization methods for solving a wide rangeof minimization problems, including second-order conic quadratic
problems. A particularly strong point of MOSEK is its state-of-theart interior-point optimizer, which allows for rapid solutions ofhigh-dimensional continuous linear, quadratic and conic optimization problems. Note that it is numerically convenient to use f h in theequilibrium condition as this simplifies considerably our calculations of the SF. Please refer to [38] for a detailed description of
the NLA method used in our study.
A.3. Boundary conditions
Once the collapse load optimization problem has been definedmathematically, one still has to define the boundary conditionsof the hillslope being studied. The boundary conditions appear in
the equilibrium matrix, G, and the vectors f and f h for the bedrock,lateral, downslope, and upslope faces, respectively. As detailed inFig. 2b, we used zero velocities in the u; v and w directions.

Appendix B. Characterization of bedrock topographic
uncertainty
Depth to bedrock is a key variable in slope stability studies but
difficult to measure directly in practice as hillslope interiors arepoorly accessible. Instead, we generated many different plausible
relizations of the bedrock surface using the calibrated depth-tobedrock (DTB) model of [33]. The approach we describe below
accounts explicitly for DTB parameters and model uncertainty,
and delivers on the promise of [33] that (page 3085) ”The posterior
prediction uncertainty of the DTB model can be propagated forward
through hydromechanical models to derive probabilistic estimates offactors of safety.” Indeed, one can treat depth to bedrock as a probabilistic variable and use Monte Carlo simulation with many different bedrock topographies to quantify uncertainty in the estimatesof slope stability.
To clarify our approach, we write the DTB model, F ðÞ as a simple nonlinear regression function

H ¼ F ðdÞðB:1Þwhich returns a K  L matrix of simulated bedrock depths, H, forparameter values, d ¼ fd1 ; . . . ; dd g. The K rows and L columns signify
different dimensions of H, i.e., the x and y coordinates of the hillslope. We can now estimate the values of the d parameters of theregression function via Bayesian inference with the DREAM algorithm [36,37]. This results in a large sample of M parameter vectors,
also called posterior solutions, whose simulated bedrock surfaces
best explain N distributed measurements of the bedrock depth.
These posterior parameter vectors are stored in a M  d matrix,
D ¼ fd1 ; . . . ; dM g, and their corresponding matrices of the bedrock , of size K  L  M, and thus
depth in a three-dimensional array, hh now stores the M differh ¼ fH1 ; . . . ; HM g. The third dimension of ent posterior maps of the bedrock surface. The expectation, EðhÞ, ofthese M realizations equates to the posterior mean bedrock depth
map simulated by the DTB model. Confidence intervals of this mean
map are easily constructed from the different posterior realizations,by sorting from low to high and for each entry of the model output
matrix the M simulated bedrock depths. The indexes 0:025M and
0:975M of each sorted vector constitute the 95% prediction rangesof the DTB model due to parameter uncertainty. One can derive
the total prediction uncertainty of the DTB model by adding to eachentry of 
h a random draw from the model error distribution. This
procedure is outlined in the main text. Fig. B.1 displays the mean
bedrock depth (red line) and its associated 95% predictions intervals(blue lines) for cross section AA’ of the hillslope using M ¼ 100(solid lines), M ¼ 200 (dashed lines) and M ¼ 1000 (dotted lines)
posterior realizations. The surface topography is indicated separately by the solid black line. The posterior mean bedrock surface
and its confidence intervals appear only minimally affected by the
choice of M. To minimize the computational requirements of theMonte Carlo framework, we therefore did set M at 100.

240

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241

M =100
M =200
M =1000
M =100
M =200
M =1000

vertical distance (m)

60

45

30

↑

mean

↑

surface↓

15095% total uncertainty

20

40

60

80

100

120

horizontal distance (m)Fig. B.1. Two-dimensional plot of the posterior mean bedrock topography (red line)
and the 95% confidence intervals due to DTB model and parameter uncertainty(blue lines) of transect AA’ of the hillslope using M ¼ 100 (solid lines), M ¼ 200(dashed lines) and M ¼ 1000 (dotted lines) posterior solutions. The solid black line
depicts the surface topography. (For interpretation of the references to colour inthis figure legend, the reader is referred to the web version of this article.)

References[1] Rempe DM, Dietrich WE. A bottom-up control on fresh-bedrock topography
under landscapes. Proc Natl Acad Sci 2014;111(18):6576–81. http://dx.doi.org/
10.1073/pnas.1404763111.
[2] Zhang Y, Chen G, Zheng L, Li Y, Zhuang X. Effects of geometries on threedimensional slope stability. Can Geotech J 2013;50(3):233–49. http://dx.doi.
org/10.1139/cgj-2012-0279.
[3] Vargas Jr EA, Velloso RC, de Campos TMP, Filho LMC. Saturated-unsaturated
analysis of water flow in slopes of Rio de Janeiro, Brazil. Comput Geotech
1990;10(3):247–61. http://dx.doi.org/10.1016/0266-352X(90)90038-W.
[4] Mukhlisin M, Taha M, Kosugi K. Numerical analysis of effective soil porosity
and soil thickness effects on slope stability at a hillslope of weathered granitic
soil formation. Geosci J 2008;12(4):401–10. http://dx.doi.org/10.1007/s12303008-0039-0.
[5] Catani F, Segoni S, Falorni G. An empirical geomorphology-based approach tothe spatial prediction of soil thickness at catchment scale. Water Resour Res
2010;46(5). http://dx.doi.org/10.1029/2008WR007450.
[6] Fan L, Lehmann P, Or D. Effects of soil spatial variability at the hillslope andcatchment scales on characteristics of rainfall-induced landslides. Water
Resour Res 2016:840–57. http://dx.doi.org/10.1002/2015WR017758.
[7] Ng CWW, Shi Q. A numerical investigation of the stability of unsaturated soil
slopes subjected to transient seepage. Comput Geotech 1998;22(1):1–28.
http://dx.doi.org/10.1016/S0266-352X(97)00036-0.
[8] Zhang LL, Zuo ZB, Ye GL, Jeng DS, Wang JH. Probabilistic parameter estimationand predictive uncertainty based on field measurements for unsaturated soil
slope.
ComputGeotech2013;48:72–81.
http://dx.doi.org/10.1016/
j.compgeo.2012.09.011.
[9] El-Ramly H, Morgenstern NR, Cruden DM. Probabilistic assessment of stabilityof a cut slope in residual soil. Géotechnique 2005;55(1):77–84. http://dx.doi.
org/10.1680/geot.2005.55.1.77.
[10] François B, Tacher L, Bonnard C, Laloui L, Triguero V. Numerical modelling of
the hydrogeological and geomechanical behaviour of a large slope movement:
the Triesenberg landslide (Liechtenstein). Can Geotech J 2007;44(7):840–57.
http://dx.doi.org/10.1139/t07-028.
[11] Griffiths DV, Huang J, Fenton GA. Influence of spatial variability on slope
reliability using 2-D random fields. J Geotech Geoenviron Eng 2009;135(10):1367–78. http://dx.doi.org/10.1061/(ASCE)GT.1943-5606.0000099.
[12] Cho SE, Lee SR. Instability of unsaturated soil slopes due to infiltration. Comput
Geotech 2001;28(3):185–208. http://dx.doi.org/10.1016/S0266-352X(00)00027-6.
[13] Tsaparas I, Rahardjo H, Toll DG, Leong EC. Controlling parameters for rainfallinduced landslides. Comput Geotech 2002;29(1):1–27. http://dx.doi.org/
10.1016/S0266-352X(01)00019-2.
[14] Duncan JM, Wright SG. Soil strength and slope stability. New York: John Wiley& Sons; 2005.
[15] Christian JT, Ladd CC, Baecher GB. Reliability applied to slope stability analysis.
J Geotech Eng-ASCE 1994;120(12):2180–207. http://dx.doi.org/10.1061/(ASCE)0733-9410.
[16] Kang F, Han S, Salgado R, Li J. System probabilistic stability analysis of soil
slopes using Gaussian process regression with Latin hypercube sampling.
ComputGeotech2015;63:13–25.
http://dx.doi.org/10.1016/
j.compgeo.2014.08.010.
[17] Malkawi AIH, Hassan WF, Abdulla FA. Uncertainty and reliability analysis
applied to slope stability. Struct Saf 2000;22:161–87. http://dx.doi.org/
10.1016/S0167-4730(00)00006-0.

[18] Phoon K, Santoso A, Quek S. Probabilistic analysis of soil-water characteristic
curves. J Geotech Geoenviron Eng 2010;136(3):445–55. http://dx.doi.org/
10.1061/(ASCE)GT.1943-5606.0000222.
[19] Wang Y, Cao Z, Au SK. Efficient Monte Carlo simulation of parameter
sensitivity in probabilistic slope stability analysis. Comput Geotech
2010;37:1015–22. http://dx.doi.org/10.1016/j.compgeo.2010.08.010.
[20] Zhang LL, Zhang LM, Tang WH. Rainfall-induced slope failure considering
variability of soil properties. Géotechnique 2005;55(2):183–8. http://dx.doi.
org/10.1680/geot.2005.55.2.183.
[21] Vrugt JA, Diks CGH, Gupta HV, Bouten W, Verstraten JM. Improved treatmentof uncertainty in hydrologic modeling: combining the strengths of global
optimization and data assimilation. Water Resour Res 2005;41(1). http://dx.
doi.org/10.1029/2004WR003059.
[22] Scharnagl B, Vrugt JA, Vereecken H, Herbst M. Inverse modelling of in situ soil
water dynamics: investigating the effect of different prior distributions of thesoil hydraulic parameters. Hydrol Earth Syst Sci 2011;15(10):3043–59. http://dx.doi.org/10.5194/hess-15-3043-2011.
[23] Liang WL, Uchida T. Effects of topography and soil depth on saturated-zone
dynamics in steep hillslopes explored using the three-dimensional Richards’
equation.
J
Hydrol2014;510:124–36.
http://dx.doi.org/10.1016/j.
jhydrol.2013.12.029.
[24] Reid M. Slope instability caused by small variations in hydraulic conductivity. J
Geotech Geoenviron Eng 1997;123(8):717–25. http://dx.doi.org/10.1061/(ASCE)1090-0241(1997)123:8(717).
[25] Cai F, Ugai K, Wakai A, Li Q. Effects of horizontal drains on slope stability underrainfall by three-dimensional finite element analysis. Comput Geotech
1998;23(4):255–75. http://dx.doi.org/10.1016/S0266-352X(98)00021-4.
[26] Dou H, Han T, Gong X, Zhang J. Probabilistic slope stability analysis considering
the variability of hydraulic conductivity under rainfall infiltrationredistribution conditions. Eng Geol 2014;183:1–13. http://dx.doi.org/
10.1016/j.enggeo.2014.09.005.
[27] Gerscovich DMS, Vargas Jr EA, de Campos TMP. On the evaluation of
unsaturated flow in a natural slope in Rio de Janeiro, Brazil. Eng Geol
2006;88(12):23–40. http://dx.doi.org/10.1016/j.enggeo.2006.07.008.
[28] Lanni C, Borga M, Rigon R, Tarolli P. Modelling shallow landslide susceptibility
by means of a subsurface flow path connectivity index and estimates of soil
depth spatial distribution. Hydrol Earth Syst Sci 2012;16(11):3959–71. http://dx.doi.org/10.5194/hess-16-3959-2012.
[29] Hopp L, McDonnell JJ. Connectivity at hillslope scale: identifying interactions
between storm size, bedrock permeability, slope angle and soil depth. J Hydrol
2009;376:378–91. http://dx.doi.org/10.1016/j.jhydrol.2009.07.047.
[30] Montgomery DR, Dietrich WE, Hefner JT. Piezometric response in shallow
bedrock at CB1: implications for runoff generation and shallow landsliding.
Water Resour Res 2002;38(12). http://dx.doi.org/10.10029/2001WR001429.
[31] Wang GH, Sassa K. Pore-pressure generation and movement of rainfallinduced landslides: effects of grain size and fine-particle content. Eng Geol
2003;69(1-2):109–25. http://dx.doi.org/10.1016/S0013-7952(02)00268-5.
[32] Kim MS, Onda Y, Kim JK, Kim SW. Effect of topography and soil
parameterisation representing soil thicknesses on shallow landslide
modelling. Quatern Int 2015;384:91–106. http://dx.doi.org/10.1016/j.
quaint.2015.03.057.
[33] Gomes GJC, Vrugt JA, Vargas Jr EA. Toward improved prediction of the bedrock
depth underneath hillslopes: Bayesian inference of the bottom-up control
hypothesis using high-resolution topographic data. Water Resour Res 2016:52.
http://dx.doi.org/10.1002/2015WR018147.
[34] Lee KT, Ho JY. Prediction of landslide occurence based on slope-instability
analysis and hydrological model simulation. J Hydrol 2009;375:489–97.
http://dx.doi.org/10.1016/j.jhydrol.2009.06.053.
[35] Ho JY, Lee KT, Chang TC, Wang ZY, Liao YH. Influences of spatial distribution ofsoil thickness on shallow landslide prediction. Eng Geol 2012;124:38–46.
http://dx.doi.org/10.1016/j.enggeo.2011.09.013.
[36] Vrugt JA, ter Braak CJF, Clark MP, Hyman JM, Robinson BA. Treatment of input
uncertainty in hydrologic modeling: doing hydrology backward with Markov
chain Monte Carlo simulation. Water Resour Res 2008;44(12). http://dx.doi.
org/10.1029/2007WR006720.
[37] Vrugt JA, ter Braak CJF, Diks CGH, Robinson BA, Hyman JM, Higdon D.
Accelerating Markov chain Monte Carlo simulation by differential evolution
with self-adaptive randomized subspace sampling. Int J Nonlin Sci Num
2009;10(3):273–90. http://dx.doi.org/10.1002/esp.3423.
[38] Camargo JT, Velloso RQ, Vargas Jr EA. Numerical limit analysis of threedimensional slope stability problems at catchment scale. Acta Geotech
2016:1–15. http://dx.doi.org/10.1007/s11440-016-0459-3.
[39] Schaap MG, Leij FJ, van Genuchten MTh. ROSETTA: a computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions.
J Hydrol 2001;251:163–76. http://dx.doi.org/10.1016/S0022-1694(01)004668.
[40] Matthies H, Strang G. The solution of nonlinear finite element equations. Int J
Num
Meth
Eng
1979;14(11):1613–26.
http://dx.doi.org/10.1002/
nme.1620141104.
[41] Bathe KJ, Cimento AP. Some practical procedures for the solution of nonlinear
finite element equations. Comput Methods Appl Mech Eng 1980;22(1):59–85.
http://dx.doi.org/10.1016/0045-7825(80)90051-1.
[42] Gomes GJC, Araujo JP, Vargas Jr EA, Fernandes NF. In-situ regolith depth
measurements by DPL tests. In: Proceedings of the XVIII Brazilian conferenceon soil mechanics and geotechnical engeneering, vol. 1. Belo Horizonte, Brazil;
2016. p. 1–8. doi:http://dx.doi.org/10.20906/CPS/GJ-01-0035.

G.J.C. Gomes et al. / Computers and Geotechnics 88 (2017) 222–241[43] Paniconi C, Aldama AA, Wood EF. Numerical evaluation of iterative andnoniterative methods for the solution of the nonlinear Richards equation.
Water
ResourRes
1991;27(6):1147–63.
http://dx.doi.org/10.1029/
91WR00334.
[44] Šimunek, J, Vogel, T, van Genuchten, MTh. The SWMS 2D code for simulatingwater flow and solute transport in two-dimensional variably saturated media– version 1.21. Research report 132, US Salinity Laboratory, Riverside,
California; 1994.
[45] Anderson MP, Woessner WW, Hunt RJ. Applied groundwater modeling. San
Diego: Academic Press; 2015.
[46] Rahman MM, Lu M, Kyi KH. Seasonality of hydrological model spin-up time: acase study using the Xinanjiang model. Hydrol Earth Syst Sci Discuss
2016:1–22. http://dx.doi.org/10.5194/hess-2016-316.
[47] van Genuchten MTh. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci Soc Am J 1980;44(5):892–8. http://dx.doi.org/10.2136/sssaj1980.03615995004400050002x.
[48] Mualem Y. A new model for predicting the hydraulic conductivity ofunsaturated porous media. Water Resour Res 1976;12(3):513–22. http://dx.
doi.org/10.1029/WR012i003p00513.
[49] Sloan S. Geotechnical stability analysis. Géotechnique 2013;63(7):531–72.
http://dx.doi.org/10.1680/geot.12.RL.001.
[50] Durand AF, Vargas Jr EA, Vaz LE. Applications of numerical limit analysis (NLA)
to stability problems of rock and soil masses. Int J Rock Mech Min Sci 2006;43(3):408–25. http://dx.doi.org/10.1016/j.ijrmms.2005.07.010.
[51] Makrodimopoulos A, Martin CM. Lower bound limit analysis of cohesivefrictional materials using second-order cone programming. Int J Num Meth
Eng 2006;66(4):604–34. http://dx.doi.org/10.1002/nme.1567.

241[52] Li AJ, Merifield RS, Lyamin AV. Three-dimensional stability charts for slopes
based on limit analysis methods. Can Geotech J 2010;47(12):1316–34. http://dx.doi.org/10.1139/T10-030.
[53] Vieira BC, Fernandes NF. Landslides in Rio de Janeiro: the role played byvariations in soil hydraulic conductivity. Hydrol Process 2004;18(4):791–805.
http://dx.doi.org/10.1002/hyp.1363.
[54] Vrugt JA, Bouten W, Gupta HV, Hopmans JW. Toward improved identifiabilityof soil hydraulic parameters: on the selection of a suitable parametric model.
Vadose Zone J 2003;2(1):98–113. http://dx.doi.org/10.2136/vzj2003.9800.
[55] Lanni C, McDonnell J, Hopp L, Rigon R. Simulated effect of soil depth andbedrock topography on near-surface hydrologic response and slope stability.
Earth Surf Proc Land 2013;38(2):146–59. http://dx.doi.org/10.1002/esp.3267.
[56] Riley SJ, DeGloria SD, Elliot R. A terrain ruggedness index that quantifies
topographic heterogeneity. Intermt J Sci 1999;5(1–4):23–7.
[57] Baecher GB, Christian JT. Reliability and statistics in geotechnicalengineering. England: John Wiley & Sons Ltd; 2003.
[58] Leach B, Herbert R. The genesis of a numerical model for the study of thehydrogeology of a steep hillside in Hong Kong. Q J Eng Geol Hydroge 1982;15(3):243–59. http://dx.doi.org/10.1144/GSL.QJEG.1982.015.03.09.
[59] Drucker DC, Prager W, Greenberg HJ. Extended limit design theorems forcontinuous media. Quart Appl Math 1952;9:381–9.
[60] Drucker DC, Prager W. Soil mechanics and plastic analysis for limit design.
Quart Appl Math 1952;10:157–65.
[61] Fredlund DG, Rahardjo H. Soil mechanics for unsaturated soils. New York: John
Wiley & Sons; 1993.
[62] MOSEK ApS. The MOSEK optimization toolbox for MATLAB manual. Version
7.1 (Revision 28); 2015. <http://docs.mosek.com/7.1/toolbox/index.html>.