s00603 023 03694 5

Rock Mechanics and Rock Engineering (2024) 57:2421–2443
https://doi.org/10.1007/s00603-023-03694-5
ORIGINAL PAPER
A Computationally Efficient Method to Determine the Probability

of Rainfall‑Triggered Cut Slope Failure Accounting for Upslope
Hydrological Conditions
Ellen Robson1,2 · David Milledge2 · Stefano Utili2 · Giuseppe Dattola3
Received: 15 May 2023 / Accepted: 20 November 2023 / Published online: 24 December 2023
© The Author(s) 2023
Abstract
We present a new computationally efficient methodology to estimate the probability of rainfall-induced slope failure based
on mechanical probabilistic slope stability analyses coupled with a hydrogeological model of the upslope area. The model
accounts for: (1) uncertainty of geotechnical and hydrogeological parameters; (2) rainfall precipitation recorded over a period
of time; and (3) the effect of upslope topography. The methodology provides two key outputs: (1) time-varying conditional
probability of slope failure; and (2) an estimate of the absolute frequency of slope failure over any time period of interest. The
methodology consists of the following steps: first, characterising the uncertainty of the slope geomaterial strength parameters;
second, performing limit equilibrium method stability analyses for the realisations of the geomaterial strength parameters
required to calculate the slope probability of failure by a Monte Carlo Simulation. The stability analyses are performed for
various phreatic surface heights. These phreatic surfaces are then matched to a phreatic surface time series obtained from
the 1D Hillslope-Storage Boussinesq model run for the upslope area to generate Factor of Safety (FoS) time series. A time-
varying conditional probability of failure and an absolute frequency of slope failure can then be estimated from these FoS time
series. We demonstrate this methodology on a road slope cutting in Nepal where geotechnical tests are not readily conducted.
We believe this methodology improves the reliability of slope safety estimates where site investigation is not possible. Also,
the methodology enables practitioners to avoid making unrealistic assumptions on the hydrological input. Finally, we find
that the time-varying failure probability shows marked variations over time as a result of the monsoon wet–dry weather.
Highlights
• Probabilistic slope stability analyses are coupled with a hydrogeological hillslope model to estimate the probability and
frequency of rainfall-induced slope failure.
• The model accounts for the uncertainty about rainfall using a time-dependent method, and for uncertainty relative to the
geomaterial properties.
• The model is tested on a road cut slope in Nepal (mountainous area subject to a monsoon season) finding that the cut slope
will fail every other year.
• Time-varying failure probability shows marked variations over time as a result of the monsoon wet–dry weather.
• The findings indicate that it is important to use a time-dependent system to represent rainfall variability for slope failure
probability analysis.
* Ellen Robson
ellen.robson@durham.ac.uk
1
Institute of Hazard, Risk and Resilience, Durham University,
Durham, UK
2
School of Engineering, Newcastle University,
Newcastle upon Tyne, UK
3
Dipartimento di Scienze dell’Ambiente e della Terra,
Universitá degli Studi di Milano, Milan, Italy
13
Vol.:(0123456789)
2422 E. Robson et al.
Keywords Slope stability · Limit equilibrium method · Monte Carlo simulation · Probabilistic stability analyses · Hillslope-
Storage Boussinesq · Nepal
1 Introduction average slope conditions over a region. Statistical models

