Computers and Geotechnics 188 (2025) 107558

Contents lists available at ScienceDirect

Computers and Geotechnics

journal homepage: [Link]/locate/compgeo

Research Paper

Analysis of the time-varying reliability of earth dam slopes considering

multiparameter spatial variability
Zhi-Ping Deng a,b , Yi Zou c , Wenqi Du a , Dian-Qing Li a,*, Mao-Xin Wang d ,
Shui-Hua Jiang c
a
State Key Laboratory of Water Resources Engineering and Management, Institute of Engineering Risk and Disaster Prevention, Wuhan University, Wuhan 430072,
China
b
College of Water Conservancy, Jiangxi University of Water Resources and Electric Power, Nanchang 330099, China
c
School of Infrastructure Engineering, Nanchang University, Nanchang 330031, China
d
Department of Civil and Environmental Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, China

A R T I C L E I N F O A B S T R A C T

Keywords: Reliability analysis is an effective method for comprehensively assessing earth-rock dam safety, considering both
Slope stability soil spatial variability and reservoir water level fluctuations. However, previous studies primarily focused on
Spatial variability static analyses at specific moments, and were not applicable to problems involving multiple high-dimensional
Time-varying water level
random variables. Therefore, a probabilistic framework integrating slice inverse regression (SIR), Gaussian
Parameter sensitivity
Surrogate model
process regression (GPR), and Monte Carlo simulation (MCS) is proposed for assessing the time-varying reli­
ability of earth dam slopes. This framework consists of three key components: (1) The Karhunen–Loève (KL)
expansion is adopted to discretize the random field (RF) of high-dimensional multivariate soil parameters, and
the dimension of the random field variables is further reduced utilizing a SIR method; (2) An efficient GPR
surrogate model is constructed based on the reduced-dimension variables; (3) The MCS method is employed to
estimate failure probability (Pf) of earth dam slopes. The proposed framework is illustrated using the Ashigong
earth dam case study. Results demonstrate that the SIR-GPR-MCS framework effectively characterizes multi-
parameter spatial variability (SV) and estimates time-dependent Pf with acceptable computational efficiency.
Parametric analyses investigating soil parameter statistics and the time-varying statistical distribution of safety
factors (FS) are also conducted. These analyses provide valuable insights for the design and maintenance of dam
slope stability under fluctuating reservoir water levels.

1. Introduction et al., 2023). Traditional deterministic analysis methods cannot suffi­