also have restrictions in terms of planning slope stabilisa-
When roads and railway lines are constructed, they often tion measures as they are necessarily based on a database of
require slope cuttings (also known as cut slopes), especially historic failures. These failures encompass not only a wide
in hilly topography. When cuttings fail, there can be huge range of site-specific conditions (e.g. geometry, hydrology
economic and social consequences, especially when debris and material properties) but also of stabilisation measures.
collides with infrastructure or vehicles. Cutting failures are Mechanistic models have been widely used to account
widespread along roads in low- to lower–middle-income for uncertainty in geomaterial mechanical properties and
countries (LIC/LMIC) where natural susceptibility to slope their spatial variability, in addition to uncertainty in geo-
failure from heavy rain and hilly topography is exacerbated metric parameters, e.g. position of contact surfaces in a
by the inadequate design of cuttings and their stabilisation multilayered slope. A large number of mechanistic mod-
measures. To carry out planning, design and costing of sta- els have been developed to predict landslides across land-
bilisation measures along transportation infrastructure, it is scapes by coupling hydrological models of varying levels
useful for practitioners to have a measure for the stability of of complexity (e.g. SHALSTAB, Montgomery and Dietrich
a cutting (Liu and Wang 2021; Corominas and Moya 2008). (1994), TRIGRS, Baum et al. (2008), GEOtop-FS, Formetta
Factor of Safety (FoS) values are approximate indicators of et al. (2014)) with the infinite slope model. Even those who
safety based on a single best estimate of the slope conditions treat the slope stability component more completely make
with the FoS expressing a margin of safety, the value of the assumption that failure occurs within the soil or at the
which is highly dependent on the problem considered. For soil–bedrock boundary (Hess et al. 2017; van Zadelhoff et al.
example, most geotechnical standards require FoS = 1.3 for 2022; Bellugi et al. 2015). These models are not appropri-
a slope whilst they typically require FoS = 4 against piping ate for cutting failures because of their representation of the
failure. One may think that the safety margin against piping failure surface. Mechanistic models that better capture the
is much higher but actually the differences in the threshold failure surface currently do not deal with complex upslope
values largely stem from the model and parameter uncertain- topography extending to the watershed for that slope (e.g.
ties. Also, two slopes exhibiting the same FoS can have a Manning et al. (2008); Rouainia et al. (2009); Holcombe
vastly different probability of failure. Instead, probability of et al. (2012)). Therefore, we suggest that there is a gap in
failure (P(f)) and its complementary probability of survival literature: complete models that represent upslope hydro-
provide a better metric to estimate a slope’s (lack of) safety logical conditions and stability of cuttings both in a high
for practitioners and public authorities since the meaning of level of detail.
such metrics is well understood outside the boundaries of Probability of failure given some trigger event can be esti-
civil engineering. mated using mechanistic models by conducting a probabil-
Landslide hazard and risk assessment studies often esti- istic stability analysis (Dou et al. 2014; Zhang et al. 2014).
mate an annual probability of slope failure, often in map Thus, the resultant probability of failure is always condi-
form (Corominas et al. 2014; Dai et al. 2002; Tang et al. tional, though this is not always explicitly recognised. The
2018). Cost assessment studies require a frequency within most common triggers of slope failure are pore water pres-
a time window. Annual probability can also be thought of sure and seismic shaking, which are driven by rainfall and
as a frequency; a probability set in a reference time frame earthquakes, respectively. These triggers are not binary, but
(Tang et al. 2018). These are typically determined either they have characteristics (e.g. pore pressure/seismic accel-
through statistical models which involve analysis of land- eration) and these characteristics both: (1) affect the prob-
slide databases (e.g. Lu et al. (2020); Guzzetti et al. (2008)) ability of failure (i.e. the probability of failure is conditional
or mechanistic models which involve slope stability analyses on the trigger characteristics); and (2) vary in space and time
(e.g. Lu et al. (2022a)). A drawback of statistical models is such that their future values at a given location are highly
that they require a large volume of reliable historical data uncertain. Therefore, estimating probability of failure within
which can be difficult (and costly) to attain, especially in the some time window (e.g. annual probability) for a given loca-
context of LIC/LMICs that are often data-poor (particularly tion should account for not only the probability of failure
in terms of data on smaller scale road cuttings). In addition, conditional on a particular trigger event (P(f|T)) but also the
statistical approaches often do not consider the local geology probability of that trigger event (or those trigger conditions)
and geomaterial conditions and, therefore, only represent the occurring (P(T)). Further, since multiple trigger events result
in non-zero failure probability (i.e. P(f|Ti)), occurrence
13
A Computationally Efficient Method to Determine the Probability of Rainfall‑Triggered Cut… 2423
probabilities are required for multiple trigger conditions assumed to be stochastic variables, with the soil properties
P(Ti), for example in the form of a probability distribution assumed to be lognormally distributed and the rainfall inten-
of trigger intensities (Mori et al. 2020; Lu et al. 2022a, b). sity is assumed to be an exponential variable (the correla-
Yet many studies do not account for variability in trigger tion between intensity and duration is taken into account).
intensities and only focus on the variability of the geomate- Mori et al. (2020) also determines the probability of slope
rial properties. In these cases, failure probabilities are really failure accounting for uncertainty in material properties (soil
conditional probabilities (i.e. P(f|T)) and it is essential that strength and permeability) and trigger variability (rainfall
they are reported as such, with the trigger conditions on intensity from an IDF curve) using a MCS framework. They
which they are conditional clearly reported, and that they use a smooth particle hydrodynamics stability model with
are not interpreted as true probabilities associated with some random field modelling to capture the spatial correlation of
time window (e.g. annual probabilities) from which esti- the soil properties. They show that variability in saturated
mates of hazard or risk can be made. Annual probability of hydraulic conductivity can have a very strong influence on
a rainfall-triggered slope failure is often assumed to be con- the probability of failure.
trolled by the most critical rainfall (the rainfall under which A key limitation of IDF curves is that there is often an
the failure probability of the slope is maximum amongst all assumption of constant rainfall over the entire duration
rainfall events in a year) (Lu et al. 2022a; Tang et al. 2018). of the storm which results in the shape of the storm time
The use of a critical rainfall to represent annual probability series being lost (the method used by Tang et al. (2018)
of failure raises additional concerns regarding its derivation: relaxes this assumption). In addition, studies that determine
(1) how do we know that this is the most critical rainfall (i.e. rainfall-induced slope failure probability using IDF curves
how do we know that it triggers max P(f|T))?; (2) how do assume that rainfall is time independent and, therefore, do
we decide what year to use (i.e. large rainfall events could not account for antecedent rainfall conditions. However, this
be a one in 10-year event, a one in 50-year event or a one in is rarely the case. In fact, many processes associated with
100-year event)? landslide triggering are rarely truly time independent: earth-
When the variability in trigger characteristics is quakes trigger aftershocks (e.g. Parsons (2002)), and damage
accounted for, this is often done assuming that trigger rock altering its material properties (e.g. Jones et al. (2021));
probability is time independent, e.g. modelling pore pres- rainfall events cluster over a range of timescales from days
sure using rainfall events drawn from an Intensity–Dura- to years (e.g. Wang et al. (2005)); and phreatic surfaces are
tion–Frequency (IDF) curve (Holcombe et al. 2012; Lu influenced by previous rainfall sometimes over weeks or
et al. 2022b). IDF curves are joint probability distributions months (e.g. (Iverson 2000)). Thus, an alternative to time-
of the intensity and duration of storm rainfall. Intensity and independent sampling becomes necessary where triggers
duration of rainfall are widely considered as the primary display consistent and/or significant time-dependence (e.g.
rainfall characteristics responsible for generating landslide- areas that experience dry and wet seasons). Some studies
triggering pore pressure distributions within slopes (Caine have sought to address this within an IDF-based framework
1980; Guzzetti et al. 2008). IDF curves are generally based by sampling different antecedent conditions (e.g. Frattini
on very long records of rainfall and have been used in both et al. (2009)), others generate synthetic time series sampled
statistical (e.g. Guzzetti et al. (2008)) and mechanistic (e.g. from IDF curves (e.g. Tang et al. (2018)). However, there
Tang et al. (2018)) rainfall-induced landslide models. Tang have been surprisingly few attempts to account for the time-
et al. (2018) determine a conditional slope failure probabil- dependent probability due to rainfall sequencing through
ity and annual failure probability (alongside a deterministic direct simulation of rainfall time series (see Jones et al.
FoS) of a partially saturated soil slope triggered by rainfall, (2021); Ozturk (2022) for earthquake-triggered equivalents).
using ‘random rainfall patterns’ (rainfall intensity over time) In accounting for rainfall variability, a hydrological model
simulated by a random cascade model characterised by an is required for the derivation of time-dependent pore water
IDF curve. They use a probabilistic framework (a Monte pressure conditions from the rainfall record. For cuttings,
Carlo Simulation, MCS) to determine the conditional prob- this has focussed on representing rainfall infiltration: De
ability of failure, varying the random rainfall patterns. They Leon and Garduño (2020) use a 2D Richard’s Equation (RE)
determine the annual failure probability by multiplying the solver to model rainfall infiltration on a 2D slope, they do
failure probability conditional on a set of rainfall intensities not take into account the whole slope domain; Tang et al.
(obtained from MCS) with the occurrence probabilities of (2018) perform rainfall infiltration analyses using 2D seep-
those rainfall intensities (obtained from an IDF curve). De age analysis characterised by the soil–water characteristic
Leon and Garduño (2020) account for uncertainty in geo- curve and permeability function curve; Mori et al. (2020)
material properties and rainfall in determining an annual simulate infiltration into the slope with rainfall intensity
failure probability for soil slopes under heavy rainfall. Soil specified through the Gumbel distribution (where the rate
density and strength, rainfall intensity and duration are of infiltration is dependent on permeability and on pressure
13
head differences of different soil layers). They use transient complex hillslopes. The continuity and Darcy equations of
seepage analysis following Darcy’s law. Holcombe et al. the Boussinesq equation are reformulated in terms of storage
(2012) are an exception in that they use a 1D RE solver to along the hillslope to derive the so-called Hillslope-Storage
model vertical infiltration but also represent groundwater Boussinesq (HSB) equation (see Eq. 1 in Sect. 2.2). By
flow below the phreatic surface with an explicit solution to introducing soil moisture storage, the 3D groundwater flow
the Darcy equation. A drawback of all these methods is that problem recast as a 1D flow problem for a 1D slope of vari-
they focus on the cutting and do not account for the influence able planform width with an inclined planar impermeable
of the surrounding topography on the groundwater regime, boundary. The HSB model can be applied to any slopes that
i.e. the effect of hydrogeology and antecedent rainfall on the can be considered as continuous media, such as soil slopes
phreatic surface. Current approaches to modelling rainfall and rock slopes with no dominant fracture orientation (as
infiltration and seepage in cuttings have many limitations: dominant fracture orientation would affect the flow motion).
(1) they are computationally expensive; (2) hydrological In these cases, Darcy’s equation and the continuity equation
parameters (e.g. hydraulic conductivity) and their variation are valid.
in space are uncertain and poorly constrained, but have a Further simplifications enable faster numerical or even
substantial influence on the model output; (3) lateral inputs analytical solutions (e.g. Troch et al. (2004); Hilberts et al.
(i.e from upslope) are rarely taken into account. Models that (2004); Talebi et al. (2008)) but these are unnecessary in
predict landslides across landscapes have used a wider range this case since the Troch et al. (2003) model is already suf-
of hydrological treatments, some only consider infiltration ficiently fast to run 1000 s of simulations within hours. Thus,
(e.g. Iverson (2000)), but many highlight the importance Troch et al. (2003) HSB model is used in this study to simu-
of lateral subsurface flow, representing it either in isolation late dynamic phreatic surface changes in response to rainfall
(e.g. Talebi et al. (2008); Montgomery and Dietrich (1994)), time series and to capture the influence of upslope hillslope
or in combination with infiltration (though these are compu- geometry. Paniconi et al. (2003) found that the HSB is able
tationally expensive, e.g. Formetta et al. (2014)). However, to capture the broad shapes of the storage and the outflow
these treatments have not been applied to the problem of profiles for all hillslope profiles.
cutting stability. Here we present a computationally efficient mechanis-
In this research, we overcome the aforementioned pit- tic model for rainfall-triggered slope failure probability
falls by choosing a hydrological model that captures some accounting for trigger variability using rainfall time series
of the hillslope hydrological dynamics, but that runs in and upslope influence by coupling to a hillslope hydrologi-
minutes rather than hours. This model uses Boussinesq’s cal model. We choose a mechanistic rather than a statisti-
groundwater theory to characterise groundwater flow (solv- cal model so that its application is not limited to data-rich
ing Darcy’s equation coupled with the conservation of mass regions. We believe that this method is novel for establishing
equation). Boussinesq’s theory provides a diffusion wave the frequency of a cutting failure as we develop a model
solution that is less restrictive than other solutions for the that treats both the complexity of the failure mechanics
governing groundwater flow equations (e.g. kinematic wave, and the hillslope hydrology in a computationally efficient
Fan and Bras (1998); Troch et al. (2002), or by simulating way enabling efficient large-scale application. Conditional
steady-state rather than dynamic flow, Beven and Kirkby probability based on a single representative storm ignores
(1979); Talebi et al. (2008)). Note, the Boussinesq equation known variability/uncertainty in trigger conditions. Our
neglects the effect of capillary rise above the groundwater model is also novel as it accounts for the uncertainty about
table and follows the Dupuit–Forcheimer approximation rainfall using a time-dependent method, as well as account-
(flow through an unconfined aquifer) that the streamlines ing for uncertainty in geomaterial properties. Where previ-
are approximately parallel to the impermeable boundary ous models have sampled from IDFs to generate triggering
(Boussinesq 1877). The original 1D Boussinesq theory can- rainfall, we drive our model from real rainfall time series
not account for the 3D geometry of a slope (e.g. convergence to correctly capture failure probability in situations where
or divergence, convexity or concavity). But slope geometry phreatic surface geometry and, thus, failure probability may
is known to significantly influence the hydrologic response be strongly dependent on antecedent conditions and rainfall
(Troch et al. 2003, Eq. 6). sequencing. This scenario is particularly prevalent in areas
Fan and Bras (1998) and Troch et al. (2002) have devel- that exhibit seasonal weather patterns (e.g. a dry and wet
oped simple models of groundwater flow for more complex season, monsoons). The hydrological model we use accounts
geometries. Building on Fan and Bras (1998) and Troch for lateral inputs by modelling the whole hillslope domain,
et al. (2002), Troch et al. (2003) reformulated the equa- whilst remaining computationally efficient. The model is
tions of Boussinesq in terms of storage to develop a more developed by coupling a probabilistic slope stability analysis
complete hillslope-storage equation (relative to the previous with a hillslope hydrological model to predict time-varying
kinematic wave approximations) that is applicable to more phreatic surface conditions resulting from input rainfall over
13
a long time window of >10 years and accounting for the conducted either due to cost or difficult terrain (i.e. moun-
entire hillslope topography. tainous, remote areas). The questions raised over the use a a
We used our model to estimate the absolute frequency of time-dependent model to represent trigger uncertainty, and
a slope failure over 11 years ( Ff ), which is closely related to what slope failure probability failure means and how it can
annual failure probability, for a road slope in Nepal triggered be used are discussed in Sect. 5.2. The uncertainties and
by rainfall. We also explore the limitations of conditional assumptions of the model are discussed in Sects. 5.1 and
probability of failure for a specific pore pressure scenario by 5.3, respectively.
examining how conditional probability, P(f|T), varies over An overview of the methodology is presented as a flow-
time for a road cut slope in Nepal. chart in Fig. 1. Table 1 outlines all the parameter symbols
used in this paper, categorised by the geotechnical, hydro-
logical and failure frequency model.
2 Methodology
Our methodology can be divided into three key sections: 2.1 Monte Carlo Simulations
1. Using MCS to determine a FoS distribution with FoS It is well known that model parameters are characterised
as a function of the phreatic surface level and the varied by two sources of uncertainty: (1) the uncertainty due to
geomaterial parameters (Sect. 2.1). intrinsic material conditions; and (2) spatial variability in the
2. Using the HSB model (Troch et al. 2003) to generate analysed domain. For the spatial variability, plenty of experi-
phreatic surface time series (Sect. 2.2). mental evidence shows that all rocks and soils are nonhomo-
3. Combining the outputs of the MCS and HSB model to geneous (Phoon and Kulhawy 1999a, b). This problem could
determine time-varying P(f|T) and a Ff (Sect. 2.3). be modelled using random field theory (e.g. Griffiths and
Fenton (2000); Dyson and Tolooiyan (2019); Gravanis et al.
First, a MCS is conducted to determine multiple FoS val- (2014)), but this necessitates the definition of length scales
ues that capture the aleatoric uncertainty in the geomate- associated with the variability (spatial correlation lengths).
rial strength properties, for several seepage scenarios each These length scales are extremely poorly constrained (Shokri
associated with a different phreatic surface height at the et al. 2019) for the most intensively monitored sites glob-
upslope boundary (explained in Sect. 2.1). The slope geoma- ally; thus, there is a lack of general guidance on reasonable
terial strength is characterised according to the Generalised length scales of variability for different material properties
Hoek–Brown (G-H–B) failure criterion which is today the in different settings (e.g. rock types). Furthermore, intensive
most popular criterion to characterise rock mass strength field testing is necessary to estimate these length scales and
(Wyllie and Mah 2017). Also well-constrained estimates of this would be prohibitively expensive. For these reasons,
some G-H–B parameters can be reliably established from we neglect the spatial arrangement of variability and we
geological and geomorphological field observations (Robson concentrate solely on the intrinsic conditions.
et al. 2022). Second, relative to slope strength and hydraulic To account for aleatory uncertainty in material properties,
conductivity, meteorological information is easier to acquire, we represent them as distribution density functions. For an
though in many countries, the network of gauges remains assigned phreatic surface level, we perform Monte Carlo
sparse. We assume that local daily rainfall is expressed by Simulations (MCS) to propagate the material uncertainties
a time series and these data are available in the following to the probability of failure of the slope. The MCS propa-
analysis. The so-called hillslope-storage Boussinesq (HSB) gate parameter variability captured in probability density
equation is solved to generate a phreatic surface time series functions through the stability model to find the conditional
for the slope to account for local hydrological conditions in probability of failure (i.e. the number of realisations which
response to rainfall (discussed in Sect. 2.2). By associating failed divided by the total number of realisations) (Fenton
the phreatic surface level time series from the HSB model and Griffiths 2008).
with those assessed in the MCS, a FoS time series is gener- In the MCS, not all material parameters are modified, but
ated for each parameter realisation. The time-varying P(f|T) only those that are significant to stability conditions. For this
and Ff can then be calculated from the FoS time series (dis- purpose, we conducted a one-at-a-time sensitivity analysis
cussed in Sect. 2.3). to determine which G-H–B parameters the model is most
The proposed methodology requires the following input sensitive to and, therefore, which G-H–B parameters should
data: rainfall time series, G-H–B parameters character- be varied as part of the MCS. To do so, literature-based
ising the slope geomaterial, estimates of hydrological estimates for the most likely, upper limit and lower limit for
parameters and slope geometry. Our approach is especially each parameter were input into the model, and the standard
advantageous when geotechnical investigation cannot be deviation (STDEV) in output FoS was computed for each
13
Fig. 1 Flowchart outlining the

key steps of methodology to
estimate the frequency of failure
and the conditional probability
of failure. Nr (number of reali-
sations of G-H–B parameters
varied in the Monte Carlo
Simulation), Nz (number of
phreatic surfaces tested in the
seepage analysis) and Nk (num-
ber of realisations of hydraulic
conductivity) were deter-
mined through convergence
tests. Abbreviations: DEM =
Digital Elevation Model, FE =
Finite Element, FoS = Factor
of Safety, G-H–B = Gener-
alised Hoek–Brown, HSB =
Hillslope-Storage Boussinesq, k
= hydraulic conductivity, LEM
= Limit Equilibrium Method,
MCS = Monte Carlo Simula-
tion, and Z = total head bound-
ary condition
parameter tested. In our case, we found that the FoS is most The parameters chosen to be varied in the MCS were
sensitive to the Geological Strength Index (GSI). From a sampled from a lognormal distribution as there is general
physical point of view, a reduction of this quantity represents agreement in studies that examine or use rock strength
the physical degradation of the exposed rock mass, taking distributions that the distribution should be uni-modal
into account the blockiness of the mass and its surface condi- with many researchers suggesting a lognormal distribu-
tions (Marinos and Hoek 2000). tion as a reasonable model for the physical properties of
rock (Hoek 1998; Parkin and Robinson 1992; Nour et al.
13
Table 1 Table of all parameter Parameter symbol Description

symbols used in this paper
Geotechnical model
D Disturbance
GSI Geological Strength Index
mi Material constant for intact rock
𝛾 Unit weight
𝜎ci Unconfined compressive strength
Hydrological model
Bp Thickness of the permeable layer
Bw Phreatic surface height
k Hydraulic conductivity
L Length of the hillslope (channel to the ridge)
nf Geomaterial porosity
q Flux
r Daily rate of rainfall
ru Water pressure ratio
S Groundwater storage
t Time
w Width of the hillslope
w0 Width of the hillslope at x = 0
x Space
Z Total head boundary condition
𝛼 Inclination of the impermeable boundary
𝛽 Degree of convergence or divergence towards the toe of the slope
Failure frequency model
Ff Absolute frequency of a slope failure
Nf Number of failures considering all of the FoS time series
Nf (t) Number of failures for a fixed time step for all of the FoS time series
Nk Number of realisations of hydraulic conductivity
Nr Number of realisations of parameters varied in the Monte Carlo Simulation
N(t*) Number of failures in a time period (t*)
Nz Number of phreatic surfaces tested in the seepage analysis
P(f) Probability of failure
P(f|T) Probability of failure conditional on a particular trigger event
P(f|Ti) Probability of failure conditional on multiple trigger events
P(T) Probability of a trigger event
𝜆 Rate of landslide occurrence
2002; Massih et al. 2008; Griffiths and Fenton 2009; Pan- slighty reduced. However, Pan et al. (2017) conclude
dit et al. 2019), particularly amongst those with a strong that the influence that these correlations have on the P(f)
physical constraint on the distributions lower bound (e.g. is ‘rather small’ (p. 6). Correlation relationships of the
values cannot be negative). Furthermore, a lognormal dis- G-H–B parameters are neglected in a number of proba-
tribution ensures that the random value of the parameter bilistic slope stability studies (Li et al. 2008, 2012; Pandit
will always be positive. In a probabilistic analysis, any et al. 2019; Farichah and Hutama 2020). Given the uncer-
correlation between parameters should be accounted for tainty in the correlation coefficients between the G-H–B
(Zeng et al. 2014). However, the correlation coefficients parameters, the correlation relationships are also neglected
between the G-H–B parameters are poorly constrained. here. Nr random realisations of GSI were generated from
Pan et al. (2017) plot failure probabilities as a function the lognormal distribution (with Nr based on a conver-
of the normalised slope height, corresponding to corre- gence analysis).
lated and independent variables and find that when the We performed Nr deterministic stability analyses, rep-
G-H–B parameter correlations are considered, P(f) is resenting each value of GSI, using the Morgenstern-Price
13
Fig. 2 Schematic diagram

highlighting the key elements of
the probabilistic (system high-
lighted in blue) and hydrologi-
cal (system highlighted in red)
models. The yellow area of the
slope shows the stability model
domain (the external bounda-
ries of the stability model were
chosen such that they have no
effect on the resulting FoS).
Nr Limit Equilibrium Method
(LEM) analyses are conducted
representing each realisation
of GSI for Nz phreatic surfaces
generated by Finite Element
(FE) method seepage analyses
from total head nodes (Z) at
the upslope boundary of the
stability model. Nk phreatic sur-
face time series are generated
through the Hillslope-Storage
Boussinesq (HSB) model repre-
senting each value of hydraulic
conductivity (k) (Color figure
online)
(M-P) Limit Equilibrium Method (LEM) using Rocsci- MCS ( Nr deterministic LEM analyses) were performed for
ence, Slide2. The M-P method was employed here as the slope cutting for Nz phreatic surfaces imposed at various
it satisfies all equations of equilibrium (rigorous equilibrium heights to capture phreatic surface variability. The phreatic
method) and is found to be most accurate across different surfaces were generated using Finite Element (FE) steady-state
slope conditions (Duncan 1996). We used the default half- seepage analyses. Nz total head boundary conditions at differ-
sine interslice force function; however, the influence of the ent heights, Z, were imposed on the upslope boundary of the
choice of interslice force is negligible for homogeneous model. The elevation range in Z is from the ground surface to
slopes (Fredlund and Krahn 1977). A convergence test was the elevation below which the phreatic surface does not influ-
first carried out to determine the optimal number of slices to ence the failure surface and thus the FoS. Z are equally spaced
be used in the LEM. The external boundaries of the model between these two elevation limits with the spacing set by the
were chosen such that they have no effect on the resulting mesh element size, since the phreatic surface is insensitive to
FoS. spacing more granular than this. The number of mesh elements
13
for the FE analysis was optimised using convergence testing from the rainfall time series at the location of the upslope
on the resultant FoS. We assumed homogeneous, isotropic boundary of the stability model domain in the hillslope.
hydraulic conductivity and no vertical inputs (these assump- Table 1 outlines the input parameters required for the HSB
tions are taken to simplify our model into 1D, allowing for model (hydrological model). The boundary conditions are
proof of concept for the overall methodology). Unsaturated presented in Table 2.
zones in the slope were neglected. The value of hydraulic con- All variables depend on time (t) and space (x), with x being
ductivity (k) prescribed to the model does not affect the result- a distance measured from the channel to the ridge. As per
ing phreatic surface, nor the FoS since the seepage analysis is the requirements of the HSB model presented by Troch et al.
at steady state. Thus, the analyses did not need to be performed (2003), Bp and 𝛼 are constant. Bw is constrained to be in the
at different k realisations. range between zero (no water) and Bp. This means that no flux
In summary, we performed Nr × Nz LEM stability analyses at the slope surface is modelled.
and we obtain FoS as a function of the phreatic surface level By combining the continuity equation with Darcy’s
Z and the GSI value. For a fixed phreatic surface level, the equation, reformulated in terms of the storage equation
FoS is a monotonically increasing curve with increasing GSI ( S(t, x) = nf wBw ), a version of the HSB partial differential
and for each Z, we obtain different curves in which the FoS equation (Troch et al. 2003) (Eq. 6) of second order in space
reduces with increasing Z. The key steps of the MCS used and first order in time is obtained:
in our model are highlighted in blue in a schematic diagram [ ( )]
𝜕S k cos 𝛼 𝜕 S 𝜕S S 𝜕w
presented in Fig. 2. In the subsequent section, we illustrate the nf = −
𝜕t nf 𝜕x w 𝜕x w 𝜕x
model employed to evaluate the phreatic surface level time (1)
series under unsteady-state conditions. + k sin 𝛼
𝜕S
+ nf rw
𝜕x
2.2 Phreatic Surface Time Series For many hillslopes, the relationship between hillslope width
(i.e the width of the flow strip) and distance upslope can be
In the methodology presented here, we aim to determine the approximated as an exponential function of the form (Troch
frequency with which the slope experiences each phreatic sur- et al. 2003, p. 7):
face level considered in the previous MCS. For this purpose,
we use a simple 1D hydrological model to simulate the phre- w(x) = w0 exp (𝛽x) (2)
atic surface height (and thus total head used in the FE seepage
where w0 is the hillslope width at x = 0 and 𝛽 is a shape fac-
analysis) at the upslope boundary of the slope stability model.
tor controlling the variation of the width along the x-axis. If
As an input in this model, we considered a rainfall time series
𝛽 = 0, a constant width profile is obtained (i.e. w(x) = w0 ),
rather than sampling storms affecting the slope of interest from
negative values of 𝛽 reflect divergent hillslope topography
an IDF distribution, to account for the correct timing and shape
with respect to the toe of the slope and positive values reflect
of storms and build up of pore water pressure in the slope (a
convergent topography.
time-dependent system). However, the structure of our model
Equation 1 is numerically integrated after the initial
is such that a storm-based slope stability analysis could easily
condition and the boundary conditions are imposed. The
be performed as an alternative if a suitable IDF were available
equations are solved using a finite difference method imple-
for the site.
mented in MATLAB and the main variable computed is the
The HSB model of Troch et al. (2003) is used to generate a
storage time series. The phreatic surface time series is then
phreatic surface time series (evolution of the phreatic surface)
computed from the storage.
In situ, k is strongly dependent on the characteristics of
the geomaterials (e.g. degree of weathering, fragmentation,
and grain size heterogeneity) so that a deterministic value
Table 2 Table of initial and boundary conditions for the hillslope is not physically reasonable. However, spatial variability
hydrological model (Hillslope-Storage Boussinesq model) in k would be prohibitively computationally expensive to
Expression Description constrain. To account for the variability of k, a statistical
approach was used so that multiple simulations of the model
Bw (t = 0, x)(x) Initial phreatic surface height were conducted, each with a unique value of k drawn from
(measured perpendicular to a lognormal distribution of k. A lognormal distribution is
the impermeable boundary)
chosen to prevent negative values. The lognormal distribu-
q(t, x = L) = 0 Impermeable flux at the ridge
tion is characterised by the 1st and 99th percentiles which
Bw (t, x = 0) = 0 River level intersects the
ground surface at the are the lower end and upper end of the k values for the slope,
channel allowing for occasional occurrences of values outside the
13
Fig. 4 Illustration of the Hillslope-Storage Boussinesq slope pro-

file highlighting the geometry of the impermeable boundary. The
Slide2 stability model domain is shown in yellow (the external
boundaries were chosen such that they have no effect on the resulting
FoS). The hillslope length (L) is measured from the channel to the
ridge. The impermeable boundary is highlighted in brown, where 𝛼
denotes its inclination. The thickness of the permeable layer ( Bp) is
measured at the upslope boundary of the Slide2 model so that the
impermeable boundary is at the level of the river channel at the river
channel. A potential phreatic surface is shown in blue, with the phre-
atic surface height above the impermeable boundary denoted by Bw
Fig. 3 Hillslope measurements to define catchment area. The extent (Color figure online)
of the slope widths is shown in light green (defined as where the
hillslope switches from convex to concave form at a gulley), measure-
ments of the slope width (w0, w1 and w2) are shown in pale pink and
by measuring the width of the hillslope (distance between
the length of the slope in white. w0, w1 and w2 are used to estimate
the shape factor (𝛽 ). Image sourced from google earth (Color figure the two slope width extents) at three points along the slope
online) length ( w0 , w1 and w2 in Fig. 3) and taking the best-fitting
exponential of these three widths.
The depth to the impermeable boundary and its incli-
range. A convergence analysis was conducted to determine nation were taken based on assumptions as there are few
the minimum number of realisations of k required for the constraints that can be drawn from observations at this site
model ( Nk ). Phreatic surface time series are generated for and there is no accepted theory on the depth of the perme-
each value of k, using the HSB model. We compare the able region in a hillslope (discussed as the critical zone) in
resulting time-varying phreatic surface fluctuations to our the scientific community (Anderson et al. 2019; Grant and
expectations based on qualitative field observations at the Dietrich 2017; Flinchum et al. 2018; Clair et al. 2015). In
site to constrain the values for the 1st and 99th percentiles literature, there is general agreement that the impermeable
of k. As with Troch et al. (2003), it is assumed that the vari- boundary tapers in inclination towards the ridge, meaning
ation in k with depth can be approximated as a step function, that the permeable layer is thicker at the ridge (Anderson
with uniform k geomaterial above an impermeable boundary et al. 2019; Grant and Dietrich 2017; Flinchum et al. 2018;
(discussed further in 5.3). Clair et al. 2015; Anderson et al. 2013; Medwedeff et al.
The value of geomaterial porosity ( nf ) used in the model 2022). Thereby, we assumed that the inclination of the
was assumed from an average of literature values for the rock impermeable boundary (𝛼 ) is 10◦ lower than the inclination
type at the site as it is not as variable as k (Singhal and Gupta of the ground surface. The sensitivity of the model to this
2010) (further discussed in Sect. 5.3). The catchment area assumption is explored in Sect. 5.3.
upslope of the cutting (i.e. the area that would drain through The downslope boundary condition of the HSB model
it) was defined using a Digital Elevation Model (DEM) in implies that the impermeable boundary intersects the ground
ArcGIS (see Fig. 3). In doing so, the watershed was defined surface at this boundary, in most cases a river at the base
by projecting a line upslope from the cutting along the line of of the hillslope with the impermeable boundary at the bed
steepest inclination to the ridge (defined at the point where of the river and a wedge of permeable material extending
inclination goes to zero). This line is the hillslope length upslope. The impermeable boundary was taken such that
(L), as depicted in white in Fig. 3. The extent of the slope it intersects the ground surface at the river (the downslope
width is defined as where the hillslope switches from convex boundary of the hillslope). Figure 4 displays a 2D section of
to concave form at a gulley (these slope extents are shown the hillslope in profile, highlighting the geometrical proper-
in lime green in Fig. 3). The shape factor ( 𝛽 ) is estimated ties of the model.
13
Fig. 5 A schematic diagram