ciently account for these potential uncertainties in geotechnical stability
Earth dams represent a pivotal dam type in contemporary hydraulic assessments (Deng et al., 2022; Jiang et al., 2022). Reliability analysis
project construction and play a vital role in water management and methods, which consider uncertainties and can be used to evaluate
control (Su et al., 2022). The Ashigong earth dam, located on the project safety and viability through a probabilistic manner (Cho, 2010;
mainstream of the Yellow River, is characterized by challenging condi­ Ji et al., 2021a; Griffiths and Fenton, 2004), offer a robust framework for
tions because of the influence of a fine sand layer and a gravel layer quantitatively analysing the impact of uncertain factors on engineering
within its foundation. Geotechnical survey data indicate that the dam safety. In recent years, reliability analyses of earth dam slopes and
construction conditions were nonoptimal. In the event of an instability corresponding disaster risk control methods have become popular
disaster, the repercussions could be severe (Gu et al., 2023). There are research topics in the geotechnical engineering field; notably, it is
many inevitable uncertainties in engineering practice, such as strati­ essential to propose disaster mitigation measures and prevent instability
graphic variability, spatial variability of soil properties (SV), and mea­ catastrophes associated with hydraulic projects (Cui et al., 2022).
surement uncertainty (Cao et al., 2016; Phoon and Kulhawy, 1999; Deng Numerous scholars have conducted valuable research for

* Corresponding author.
E-mail addresses: dengzhiping@[Link] (Z.-P. Deng), zouyi@[Link] (Y. Zou), wqdu309@[Link] (W. Du), dianqing@[Link] (D.-Q. Li),
maoxinwang@[Link] (M.-X. Wang), sjiangaa@[Link] (S.-H. Jiang).

[Link]
Received 20 March 2025; Received in revised form 14 July 2025; Accepted 5 August 2025
Available online 11 August 2025
0266-352X/© 2025 Elsevier Ltd. All rights are reserved, including those for text and data mining, AI training, and similar technologies.
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

characterizing soil parameter uncertainty (Jiang et al., 2020a), exam­ to predict the compressive strength of stabilized pond ash samples.
ining reliability under complex load conditions (Li et al., 2022), and Among these methods, the Gaussian process regression (GPR) technique
developing efficient reliability analysis methods (Liu and Cheng, 2016), stands out as particularly pervasive in engineering applications, span­
providing insights and tools for unsaturated slope reliability analyses. ning structural reliability analysis (Li et al., 2021; Ray et al., 2021; Zhu
However, the existing studies generally concentrated solely on the static et al., 2020), rock bearing capacity prediction (Fakhri et al., 2023; Chen
reliability of unsaturated slope stability at particular time points. To and Zhang, 2022), geosynthetic settlement prediction (Raja et al.,
address this gap, Liao and Ji (2021) discretized the rainfall process to 2023), the prediction of lateral displacement due to seismic liquefaction
formulate an instantaneous limit state function for a slope at any (Ahmad et al., 2022), and various other geotechnical engineering tasks
moment in time and combined it with a reliability method to assess time- (Fathipour-Azar, 2021; Mahmoodzadeh et al., 2021; Shuku and Phoon,
varying reliability. Song et al. (2019), leveraging generalized probabil­ 2023; He et al., 2021). The widespread geotechnical applications of
ity density evolution techniques and seismic dynamic stability analysis, machine learning methods can be attributed to its ease of implementa­
developed a time-varying reliability analysis approach for three- tion, flexible nonparametric inference method, and adaptive determi­
dimensional slopes. Hu et al. (2023) combined the Monte Carlo simu­ nation of hyperparameters. However, to the authors’ knowledge, the
lation (MCS) method with a modified Newmark method to assess both application of the GPR method in the time-varying reliability analyses of
the reliability and displacement of a slope under the influence of sto­ unsaturated earth dam slopes considering multiparameter SV is rela­
chastic earthquake dynamics. Considering soil strength degradation, Ji tively rare. Moreover, the GPR model and traditional surrogate models
et al. (2021b and 2023) employed dynamic stability analysis to assess may be influence by the curse of dimensionality (Guo et al., 2019),
the displacement and dynamic steady state of slopes subjected to the which means that the computational cost and the difficulty of statistical
combined effects of rainfall infiltration and seismic shaking. While the analyses increase exponentially as the dimensionality of the problem
previous works extensively investigated the seismic response and reli­ increases. This not only drastically increases the computational cost but
ability of embankments and dams during earthquakes and rainfall also affects the accuracy of results. To address this issue, Duan et al.
events, limited attention has been given to exploring the time-varying (1991) introduced the slice inverse regression (SIR) method for dimen­
reliability of earth dam slopes subjected to fluctuations in reservoir sionality reduction. This methodology posits that the essential infor­
levels. mation of the response variable resides in directional variations within
Manipulating the flow rate to discharge water in advance of flood the predictor space resulting from slice transformations. By identifying
waves can reduce the water level in reservoirs, thereby altering both the and extracting these maximally informative directions, high-
internal status (pore water pressure and permeability) and external dimensional data can be projected into a low-dimensional subspace −
loadings (hydrostatic pressure and hydrodynamic effects) in earthen the sufficient dimension reduction (SDR) subspace − spanned by a basis
dam system assessments. The disparity between internal and external of statistically significant directions. The subspace preserves all suffi­
pressures can lead to changes in dam reliability. Siacara et al. (2020) cient structural information relevant to the response variable from the
investigated the time-varying reliability of earth dams by setting varying original space, enabling both efficient data compression and targeted
rates of water level decline and employing a first-order reliability information extraction. The SIR method has gained widespread accep­
method. Liao et al. (2020) conducted dry–wet cycle tests to simulate tance because of its robust reliability and ease of implementation.
reservoir operation-induced water level fluctuations and then analysed However, despite the proven efficacy of SIR in handling high-
changes in the seepage field and reliability in the context of the Majia­ dimensional problems, its application in the specific context of time-
gou landslide. Guardiani et al. (2021) developed two surrogate models varying reliability analyses of unsaturated earth dam slopes with
and integrated them with the MCS method to assess the reliability of multiparameter SV remains a research area that warrants exploration.
reservoir bank slopes under various water-level decline rates. Although This study develops a SIR-GPR-MCS method for efficient time-
the above studies evaluated the time-varying reliability of slopes in varying reliability analysis of earth dam slopes considering multipa­
response to reservoir water level variations, the random variable model rameter SV. Initially, an RF is used to simulate the SV of soil, and the
without representing the possible spatial variation of soil parameters Karhunen–Loève (KL) expansion method is then applied for RF dis­
was often adopted. cretization. The high-dimensional random variables resulting from RF
Previous studies have demonstrated that the SV of soil parameters discretization are then downscaled with the SIR method introduced by
have significant effects on failure mechanisms and risk of geotechnical Duan et al. (1991). The GPR model is subsequently established, and the
systems (Li et al., 2015; Griffiths and Fenton, 2004). Disregarding the SV MCS method is employed to analyse the reliability of earth dam slopes. A
of soil parameters may lead to inaccurate estimates of geotechnical real engineering example, namely the Ashigong earth dam, is used to
structure reliability. As a well recognized approach to describe the SV, illustrate the proposed SIR-GPR-MCS approach. Meanwhile, the impact
random field (RF) theory has been extensively applied in geotechnical of the statistical characteristics of soil parameters on failure probability
reliability and risk assessments (Nguyen et al., 2022; Jiang et al. 2022; (Pf) of the Ashigong earth dam is investigated.
Zhang et al. 2022). Thus, it is imperative to introduce RF theory to ac­
count for the SV of soil parameters in time-varying reliability analyses of 2. Deterministic analysis
earth dam slopes. Notably, the time-varying reliability issue can be
viewed as a series of static problems at distinct time intervals (Wang 2.1. Seepage analysis
et al., 2020), entailing greater complexity than conventional static
reliability analysis. In addition, due to the generally unsaturated natural Seepage impacts the regular operation of earth dams invariably. A
state of soils, incorporating the SV of shear strength and hydraulic pa­ profound understanding of key seepage parameters (e.g., seepage vol­
rameters simultaneously introduces a large number of random variables, ume, flow rate, infiltration line, and pore water pressure) is essential to
further increases the complexity of the earth dam slope reliability ensure the stable operation of earth dams and to achieve effective flood
regression analysis (Tan et al., 2017; Jiang et al., 2020b). Therefore, it is prevention and mitigation.
necessary to develop a new time-varying reliability analysis method for Compared with saturated seepage in conventional scenarios, seepage
the unsaturated earth dam slope considering multiparameter SV. in earth dam cases with time-varying water levels is a more intricate
In recent years, machine learning techniques have garnered challenge, involving both saturated and unsaturated states. Richards
increasing attention and been frequently applied in geotechnical engi­ (1931) formulated a two-dimensional differential equation, thereby
neering. For instance, Qi et al. (2021) employed the multivariate establishing the foundation for the analysis of unsaturated seepage:
adaptive regression spline method to predict stratigraphic structural
surfaces. Suthar et al. (2019) utilized various machine learning methods

2
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

( ) ( )
∂ ∂H ∂ ∂H ∂θ ∂H (4), the outcomes of seepage analysis are of paramount importance for
k + k = γw (1)
∂x x ∂x ∂y y ∂y ∂ψ ∂t stability analysis. In this study, transient seepage and stability analyses
were performed via GEO-STUDIO software (GEO-SLOPE International
where kx and ky denote the horizontal and vertical permeability co­ Ltd., 2021). The time step was strategically set to record the time-
efficients, respectively; H is the total head; ψ is the matric suction of soil, varying FS. These deterministic analyses serve as the foundation for
which is obtained by subtracting the water pressure uw from the air subsequent reliability calculations. Additionally, the uncertainty prop­
pressure ua; γ w is the weight of water; t denotes the time; and ∂H/∂t=0 erties of the soil (e.g., the SV of the parameters) also need to be
indicates steady seepage. Usually, k and θ for saturated soil are constant reasonably characterized, as discussed in the next section.
(i.e., ks and θs, respectively), whereas the seepage capacity of unsatu­
rated soils usually varies with the magnitude of the pore water pressure 3. SV of the soil parameters simulated via an RF
(PWP). The relationship between θ and ψ in unsaturated soils is char­
acterized by the soil–water characteristic curve (SWCC), and various The spatial variability of soil, manifested as an inherent character­
SWCC models have been proposed by scholars over the past few decades. istic arising from depositional process, underlies the evolution of engi­
In this study, the van Genuchten (VG) model was employed (Wang et al., neering risks and failure mechanisms. Neglecting this variability may
2018), with the following expression: result in misidentification of both the critical slip surface and probability
[ ]m of failure (Li et al., 2015; Wang et al., 2020; Liao and Ji, 2021). The RF
θ − θr 1
Se = = n (m = 1 − 1/n, n > 1) (2) model proposed by Vanmarcke et al. (1986) has garnered widespread
θ s − θr 1 + (ψ /a)
attention owing to its ability to characterize the SV of soil parameters.
where Se is the degree of saturation; θs and θr are the saturated and re­ First, random variables are generated on the basis of the statistical in­
sidual θ values, respectively; and a, n, and m are the fitting parameters of formation of soil parameters. A set of ordered random variables is then
the VG model. a is the inlet value, where a small soil particle size and integrated following a specific rule to construct an RF. The RF is sub­
pore size correspond to a large a; n is a parameter related to the outflow sequently mapped onto an FEM mesh to determine the soil parameters
rate function of the soil when ψ reaches the inlet value of the soil, and n within each element of the geotechnical structural model. This entire
increases as the pore size distribution of the soil becomes increasingly process is known as RF simulation, with RF discretization commonly
uniform. m is related to the soil parameter θr, which reflects the overall implemented via methods such as midpoint and series expansion.
symmetry of the SWCC. The relationship between k and ψ for unsatu­ In this study, the KL expansion method is selected for RF dis­
rated soil is expressed by the permeability coefficient function (PCF), cretization because of its conceptual simplicity, clear principles, and
which corresponds to the VG model with the following functional widespread applicability. The autocorrelation function is used to
expression: quantify the correlations among arbitrary soil parameters, whereas the
autocorrelation distance characterizes the degree of spatial variation in
k(ψ ) [ ( )m ]2 soil parameters. A higher autocorrelation distance signifies stronger
kr (ψ ) = = S1/2 1 − 1 − S1/m (3)
ks e e
spatial similarity between parameters at two positions. However, owing
to the limited availability of geotechnical data in engineering practice,
where kr is the relative hydraulic conductivity. The unsaturated prop­ empirical autocorrelation functions with site-specific statistics are usu­
erties of soil are jointly influenced by several parameters in the SWCC ally used instead of accurately determining the spatial autocorrelation
and PCF cases, so accurately characterizing the uncertainties of these structure of the parameters. Li et al. (2015) summarized some
parameters is fundamental for seepage and stability analyses of earth commonly used autocorrelation functions and ranges of autocorrelation
dams. distances in geotechnical engineering reliability analysis. Among them,
Given the highly nonlinear characteristics of the aforementioned the Gaussian autocorrelation function is employed in this study:
control functions, direct solutions via analytical methods are chal­ ( [( )2 ( )2 ] )
lenging to obtain. Various numerical methods have been employed for |x − x ʹ| |y − y ʹ|
ρ(x, y) = exp − + (5)
saturated–unsaturated seepage analysis, with the finite element model lx ly
(FEM) method being one notable example. The PWP and θ values ob­
tained via seepage analysis are subsequently incorporated into slope where (x, y) and (x’, y’) are the coordinates of any two points in the RF
stability analyses, as detailed in the following subsection. space and lx and ly are the horizontal and vertical autocorrelation dis­
tances, respectively. The analysis area is discretized into multiple RF
elements to construct a correlation matrix. The matrix is subjected to
2.2. Stability analysis Cholesky decomposition, yielding a lower triangular matrix. This matrix
is then multiplied with a random sample matrix to generate the corre­
In geotechnical engineering practice, the safety factor (FS) is a lated standard Gaussian RF. A correlated non-Gaussian RF is finally
commonly utilized metric for evaluating the stability of structures or simulated through an equiprobability transformation, and the detailed
systems. The Morgenstern-Price method (Morgenstern and Price, 1965) procedure can be referred to the study by Jiang et al. (2014).
is used to calculate the FS in this study, belonging to the range of limit This section briefly introduces the implementation steps of the KL
equilibrium method (LEM). The shear strength of the soil directly affects expansion method. As the KL method has been extensively discussed in
the stability of the earth dam. Integrating dual-stress-state variables with the literature (Tan et al., 2017; Deng et al., 2021), the details of this
the Mohr–Coulomb strength criterion, Fredlund et al. (1978) derived the method are not given in this paper. In this study, φ’, ks, a, and n for an
shear strength formula for unsaturated soil, expressed as: earth dam were simulated as RFs and subsequently integrated into the
τ = c ʹ + (σ − ua )tanφ ʹ + (ua − uw )Se tanφ ʹ (4) reliability analysis.

where τ represents the shear strength, σ denotes the normal stress, and c’ 4. The SIR-GPR-MCS reliability analysis method
and φ’ correspond to the effective cohesion and internal friction angle,
respectively. Mechanical stress within soils is defined by the net normal Reliability analyses of unsaturated earth dam slopes involve three
stress variable (σ − ua), whereas internal stress induced by capillary primary challenges: (1) the reliability of the earth dam slope under
action and adsorption effects is characterized through the matrix suction seepage is affected by the SV of φ’, ks, a, and n, so the discretization of
variable (ua − uw). Hence, gas–liquid interactions in unsaturated soils the RFs typically leads to a high-dimensional problem; (2) the direct
may be rigorously quantified by this strength model. As delineated in Eq. application of the MCS method for reliability assessment entails

3
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

numerous intricate FEM analyses, resulting in a substantial computa­ (6) Sort the eigenvalues ̂ λ i , and the eigenvector βi (i = 1…, s) cor­
tional burden, whereas traditional surrogate models exhibit poor accu­ responding to the first sth largest eigenvalues is identified as the pro­
racy when solving highly nonlinear problems of this type; and (3) jection vector in the effective reduced-dimension space;
existing reliability analyses of earth dam slopes commonly focus on only (7) Calculate the projection value Xi = [Xi1, Xi2…, Xis] (i = 1…, N) for
the static reliability problem at a certain moment and lack thorough each slice of ξi in the effectively reduced dimension space. These
consideration on the time-varying reliability problem. To address the correspond to the low-dimensional random variables of the soil pa­
abovementioned problems, a SIR-GPR-based surrogate model is pro­ rameters after dimensionality reduction with the SIR method.
posed as an approximate substitute for the real function associated with The SIR method is founded on the concept of inverse regression,
the FS of earth dams. Notably, this approach helps avoid ’dimensional which inverts the traditional approach by regressing x on y instead of the
disasters’. In combination with the MCS method, a time-varying reli­ more common y on x. The primary advantage of this switch is that it
ability analysis of earth–rock dams is conducted, and the changes in circumvents the issue of high-dimensionality. The inverse model-based
earth–rock dam slopes Pf under different water level drop conditions are regression enables the regression of each component of x on y, effec­
investigated. tively transforming a complex, high-dimensional regression problem
into a series of simpler, one-dimensional regression problems. SIR delves
into the inverse regression curve of the conditional expectation E (x| y)
4.1. SIR method for dimensionality reduction preprocessing
to understand how x changes in response to y. Throughout the SIR
process, the response variables y are treated as independent variables,
Characterized by its model-free nature and robustness to error dis­
while the input variable x is considered the dependent variable. The SIR
tributions, the SIR method demonstrates strong universality while
method is theoretically robust, and a comprehensive explanation of the
exerting minimal impact on the accuracy of integrated surrogate
’inverse’ process, along with further proofs of the method, can be found
models. It is a technique that projects independent variables into an
in the relevant literature (Duan et al, 1991).
effective dimensionality reduction space, yielding a new set of mutually
orthogonal composite variables. The new composite variables are sub­
4.2. Gaussian process regression method
sequently regressed on the basis of the dependent variable. This
approach preserves the main features of the data, reduces dimension­
The GPR (Rasmussen and Williams, 2006) method was developed on
ality, and enhances the accuracy of regression. In the context of earth
the basis of Bayesian theory, offering a key advantage in providing
dam slope reliability, we assume that the original data model is as fol­
regression results with probabilistic significance. This approach excels
lows:
in addressing small-sample-size problems. Furthermore, its hyper­
Y = f(ξ) (6) parameters are adaptively optimized by maximum likelihood estimation
(MLE), eliminating the need for repetitive cross-validation across
Here, ξ represents high-dimensional random variables obtained from the hyperparameter combinations, thereby conserving computational re­
RF discretization of the soil parameters, given by ξ = [ξ1, ξ2…, ξN]T, sources and time. The construction of the GPR model is based on the
where N is the FEM capacity. Y is the matrix of the earth dam FS, matrix X of the soil parameters after dimensionality reduction and the
expressed as Y = [Y1, Y2…, YN]T. Additionally, ξi = [ξi1, ξi2…, ξir] (i = matrix Y of the safety factors. A training set D= (X, Y) = [(Xi, Yi)] with N
1…, N), and r is the overall dimension of the variables after discretiza­ sets of observations, where (i = 1…, N), is assumed. For a new input X*,
tion of the RF of the soil parameters. Following SIR dimension reduction, the objective of the GPR model is to predict the output value Y* corre­
the data model is transformed to: sponding to X* on the basis of prior knowledge. Assuming the presence
( ) of noise in a training set that adheres to a standard Gaussian distribu­
Y ≈ g(X) = f βT1 ξ, βT2 ξ, ..., βTs ξ, ε (7)
tion, the regression model under these assumptions is established as:
where X = [X1, X2…, XN]T, Xi = [Xi1, Xi2…, Xis] (i = 1…, N) is the
( )
Y = g(X) + ε, ε ∼ N 0, σ2N IN (8)
mapping variable after dimensionality reduction, s is the reduced
dimension, βi (i = 1…, s) is the projection vector and ε is a random where Y = [Y1, Y2, …, YN]T is the N × 1-dimensional output variable; X
perturbation. Eq. (7) clearly shows that through the SIR method of = [X1,s, X2,s, …, XN,s]T is the N × s-dimensional input random variable,
dimensionality reduction, the original r-dimensional independent vari­ with s representing the dimensionality of the random variable after
able is projected by the projection vector βi into the s-dimensional dimensionality reduction; ε is the observation noise vector, character­
effective dimensionality reduction space, where s < r. The information ized by a mean of zero and a variance of σN2; and IN is the unit matrix.
in the reduced dimensional space is more concise and effective, making The GPR function g(X) is defined as:
the process of regression analysis simpler. The objective of SIR analysis
is to discover the projection vector that forms the effective dimension­ g(X) ∼ N(m(X), K(X, X) ) (9)
ality reduction space, achieving the purpose of dimensionality reduc­
m(X) = E[g(X) ] (10)
tion. The specific operation is as follows.
(1) Standardize ξi; [ ]
K(X, X) = E (g(X) − m(X) )(g(X) − m(X) )T (11)
(2) Sort Yi and divide the results into D nonoverlapping interval slices
Id (d = 1…, D). Sort ξi into slices as well. Notably, the SIR method
Here, m(X) represents the mean function of the input sample X, which is
converges and exhibits asymptotic normality when the number of slices
√̅̅̅̅ typically assumed to be zero, and K(X, X) is the covariance function of
is in the range of N to N/2. Duan et al. (1991) also highlighted that the
the input sample X. K(X, X) is a N × N-dimensional matrix used to
number of slices has minimal impact on the dimensionality reduction
quantify the correlation among the variables. Combining Eq. (8) and Eq.
outcome of the SIR method;
(9), the prior distribution of the observed target value Y is
(3) Calculate the mean value ξd (d = 1…, D) of each slice of ξi and the ( )
Y ∼ N 0, K(X, X) + σ 2N IN .
mean value Yd (d = 1…, D) of each slice of Yi after sorting; The joint Gaussian prior distribution, formed on the basis of the
(4) Calculate the weighted covariance matrix, accounting for po­ observed values Y of the N training samples and the output values Y*
tential differences in sample sizes among slices; (also known as the predicted values) of the M test samples, is given by:
(5) Conduct principal component analysis (PCA) on the mean of each
slice ξd , and derive the eigenvalues ̂ λ i and corresponding eigenvectors
̂ i (i = 1…, r) on the basis of the weighted covariance matrix;
α

4
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

{ [ ]}
[
Y
]
K(X, X) + σ 2N IN K(X, X* ) dam slope at moment t; Xi=(φ’, ks, a, n) is the realized value of the soil
∼ N 0, (12) parameter RF for the ith FEM sample; Pf (t) is the earth dam slope Pf at
Y* K(X, X )* T
K(X , X )
* *
moment t; and I{⋅} is the indicator function, with I{⋅} = 1 for FS (Xi, t) <
1 and zero otherwise.
Here, K(X, X*) is the N × M− dimensional covariance matrix, encom­
passing all inputs values X of the training set and the input values X* of
5. Implementation steps
the test sample, and K(X*, X*) represents the covariance matrix of the
test sample. The covariance function type chosen in this paper is the
To help engineers understand and apply the SIR-GPR-MCS reliability
commonly used squared exponential type, expressed as:
analysis methodology, an analysis flowchart is provided (Fig. 1), and it
( )
( ) 1( )2 encompasses four key components: deterministic analysis, RF dis­
K Xe , Xf = σ2f exp − Xe − Xf + σ 2N δef (13)
2l2 cretization, SIR-GPR model construction and calibration, and the pre­
diction of Pf. The initial step involves providing essential information,
where Xe and Xf are two arbitrary elements of a sample input value, σ f2 is such as boundary conditions and statistical characterizations of soil
the signal variance, l is a length scale parameter, σ N2 is the noise vari­ parameters. GEO-STUDIO software is employed to model unsaturated
ance, and δef is the Kronecker function. The parameters l, σf2, and σ N2 are earth dams, and coupled flow-solidity analysis is performed to ascertain
collectively referred to as the three hyperparameters of the squared the slope’s stability state. The FEM calculation file is subsequently
exponential covariance function, and the optimal solution can be exported. The SV of the soil parameters is characterized by the RF and
adaptively obtained through the maximum likelihood estimation discretized via the KL method. The parameters in the original FEM file
method: are rewritten with a MATLAB program, generating a new input file. A
1 T( )− 1 1 ⃒ ⃒ N WINDOWS BAT batch script is used to invoke the GEO-STUDIO software
L = lgp(Y|X) = − Y K + σ 2N IN Y − lg⃒K + σ 2N IN ⃒ − lg(2π ) (14) kernel to perform analyses of the new input file, with each multipa­
2 2 2
rameter random simulation of the soil corresponding to a time-varying
The training set is learned, and the probability distribution of the pre­ sequence of FSs. At each time step, the random variables associated
dicted values is continuously updated with new observations according with soil parameters and the corresponding FS values are subsequently
to the Bayesian method. The maximum possible predictive distribution separated into a training set and a test set. These sets are employed to
of Y* is inferred from the new X*, X, and Y values: develop the SIR-GPR model, which predicts the time-varying FS of the
earth dam slope. Two error assessment metrics are used to evaluate the
p(Y * |X* , X, Y) ∼ N(m(Y * ), K(Y * , Y * ) ) (15) predictive performance of the model at each moment. Finally, the MCS
method is applied to calculate the corresponding Pf at each moment, and
where the specific expressions for the mean m(Y*) and variance K(Y*,
the time-varying reliability trend of earth dam slopes is subsequently
Y*) of Y* are as follows:
fitted.
⎧ [ ]− 1 To enhance the practical applicability of the suggested SIR-GPR-MCS
⎨ m(Y * ) = K(X, X* )T K(X, X) + σ 2N IN Y
(16) method, the aforementioned execution steps are encapsulated as func­
⎩ K(Y * , Y * ) = K(X* , X* ) − K(X, X* )T [K(X, X) + σ 2 I ]− 1 K(X, X* )
N N tions and toolboxes. These can be treated as black boxes, enabling
geotechnical engineers to seamlessly execute reliability analysis tasks in
4.3. Error assessment future projects.

Although replacing the original time-consuming FEM with the sur­ 6. Application of the proposed method to the Ashigong earth
rogate model significantly reduces computational costs, it is worth dam
emphasizing that the surrogate model is an approximation and that
there will be deviations from real FEM simulations. To quantitatively 6.1. Overview and deterministic analysis of the project
assess the predictive ability of the SIR-GPR model, two evaluation
indices, namely, the root mean square error (RMSE) and coefficient of In this section, the SIR-GPR-MCS method is utilized for conducting a
determination (R2), are employed in this study. time-varying reliability analysis of the Ashigong earth dam, which is
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ located in the upstream section of the Yellow River mainstream in Guide
√N
√∑ County, Qinghai Province, China. It represents the most downstream
√ (Yi − yi )2
√ stage in the recent four-stage navigation hub development program for
RMSE = i=1 (17)
N the Guide section of the Yellow River. With a total installed capacity of
180 MW, a normal storage level of 2198 m, a corresponding reservoir
∑N
1
i=1 (Yi − Y)2 capacity of 0.6836 billion m3, and an average annual power generation
R2 = 1 − N
(18)
1
∑N
− y)2 totalling 698.52 million kW⋅h, the power station falls into the category
N i=1 (yi
of a medium-sized level-III project. The dam body, which is filled with
where yi is the true FS calculated by the ith FEM and Yi is the FS predicted sand and gravel, boasts a crest elevation of 2202 m and a bottom
by the SIR-GPR model for the ith earth dam slope. elevation of 2179 m. The main structures of the hub, which are arranged
from left to right, include a sand and gravel dam with clay-inclined
4.4. MCS method for assessing Pf walls, a lock, a spillway gate, a concrete connecting section, a power
plant, and a connecting section to the earth dam on the right bank. The
A time-varying reliability analysis is conducted by integrating the dam line stretches for 841.5 m and is arranged in a zigzag pattern.
MCS method into the above-described approach. The earth dam slope Pf Geological exploration data indicate a substantial depth for dam
at the tth moment under the action of a time-varying water level is given bedrock. Considering the 50 m thickness of the dam foundation
by: (elevation range of 2129 ~ 2179 m), the dam foundation exhibits a
three-layer structure: a top layer of sandy gravel, a central part with a
Nsim
1 ∑ general distribution of sand-containing silt and fine-sand layers, and
Pf (t) = I{FS(Xi , t) < 1 } (19)
Nsim i=1 silty clay layers distributed in the form of interlayers. The bottom fea­
tures a continuous layer of sand and gravel. The stratigraphic structure
where Nsim is the number of MCS samples; FS (Xi, t) is the FS of the earth and distribution of the dam are shown in Fig. 2.

5
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 1. Flowchart of the SIR-GPR-MCS reliability analysis method.

Fig. 2. Ashigong earth dam profile.

To simulate the effect of discharged floodwater on the stability of the decrease in the total head is considered upstream; a constant total head
earth dam, four upstream water level drawdown scenarios were artifi­ located 2 m above the dam foundation (elevation of 2181 m) is
cially defined at rates of 0.5, 1, 1.8, and 3 m/d (Fig. 3). At the initial considered downstream; and the left, right, and bottom of the dam
moment, the upstream water level is 19 m above the dam foundation foundation are set as zero-flow boundaries.
(elevation of 2198 m), and it is linearly decreased at a specific rate until The soil parameters for this study were partially obtained from
reaching the stabilized water level 10 m above the dam foundation existing studies (Wang et al., 2020) and are presented in Table 1 and
(elevation of 2189 m) on days 18, 9, 5, and 3. Transient seepage analyses Table 2. γ is the saturated gravity of the soil. The soil parameters of earth
are performed with preset hydraulic boundary conditions: a linear dam body follow a lognormal distribution. The θs of the soil body can be

6
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

obtained directly from tests, but θr is usually difficult to estimate and is

generally considered zero. The SWCC and PCF of the unsaturated soils of
the dam body are depicted in Fig. 4(a)-(b), where the dam foundation is
saturated for an extended period. An analysis model of the Ashigong
earth dam (Fig. 5) was constructed, and most of the RF cells of the dam
body were rectangular cells with a length of 2 ~ 3 m and a width of 1 m.
The finite element analysis grid near the seepage control wall was
properly encrypted, and the whole model was divided into 1935 cells.
The seepage calculation results are shown in Fig. 6(a)-(d), high­
lighting the effectiveness of the clay wall and concrete cut-off wall in
reducing the PWP within the dam body and significantly decreasing the
position of the infiltration line. The stability calculations were carried
out separately for the upstream and downstream slopes of the Ashigong
dam, and the results are shown in Fig. 7(a)-(b). When the upstream
water level of the reservoir was decreased, the stability of the upstream
dam slope initially decreased and then increased over time, with the
Fig. 3. Variation in the upstream water level under different calcula­ highest rate of water level decrease resulting in the earliest point of
tion conditions. maximum instability. This is because there existed a large hydraulic
gradient inside the dam slope during the phase of water level descent,
and the hydrostatic pressure acting on the upstream slope decreased.
Table 1 The water pressure and permeability were oriented outward, generating
Soil parameters of earth dam body. an outward pulling force that resulted in a decrease in the FS. Later,
Soil type Parameter Mean COV Distribution Correlation during the phase of water level stabilization, the infiltration line within
distance the dam continued to drop due to the hysteresis effect, leading to a
Sandy gravel ks (m/s) 5 × 10- 0.5 Lognormal lx = 50 m redistribution of stress. Consequently, the FS for the upstream dam slope
5
filler ly = 10 m exhibited a slight increase. In comparison, the stability of the down­
φ’ (◦ ) 29 0.2
stream embankment increased over time and then tended to plateau. For
c’ (kPa) 0 −
γ (kN/m3) 22 − example, at a drawdown rate of 1.8 m/d, the FS of the upstream
a (kPa) 5 0.4 embankment decreased from 1.239 to 1.180 and then increased to
n 2.68 0.2 1.197, whereas the FS of the downstream embankment increased from
θs (m3/ 0.39 − 1.220 to 1.372.
m3)
θr (m3/ 0 −
For the research case in this paper, the deterministic analysis of slope
m3) reveals that the FS for the upstream dam slope is lower than that of the
downstream slope, indicating that the downstream slope exhibits
greater stability (as shown in Fig. 7(a)-(b)). However, when accounting
Table 2
for the SV of soil parameters in reliability assessments, we find that the
Soil parameters of other parts of earth dam. Pf for the upstream and downstream slopes are 10-4 and 10-2 respec­
tively, indicating that the downstream slope is more susceptible to
Soil type ks (m/s) φ’ (◦ ) c’ (kPa) γ (kN/m3)
instability. Previous studies (Duncan, 2000) have also shown that results
-5
Sandy gravel layer 3.2 × 10 27 0 19.9 for stability analyses obtained using two metrics (the FS and the Pf) may
Silty clay layer 5 × 10-7 11 16 18
be different, it is related to the SV of soil parameters. Consequently, this
Silty fine-sand interlayer /layer 5 × 10-6 21 0 18.8
Gravel layer 3.2 × 10-5 30 0 21.2 study concentrates on the time-varying reliability analysis of the
Argillaceous siltstone layer 5 × 10-7 21 0 18.8 downstream slope. When the upstream water level of the reservoir
Slope protection 5 × 10-7 30 0 21.2 decreased, the PWP of the Ashigong dam decreased, leading to an
Loam layer 3.2 × 10-5 11 16 18 improvement in the stability of the downstream dam slope. Specifically,
Filter layer 5 × 10-7 27 0 19.9
Concrete cut-off wall 5 × 10-7 45 1 × 103 2.36 × 103
the reduction in the upstream water level reduced the hydrostatic load
on the dam body and the dynamic water pressure within the earth-rock
embankment, improving the stability of the downstream embankment.

Fig. 4. SWCC and PCF: (a) SWCC; and (b) PCF.

7
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 5. Finite element model.

6.2. Simulation of the multiparameter SV of the soil body results in a trade-off with model accuracy. In this study, when the
dimension reduction target is set to 15, the RMSE of the SIR-GPR model
In this work, the parameters φ’, ks, a, and n of the dam body are is approximately 0.02, and the R2 is approximately 0.94 across various
regarded as RF parameters to assess their SV. Determining the co­ computational working conditions. This demonstrates a high degree of
efficients of variation (COVs), probability distributions, autocorrelation modelling accuracy. A further reduction in dimensionality leads to a
functions and distances required for RF simulations is not a simple task more rapid increase in model error, suggesting that the optimal
because of the absence of data from actual surveys and experimental dimensionality reduction target is 15.
SWCCs. In this research, we refer to the data available in the literature To demonstrate the accuracy of the SIR-GPR model constructed for
(Wang et al., 2020), and the values are presented in Table 1. Since Wang predicting the stability of earth dams, the LEM analysis method and the
et al. (2020) did not address the SV of hydraulic parameters a and n, the SIR-GPR surrogate model prediction method are employed to calculate
COVs for a and n in this study are set according to the literature the FS for 2000 sets of training samples at the initial moment. The results
(Guardiani et al., 2021). Moreover, a negative correlation between a and are compared in Fig. 10. The FSs computed by the two methods are fitted
n is existed in the natural state, assumed a correlation coefficient of ρa,n to the 45◦ line, indicating a good prediction performance of the SIR-GPR
= -0.25 in this study. model. The RMSE for 2000 groups of samples is 0.029, and the R2 is
Given that the parameters φ’, ks, a, and n are all positive values, they 0.928. In addition, the constructed SIR-GPR model is used to predict the
are modelled as lognormal-type RFs. The KL method is employed for the FSs of 100 sets of test samples, as illustrated in Fig. 11. The extremely
discretization of the RFs, ensuring that the discretization error of each high percentage of predicted values within the 95 % confidence interval
RF meets accuracy requirements when the discretization term is set to M = [Lower 95 %, Upper 95 %] further validates the accuracy of the
10. A total of 40 random variables are obtained after discretizing the RF model predictions.
simulations for the earth dam slope parameters φ’, ks, a, and n. The
model is divided into 1935 RF elements on the basis of the actual size of
the earth dam slope. A typical realization of the RF for φ’, ks, a, and n is 6.4. Time-varying reliability analysis
illustrated in Fig. 8(a)-(d), displaying good smoothness and continuity.
Notably, the negative correlation between the RFs of a and n is obvious In the above subsection, the SIR-GPR model constructed from 2000
from Fig. 8(c)-(d). training samples is found to be effective for predicting the FS of the
A series of RF samples for soil parameters provides inputs for earth Ashigong embankment slope, with the initial moment used as an
dam stability modelling. The FS corresponding to each FEM sample is example. In this section, the SIR-GPR model is further implemented to
evaluated by repeating the transient seepage and stability analyses. On predict the time-varying FS of the dam under varying water level drop
the basis of the generated RF samples and corresponding FSs, SIR-GPR conditions. First, 1 × 104 RF simulations are conducted under four
method is applied to accurately obtain the connection among the soil different water level drop conditions. Each FEM is configured with a
parameters and the FS of the Ashigong dam. one-day time step to capture the changes in dam stability over time, and
the daily FS of the RF samples is predicted on the basis of the established
SIR-GPR model. The mean FS of the 1 × 104 samples predicted with the
6.3. Construction of the SIR-GPR model SIR-GPR model is displayed in Fig. 12, and the FS of the dam is nega­
tively correlated with the upstream water level. Specifically, as the
Sufficient training samples are vital for constructing the SIR-GPR water level decreases, the FS of the dam synchronously and steadily
model. Nevertheless, for earth dam stability problems, preparing these increases from the initial 1.25 to 1.346, and the rate of FS rise is faster
training samples in advance requires repeating numerous deterministic with a more rapid water level decline; when the water level decreases to
analyses, which is a highly time-consuming process. For this research, the lowest level, the FS stabilizes at approximately 1.346. The proba­
the number of training samples is set to 50 times the number of random bility density function (PDF) and cumulative distribution function (CDF)
variables (N = 2000) for various water level drop conditions based on of the FS are plotted in Fig. 13(a)-(b) for the case of a 1.8 m/d water level
our empirical observations. The RF samples and the corresponding FS drawdown rate. The low water level corresponds to PDF with a wide
values are selected through Latin hypercube sampling (LHS) to form the distribution range, and PDF kurtosis decreases as skewness increases. As
training sample set used to construct the SIR-GPR model. the water level decreases, both the PDF and CDF horizontally shift to­
The discretization of RF for a single finite element sample results in wards the region of high FS values; when the level decreases to the
the generation of 40 random variables, resulting in a relatively high stabilization level, the PDF and CDF remain basically unchanged.
dimension for the problem. The computational complexity in high- On this basis, 1 × 106 MCS are conducted to evaluate the Pf, and the
dimensional space is generally high, and the number of samples results are shown in Fig. 14(a). Pf displays the opposite trend as observed
required can increase exponentially with the increasing dimension, for FS. As the water level decreases, Pf generally tends to decrease,
necessitating the use of the SIR method for dimensionality reduction. To gradually decreasing from 0.02 at the initial moment to approximately
determine the optimal dimensionality reduction with the SIR-GPR 0.01. Moreover, the faster the water level decreases, the faster Pf de­
model, with the initial moment as an example, the model errors in creases; once the water level decreases to the lowest level, Pf stabilizes at
different dimensions (s = 1, 2, …, 40) are calculated, as depicted in approximately 0.01. This suggests that for the practical case of
Fig. 9(a)-(b). Reducing the dimensionality of the variables inevitably embankment dams discussed herein, properly lowering the upstream

8
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 6. PWP distribution and sliding surface at different times: (a) initial moment; (b) day 3; (c) day 15; and (d) day 30.

water level would be beneficial to the stability of the dam slope. results are shown in Fig. 14(b). Fig. 14 (a)-(b) show that the time-
To confirm the accuracy of the time-varying Pf obtained via the SIR- varying patterns of the Ashigong earth dam slope Pf calculated via two
GPR-MCS method, the reliability of the Ashigong earth dam slope is reliability analysis methods are consistent. This finding indicates that
further assessed via the direct LHS method with 1 × 104 samples, and the the proposed SIR-GPR-MCS method is highly effective for performing

9
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 7. Stability of upstream and downstream slopes: (a) upstream; and (b) downstream.

Fig. 8. Realization of RF for different soil parameters: (a) the effective internal friction angle φ’; (b) the saturated permeability coefficient ks; (c) the hydraulic
parameter a; and (d) the hydraulic parameter n.

time-varying reliability analyses of earth dam slopes considering time-varying FS prediction for the Ashigong earth dam are presented in
multiparameter SV. Fig. 15(a)-(b). The RMSE of the model ranges from 0.025 to 0.04, and
The errors of the SIR-GPR surrogate model for the whole process of the R2 values ranges from 0.92 to 0.94. Both of these error metrics

10
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 9. SIR-GPR model errors corresponding to different downscaled dimensions: (a) RMSE; and (b) R2.

Fig. 10. Comparison between the real and predicted values of 2000 sets of
training samples. Fig. 12. The average FS changes over time.

With respect to computational efficiency, the most time-consuming

part of the whole process of reliability analysis is the numerical
computation of various RF finite element samples. To address the reli­
ability problem of an earth dam with Pf on the order of 10-2, similar to
the case in this study, the direct LHS method requires at least 1 × 104
FEM samples for seepage and stability analyses. In contrast, the SIR-
GPR-MCS reliability analysis method utilized in this study requires
only 2 × 103 FEM samples. This not only effectively reduces the
execution time but also reasonably mitigates computational costs. Thus,
for determining the time-varying reliability of other complex earth dam
slope systems and addressing small Pf problems, the computational ef­
ficiency of the presented approach is highly valuable.
In this paper, the time-varying reliability of the Ashigong earth dam
considering multiparameter SV is initially investigated via the SIR-GPR-
MCS reliability analysis method. This method can be readily extended to
analyses of other earth dam slopes, contingent upon obtaining specific
information regarding hydrological and geological conditions,
Fig. 11. FS of 100 sets of test samples.
geotechnical parameters, etc. However, it is crucial to note that the SIR-
GPR model calibration dataset will vary, as each earth dam slope re­
quires site-specific data. Therefore, geotechnical investigators and op­
indicate that the SIR-GPR model’s accuracy is satisfactory. Notably, as
erators intending to employ the SIR-GPR-MCS approach in practical
the Pf of the Ashigong dam decreases with decreasing water level, the
scenarios will need to reprogram the model on the basis of new datasets
model accuracy based on the machine learning method decreases when
generated for the target earth dam slopes.
dealing with the small failure probability problem. The RMSE of the
In the case of the Ashigong earth dam, the proposed SIR-GPR-MCS
surrogate model increases slightly, from 0.025 at the initial moment to
reliability analysis method can accurately depict the time-varying
almost 0.037; the R2 values of the surrogate model decreases slightly
pattern of Pf. Nonetheless, reservoir levels often undergo periodic fluc­
from 0.91 at the initial moment to almost 0.89. This phenomenon may
tuations, and long-term dynamic actions may result in the deterioration
be ameliorated by appropriately extending the number of training
of geotechnical properties. This deterioration significantly impacts the
samples to further optimize the SIR-GPR model.
reliability of the earth dam slope but is not accounted for in this study.

11
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 13. The PDF and CDF of FS changes over time: (a) the PDF of FS; and (b) the CDF of FS.

Fig. 14. Time-varying Pf trend of the Ashigong earth dam under different water level drop conditions: (a) Pf calculated via the SIR-GPR-MCS method; and (b) Pf
calculated via the LHS method.

Fig. 15. Error analysis of the SIR-GPR model: (a) RMSE over time; and (b) R2 over time.

Therefore, future studies need to consider the deterioration properties of 6.5. Effects of the statistical characterization of soil parameters on Pf
soil in time-varying reliability analyses involving practical applications.
The reliability of earth embankments is closely related to the SV of
soil parameters and associated correlations. The effects of the statistical

12
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

characteristics of the soil parameters on the Pf of the Ashigong earth dam from 0.791 ~ 1.4717. This finding suggests that large lx and ly values
can be quantitatively analysed from three aspects: (1) the impact of the correspond to a PDF with a wide range of intervals, and the CDF shifts
COV of the soil parameters on Pf; (2) the influence of lx and ly of the soil towards the region of small FS values. Consequently, Pf increases.
on Pf; and (3) the consequences of neglecting the effect of the multipa­ Moreover, the effect of ly on Pf for the Ashigong earth dam is greater than
rameter SV of the soil on Pf. the effect of lx for the same increment.
First, the impact of the COV of the soil parameters on Pf is investi­ Neglecting spatial variation of multiple soil parameters in reliability
gated. With the variability of φ’ as an example, the Pf of the earth analysis impacts the Pf assessment of the Ashigong earth dam. lx and ly
embankment is calculated with a coefficient of variation, COVφ’, varying can be set as (lx = ly = 1000 m) if the SV of soil parameters of the earth
from 0.05 to 0.3 under a rate of water level drop of 1.0 m/d. The results dam slope is not considered. As depicted in Fig. 18(a), for the case in this
are depicted in Fig. 16(a), indicating that a large COVφ’ corresponds to a paper under the working condition of a 1.0 m/d water level drop rate,
large Pf. To elucidate this phenomenon, the distributions of the FS cor­ the Pf results obtained by considering and ignoring the SV of the soil
responding to 2000 sets of RF realizations with different COVφ’ values parameters vary by a factor of 2 to 3. Ignoring the SV of soil parameters
are obtained, and the PDF and the CDF of the FS are fitted accordingly, leads to the overestimation of the Pf of the earth dam slope, leading to a
as shown in Fig. 16(b)-(c), taking day 30 as an example. The diverse conservative design result, which may result in an uneconomical engi­
PDFs and CDFs exhibit significant variations. For an increase in COVφ’ neering design solution. The PDF and CDF obtained for the FS are also
from 0.05 to 0.3, the FS ranges from 1.345 ~ 1.432, 1.149 ~ 1.468, inconsistent, as displayed in Fig. 18(b)-(c), taking day 30 as an example.
0.978 ~ 1.471, 0.791 ~ 1.472, 0.682 ~ 1.4721, and 0.567 ~ 1.472, Therefore, the SV of multiple soil parameters needs to be considered in
with kurtosis values of 9.419, 8.088, 7.109, 5.350, 3.918, and 3.214, Pf assessments of earth dam slopes.
respectively. A large COVφ’ corresponds to a PDF with low kurtosis, and
the CDF shifts towards the region of smaller FS values, leading to an 7. Conclusions
increase in Pf.
The impact of the soil correlation distance on Pf is subsequently In this study, a SIR-GPR-MCS method is developed for the time-
investigated. Three control groups of lx and ly were established under the varying reliability analysis of earth dam slope stability considering
working condition of a 1.0 m/d water level drop rate: (lx = 50 m, ly = 5 multiparameter SV. An important output of this method is the dynamic
m), (lx = 100 m, ly = 5 m), and (lx = 50 m, ly = 10 m). The results are probability of dam slope failure over time. The proposed method is
depicted in Fig. 17 (a), where large lx and ly values correspond to a large illustrated through a real example of the Ashigong earth dam subjected
Pf. This phenomenon is further elucidated by the PDF and CDF of FS, to fluctuations in upstream reservoir levels, providing insights into the
taking day 30 as an example, as illustrated in Fig. 17 (b)-(c). With (lx = practical engineering design. The key findings are summarized below.
50 m, ly = 5 m) as the base case, doubling lx expands the distribution of (1) The proposed SIR-GPR-MCS reliability analysis method can
the FS from the initial 0.907 ~ 1.4716 to 0.841 ~ 1.4717. Similarly, effectively depict the trend of the time-varying Pf of the Ashigong earth
doubling ly results in an even larger distribution of FS values ranging dam. Combining the SIR-GPR-MCS method with time sequences for

Fig. 16. Influence of the variation coefficient of φ’ on Pf: (a) effects on Pf; (b) effects on the PDF of FS; and (c) effects on the CDF of FS.

13
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Fig. 17. Influence of the correlation distance on Pf: (a) effects on Pf; (b) effects on the PDF of FS and (c) effects on the CDF of FS.

Fig. 18. Influence of the SV of soil parameters on Pf: (a) effects on Pf; (b) effects on the PDF of FS; and (c) effects on the CDF of FS.

14
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

reliability assessment not only alleviates the mathematical burden for References
time-varying reliability analysis of embankment slopes considering
multiparameter SV but also provides information to support geotech­ Ahmad, M., Amjad, M., Al-Mansob, R.A., Kamiński, P., Olczak, P., Khan, B.J., Alguno, A.
C., 2022. Prediction of liquefaction-induced lateral displacements using Gaussian
nical engineers and managers to further characterize the stability con­ process regression. Appl. Sci. 12 (4), 1977.
ditions of earth dam slopes. This, in turn, can aid in the design of Cao, Z.J., Wang, Y., Li, D.Q., 2016. Site-specific characterization of soil properties using
reasonable measures for the prevention of dam instability disasters. multiple measurements from different test procedures at different locations – a
Bayesian sequential updating approach. Eng. Geol. 211, 150–161.
(2) The statistical characteristics of multiple soil parameters collec­ Chen, H.H., Zhang, L.Y., 2022. A machine learning-based method for predicting end-
tively influence the statistical distribution of FS and the associated Pf in bearing capacity of rock-socketed shafts. Rock Mech. Rock Eng. 55 (3), 1743–1757.
earth dam slope stability assessment. Therefore, the SV of diverse soil Cho, S.E., 2010. Probabilistic assessment of slope stability that considers the spatial
variability of soil properties. J. Geotech. Geoenviron. Eng. 136, 975–984.
parameters, including hydraulic parameters and shear strength param­ Cui, Q., Zhang, L.L., Chen, X.Y., Cao, Z.J., Wei, X., Zhang, J., Xu, J.B., Liu, D.S., Du, C.L.,
eters, should be considered when evaluating the reliability of reservoir 2022. Quantitative risk assessment of landslides with direct simulation of pre-failure
slopes and earth dam slopes. to post-failure behaviors. Acta Geotech. 17, 4497–4514.
Deng, Z.P., Pan, M., Niu, J.T., Jiang, S.H., 2022. Full probability design of soil slopes
(3) In engineering practices related to the slopes of earth and rock
considering both stratigraphic uncertainty and spatial variability of soil properties.
dams, disregarding the SV of various soil parameters can result in an Bull. Eng. Geol. Environ. 81 (5), 195.
overestimation of the slope’s Pf. This overestimation, consequently, Deng, Z.P., Pan, M., Niu, J.T., Jiang, S.H., Qian, W.W., 2021. Slope reliability analysis in
leads to the development of design schemes that are not only uneco­ spatially variable soils using sliced inverse regression-based multivariate adaptive
regression spline. Bull. Eng. Geol. Environ. 80, 7213–7226.
nomical but also inefficient. Deng, Z.P., Pan, M., Niu, J.T., Jiang, S.H., Wu, B.B., Li, S.L., 2023. Spatial prediction of
(4) The SIR-GPR-MCS method provides a new tool for solving other rockhead profile using the Gaussian process regression method. Can. Geotech. J. 60,
complex earth dam slope system reliability problems, multiparameter 1849–1860.
Duncan, J.M., 2000. Factors of safety and reliability in geotechnical engineering.
high-dimensional reliability problems, and time-varying reliability J. Geotech. Geoenviron. 126 (4), 307–316.
problems at low failure probability levels. Duan, N., Li, K.C., 1991. Slicing regression: a link-free regression method. Ann. Stat. 19,
Notably, in the case of the Ashigong earth dam, the surrogate model 505–530.
Fakhri, D., Mahmoodzadeh, A., Mohammed, A.H., Khodayari, A., Ibrahim, H.H.,
predicts only the reliability at a specific time step, which is then inte­ Rashidi, S., Karim, S.H.T., 2023. Forecasting failure load of sandstone under
grated to obtain the time-varying pattern of dam reliability over a different freezing-thawing cycles using Gaussian process regression method and grey
selected period. However, the temporal information embedded in the wolf optimization algorithm. Theor. Appl. Fract. Mech. 125, 103876.
Fathipour-Azar, H., 2021. Machine learning-assisted distinct element model calibration:
time sequence data are not captured or used for learning. Therefore, ANFIS, SVM, GPR, and MARS approaches. Acta Geotech. 17 (4), 1207–1217.
future work should focus on developing a reliability growth model to Fredlund, D.G., Morgenstern, N.R., Widger, R.A., 1978. The shear strength of
further predict the reliability of earth dam slopes on future time scales in unsaturated soils. Can. Geotech. J. 15 (3), 313–321.
GEO-SLOPE International Ltd, 2021. Heet and Mass transfer Modeling with GeoStudio.
a data-driven manner. Furthermore, the hysteresis effect present in the
GEO-SLOPE International Ltd., Calgary, Canada.
soil–water characteristic curve significantly influences the shear Griffiths, D.V., Fenton, G.A., 2004. Probabilistic slope stability analysis by finite
strength characteristics and mechanical behavior of unsaturated soils. elements. J. Geotech. Geoenviron. Eng. 130, 507–518.
Consequently, in further research, it is essential to conduct a series of Gu, X., Wang, L., Ou, Q., Zhang, W.G., 2023. Efficient stochastic analysis of unsaturated
slopes subjected to various rainfall intensities and patterns. Geosci. Front. 14 (1),
specific tests to gain a more precise understanding of the soil conditions 101490.
affecting the stability of slopes in earth and rock dams. Guardiani, C., Soranzo, E., Wu, W., 2021. Time-dependent reliability analysis of
unsaturated slopes under rapid drawdown with intelligent surrogate models. Acta
Geotech. 17, 1071–1096.
CRediT authorship contribution statement Guo, X.F., Dias, D., Carvajal, C., Peyras, L., Breul, P., 2019. A comparative study of
different reliability methods for high dimensional stochastic problems related to
Zhi-Ping Deng: Writing – original draft, Software, Methodology, earth dam stability analyses. Eng. Struct. 188, 591–602.
He, P., Wang, G., Xu, F., Sun, S.Q., 2021. GPR-MCS model of reliability analysis of key
Formal analysis. Yi Zou: Writing – review & editing, Investigation. blocks and its engineering application. Int. J. Numer. Anal. Meth. Geomech. 45 (12),
Wenqi Du: Writing – review & editing, Supervision, Funding acquisi­ 1739–1755.
tion. Dian-Qing Li: Supervision, Conceptualization. Mao-Xin Wang: Hu, C., Lei, R.D., 2023. A probabilistic framework for the stability analysis of slopes
considering the coupling influence of random ground motion and conditional
Visualization, Validation, Data curation. Shui-Hua Jiang: Visualization, random field. Soil Dyn. Earthq. Eng. 164, 107632.
Investigation. Ji, J., Wang, C.W., Gao, Y.F., Zhang, L.M., 2021a. Probabilistic investigation of the
seismic displacement of earth slopes under stochastic ground motion: a rotational
sliding block analysis. Can. Geotech. J. 58 (7), 952–968.
Declaration of competing interest
Ji, J., Wang, C.W., Cui, H.Z., Li, X.Y., Song, J., Gao, Y.F., 2021b. A simplified nonlinear
coupled Newmark displacement model with degrading yield acceleration for seismic
The authors declare that they have no known competing financial slope stability analysis. Int. J. Numer. Anal. Meth. Geomech. 45 (10), 1303–1322.
interests or personal relationships that could have appeared to influence Ji, J., Zhang, T., Cui, H.Z., Yin, X., Zhang, W.J., 2023. Numerical investigation of post-
earthquake rainfall-induced slope instability considering strain-softening effect of
the work reported in this paper. soils. Soil Dyn. Earthq. Eng. 171, 107938.
Jiang, S.H., Li, D.Q., Zhang, L.M., Zhou, C.B., 2014. Slope reliability analysis considering
Acknowledgments spatially variable shear strength parameters using a non-intrusive stochastic finite
element method. Eng. Geol. 168, 120–128.
Jiang, S.H., Huang, J.S., Qi, X.H., Zhou, C.B., 2020a. Efficient probabilistic back analysis
This study has been supported by the Visiting Researcher Fund of spatially varying soil parameters for slope reliability assessment. Eng. Geol. 271,
Program of State Key Laboratory of Water Resources Engineering and 105597.
Jiang, S.H., Liu, X., Huang, J.S., 2020b. Non-intrusive reliability analysis of unsaturated
Management (No. 2022SGG03), the National Natural Science Founda­ embankment slopes accounting for spatial variabilities of soil hydraulic and shear
tion of China (No. 52378344), the Young Elite Scientists Sponsorship strength parameters. Eng. Comput. 38, 1–14.
Program by JXAST (No. 2023QT08), and the Natural Science Founda­ Jiang, S.H., Liu, X., Huang, J.S., Zhou, C.B., 2022. Efficient reliability-based design of
slope angles in spatially variable soils with field data. Int. J. Numer. Anal. Meth.
tion of Jiangxi Province (Project No. 20242BAB23045). The financial Geomech. 46 (13), 2461–2490.
support is gratefully acknowledged. Li, D.Q., Jiang, S.H., Cao, Z.J., Zhou, W., Zhou, C.B., Zhang, L.M., 2015. A multiple
response-surface method for slope reliability analysis considering spatial variability
of soil properties. Eng. Geol. 187, 60–72.
Data availability
Li, M.Y., Wang, G., Qian, L., Li, X.P., Ma, Z.Y., 2021. ASS-GPR: Adaptive sequential
sampling method based on gaussian process regression for reliability analysis of
Data generated or analysed during this study are available from the complex geotechnical engineering. Int. J. Geomech. 21 (10), 04021192.
corresponding author upon reasonable request. Li, Y.L., Tang, W., W, L.F., Wang, J., 2022. Study on seismic failure probability of high
earth-rock dam considering dam body deformation and slope stability. Eur. J.
Environ. Civ. Eng. 26 (9), 3981–3995.

15
Z.-P. Deng et al. Computers and Geotechnics 188 (2025) 107558

Liao, K., Wu, Y.P., Miao, F.S., Li, L.W., Xue, Y., 2020. Time-varying reliability analysis of Shuku, T., Phoon, K.K., 2023. Comparison of data-driven site characterization methods
Majiagou landslide based on weakening of hydro-fluctuation belt under wetting- through benchmarking: methodological and application aspects. ASCE-ASME J. Risk
drying cycles. Landslides 18, 267–280. Uncertain. Eng. Syst. a. Civ. Eng. 9 (2), 04023006.
Liao, W.W., Ji, J., 2021. Time-dependent reliability analysis of rainfall-induced shallow Siacara, A.T., Beck, A.T., Futai, M.M., 2020. Reliability analysis of rapid drawdown of an
landslides considering spatial variability of soil permeability. Comput. Geotech. 129, earth dam using direct coupling. Comput. Geotech. 118, 103336.
103903. Song, L.F., Xu, B., Kong, X.J., Zou, D.G., Pang, R., Yu, X., Zhang, Z.Y., 2019. Three-
Liu, L.L., Cheng, Y.M., 2016. Efficient system reliability analysis of soil slopes using dimensional slope dynamic stability reliability assessment based on the probability
multivariate adaptive regression splines-based Monte Carlo simulation. Comput. density evolution method. Soil Dyn. Earthq. Eng. 120, 360–368.
Geotech. 79, 41–54. Su, Z.Y., Zhang, K., Liu, C.D., 2022. Dynamic risk assessment of slope stability of
Mahmoodzadeh, A., Mohammadi, M., Ali, H.F.H., Abdulhamid, S.N., Ibrahim, H.H., homogeneous earth-rock dam under action of multiple hazards. SIMULATION 98
Noori, K.M.G., 2021. Dynamic prediction models of rock quality designation in (8), 699–710.
tunneling projects. Transp. Geotech. 27, 100497. Suthar, M., 2019. Applying several machine learning approaches for prediction of
Morgenstern, N.R., Price, V.E., 1965. The analysis of the stability of general slip surface. unconfined compressive strength of stabilized pond ashes. Neural Comput. & Applic.
Géotechnique 15 (1), 79–93. 32, 9019–9028.
Nguyen, T.S., Likitlersuang, S., Tanapalungkorn, W., Phan, T.N., Keawsawasvong, S., Tan, X.H., Wang, X., Khoshnevisan, S., Hou, X.L., Zha, F.Z., 2017. Seepage analysis of
2022. Influence of copula approaches on reliability analysis of slope stability using earth dams considering spatial variability of hydraulic parameters. Eng. Geol. 228,
random adaptive finite element limit analysis. Int. J. Numer. Anal. Meth. Geomech. 260–269.
46 (12), 2211–2232. Vanmarcke, E., Shinozuka, M., Nakagiri, S., Schuëller, G.I., Grigoriu, M., 1986. Random
Phoon, K.K., Kulhawy, F.H., 1999. Characterization of geotechnical variability. Can. fields and stochastic finite elements. Struct. Saf. 3 (3–4), 143–166.
Geotech. J. 36, 612–624. Wang, L., Cao, Z.J., Li, D.Q., Phoon, K.K., Au, S.K., 2018. Determination of site-specific
Qi, X.H., Wang, H., Pan, X.H., Chu, J., Chiam, K., 2021. Prediction of interfaces of soil-water characteristic curve from a limited number of test data – a Bayesian
geological formations using the multivariate adaptive regression spline method. perspective. Geosci. Front. 9 (6), 1665–1677.
Undergr. Space. 6 (3), 252–266. Wang, L., Wu, C.Z., Tang, L.B., Zhang, W.G., Lacasse, S., Liu, H.L., Gao, L., 2020. Efficient
Raja, M.N.A., Jaffar, S.T.A., Bardhan, A., Shukla, S.K., 2023. Predicting and validating reliability analysis of earth dam slope stability using extreme gradient boosting
the load-settlement behavior of large-scale geosynthetic-reinforced soil abutments method. Acta Geotech. 15, 3135–3150.
using hybrid intelligent modeling. J. Rock Mech. Geotech. Eng. 15 (3), 773–788. Zhang, W.G., Wu, J.H., Gu, X., Han, L., Wang, L., 2022. Probabilistic stability analysis of
Rasmussen, C.E., Williams, C.K.I., 2006. Gaussian processes for machine learning. MIT embankment slopes considering the spatial variability of soil properties and seismic
press, Cambridge, MA. randomness. J. Mt. Sci. 19 (5), 1464–1474.
Ray, R., Choudhary, S.S., Roy, L.B., 2021. Reliability analysis of soil slope stability using Zhu, B., Hiraishi, T., Pei, H.F., Yang, Q., 2020. Efficient reliability analysis of slopes
MARS, GPR and FN soft computing techniques. Model. Earth Syst. Environ. 8, integrating the random field method and a Gaussian process regression-based
2347–2357. surrogate model. Int. J. Numer. Anal. Meth. Geomech. 45 (4), 478–501.
Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. J. Appl,
Phys, p. 1.

16