highlighting how the outputs of
the Monte Carlo Simulations
( Nz FoS distributions) and the
Hillslope-Storage Boussinesq
model ( Nk phreatic surface time
series) are combined to deter-
mine a time-varying conditional
probability of failure (P(f|T))
and a frequency of failure ( Ff )
We obtained an 11-year daily rainfall time series from 2.3 Combining the Models
a local meteorological station and duplicated this resulting
in a 22-year time series to avoid the effects of the condition In this section, we combine the results obtained from
occurred before our analysis on the subsurface flux. If there Sects. 2.1 and 2.2 to obtain the probability of failure and
is no information on the phreatic surface, the initial condi- the frequency of failure (presented as a schematic dia-
tion of the phreatic surface is set to a height of zero (i.e. the gram in Fig. 5). For an assigned value of k, we discretised
level of the impermeable boundary). The model should be the phreatic surface time series according to Z values that
spun up to minimise the influence that this assumption could were used to generate the phreatic surfaces for the MCS.
have on the model. Prior to doing so, the reference frame of the output of the
The HSB model is written in MATLAB. Nk model simula- HSB model (total head time series) has to be converted to
tions are conducted each with a unique value of k, resulting match the horizontal/vertical Cartesian reference frame of
in Nk phreatic surface time series. The key steps of the HSB the phreatic surface total head boundary condition in the
model are highlighted in red in a schematic diagram pre- stability model (Z). This is done by dividing the phreatic
sented in Fig. 2. The uncertainties and assumptions made in surface from the HSB model by the cosine of the incli-
the HSB model are discussed in Sect. 5. nation of the impermeable boundary. Once the reference
frame conversion is made, the ground level in the phreatic
surface time series is made equal to the ground level of
the stability model. By discretising each phreatic surface
time series according to Z, FoS can be expressed as a
13
Fig. 6 Map showing the loca-

tion of the case study site along
the Narayanghat-Mugling high-
way in the Chitawan District of
Nepal
discretised function of the phreatic surface for one GSI from being counted) for a number of remediation days (the
realisation and k realisation; therefore, we obtain a unique choice of remediation days is discussed further in Sect. 3.3
discretised FoS time series for each GSI realisation and k whilst the sensitivity of the model output to this parameter is
realisation. tested in Sect. 5.3) to allow time for debris to be cleared and
Combining the results of the probabilistic stability anal- the cutting to be reinstated as it was. If we denote Nf as the
yses and the hydrological model in this way significantly number of failures considering all of the FoS time series, then
reduces the computational time that would otherwise be the number of landslides per the timescale of the rainfall time
required to determine a FoS for every water level in this series (Ff ) is given by
time series for every rock parameter realisation and every
value of k.
Ff = Nf ∕(Nr × Nk ) (4)
The FoS time series for each realisation is converted to a
binary ‘failure’ time series with failure for FoS<1 and stabil- 2.4 Case Study
ity for FoS>1.
We demonstrated this methodology on a road cutting on
2.3.1 Method to Determine Time‑Varying Conditional the Narayanghat-Mugling road in Chitawan, Nepal (see
Failure Probability Fig. 6 for map). Nepal has an elevation range of less than
100 m in the Terai region in the south to up to 8000 m in the
After doing the previous combination, we obtain Nr × Nk dis- Himalayan mountains. The region is tectonically active and
cretised FoS time series. We consider slope failure when the experiences a summer monsoon season for four months of
FoS is <1. If we denote Nf (t) as the number of failures for a the year during which time 80% of Nepal’s annual rainfall
fixed time step (t) considering all of the FoS time series then occurs (Shakya and Nirula 2008). The case study site is at an
P(f|T) (probability accounting for the phreatic surface) for each elevation of around 240 m. It is situated in an area of sedi-
time step P(f|T(t)) is given by mentary to low-grade metamorphic rocks of Proterozic age
P(f |T(t)) = Nf (t)∕(Nr × Nk ) (3) (2500–539 Ma) in the Lesser Himalaya geological region.
At this site, there is an above-road cut slope around 25 m in
height and 70◦ inclination made up of weathered phyllite
2.3.2 Method to Estimate the Frequency of Failure (see Fig. 7 for image of cutting). This cut slope exists in a
valley side hillslope that is inclined at c. 25◦ above the cut-
Assuming that the slope is returned to its pre-failure condi- ting. The road was originally excavated by blasting around
tion after each failure (this assumption is discussed further 40 years ago (personal communications with consultant in
in Sect. 5.3), each ‘failure’ time series is worked through November 2019). A 2 m tall gabion wall constructed along
chronologically, and when failure occurs the following steps the cutting collapsed due to a minor rockfall during the 2019
are taken: (1) iterate a failure count for that failure time series; monsoon season. Below the road, there is a 15 m long slope
and (2) zero the failure time series (preventing further failures descending into the Trishuli River.
13
Table 3 Most likely, upper and lower limit estimates for values of
G-H–B parameters and unit weight (𝛾 ) for a slope along the Naray-
anghat-Mugling road based on geological and geomorphological field
observations
𝛾 GSI mi 𝜎ci
[kN/m ]3 [−] [−] [MPa]
Most likely 24 22 7 25
Lower limit 23 15 4 15
Upper limit 25 30 10 35
The Hoek–Brown failure criterion include: Geological Strength Index

(GSI), material constant for intact rock (m i ) and, unconfined com-
pressive strength (𝜎ci)
Fig. 7 Image of road cutting case study along the Narayanghat-Mug-

ling road in Chitawan, Nepal. Image taken in November 2019
The slope cutting of exposed rock mass is characterised

according to G-H–B failure criterion. For all the G-H–B
parameters and unit weight ( 𝛾 ), the values for the most
likely, upper and lower limits of the parameter were esti-
mated based on field observations of the slope and values
from the literature (outlined in Table 3). GSI is a measure
of the physical condition of the exposed rock mass (Marinos
Fig. 8 Daily rainfall (m/s) from 2010–2020 recorded at Sakhar mete-
and Hoek 2000), whilst the unconfined compressive strength orological station in the district of Tanahu acquired from the Govern-
(𝜎ci ) accounts for the strength of the rock (Hoek 1998). GSI ment of Nepal Department of Hydrology and Meteorology. Data for
and 𝜎ci were estimated according to field observations (the 2013 were missing and have been replaced with a copy of the data
rock, phyllite, is disintegrated with a highly weathered sur- from 2012 (highlighted in orange) (Color figure online)
face), using Marinos and Hoek (2000) as guidance. The
values for the G-H–B material constant, m i and 𝛾 were esti- simulation since a break in the simulation would have neces-
mated based on expected ranges for the rock type (weathered sitated re-defining the model initial conditions.
phyllite) taken from literature ( 𝛾 from Fine (2021) and mi The application of our methodology to this case study is
from Marinos and Hoek (2000)). discussed in the next section (Sect. 3).
The Disturbance Factor (D) in the G-H–B criterion
accounts for the blast damage and stress relaxation of a rock
mass. Despite this cutting being excavated by historical
blasting, which is likely to have caused disturbance, D was 3 Site‑Specific Methodology
taken as zero across the slope for the sake of simplicity. This
choice is discussed further in Sect. 5.1. 3.1 Site‑Specific Probabilistic Slope Stability Model
Daily rainfall data spanning an 11-year period from Janu-
ary 2010 to December 2020 recorded at the Sakhar mete- 2D stability analyses were performed using M-P LEM in
orological station in the Tanahu district (Gandaki Province) Rocscience, Slide2. Literature-based estimates
were acquired from the Government of Nepal Department for most likely, upper limit and lower limit values of each
of Hydrology and Meteorology (Fig. 8). Sakhar meteoro- parameter (outlined in Sect. 2.4) were input to the model,
logical station is around 13 km from the research site and and the STDEV in output FoS was computed for each
lies at an elevation of around 60 m higher than the road by parameter tested. For 𝛾 , the STDEV was found to be negligi-
the research site. It should also be noted that no data were ble at 0.01, and, therefore, a single value of 𝛾 was employed.
available from the station for 2013, thus missing data were The STDEV of the FoS for GSI, m i and 𝜎ci were found to be
replaced by repeating the 2014 record to enable continuous 0.5, 0.1 and 0.1, respectively.
13
Table 4 Table displaying values for Generalised Hoek–Brown param-

eter (intact rock parameters, m i , unconfined compressive strength,
𝜎ci , and disturbance, D) and unit weight (𝛾 ) used in each deterministic
analysis
𝛾 mi 𝜎ci D GSI
kN/m 3 – MPa – –
24 7 25 0 15→30
Constant Constant Constant Constant Varied
One thousand values of Geological Strength Index (GSI) are sampled

from a lognormal distribution characterised in terms of its 1st and
99th percentiles being equal to 15 and 30, respectively, to be used in
the MCS. The values of 𝛾 , m i , 𝜎ci and D were taken as constant in the
MCS
Fig. 9 Lognormal distribution of the Geological Strength Index cat-

egorised in terms of the 1st and 99th percentiles being the lower and
upper limits of the reasonable range for the cutting made up of phyl-
lite
Given that the model exhibits high sensitivity to the GSI

compared to the other parameters (Pandit et al. 2019; Chen
and Lin 2019; Hoek 1998; Pan et al. 2017), it was varied
as part of the MCS. Instead, 𝜎ci and m i were kept constant.
The rock mass of the cutting examined here was identi-
fied through geological observations to be phyllite. Field
observations at the cutting suggested that GSI ranged from
15 to 30 based on the GSI charts of Marinos and Hoek Fig. 10 Stability analysis model set up in Rocscience, Slide2.
(2000) as guidance. These values were used to parameterise Cutting is 25 m in height inclined at 70◦ . The topography upslope
of the cutting is inclined at 25◦ . Total head boundary condition (Z)
a lognormal distribution assuming that the lower and upper
equally spaced between 45 and 92.48 ms imposed at the ridge side of
limits reflected 1st and 99th percentiles of the distribution, the slope to generate phreatic surfaces
respectively. This is to reflect its expected uni-modal form
and allow for occasional occurrence of values outside the
typical range (1st and 99th percentiles chosen so that only a 92.48 m above river level (when the phreatic surface is at
very small proportion fall outside this range). The lognormal the ground surface) to 45 m above river level (below which
distribution of GSI is shown in Fig. 9. the phreatic surface does not influence the failure surface
Deterministic stability analyses were performed using and thus the FoS). Seepage analyses were performed for 25
the M-P 2D LEM in Rocscience, Slide2. A conver- total head values ( Nz = 25) each with a different Z equally
gence test was first carried out to determine that the optimal spaced from 45 to 92.48 m (see Fig. 10 for model set up).
number of slices to be used in the LEM slope stabilisation One thousand deterministic slope stability analyses (vary-
analyses is 50. The external boundaries of the model were ing the GSI according to the lognormal distribution) were
chosen such that they have no effect on the resulting FoS. carried out for twenty-five phreatic surface scenarios (i.e.
By conducting a convergence analysis, we find that the opti- 1000 × 25 LEM stability analyses were conducted). The
mal number of realisations of GSI for the MCS is 1000 (i.e. predominant failure mechanisms observed in these analy-
Nr = 1000). Further realisations result in <1% change in the ses were shallow in terms of aspect ratio and constrained to
output of the model. See Table 4 for parameter values used the cut slope itself. Figure 11 displays probability density
in each MCS analysis. functions for the 1000 factors of safety determined through
FE steady-state seepage analyses were performed in these deterministic slope stability analyses for every phreatic
Slide2, solving for Darcy’s equation and the continuity surface. As expected, the distribution of FoS shifts towards
equation. The mesh for the FE analyses was made up of lower values as Z increases and the fraction of runs with FoS
10,000 six-noded triangles (optimised based on convergence <1, i.e. the conditional probability of failure for that phreatic
testing). Total head Z at the upslope boundary ranges from surface scenario increases.
13
Table 5 Case study input parameters for the HSB model
Depth to Inclina- Shape Initial Length of Material

imper- tion of factor width slope porosity
meable imper-
boundary meable
boundary
Bp 𝛼 𝛽 w0 L nf
[m] [] ◦
– [m] [m] [%]
51.28 15 − 0.0008 1619.9 1457.0 7.5
The inclination of impermeable boundary is assumed to be less than

the inclination of the hillslope. The shape factor, initial width and
length of slope are estimated from a DEM in ArcGIS. Porosity is
taken as an average of crystalline metamorphic rock from literature
Fig. 11 Probability density function of the 1000 factors of safety
values determined for the cutting for 25 phreatic surface scenarios.
Phreatic surfaces are determined by carrying out Finite Element (FE)
method seepage analysis with a boundary condition of total head (Z) Based on the observation that there is no evidence of the
on the ridge side of the model varying from 45 m in height to 92.48 water table daylighting at the ground surface above the cut-
m ting on the hillslope, it is determined that the impermeable
boundary is at a depth below the ground surface across the
3.2 Site‑Specific Hillslope Hydrological Model length of the hillslope. Based on the literature, we assume
that the impermeable boundary tapers in inclination with
The cutting examined here is made up of weathered phyllite. distance towards the ridge (Anderson et al. 2019; Grant
Phyllite is a crystalline, foliated, low-grade metamorphic and Dietrich 2017; Flinchum et al. 2018; Clair et al. 2015;
rock. Such rocks generally have comparatively high porosity Anderson et al. 2013; Medwedeff et al. 2022). Thereby, it is
but low permeability (Singhal and Gupta 2010). Values of k assumed that the impermeable boundary is at an inclination
for fractured crystalline metamorphic rock were taken from of 15◦ , 10◦ lower than the ground surface. The slope profile
Singhal and Gupta (2010) (10−9 to 10−5 m/s) and applied in is composed of a planar 25◦ slope with a steeper 70◦ cutting
the model to predict time-varying phreatic surface fluctua- near its toe. The impermeable boundary is configured such
tions. We observed no evidence of surface run off at this that it intersects the ground surface at the river resulting in
cutting and inferred that surface saturation is rare within or a permeable layer thickness of c. 51.28 m at the upslope
upslope of the cutting. Therefore, we imposed a rule that the boundary of the Slide2 model. The key input parameters
1st percentile value of k cannot result in a phreatic surface for the HSB model for this case study are outlined in Table 5.
at the ground level for a sustained amount of time (more The rainfall time series is shown in Fig. 8. This is an
than a day). This considerably increased the 1st percentile 11-year record of daily rainfall. Given that we have no infor-
of k to 1.7 × 10−6 m/s from the minimum value found in the mation on the phreatic surface, the initial condition of the
literature (10−9) which is perhaps indicative of a relatively phreatic surface was set to a height of zero. To minimise the
high degree of weathering on the slope. The 99th percentile influence that this assumption could have on the model, the
is given as the maximum values for k in metamorphosed model was spun up using a duplication of the 11-year rainfall
crystalline rock (1 × 10−5 m/s). The 1st and 99th percentiles record. We believe that this duration (an additional 11 years)
are use to characterise a lognormal distribution of k. Two is sufficient for the subsurface flux to avoid the effects of the
thousand realisations of k ( Nk = 2000 ) are then generated conditions that occurred before our analysis.
from the lognormal distribution. Two thousand is the mini- MATLAB was used to run the HSB model according to
mum number of realisations needed in the model determined the method outlined in Sect. 2.2 for 2000 realisations of k
by a convergence analysis. taken from the lognormal distribution. The model output
The value of nf used in the model was taken from an was 2000 phreatic surface time series. Figure 12 displays
average of fractured crystalline metamorphic rocks (Singhal the model output for 100 realisations of k varying from the
and Gupta 2010). The catchment area upslope of the cut- minimum to maximum k in the distribution over 11 years.
ting is defined using the Shuttle Radar Topography Mis- The lower the value of k, the higher the phreatic level time
sion (SRTM) DEM of 1-arc second resolution in ArcGIS. series. The effect of the annual variation in rainfall can be
The slope width is measured at three places along the slope seen in this figure; during the monsoon season, the phreatic
length (at x = 0 plus two random locations further towards surface in the slope increases, and then decreases again dur-
the ridge, see Fig. 3) to estimate 𝛽 (by taking the best-fitting ing the dry season.
exponential of these three widths).
13
3.3 Site‑Specific Model Coupling
To couple the outputs from the probabilistic (the MCS)

and hillslope hydrological (the HSB) models, the phreatic
surface time series (output from the HSB) were discretised
according to 25 Z values ( Nz = 25), that were used to gener-
ate the phreatic surfaces for the MCS. Prior to doing so, a
reference frame convergence is made by dividing the phre-
atic surface from the HSB model by cosine of the inclination
of the impermeable boundary (15◦ ). By discretising each
phreatic surface time series according to Z, a unique FoS
time series can be generated for each GSI realisation and
k to generate 1000 × 2000 ( Nr × Nk ) FoS time series. The
FoS time series for each realisation is converted to a binary
‘failure’ time series with failure for FoS<1 and stability for
FoS>1.
Fig. 12 One hundred phreatic surface time series derived using the The time-varying conditional probability of failure
HSB model based on one hundred realisations of k. 50 m represents
the ground surface and 0 m is the impermeable boundary (P(f|T)) was determined by dividing the number of failures
(where FoS < 1) for a fixed time step considering all FoS
time series by the number of realisations of GSI ( Nr = 1000)
and the number of realisations of k ( Nk = 2000), according
to Eq. 3. Thereby, determining a P(f|T) for each time step.
To determine a frequency of failure ( Ff ), each ‘failure’
time series was worked through chronologically, and when
failure occurs, there is an iteration for the failure count for
that failure time series and the failure time series is zeroed
for 90 remediation days preventing further failures from
being counted in a 90-day window. The timescale for reme-
diation was assumed to be 90 days for Nepal because the
monsoon season lasts 3 to 4 months and work cannot be car-
ried out during this time. The total number of landslides is
then summed across all ‘failure’ time series and normalised
Fig. 13 Probability density distributions of phreatic surface height by 1000 realisations ( Nr = 1000 ) and 2000 realisations of
above the impermeable surface (across all values of k) at four time
steps: (1) time = 0.74 years (peak, higher peak values); (2) time = k ( Nk = 2000 ) to determine a number of landslides per 11
1.33 years (trough, higher peak values); (3) time = 6.80 years (peak, years (according to Eq. 4).
lower peak values) and (4) 7.46 years (trough, lower peak values).
The peaks are shown in a solid line, whereas the troughs are in a
dashed line
4 Results
Figure 13 displays probability density distributions of 4.1 Time‑Varying Conditional Probability of Failure

the phreatic surfaces (across all realisations of k) at four
time steps to showcase the variation in trends: (1) time = Figure. 14 displays the time-varying conditional probability
0.74 years (peak, higher peak values); (2) time = 1.33 years of failure determined by following the method outlined in
(trough, higher peak values); (3) time = 6.80 years (peak, Sect. 2.3.1. The figure exhibits similar general trends that
lower peak values) and (4) 7.46 years (trough, lower peak can be observed in the phreatic surface time series plot
values). As expected, the spread of the phreatic surfaces at (Fig. 12). Figure 12 displays annual cyclicity, with peaks
the peaks are much higher than those at the troughs. Another in the phreatic surfaces in correspondence of the monsoon
finding is that at the times where the phreatic surface peaks season and troughs in the dry season. The same trend can
at the highest values, they have a wider spread of values (and also be observed in Fig. 14, with peaks in the P(f|T) during
in the trough following the high peak). Where the peak is the monsoon season and significant dips during the dry sea-
lower, the spread of values is more constrained (also in the son. In Fig. 12, it can be seen that the peaks of the phreatic
trough following the low peak). surfaces during the monsoon seasons for the first 2 years
are the highest, and then the peaks of the phreatic surfaces
13
Fig. 15 Histogram of the number of landslides in 11 years for the

case study site binned according to the discretised realisations of
Fig. 14 Time-varying conditional probability of failure (P(f|T)) over hydraulic conductivity. There are 40 realisations of k in each bin
11 years for the case study site. P(f|T) calculated for every time step
by dividing the total number of failures across all ‘failure’ time series
by the number of realisations and number of realisations of hydraulic
conductivity
decline every year to the lowest in the 7th year. After the 7th
year, the peaks of the phreatic surfaces during the monsoon
season start to increase, and then drop off again in the 10th
year, followed by a sharp increase for the 11th year.
4.2 Frequency of Failure
By summing all the landslides from each FoS time series,

and normalising by the number of k realisations ( Nk = 2000)
and the number of realisation ( Nr = 1000 ), it is estimated
that 5.10 failures of this cutting will occur every 11 years. Fig. 16 Histogram of the number of landslides in 11 years for the
As shown by Ross (1972), the rate of landslide occur- case study site binned according to realisations of the Geological
Strength Index (GSI). There are 40 realisations of GSI in each bin
rence (𝜆) can be estimated as:
N(t∗ )
𝜆= (5) According to the Poisson model, the probability of one or
t∗ more landslides during a time t is:
where N(t∗ ) is the number of failures in a time period t∗.
In this case, 𝜆 = 5.10∕11 = 0.46 , meaning that failure is
P{N(t) ≥ 1} = 1 − exp −𝜆t (6)
expected at a rate of 0.46 per year (annual Ff ). This can be Thereby, it is estimated that the annual probability of failure
equated to a failure approximately every other year. Infor- for this cutting is 0.37 (1 − exp −0.46×1 = 0.37). The prob-
mal comments from consultants on the Narayanghat-Mug- ability of one or more cutting failures occurring in 11 years
ling road suggest that this frequency of failures has been is 0.99 (1 − exp −0.46×10 = 0.99) and the probability of one or
observed in the past in the area. This methodology devel- more cutting failures occurring in 100 years is close to one
oped using an 11-year rainfall time series, assumes that this (1 − exp −0.46×100 = 1).
Ff is still valid for years into the future (i.e. in 20 years time). Figure 15 displays a histogram of the number of slope
Probability of failure can be determined from the fre- failures over 11 years into 100 bins of hydraulic conductiv-
quency of failure using simple statistics, following a Pois- ity values. The number of failures were normalised by the
son model which is a continuous-time model including the number of realisations ( Nr = 1000 ) and the size of the bin
occurrence of random point-events, or in this case landslides (40 realisations of k). This figure shows that the lower the k
(Crovelli 2000), the probability of failure can be determined. of the slope, the greater the Ff .
13
Figure 16 displays a histogram of the number of slope The uncertainties in the rock parameters used in the LEM
failures per 11 years into 100 bins of GSI. The number analysis are partially accounted for (the parameter values
of failures were normalised by the number of realisations used for the probabilistic analysis are themselves uncertain)
( Nr = 1000) and the size of the bin (40 realisations of GSI). by carrying out probabilistic stability analyses. The assump-
This figure shows that the lower the GSI of the slope, the tion of infinite correlation length is also made by neglecting
greater the Ff as it is reasonable to expect. local spatial variability in the cutting.
As previously stated, the predominant failure mechanisms
observed in the LEM stability analyses were shallow in
terms of aspect ratio and constrained to the cut slope itself.
5 Discussion Shallow failure of the road cut slope itself is commonly
observed along this road, as well as multiple other roads in
Following the methodology presented in this paper, the time- Nepal, due to over-steepened cut slopes.
varying conditional probability of failure and frequency of The Narayanghat-Mugling road was historically blasted
failure of a slope cutting triggered by rainfall can be esti- when it was initially excavated. However, disturbance caused
mated. In this section, the uncertainties of the method and by blasting is not accounted for. This was decided to keep
the sensitivity of the model to key assumptions made are the slope model simple, as this case study is being used
discussed. to demonstrate the methodology presented in this paper. In
addition, there is no data on the intensity and extent of dam-
5.1 Uncertainties age that was caused by blasting. If the blasting was carried
out in an uncontrolled manner, the disturbance towards the
The HSB model contains a large number of assumptions face of the cutting may be quite high which can significantly
necessary to enable a fast and simple solution to the RE. reduce the stability of the cutting.
Paniconi et al. (2003) carried out a comparison between the
Hillslope-Storage Boussinesq equation and RE examining 5.2 Probability of Failure
various hillslope geometries. They determined that the two
models have closer matches in outcomes for convergent than The similarities in the trends of the phreatic surface time
divergent hillslopes, and under drainage conditions than series (Fig. 12) and the P(f|T) time series (Fig. 14) highlight
recharge conditions. They state that there are “remarkably the influence of the phreatic surface level on the P(f|T). The
good matches of the diversity of shapes, including peaks and phreatic level peaks in the first two monsoon seasons, as the
spreads, that characterise the storage and outflow dynamics rainfall in the previous year during the spin up is very high
of the different hillslopes" (Paniconi et al. 2003) (p. 9). The (a copy of the daily rainfall during 2020, see Fig. 8). Based
reason for the difference in the outcomes of the recharge on this finding, it can be said that daily rainfall is not the
scenarios is attributed to the role of the unsaturated zone in key driver of failure, but prolonged heavy rain instead. This
the RE, slowing the vertical transmission of rainfall through suggests that it is important to account for rainfall variability
the hillslope soil. Due to the influence of the unsaturated in time using a time-dependent system (e.g. a rainfall time
zone, the storage profile simulated by RE is lower than that series) for an area that hosts a monsoon season, rather than
from the HSB model. This relationship is exhibited in the using a time-independent system (e.g. using rainfall events
hydrographs they present, in that there is a closer match drawn from an IDF curve).
between the HSB and RE hydrographs for a convergent The annual cyclicity and variability of the P(f|T) observed
slope than a divergent slope, as the convergent slope drains in Fig. 14 also suggest that annual probability may not be
slower, remaining more saturated, meaning the unsaturated representative for an area that hosts a monsoon season and
zone plays less of an important role. On the other hand, the has a long-term memory system. If an instantaneous failure
HSB and RE hydrographs for a divergent hillslope exhibit probability is used in stability analysis and design, it could
greater differences given that divergent slopes are faster result in dramatically over-conservative or under-conserva-
draining and, therefore, the unsaturated zone plays a more tive results. For example, we estimated an annual probability
important role. Thus, whilst the assumptions of the HSB of failure of 0.37 for this cutting. Looking at the time-var-
are known to introduce some additional model uncertainty, ying probability of failure (Fig. 14), this would be a gross
they are necessary to render the approach tractable and this overestimate of the probability of failure for this cutting for
additional uncertainty is small in the context of the very the majority of the 11 years analysed.
large uncertainty in material properties (e.g. permeability, As discussed, Fig. 14 can be used to observe P(f|T) over
porosity, and geometry of impermeable layer) for hillslopes time and how this reflects fluctuations in the phreatic sur-
in general (this cutting is no exception). face time series. However, it is difficult to estimate a rate of
slope failure occurrence from this plot as the time-varying
13
Table 6 Sensitivity test on how the number of remediation days (days Table 7 Sensitivity of the frequency of failure ( Ff ) model to the 99th
after a landslide) affects the frequency of slope failures over 11 years percentile of the lognormal distribution defining hydraulic conductiv-
ity
Remediation days Ff Failure/11 years Percentage
difference 99th percentile of k m/s Ff Failure/11 years Percentage
% difference
%
7 61.97 1115
14 27.24 434 6.0 × 10−6 6.82 34
30 13.29 161 1.0 × 10−5 5.10 –
60 7.15 40 1.4 × 10−5 3.26 − 36
90 5.10 – 1.8 × 10−5 2.60 − 49
180 3.10 − 39 2.2 × 10 −5 2.45 − 52
270 2.52 − 51
The percentage difference is calculated using the value used in the
365 2.21 57
model (1.0 × 10−5 m/s) as the baseline
The percentage difference is calculated using the number of remedia-
tion days chosen for the model (90 days) as the baseline
working on a project who is likely to have a good knowledge
of the expected number of days until remediation occurs. It
conditional probability is not bound by time, and, as dis- is also assumed that the cutting is reinstated to its pre-fail-
cussed, a rate of occurrence may not be representative of this ure condition, during each remediation following a failure.
plot. It is not straightforward how to use time-varying P(f|T) Although it is unlikely that the slope would be returned to
to inform slope stability design. Moreover, there are not yet exactly the same conditions, it may not be dissimilar.
slope design standards at present expressing threshold values The 1st and 99th percentiles of the lognormal distribution
in terms of failure probability. Conversely, the Ff is a much of k were initially defined through evaluating values of k for
more intelligible value, which could be hugely beneficial to fractured crystalline metamorphic rock from Singhal and
stability design. For example, Ff can be used in a cost–ben- Gupta (2010) (10−9 to 10−5 m/s). The value for the 1st per-
efit analysis where stabilisation measures can be compared centile was then adjusted based on predictions from the HSB
according to their Ff value and cost. model evaluated against the observational constraint that the
cutting had no evidence of overflow and, therefore, should
5.3 Assumptions not experience prolonged periods of surface saturation. For
the case of the 99th percentile, there is no observational
A significant assumption in the HSB model is that flow is constraint and, therefore, it is assumed that this value is the
oriented parallel to the bed slope. This differs from the RE higher end of values found in literature ( k = 1 × 10−5 m/s).
model where flow can be resolved in any direction and flow The sensitivity of the model to this value is tested (results
direction emerges from the solution (Paniconi et al. 2003). displayed in Table 7). Values tested include 6 × 10−6 to
Another assumption of the HSB model is that it does not 2.2 × 10−5 m/s in increments of 4 × 10−6 m/s. Values lower
account for vertical infiltration meaning that there will be than 6 × 10−6 m/s were not tested as these would be too
an increase in water flux at the toe of the slope as compared close to the 1st percentile (1.7 × 10−6 m/s). Values higher
to a model which accounts for vertical infiltration (e.g. RE) than 2.2 × 10−5 m/s were not tested as the phreatic surface
(Paniconi et al. 2003). was levelling out at the impermeable boundary. This test
The model’s sensitivity to the remediation time required shows that the greater the value of the 99th percentile of k,
to return the slope to its pre-failure condition (during which the lower the Ff . There is a steep reduction in the Ff from
FoS<1 does not result in a countable landslide) was tested 6 × 10−6 to 1.4 × 10−5 m/s, with little change thereafter. The
using values of 7–365 days (Table 6). There is a steep range of values of k from literature for fractured crystalline
decline in the number of failures in increasing the number of rock is wide (10−9 to 10−5 m/s), but based on observational
remediation days from 7 to 90 days, this then plateaus after constraints, the potential range for the cutting is reduced to
90 days. This trend demonstrates that reducing remediation 1.7 × 10−6 to 10−5 m/s. Given that this range is at the higher
time can considerably increase estimated frequency of fail- end for fractured crystalline rock, it can be said that the cut-
ure. It also shows that the cutting used to demonstrate this ting is highly fractured.
methodology is very unstable, and that if remediation sim- The HSB model assumes that the variation in k with
ply returns it to its pre-failure condition, it will quickly fail depth can be approximated as a step function, with uniform
again. This indicates that slope stabilisation measures are k material above an impermeable boundary. This is clearly
necessary. Although the Ff is strongly sensitive to the reme- an approximation, but underpinned by theoretical and obser-
diation time, this assumption can be refined by a practitioner vational studies on the permeable layer sometimes referred
13
Table 8 Sensitivity test on how the effective porosity in the HSB Table 9 Sensitivity test on how the inclination of the impermeable
affects frequency of slope failures ( Ff ) over 11 years boundary in the HSB effects frequency of slope failures over 11 years
Effective porosity % Ff Failure/11 years Percentage Inclination of impermeable Ff Failures/11 years Percentage
difference boundary ◦ difference
% %
5 5.08 0 15 5.10 –
7.5 5.10 – 20 13.66 168
10 4.51 − 12 25 35.34 593
The percentage difference is calculated using the effective porosity The percentage difference is calculated using the inclination of the
chosen for the model (7.5%) as the baseline impermeable boundary chosen for the model (25◦ ) as the baseline
to as the critical zone (the weathered mantle between fresh layer is thicker towards the ridge based on evidence from lit-
bedrock and the atmosphere) (Anderson et al. 2019; Grant erature (Clair et al. 2015). As with the boundary conditions
and Dietrich 2017; Flinchum et al. 2018; Clair et al. 2015). imposed by Troch et al. (2003), the impermeable boundary
Many of these studies note a decline in permeability, fracture at the river channel is fixed to the depth of the river channel.
density or openness, and/or degree of weathering with depth We tested the model’s sensitivity to the inclination of the
and also note that this decline is typically nonlinear (with impermeable boundary (see Table 9). Inclinations steeper
much of the permeability concentrated near the surface; e.g. than parallel to the ground surface had not been tested as
Jiang et al. 2009; Ameli et al. 2016) and/or characterised there is no evidence of this subsurface architecture in the lit-
by a sharp boundary between disturbed and fresh bedrock erature (Clair et al. 2015; Riebe et al. 2017; Anderson et al.
(Clair et al. 2015). 2013; Lebedeva and Brantley 2013; Flinchum et al. 2018;
Hydraulic conductivity and porosity are linked as they Hayes et al. 2019; Medwedeff et al. 2022). Bp (the thickness
are both influenced by the lithology type and the density of the permeable layer) was varied with layer inclination,
of fractures in the rock, meaning that their values used in a to ensure the impermeable boundary intersects the ground
model should be correlated. However, we did not correlate surface at the river channel.
these values and the value for porosity in the HSB model A steep increase in the Ff can be observed as the inclina-
was held constant. This was done so for the sake of simplic- tion of the impermeable boundary increases. This correla-
ity to demonstrate the overall methodology. The value for tion is due to increasing the gradient of the impermeable
porosity used in the HSB model was estimated as an aver- boundary resulting in an increased height of the phreatic
age of a range of porosity values for fractured crystalline surface leading to greater instability. Increasing the gradi-
metamorphic rock (5–10%) from Singhal and Gupta (2010). ent of the impermeable boundary results in two competing
The sensitivity of the model to the porosity was tested and drivers: (1) faster and, therefore, thinner lateral subsurface
the results are displayed in Table 8. The higher the value of flow reducing the height of the phreatic surface; and (2) an
porosity, the lower the frequency of failure. If the porosity increased height of the phreatic surface as the phreatic sur-
of the cutting was actually at the lower end of the range for face is perched on the impermeable boundary (increasing the
fractured crystalline metamorphic rock, the model output gradient of the impermeable boundary increases its height
would remain the same. Alternatively, if the porosity was above datum and thus increases the height of the phreatic
at the higher end, the model would overestimate the fre- surface above datum). Based on our results (Table 9), it can
quency of failure by 12%. These differences are relatively be said that the geometric effect of the height of the bound-
small compared to the other uncertainties examined here ary on which the subsurface flow is perched increasing has
suggesting that the model is not too sensitive to porosity for a greater effect on the phreatic surface level than that of the
this rock type, and thus it is acceptable to select a determin- gradient increasing the flow velocity. The outputs of this sen-
istic value for this parameter in this case. The range of poros- sitivity analysis are slightly biassed towards being sensitive
ity values for most rock types generally varies by 5–15% given the link between k and the inclination. Different com-
(with the exception of basalt and highly weathered crystal- binations of permeability and inclination of the permeable
line rock) (Singhal and Gupta 2010). With this in mind, the layer can lead to the same phreatic surface (Beven 1996).
model is likely to be insensitive to porosity across the range The 1st percentile of k was chosen based on the observation-
of porosity variability found in most rocks and, therefore, driven constraint that it should not result in frequent surface
it is broadly accepted to neglect porosity spatial variability. saturation for the slope. This constraint will be violated by
Given the lack of in situ hydrological data, we made the some of the realisations of k from the previous distribution
assumption that the impermeable boundary is at a lower under scenarios where the bedrock inclination is increased.
inclination than the ground surface, so that the permeable Conversely, the constraint will be over-conservative (and the
13
1st percentile of k over-estimated) for scenarios where the markedly seasonal wet–dry weather. The outcome of fre-
bedrock inclination is reduced. quency of failure methodology suggested that this cutting
In conclusion, this analysis shows that the methodology will fail every other year in its current condition. The time-
is particularly sensitive to the number of remediation days varying conditional probability of failure showed annual
and the 1st percentile of k (linked with the inclination of the cyclicity with the monsoon season (peak in probability of
impermeable boundary). Given that practitioners should be failure during the monsoon season) and reflected trends dis-
able to confidently estimate the number of remediation days, played in the phreatic surface time series.
the sensitivity of this input parameter is not a concern. Our results point to the need to calculate the phreatic sur-
face to be employed in the seepage analyses of slope stability
5.4 Use by Practitioners from a hydrogeological model rather than postulating the
phreatic groundwater level as often done in current prac-
The estimation of Ff can be very useful in cost–benefit anal- tise. This is for two reasons: (1) the phreatic level is highly
yses where different slope stabilisation measures are consid- dependent on the hydrogeology of the upslope which in turn
ered. For each measure, a different cost would be incurred is a function of rainfall and geometry of the catchment area;
and a value of Ff would be calculated. (2) the phreatic level exhibits significant variation over time.
Although this methodology is not costly nor challeng- Therefore, we believe a relatively simple model like the 1D
ing in terms of ground investigation, the bulk deterministic HSB as we adopted is a good model since it captures the key
analyses in LEM can be costly in terms of computational physical drivers of hydraulic flow within the slope without
time. To carry out the stability analyses (1000 analyses for being computationally expensive.
each phreatic surface, of which there were 25), 10 virtual These findings indicate that it is important to use a time-
machines were used, which each had 32GB and 8 cores, dependent system to represent rainfall variability for slope
running in parallel. Alternative methods could be used to failure probability analysis, rather than a time-independent
conduct the stability analyses using an analytical approach system, e.g. from an Intensity–Duration–Frequency curve.
and carrying out the probabilistic analysis using a first- Our sensitivity analyses show that the frequency of failure
order second-moment (FOSM) method (e.g. Huang 2021). is very sensitive to the value chosen to represent the 99th
However, in using an analytical approach, pore water pres- percentile of k in the distribution, the inclination of the
sure can only be included as a pore water pressure ratio (ru ) impermeable boundary, the number of remediation days
which could result in an unrealistic set of trigger scenarios. between landslide counts and the inclination of the imper-
For this approach, the G-H–B parameters would have to meable boundary. The number of remediation days can be
be converted to M-C parameters. This conversion can be better constrained by knowledge of the typical remediation
done using the approach presented by Renani and Martin measures employed locally following a failure event. How-
(2020) who propose a new equation for the conversion of ever, constraining the value for hydraulic conductivity may
parameters for slope stability analysis based on elastic stress be more difficult and, therefore care should be taken in the
analysis. Conversely, the hydrological model is very fast to result if the value for hydraulic conductivity is based on
run (2000 values of k takes < five minutes). assumptions.
Acknowledgements This research was supported by NERC IAPETUS

Doctoral Training Partnership [Grant No: NE/S007431/1].
6 Conclusion
Funding Funding was provided by Natural Environment Research
We introduced a novel computationally efficient methodol- Council (Grant No. NE/S007431/1)
ogy coupling probabilistic slope stability analyses with the Data availability The data that support the findings of this study are
Hillslope-Storage Boussinesq hydrological model to esti- available from the corresponding author upon reasonable request.
mate the time-varying probability of failure and frequency of
failure of a cutting triggered by rainfall. Unlike many other Declarations
studies, this method treats both the complexity of the failure
Conflict of Interest The authors have no competing interests to declare
mechanics and the hillslope hydrology. We aimed to further that are relevant to the content of this article.
understand the implications of using a time-dependent sys-
tem to represent rainfall variability (a rainfall time series) Open Access This article is licensed under a Creative Commons Attri-
for an area that hosts a monsoon and to further explore the bution 4.0 International License, which permits use, sharing, adapta-
tion, distribution and reproduction in any medium or format, as long
implications of using an instantaneous conditional probabil- as you give appropriate credit to the original author(s) and the source,
ity of failure. The mechanistic model is suitable for use in provide a link to the Creative Commons licence, and indicate if changes
a LIC/LMIC setting. We demonstrate the methodology on were made. The images or other third party material in this article are
a road cutting in Nepal in a mountainous area subject to a included in the article’s Creative Commons licence, unless indicated
13
otherwise in a credit line to the material. If material is not included in Fan Y, Bras R (1998) Analytical solutions to hillslope subsurface
the article’s Creative Commons licence and your intended use is not storm flow and saturation overland flow. Water Resour Res
permitted by statutory regulation or exceeds the permitted use, you will 34(4):921–927
need to obtain permission directly from the copyright holder. To view a Farichah H, Hutama D (2020) Probabilistic stability and sensitivity
copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. analysis of rock slope subjected to seismic loading. J Phys Conf
Ser 1517:1–6
Fenton G, Griffiths D (2008) Risk assessment in geotechnical engineer-
ing. John Wiley & Sons Inc, Hoboken
References Fine (2021) Unit weight of rocks. https://www.finesoftware.eu/help/
geo5/en/unit-weight-of-rocks-01/. Accessed: 19 Feb 2021
Flinchum B, Steven H, Rempe D, Moon S, Riebe C, Carr B, Hayes J,
Ameli A, McDonnell J, Bishop K (2016) The exponential decline in
St. Clair J, Peters M (2018) Critical zone structure under a granite
saturated hydraulic conductivity with depth: a novel method for
ridge inferred from drilling and three-dimensional seismic refrac-
exploring its effect on water flow paths and transit time distribu-
tion data. J Geophys Res Earth Surf 123(6):1317–1343
tion. Hydrol Process 30(14):2438–2450
Formetta G, Rago V, Capparelli G, Rigon R, Muto F, Versace P (2014)
Anderson R, Rajaram H, Anderson S (2019) Climate driven coevolu-
Integrated physically based system for modeling landslide suscep-
tion of weathering profiles and hillslope topography generates
tibility. Proc Earth Planet Sci 9:74–82
dramatic differences in critical zone architecture. Hydrol Process
Frattini P, Crosta G, Sosio R (2009) Approaches for defining thresh-
33(1):4–19
olds and return periods for rainfall-triggered shallow landslides.
Anderson RS, Anderson SP, Tucker GE (2013) Rock damage and rego-
Hydrol Process Int J 23(10):1444–1460
lith transport by frost: an example of climate modulation of the
Fredlund DG, Krahn J (1977) Comparison of slope stability methods
geomorphology of the critical zone. Earth Surf Process Landf
of analysis. Can Geotech J 14(3):429–439
38(3):299–316
Grant G, Dietrich W (2017) The frontier beneath our feet. Water Res
Baum RL, Savage WZ, Godt JW (2008) TRIGRS: a Fortran program for
Res 53(4):2605–2609
transient rainfall infiltration and grid-based regional slope-stability
Gravanis E, Pantelidis L, Griffiths D (2014) An analytical solution
analysis, version 2.0. US Geological Survey Reston, VA, USA
in probabilistic rock slope stability assessment based on random
Bellugi D, Milledge DG, Dietrich WE, Perron JT, McKean J (2015)
fields. Int J Rock Mech Min Sci 71:19–24
Predicting shallow landslide size and location across a natural
Griffiths D, Fenton G (2000) Influence of soil strength spatial vari-
landscape: application of a spectral clustering search algorithm. J
ability on the stability of an undrained clay slope by finite ele-
Geophys Res Earth Surf 120(12):2552–2585
ments. Slope stability 2000. American Society of Civil Engineers,
Beven K (1996) Equifinality and uncertainty in geomorphological
Denver, pp 184–193
modelling. In: The Scientific Nature of Geomorphology: Proceed-
Griffiths D, Fenton G (2009) Influence of spatial variability on slope
ings of the 27th Binghamton Symposium in Geomorphology. John
reliability using 2-D random fields. J Geotech Geoenviron Eng
Wiley & Sons
135(10):1367–1378
Beven K, Kirkby M (1979) A physically based, variable contributing
Guzzetti F, Peruccacci S, Rossi M, Stark C (2008) The rainfall inten-
area model of basin hydrology. Hydrol Sci J 24(1):43–69
sity-duration control of shallow landslides and debris flows: an
Boussinesq J (1877) Essai sur la théorie des eaux courantes. Impr.
update. Landslides 5:3–17
nationale
Hayes JL, Riebe CS, Holbrook WS, Flinchum BA, Hartsough PC
Caine N (1980) The rainfall intensity-duration control of shallow land-
(2019) Porosity production in weathered rock: where volu-
slides and debris flows. Geografiska annaler Ser A Phys Geogr
metric strain dominates over chemical mass loss. Sci Adv
62(1–2):23–27
5(9):eaao0834
Chen Y, Lin H (2019) Consistency analysis of Hoek–Brown and equiv-
Hess DM, Leshchinsky BA, Bunn M, Benjamin Mason H, Olsen MJ
alent Mohr–Coulomb parameters in calculating slope safety factor.
(2017) A simplified three-dimensional shallow landslide suscep-
Bull Eng Geol Environ 78(6):4349–4361
tibility framework considering topography and seismicity. Land-
Corominas J, Moya J (2008) A review of assessing landslide frequency
slides 14:1677–1697
for hazard zoning purposes. Eng Geol 102(3–4):193–213
Hilberts A, van Loon E, Troch P, Paniconi C (2004) The hillslope-stor-
Corominas J, van Westen C, Frattini P, Cascini L, Malet J-P, Fotopou-
age Boussinesq model for non-constant bedrock slope. J Hydrol
lou S, Catani F, Van Den Eeckhaut M, Mavrouli O, Agliardi F et al
291(3–4):160–173
(2014) Recommendations for the quantitative analysis of landslide
Hoek E (1998) Reliability of Hoek-Brown estimates of rock mass
risk. Bull Eng Geol Environ 73(2):209–263
properties and their impact on design. Int J Rock Mech Min Sci
Crovelli R (2000) Probability models for estimation of number and
35(1):63–68
costs of landslides. US Geological Survey
Holcombe E, Smith S, Wright E, Anderson MG (2012) An integrated
Dai F, Lee C, Ngai Y (2002) Landslide risk assessment and manage-
approach for evaluating the effectiveness of landslide risk reduc-
ment: an overview. Eng Geol 64(1):65–87
tion in unplanned communities in the caribbean. Nat Hazards
De Leon D, Garduño J (2020) Time-variant failure probability of
61:351–385
critical slopes under strong rainfall hazard including mitigation
Huang W (2021) Efficient analytical approach for stability analysis of
effects. Struct Infrastruct Eng 16(10):1481–1492
infrastructure slopes. In: Proceedings of the Institution of Civil
Dou H-Q, Han T-C, Gong X-N, Zhang J (2014) Probabilistic slope
Engineers-Geotechnical Engineering, pp 1–14
stability analysis considering the variability of hydraulic conduc-
Iverson R (2000) Landslide triggering by rain infiltration. Water Resour
tivity under rainfall infiltration–redistribution conditions. Eng
Res 36(7):1897–1910
Geol 183:1–13
Jiang X-W, Wan L, Wang X-S, Ge S, Liu J (2009) Effect of exponential
Duncan J (1996) State of the art: limit equilibrium and finite-element
decay in hydraulic conductivity with depth on regional groundwa-
analysis of slopes. J Geotech Eng 122(7):577–596
ter flow. Geophys Res Lett 36(24)
Dyson AP, Tolooiyan A (2019) Prediction and classification for finite
Jones J, Boulton S, Bennett G, Stokes M, Whitworth M (2021) Tempo-
element slope stability analysis by random field comparison.
ral variations in landslide distributions following extreme events:
Comput Geotech 109:117–129
13
Implications for landslide susceptibility modeling. J Geophys Res Phoon K, Kulhawy F (1999) Evaluation of geotechnical property vari-
Earth Surf 126(7):e2021JF006067 ability. Can Geotech J 36(4):625–639
Lebedeva MI, Brantley SL (2013) Exploring geochemical controls Renani H, Martin C (2020) Slope stability analysis using equivalent
on weathering and erosion of convex hillslopes: beyond the Mohr–Coulomb and Hoek–Brown criteria. Rock Mech Rock Eng
empirical regolith production function. Earth Surf Process Landf 53(1):13–21
38(15):1793–1807 Riebe CS, Hahm WJ, Brantley SL (2017) Controls on deep critical
Li A, Cassidy M, Wang Y, Merifield R, Lyamin A (2012) Parametric zone architecture: a historical review and four testable hypotheses.
Monte Carlo studies of rock slopes based on the Hoek–Brown Earth Surf Process Landf 42(1):128–156
failure criterion. Comput Geotech 45:11–18 Robson E, Agosti A, Utili S, Milledge D (2022) A methodology for
Li A, Merifield R, Lyamin A (2008) Stability charts for rock slopes road cutting design guidelines based on field observations. Eng
based on the Hoek–Brown failure criterion. Int J Rock Mech Min Geol 307:106771
Sci 45(5):689–700 Ross S (1972) Introduction to probability models. Academic Press
Liu X, Wang Y (2021) Reliability analysis of an existing slope at a Rouainia M, Davies O, O’Brien T, Glendinning S (2009) Numerical
specific site considering rainfall triggering mechanism and its past modelling of climate effects on slope stability. In: Proceedings of
performance records. Eng Geol 288:106144 the Institution of Civil Engineers-Engineering Sustainability, vol
Lu M, Zhang J, Zhang L, Zhang L (2020) Assessing the annual risk of 162, Thomas Telford Ltd, pp 81–89
vehicles being hit by a rainfall-induced landslide: a case study on Shakya NM, Nirula DR (2008) Integration of bio-engineering tech-
Kennedy Road in Wan Chai, Hng Kong. Nat Hazards Earth Syst niques in slope stabilisation works: a cost-effective approach for
Sci 20(6):1833–1846 developing countries. International Consortium on Landslides,
Lu M, Zhang J, Zheng J, Yu Y (2022) Assessing annual probability of United Nations University, Japanese Landslide Society, Tokyo,
rainfall-induced slope failure through a mechanics-based model. In First World Landslide Forum
Acta Geotech 17(3):949–964 Shokri S, Shademan M, Rezvani M, Javankhoshdel S, Cami B, Yacoub
Lu M, Zheng J, Zhang J, Huang H (2022b) On assessing the probability T (2019) A review study about spatial correlation measurement in
of rainfall-induced slope failure during a given exposure time. rock mass. In: ISRM Congress
Acta Geotech, pp 1–13 Singhal B, Gupta R (2010) Applied hydrogeology of fractured rocks.
Manning L, Hall J, Kilsby C, Glendinning S, Anderson M (2008) Spa- Springer Science & Business Media
tial analysis of the reliability of transport networks subject to rain- St. Clair J, Moon S, Holbrook W, Perron J, Riebe C, Martel S, Carr B,
fall-induced landslides. Hydrol Process Int J 22(17):3349–3360 Harman C, Singha K, deB Richter D (2015) Geophysical imaging
Marinos P, Hoek E (2000) GSI: a geologically friendly tool for rock reveals topographic stress control of bedrock weathering. Science
mass strength estimation. Int Soc Rock Mech Rock Eng, In ISRM 350(6260):534–538
international symposium Talebi A, Troch P, Uijlenhoet R (2008) A steady-state analytical slope
Massih D, Soubra A, Low B (2008) Reliability-based analysis of strip stability model for complex hillslopes. Hydrol Process Int J
footings using response surface methodology. Int J Geomech 8:2 22(4):546–553
Medwedeff WG, Clark MK, Zekkos D, West AJ, Chamlagain D (2022) Tang G, Huang J, Sheng D, Sloan S (2018) Stability analysis of
Near-surface geomechanical properties and weathering character- unsaturated soil slopes under random rainfall patterns. Eng Geol
istics across a tectonic and climatic gradient in the central Nepal 245:322–332
Himalaya. J Geophys Res Earth Surf 127(2):e2021JF006240 Troch P, Paniconi C, van Loon E (2003) Hillslope-storage Boussinesq
Montgomery DR, Dietrich WE (1994) A physically based model for model for subsurface flow and variable source areas along com-
the topographic control on shallow landsliding. Water Resour Res plex hillslopes: 1. Formulation and characteristic response. Water
30(4):1153–1171 Resour Res 39(11)
Mori H, Chen X, Leung Y, Shimokawa D, Lo M (2020) Landslide haz- Troch P, van Loon A, Hilberts A (2004) Analytical solution of the
ard assessment by smoothed particle hydrodynamics with spatially linearized hillslope-storage Boussinesq equation for exponential
variable soil properties and statistical rainfall distribution. Can hillslope width functions. Water Resour Res 40(8)
Geotech J 57(12):1953–1969 Troch P, van Loon E, Hilberts A (2002) Analytical solutions to a
Nour A, Slimani A, Laouami N (2002) Foundation settlement statistics hillslope-storage kinematic wave equation for subsurface flow.
via finite element analysis. Comput Geotech 29(8):641–672 Adv Water Resour 25(6):637–649
Ozturk U (2022) Geohazards explained 10: Time-dependent landslide van Zadelhoff FB, Albaba A, Cohen D, Phillips C, Schaefli B, Dor-
susceptibility. Geol Today 38(3):117–120 ren L, Schwarz M (2022) Introducing slideformap: a proba-
Pan Q, Jiang Y-J, Dias D (2017) Probabilistic stability analysis of a bilistic finite slope approach for modelling shallow-landslide
three-dimensional rock slope characterized by the Hoek–Brown probability in forested situations. Nat Hazards Earth Syst Sci
failure criterion. J Comput Civ Eng 31(5):04017046 22(8):2611–2635
Pandit B, Tiwari G, Latha G, Babu S (2019) Probabilistic characteri- Wang B, Ding Q, Fu X, Kang I-S, Jin K, Shukla J, Doblas-Reyes F
zation of rock mass from limited laboratory tests and field data: (2005) Fundamental challenge in simulation and prediction of
associated reliability analysis and its interpretation. Rock Mech summer monsoon rainfall. Geophys Res Lett 32(15)
Rock Eng 52(9):2985–3001 Wyllie D, Mah C (2017) Rock slope engineering. CRC Press
Paniconi C, Troch P, van Loon E, Hilberts A (2003) Hillslope-storage Zeng P, Senent S, Jimenez R (2014) Reliability analysis of circular
Boussinesq model for subsurface flow and variable source areas tunnel face stability obeying Hoek–Brown failure criterion con-
along complex hillslopes: 2. Intercomparison with a three-dimen- sidering different distribution types and correlation structures. J
sional richards equation model. Water Resour Res 39(11) Comput Civ Eng 30(1):04014126
Parkin T, Robinson J (1992) Analysis of lognormal data. Advances in Zhang J, Huang H, Zhang L, Zhu H, Shi B (2014) Probabilistic pre-
soil science. Springer, New York, pp 193–235 diction of rainfall-induced slope failure using a mechanics-based
Parsons T (2002) Global Omori law decay of triggered earthquakes: model. Eng Geol 168:129–140
large aftershocks outside the classical aftershock zone. J Geophys
Res Solid Earth 107(B9):ESE–9 Publisher's Note Springer Nature remains neutral with regard to
Phoon K, Kulhawy F (1999) Characterization of geotechnical vari- jurisdictional claims in published maps and institutional affiliations.
ability. Can Geotech J 36(4):612–624
13

s00603 023 03694 5

Uploaded by

Copyright:

Available Formats

s00603 023 03694 5

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

s00603 023 03694 5

Uploaded by

Copyright:

Available Formats

Rock Mechanics and Rock Engineering (2024) 57:2421–2443

A Computationally Efficient Method to Determine the Probability

1 Introduction average slope conditions over a region. Statistical models

Fig. 1 Flowchart outlining the

Table 1 Table of all parameter Parameter symbol Description

Fig. 2 Schematic diagram

Fig. 4 Illustration of the Hillslope-Storage Boussinesq slope pro-

Fig. 5 A schematic diagram

Fig. 6 Map showing the loca-

The Hoek–Brown failure criterion include: Geological Strength Index

Fig. 7 Image of road cutting case study along the Narayanghat-Mug-

The slope cutting of exposed rock mass is characterised

Table 4 Table displaying values for Generalised Hoek–Brown param-

One thousand values of Geological Strength Index (GSI) are sampled

Fig. 9 Lognormal distribution of the Geological Strength Index cat-

Given that the model exhibits high sensitivity to the GSI

Table 5 Case study input parameters for the HSB model

Depth to Inclina- Shape Initial Length of Material

The inclination of impermeable boundary is assumed to be less than

3.3 Site‑Specific Model Coupling

To couple the outputs from the probabilistic (the MCS)

Figure 13 displays probability density distributions of 4.1 Time‑Varying Conditional Probability of Failure

Fig. 15 Histogram of the number of landslides in 11 years for the

By summing all the landslides from each FoS time series,

Acknowledgements This research was supported by NERC IAPETUS

You might also like

3.3 Site‑Specific Model Coupling

Figure 13 displays probability density distributions of 4.1 Time‑Varying Conditional Probability of Failure