# Reliability analysis of backward erosion piping in an embankment dam considering the spatial variability of soil properties

Fengdong Chi, Claudio Carvajal, Pierre Breul, Laurent Peyras

► To cite this version:

Fengdong Chi, Claudio Carvajal, Pierre Breul, Laurent Peyras. Reliability analysis of backward erosion piping in an embankment dam considering the spatial variability of soil properties. Scientific Reports, 2025, 15 (1), pp.23974. 10.1038/s41598-025-09175-9. hal-05460371

HAL Id: hal-05460371

https://hal.inrae.fr/hal-05460371v1

Submitted on 15 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

OPEN

# Reliability analysis of backward erosion piping in an embankment dam considering the spatial variability of soil properties

Fengdong Chi<fnref n="1" />,²✉, Claudio Carvajal<fnref n="1" />, Pierre Breul² & Laurent Peyras<fnref n="1" />

The spatial variability of soil seepage parameters is an essential component of the uncertainty of earth dams. This paper presents a method for backward erosion reliability analysis that considers the spatial variability of seepage parameters in earth dams to assess the distribution of the Factor of Safety (FoS) and the initiation probability ($P_f$). In this study, the Karhunen-Loeve method is used to generate Random Fields (RFs) of soil hydraulic properties, and Monte Carlo simulations (MCS) are used to generate realisations of these RFs, which are then introduced into the numerical model to assess the spatial variability of the seepage analysis results. The hydraulic gradient distribution is compared with the critical hydraulic gradient to assess the FoS and the $P_f$ due to backward erosion at each point of the dam. This study also presents the effect of water level $H$, permeability anisotropy coefficient $ξ$, and RFs parameters (coefficient of variation $COV_{Ks}$, correlation length $L_h$ & $L_v$) on the FoS and $P_f$ of Backward Erosion Piping (BEP). The results illustrate that the initiation probability $P_f$ increases with $H$, $COV_{Ks}$ and $ξ$. The initiation probability $P_f$ is not sensitive to the vertical correlation length, but may decrease with the horizontal correlation length. Furthermore, the study discusses the potential initiation locations of internal erosion and presents results based on criteria using both local and global hydraulic gradients.

**Keywords** Random fields, Initiation probability, Random finite element method, Reliability analysis, Backward erosion piping, Embankment dam

As one of the most common types of earth dam failure, internal erosion has been a major concern for experts and scholars worldwide. Statistics on the failures of large earth dams show that almost 50% of the failures are due to internal erosion, the rest are mainly attributed to overtopping<fnref n="1" /> and a small fraction to the sliding mechanism². Traditional seepage analysis adopts the deterministic analysis method and takes hydraulic gradient or flow velocity as the evaluation criterion. However, there are significant uncertainties associated with the geotechnical parameters used in seepage analysis and in assessing the risk of internal erosion. It is interesting to use probabilistic methods to quantify these uncertainties and to assess the reliability of seepage analysis results³,⁴. Probabilistic methods provide a systematic approach to incorporating uncertainty into the analysis by considering the probability distribution of input parameters⁵⁻⁸. The variability and uncertainty associated with these parameters can be effectively captured by assigning probability distributions to geotechnical parameters such as permeability⁹, porosity¹⁰, and soil strength¹¹. This approach allows for a more comprehensive understanding of the potential range of outcomes and the associated risks.

Backward erosion piping is one form of internal erosion that can cause the failure of an earth dam. Several scholars have studied backward erosion using analytic and deterministic analyses: Some studies¹²⁻¹⁵ have proposed models that are briefly introduced in ICOLD Bulletin 164¹⁶. Terzaghi's model¹² showed that backward erosion piping occurs when the effective stress at the downstream toe of the dam reaches zero. In this model, the critical hydraulic gradient can be estimated based on the soil's specific gravity and porosity. Sellmeijer's model¹³ was proposed based on visual observations of experiments of horizontal flow. In this model, the factor of safety (FoS) is influenced by the resistant factor $F_R$, scale factor $F_S$, and geometrical shape factor $F_G$. Schmertmann¹⁴ proposed an analytical/empirical model based on Sellmeijer's research work. This model uses an extensive set of experimental results from flume tests to estimate the safety factor against piping at any location along an anticipated progression path. Hoffmans¹⁵ described an analytical/empirical model based on Darcy's law. The equation for predicting the critical hydraulic gradient comprises the critical hydraulic gradients of

<footnote n="1">¹Aix-Marseille University, INRAE, RECOVER, Aix-en-Provence 13182, France. ²Université Clermont Auvergne, CNRS, Clermont Auvergne INP, Institut Pascal, Clermont-Ferrand 63000, France. ✉email: fengdong.chi@inrae.fr</footnote>
Shields (critical flow resistance in the pipes) and the critical hydraulic gradients of Darcy (seepage resistance in the aquifer). Fell et al.<sup>17</sup> proposed a method to approximately estimate the time for the progression of internal erosion piping and the development of a breach leading to failure in embankment dams and their foundations. Benaissa et al.<sup>18</sup> modeled water-clay mixtures flowing through a cylindrical pipe domain using the computational fluid dynamics (CFD) simulation method, and the results were validated by conducting a hole erosion test. The study showed that the fine particles are separated from the skeleton under the seepage flow, which reduces soil strength and increases porosity. El Shamy et al.<sup>19</sup> investigated the interaction between eroded soil particles, water flow, and soil skeleton. They used Navier-Stokes equations to solve problems involving continuous fluids and employed particle flow software (PFC) to simulate fine particles with internal piping erosion. Guo & Yu<sup>20</sup> coupled computational fluid dynamics and discrete element method (DEM) computational models, using CFD-DEM to simulate internal soil erosion and predict scouring depth and soil erosion rate. The study demonstrated the effectiveness of the particle flow software in calculating the degree of damage caused by internal erosion to the dam. Peng & Rice<sup>21</sup> employed an inverse analysis to interpret pore pressure data and observations during backward erosion piping (BEP) initiation and progression in sandy soils, identifying critical hydraulic conditions for stages such as channel initiation and sand fluidization.

However, the soil properties of dams are associated with uncertainty, leading to uncertainty in the factor of safety against backward erosion piping. Some scholars have conducted reliability analysis to estimate the initiation probability of the dam. Most of these studies use the random variables (RVs) method, while several studies use the random fields (RFs) method that is capable of considering the spatial correlation characteristics of the geotechnical properties in the study area. Cho<sup>22</sup> performed a reliability analysis of seepage through an embankment and its foundation to study the effects of uncertainty due to the spatial heterogeneity of permeability on seepage flow. Liu et al.<sup>23</sup> investigated the influence of five commonly used autocorrelation functions (ACFs) on the seepage flow problem based on Cho's model. The study found that the squared exponential (SQX) and single exponential (SNX) ACFs may overestimate and underestimate the seepage flow rate, respectively, while the differences between different ACFs are not significant. Tan et al.<sup>24</sup> investigated the influence of the spatial variability of hydraulic parameters on the flow in earth dams by considering soil-water characteristic curve parameters as lognormal random fields. Liang et al.<sup>25</sup> developed a numerical algorithm and conducted simulations to study backward erosion progressions in geologic media with spatially stochastic parameters, including void ratio, permeability, and the ratio of the particle contents. The study assessed the probability of failure due to seepage instability at the site investigated. Hekmatzadeh et al.<sup>26</sup> studied the performance of four cutoff wall configurations and evaluated them stochastically using the random finite element method. The study found that the probabilities of failure for different cutoff wall configurations are similar when considering isotropic permeability, but there are noticeable differences in anisotropic situations. Moreover, Chi et al.<sup>27</sup> compared the effects of three types of permeability RFs (stationary, conditional, and non-stationary RFs) on hydraulic parameters such as flow velocity, hydraulic gradient, pore water pressure, and flow rate in stochastic seepage analysis. The results showed that the different types of permeability RFs had a significant impact on flow velocity, a moderate impact on hydraulic gradient, and a minor impact on pore water pressure.

Although research on the reliability analysis of backward erosion is limited, the main objective of this study is to develop a probabilistic approach to evaluate the safety of earth dams by focusing on the backward erosion mechanism. The specific research objectives of this study are: (1) to build a hydraulic model for the probabilistic evaluation of a dam taking the internal erosion mechanism into account; (2) to analyse the spatial variability of internal flows and their impact on assessing the backward erosion mechanism in an earth dam; (3) to includes a parameter sensitivity analysis for the upstream water level and the spatial variability of permeability. The results of backward erosion are discussed and compared using the local hydraulic gradient method (Terzaghi method) and the global critical gradient method (Sellmeijer method & Hoffmans method). The results of the position of initiation of backward erosion are discussed and compared, based on the FoS and initiation probability $P_f$ assessed at each point of the dam.

The originality of this study does not lie in the development of new erosion models or novel mechanistic descriptions of backward erosion. Rather, it stems from the integration of existing backward erosion models (e.g., Terzaghi, Sellmeijer, and Hoffmans criteria) with a probabilistic framework that accounts for the spatial variability of soil properties through random fields. This approach enables a spatially-resolved assessment of the factor of safety (FoS) distribution (mean and standard deviation) and the initiation probability of internal erosion at each point of the dam, thus providing valuable insights into dam safety under uncertainty.

BEP is a progressive, time-dependent process involving detachment, transport, and enlargement phases. It should also be clarified that the scope of this study is limited to the initiation phase of backward erosion, which represents a necessary condition for failure. The initiation criterion—based on the comparison between the local hydraulic gradient and the critical gradient—is commonly used to identify high-risk zones in geotechnical assessments¹⁶.

The article is structured as follows to achieve these goals. Section 2 briefly presents the methods used for the seepage analysis of earth dams and the backward erosion mechanism. Section 3 introduces the case of an earth dam, the simulation of random fields and the procedure of a finite element method combined with Monte Carlo simulations (MCS). Section 4 presents the results of the probabilistic seepage analysis, which provides the mean and standard deviation of FoS and initiation probability $P_f$ for backward erosion piping initiation. Furthermore, global and parametric sensitivity analyses are conducted to investigate the impact of the uncertainty parameters of spatial variability on backward erosion. Section 5 discusses the initiation position of backward erosion and the initiation probability $P_f$ obtained by using a local or global hydraulic gradient. Finally, Sect. 6 presents the main conclusions.
Materials and methods for seepage analysis and backward erosion initiation analysis
Seepage analysis

Internal erosion often causes damage to hydraulic structures and foundations. According to a survey conducted by Foster et al.²⁸ on more than 11,000 large earth dams, 45% of failure modes are caused by internal erosion.

All dams have some seepage as the impounded water seeks paths of least resistance through the dam and its foundation. The seepage analysis model used in this study is based on governing equations solved by the finite element method. Richards' equation²⁹ is used to describe the flow through porous media. The generalised Darcy's law (see Eq. (1)) and steady-state are assumed in this study. The governing equation for steady-state seepage flow is written as Eq. (2).

$$v = \frac{Q}{A} = k \cdot \frac{h_1 - h_2}{L} = ki \qquad (1)$$

$$\frac{\partial}{\partial x} \left( k_x \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial z} \left( k_z \frac{\partial h}{\partial z} \right) = 0 \qquad (2)$$

Where v=flow velocity, Q=flow rate, A=cross-sectional area, k=permeability (hydraulic conductivity), h₁-h₂=elevation difference between the water levels, L length of the soil sample, i=hydraulic gradient, kₓ and kᵧ are the permeability of the soil along directions x and z, respectively, and h is the water level. The Optum-G2 software was used for the finite element analysis of seepage in this study.

This study considers a saturated-unsaturated seepage analysis in earth dams. Unsaturated seepage refers to the flow of water through the soil where the water content is less than the maximum capacity of the material. Compared with the saturated seepage calculation, the unsaturated seepage calculation requires considering the soil water retention curve (SWRC). The Van Genuchten model³⁰ is considered in this study. The volumetric water content of the soil and permeability can be expressed as Eqs. (3–1), (3–2):

$$S_e(\psi) = \frac{\theta - \theta_r}{\theta_s - \theta_r} = \frac{1}{[1 + (\psi\alpha)^{n_r}]^m} \qquad (3-1)$$

$$K = K_s S_e^{1/2} \left[ 1 - \left( 1 - S_e^{1/m} \right)^m \right]^2 \qquad (3-2)$$

where: S_e (ψ) is a function of the volumetric moisture content of the soil, ψ is the matric suction, S_e is for the effective saturation, θ is the volumetric water content, θ_s is the saturated volumetric water content, θ_r is the residual volumetric water content, arameter α represents the suction corresponding to the inflexion point when the SWRC changes from the saturated state to the unsaturated state, parameter m represents the soil parameters corresponding to the residual moisture content of the soil, m=1-1/n. K is the hydraulic conductivity function, and K_s is the saturated hydraulic conductivity.

Backward erosion initiation

As one of the main mechanisms of internal erosion, backward erosion piping commences with the detachment of particles at the exit of leakage paths through or under the embankment. The erosion process begins at a free surface on the downstream side of the dam. Backward erosion piping may occur when the hydraulic gradient (i) exceeds the critical hydraulic gradient (i_c). Piping refers to the movement and loss of fine-grained soil in the pores formed by the coarse-grained skeleton. The increase of pores in the soil eventually leads to the formation of seepage channels inside the soil, resulting in soil damage²⁸. Internal erosion is a progressive, time-dependent process involving detachment, transport, and enlargement phases. This study focuses on backward erosion initiation, a common simplification that enables efficient spatial analysis under hydraulic uncertainty. ICOLD Bulletin 164¹⁶ recommends several methods to estimate the local and global critical hydraulic gradient for assessing backward erosion initiation.

The calculation formula for the local critical hydraulic gradient, proposed by Terzaghi³¹ and based on the force balance of unit soil in the seepage field, is given by Eq. (4) where n_p is the porosity, and G_s is the specific gravity of the soil. Physically, i_c reflects the local resistance of soil to the detachment of particles under seepage. Zones with higher porosity or lower particle density exhibit lower i_c, indicating weaker structural integrity and higher susceptibility to backward erosion. The factor of safety (FoS) against backward erosion piping is defined as the ratio of the critical hydraulic gradient i_c to the hydraulic gradient i, as shown in Eq. (5).

$$i_{cr} = (G_s - 1)(1 - n_p) \qquad (4)$$

$$\text{FoS} = \frac{i_{cr}}{i} = \begin{cases} \geq 1, \text{safe} \\ < 1, \text{failure} \end{cases} \qquad (5)$$

Regarding the global critical hydraulic gradient, there are several methods used in practical engineering to evaluate dam safety, such as the Sellmeijer method³², and the Hoffmans method³³.
In the Sellmeijer model, the critical hydraulic gradient is determined as a product of three contributions: resistance factor, scale factor and geometrical shape factor (as shown in Eq. (6)):

$$ 
\left\{ 
\begin{aligned}
\frac{H_c}{L} &= \frac{1}{c} = F_r F_s F_g \\
F_r &= \eta (G_s - 1) \tan \vartheta \\
F_s &= \frac{d_{70}}{\sqrt{\kappa L_{0.28}}} \\
F_G &= 0.91 \left(\frac{D}{L}\right)^{\left(\frac{D}{L}\right)^{2.8}-1} + 0.04 \\
\kappa &= \frac{\nu}{g} k
\end{aligned}
\right. 
\qquad (6) 
$$

where $H_c$ (m) is the critical pressure head; $L$ (m) is the seepage length (base length of the embankment); $H_c/L$ is the critical hydraulic gradient; $D$ (m) is the thickness of the sand layer under the embankment; $c$ is the erosion coefficient; $F_r$ is the resistance factor; $F_s$ is the scale factor; $F_G$ is the geometrical shape factor; $d_{70}$ (m) is the soil particle diameter for which 70% by weight of the soil is finer; $\eta$ is the Whites drag coefficient; $\vartheta$ is the bedding angle of the sand; $\kappa$ (m²) is the intrinsic permeability; $\nu$ (m²/s) is the kinematic viscosity, and $g$ (m²/s) is the gravity.

In Hoffmans model, the most important parameters are permeability ($k$), particle size $d_{50}$, $d_{15}$ and the aquifer dimensions: the depth of the sand layer ($D$) and seepage length ($L$). The equations of Hoffmans backward erosion are expressed as Eq. (7):

$$ 
\left\{ 
\begin{aligned}
\frac{(H_1 - H_2)_c}{L} &= \frac{\sqrt{g(\psi_{lam,c}(G_s - 1)d_{15})^{3/2}}}{\nu\sqrt{\alpha_R}} + \left(1 - \frac{l_c}{L}\right) \frac{d_{50}\nu}{l_{Re}kD} \\
\frac{l_c}{L} &= \exp\left(-\left(\frac{\alpha_f D}{L}\right)^2 \frac{\sqrt{g(\psi_{lam,c}(G_s - 1)d_{15})^{3/2}}}{\nu\sqrt{\alpha_R}}\right) \\
\psi_{lam,c} &= 0.2(D^*)^{-0.3} \\
D^* &= d_{50}((G_s - 1)g/\nu^2)^{1/3}
\end{aligned}
\right. 
\qquad (7) 
$$

where $H_1$ is the upstream water level; $H_2$ is the water level at the soil surface (usually taken as 0 m); $(H_1 - H_2)_c/L$ is the critical hydraulic gradient; $l_c$ is the ratio of critical pipe length and seepage length; $l_{Re}$ is the geometrical length scale; $\alpha_f$ is the geometrical groundwater coefficient; $\alpha_R$ is the geometrical pipe coefficient; $\psi_{lam,c}$ is the critical Shields' parameter.

For both the Sellmeijer and Hoffmans models in the global approach, the global $FoS_{global}$ is equal to the ratio between the global critical hydraulic gradient $i_{cr_global}$ and the global hydraulic gradient $i_{global}$, see Eq. (8):

$$ FoS_{global} = \frac{i_{cr\_global}}{i_{global}} = \frac{H_c/L}{H/L} = \frac{H_c}{H} \qquad (8) $$

## Methodology for reliability analysis of backward erosion
### Case study presentation

In this study, the numerical modelling for the seepage analysis was performed within a hypothetical earth dam²² and the geometry of the dam is shown in Fig. 1. The embankment dam is 5 m high, with upstream and downstream slopes of 2 h:1 $\nu$, and is placed on a 5 m thick soil foundation. This simplified geometry, commonly used in benchmark studies²³,³⁴, facilitates the analysis of spatial variability in soil hydraulic properties and allows comparison with previous studies. The predicted flow rate and hydraulic gradient are consistent with published results, supporting the validity of the modeling approach. The research focused on the seepage analysis and backward erosion of the dam under steady-state saturated-unsaturated flow conditions.

The studied domain was discretised to 8009 triangular elements with 4154 nodes. A mesh convergence analysis was conducted to verify that the output deterministic analysis results—particularly the factor of safety

<figure><img src="image_3.png" /><figcaption>Fig. 1. The finite element model in Optum G2.</figcaption></figure>
<table><caption>**Table 1.** Mesh convergence analysis for flow rate and factor of safety (FoS).</caption><thead><tr><th>Case</th><th>Elements</th><th>Nodes</th><th>Mesh size (m)</th><th>Flow rate Q (m³/day/m)</th><th>FoS</th><th>|ε(Q)|</th><th>|ε(FoS)|</th></tr></thead><tbody><tr><td>1</td><td>2060</td><td>1106</td><td>0.5</td><td>1.487</td><td>1.773</td><td>2.04%</td><td>10.70%</td></tr><tr><td>2</td><td>4007</td><td>2110</td><td>0.38</td><td>1.518</td><td>1.601</td><td>0.98%</td><td>5.00%</td></tr><tr><td>3</td><td>6084</td><td>3171</td><td>0.33</td><td>1.533</td><td>1.525</td><td>0.65%</td><td>4.28%</td></tr><tr><td>4</td><td>8009</td><td>4154</td><td>0.25</td><td>1.543</td><td>1.462</td><td>0.52%</td><td>3.51%</td></tr><tr><td>5</td><td>10,123</td><td>5228</td><td>0.23</td><td>1.551</td><td>1.413</td><td>0.32%</td><td>3.11%</td></tr><tr><td>6</td><td>12,054</td><td>6210</td><td>0.22</td><td>1.556</td><td>1.370</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>**Table 2.** Statistical properties of soil parameters used for backward erosion analysis.</caption><thead><tr><th>Parameter</th><th>Mean</th><th>COV</th><th>Auto-correlation distance</th></tr></thead><tbody><tr><td>K<sub>s</sub></td><td>0.864 m/day (10<sup>-5</sup> m/s)</td><td>0.75</td><td rowspan="5">L<sub>h</sub>=20 m, L<sub>v</sub>=2 m</td></tr><tr><td>Critical hydraulic gradient i<sub>cr</sub></td><td>1.06</td><td>0.255</td></tr><tr><td>Specific gravity G<sub>s</sub></td><td>2.65</td><td>0.06</td></tr><tr><td>Porosity n<sub>p</sub></td><td>0.357</td><td>0.055</td></tr><tr><td>θ<sub>r</sub></td><td>0.045</td><td>-</td></tr><tr><td>θ<sub>s</sub></td><td>0.430</td><td>-</td><td>-</td></tr><tr><td>α (m<sup>-1</sup>)</td><td>1.478</td><td>-</td><td>-</td></tr><tr><td>n</td><td>2.68</td><td>-</td><td>-</td></tr><tr><td>Anisotropy coefficient ξ (K<sub>h</sub>/K<sub>v</sub>)</td><td>1</td><td>-</td><td>-</td></tr></tbody></table>

(FoS) and flow rate (Q)—stabilized with increasing mesh refinement. The results are presented in Table 1. |ε(Q)| represents the relative difference in flow rate between a given mesh and the next finer mesh, while |ε(FoS)| denotes the relative difference in the factor of safety. Both metrics are used to evaluate the sensitivity of the results to mesh refinement and to assess numerical convergence. As shown, beyond approximately 8000 elements, further refinement had a minimal effect on the estimated FoS and flow rate. In addition, the adopted mesh size of 0.25 m satisfies the criterion³⁵, which recommends that the maximum element length should be less than half of the autocorrelation distance in any direction. Therefore, the mesh with 8009 triangular elements and 4154 nodes was selected for all the subsequent probabilistic simulations, balancing accuracy and computational cost.

The hydraulic parameters (K<sub>s</sub>, θ<sub>r</sub>, θ<sub>s</sub>, α, n) adopted for the embankment and foundation can be found in Cho's study²², and are in turn derived from the soil database³⁶. The soil parameters (specific gravity G<sub>s</sub> and porosity n<sub>p</sub>) that contribute to the resistance of backward erosion piping based on Terzaghi criteria (critical gradient i<sub>cr</sub>) can be found in Mainguenaud's study³⁷ (see Table 2). These parameters are representative of moderately well-graded sandy soils. They are consistent with laboratory-tested materials from a real French dam site³⁸, and their classification aligns with the Unified Soil Classification System³⁹.

Moreover, the seepage analysis results—such as flow rate (see Table 1) and hydraulic gradient (see Fig. 9a)—are in close agreement with those reported in benchmark studies²²,²³.

By simulating random samples of G<sub>s</sub> and n<sub>p</sub> according to the probability distribution shown in Table 2, it is possible to use Eq. (4) to calculate the critical gradient i<sub>cr</sub> for each sample of G<sub>s</sub> and n<sub>p</sub>. The resulting sample of critical gradient values can be used to estimate the mean and standard deviation of the critical gradient. Figure 2a,b show the histograms and probability density function (PDF) curves for the simulated samples of G<sub>s</sub> and n<sub>p</sub> respectively. Figure 2c shows the histogram and the fitted PDF curve for samples of calculated critical gradient values.

Figure 3 presents the grain size distribution curves used to determine the particle sizes d<sub>20</sub>, d<sub>50</sub>, and d<sub>15</sub>. Figure 4 shows the soil-water retention curve used in this study, which is based on the Van Genuchten model fitting parameters used in Cho's study²². The prescribed boundary conditions for this model are as follows: the upstream water level H is assumed to be 9.5 m, and the downstream water table coincides with the surface of the dam foundation.

In this research, the Auto-Correlation Function (ACF) method, more specifically the single exponential autocorrelation function (SNX), was employed to generate random fields and represent the spatial variability of soil property parameters. The correlation length of each ACF type is associated with its Scale of Fluctuation. Regarding SNX, the Scale of Fluctuation used in this investigation is twice the correlation length. The values for L<sub>h</sub> (horizontal correlation length) and L<sub>v</sub> (vertical correlation length) were selected as 20 m and 2 m, respectively, based on assumptions made in earlier research²⁴,⁴⁰,⁴¹.

As mentioned in Sect. 2.2, besides the Terzaghi criterion based on the local hydraulic gradient, simplified methods based on the global hydraulic gradient (such as the Sellmeijer and Hoffmans methods) can also be used to assess the risk of the backward erosion of an earth dam. Compared to the Terzaghi method, which considers variable parameters such as permeability K<sub>s</sub>, specific gravity G<sub>s</sub>, and porosity n<sub>p</sub>, as shown in Table 2, the Sellmeijer and Hoffmans methods require considering additional variable parameters, as illustrated in Table 3.
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Distribution of soil parameters about backward erosion strength.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Grain size distribution for dam.</figcaption></figure>

The spatial variability of the bedding angle of sand and the soil particle diameters (e.g. *d*<sub>70</sub>, *d*<sub>50</sub>, and *d*<sub>15</sub>) are considered in this study, the suggested parameter ranges are based on some researchers' works<sup>32,33,42</sup>. Other parameters can be considered constants due to the inherent properties of the earth dam and its operating environment. For example, when the water temperature is at 20 degrees Celsius, the kinematic viscosity and drag coefficient are identified (ν = 1.0*10<sup>-6</sup> m/s, η = 0.25). The comparison and discussion of these three methods in estimating the risk for backward erosion will be presented in Sect. 5.2.

Simulation of random fields

A random field (RF) can describe the spatial correlation of material properties in different locations and represent their heterogeneity characteristics<sup>43</sup>.

The Karhunen-Loève expansion (K-L) method was selected in this study to discretise random variables and generate random fields in order to simulate the spatial variability of soil properties. This method has been widely
<figure><img src="image_6.png" /><figcaption>Fig. 4. SWRC curves of van Genuchten model.</figcaption></figure>

<table><caption>**Table 3.** The statistical properties of additional parameters used in the Sellmeijer method & Hoffmans' method.</caption><thead><tr><td>Parameters</td><td>Mean</td><td>COV</td><td>Comment</td></tr></thead><tbody><tr><td>tanθ</td><td>0.676</td><td>0.04</td><td>Normal distribution, used in Sellmeijer method</td></tr><tr><td>d<sub>70</sub> (10<sup>-3</sup> m)</td><td>0.6</td><td>1.5</td><td>Lognormal Distribution, used in Sellmeijer method</td></tr><tr><td>d<sub>50</sub> (10<sup>-3</sup> m)</td><td>0.4</td><td>1.5</td><td>Lognormal Distribution, used in Hoffmans' method</td></tr><tr><td>d<sub>15</sub> (10<sup>-3</sup> m)</td><td>0.2</td><td>1.5</td><td>Lognormal Distribution, used in Hoffmans' method</td></tr><tr><td>η</td><td>0.25</td><td>-</td><td>Used in Sellmeijer method</td></tr><tr><td>ν (m<sup>2</sup>/s)</td><td>1.0*10<sup>-6</sup></td><td>-</td><td>Used in Sellmeijer method</td></tr><tr><td>α<sub>R</sub></td><td>6</td><td>-</td><td>Used in Hoffmans' method</td></tr><tr><td>α<sub>f</sub></td><td>5</td><td>-</td><td>Used in Hoffmans' method</td></tr><tr><td>D (m)</td><td>5</td><td>-</td><td>Used in Sellmeijer &amp; Hoffmans' method</td></tr><tr><td>L (m)</td><td>22</td><td>-</td><td>Used in Sellmeijer &amp; Hoffmans' method</td></tr><tr><td>g (m/s<sup>2</sup>)</td><td>9.81</td><td>-</td><td>Used in Sellmeijer &amp; Hoffmans' method</td></tr></tbody></table>

used in geotechnical problems such as slopes<sup>44</sup>, tunnels<sup>45</sup>, and foundations<sup>46</sup>. The Gaussian random field $H_i(x,θ)$ can be expressed as Eq. (9):

$$H_i(x, \theta) = \mu + \sum_{i=1}^{\infty} \sigma \sqrt{\lambda_i \varphi_i(x_{RF})} \xi_i(\theta) \approx \mu + \sum_{i=1}^{N_{KL}} \sigma \sqrt{\lambda_i \varphi_i(x_{RF})} \xi_i(\theta) \quad (x \in \Omega) \quad (9)$$

where $\mu$ is the mean of the parameters in RFs, $\sigma$ is the standard deviation of the parameters in RFs, $\lambda_i$ and $\phi_i(x_{RF})$ are the eigenvalues and eigenfunctions of the autocovariance function for the RF, respectively.

The truncation number $N_{KL}$ is determined by evaluating the truncated error with a target threshold<sup>47</sup>. The first task is to determine the truncation term number $N_{KL}$ once the necessary probabilistic parameters (mean, COV<sub>Ks</sub>, L<sub>h</sub>, and L<sub>v</sub>) are defined. It can be achieved by evaluating the truncation error of a K-L expansion RF evaluated by Eq. (10) with a target accuracy. In this paper, the $N_{KL}$ is determined for a $\epsilon_{KL}$ lower than 2%. As shown in Eq. (10), $\epsilon_{KL}$ is the variance-based error globally estimated in the RF domain $\Omega$, $x_{RF}$ is the coordinate of an arbitrary point in the field.

$$\epsilon_{KL} = \frac{1}{\Omega} \int_{\Omega} \left[ 1 - \sum_{i=1}^{N_{KL}} \lambda_i \varphi_i^2(x_{RF}) \right] d\Omega \quad (10)$$

## Methodology based on the stochastic finite element method

The Monte Carlo simulation method, also known as the stochastic simulation or random experiment method, is a numerical technique that uses repeated random sampling to obtain numerical results<sup>48</sup>. Monte Carlo simulation has been widely used as a reliability method to evaluate structural safety in geotechnical engineering<sup>49-52</sup>. Among the many methods of reliability analysis, it is a relatively accurate method provided that a sufficiently large number of simulations are carried out. The aim of this paper is to study the influence of the spatial variability of soil parameters on the reliability or initiation probability, which has to be resolved through probabilistic seepage analysis. This study uses the Monte Carlo simulation method to simulate the distribution of the factor of safety (FoS) to initiate backward erosion.
This section presents the seepage analysis procedure using a random finite element method, as illustrated in Fig. 5. The procedure consists of the following main steps:

(1) Identify the parameters with spatial variability and determine their statistical characteristics, including means, coefficients of variation (COVs) and distributions related to the seepage analysis. Then, select an appropriate autocorrelation function and the autocorrelation distances for the 2D random-field model.

(2) Apply the Optum-G2 software to build the finite element analysis model for seepage analysis, using the mean values of input parameters. Then, obtain the coordinates $(x_p, y_i)$ of the nodes, where $i=1, 2, ..., N_p$, $N_n$ being the number of the nodes in the finite element model. The deterministic seepage analysis model files are saved as the template file named 'Template.g2x'. This file contains all the information used by OPTUM-G2, and it can also be accessed and viewed directly using a text editor.

(3) Create a standard normal sample matrix using Latin hypercube sampling points. Next, apply the K-L expansion from Sect. 3.2 to simulate non-Gaussian random fields of spatially varying variables using the sample matrix and coordinates $(x_p, y_i)$. Then, generate $N_p$ realisations of non-Gaussian random fields that represent spatially varying permeability values in the physical space. Finally, assign each realisation to corresponding nodes in the finite element model. This stratified sampling technique improves the efficiency and uniformity of sample distribution compared to standard Monte Carlo sampling, especially for moderate sample sizes.

(4) Substitute the mean values of the corresponding uncertain input parameters in the 'Template.g2x' file (created in step 2) with each pair of spatially varying variables (e.g. $K_s$ and $i_{cr}$ in Table 2) in each realisation of random fields from step (3). This process will generate $N_p$ distinct new 'Template.g2x' input files, named 1.g2x, 2.g2x, ..., $N_p$.g2x.

(5) Execute OPTUM-G2 using each of the new input files, namely 1.g2x, 2.g2x, ..., $N_p$.g2x, which were generated in step (5) to conduct seepage analysis. This procedure will yield $N_n * N_p$ distinct factors of safety for backward erosion at each node of the finite element model across $N_p$ realisations, resulting in a set of factors of safety denoted as Eq. (11). As shown in the matrix FoS: the elements in the 1st column present the FoS

<figure><img src="image_7.png" /><figcaption>Fig. 5. Flowchart of the reliability analysis.</figcaption></figure>
of all the nodes in the 1st seepage analysis realisation; the elements in the 1st row presen the FoS of the 1st node in all the seepage analysis realisations. For each realisation, the position of the dangerous (initiate backward erosion) node can be identified by looking for the minimum FoS in each column of matrix FoS. At each point of the dam (for each node $N_p$), the mean and standard deviation of the FoS can be calculated from the elements in each row of matrix FoS following the $N_p$ realisations from the MCS.

$$ FoS = \begin{cases} FoS_{1,1} & FoS_{1,2} \cdots FoS_{1,N_p} \\ FoS_{2,1} & FoS_{2,2} \cdots FoS_{2,N_p} \\ \cdots & \cdots \\ FoS_{N_n,1} & FoS_{N_n,2} \cdots FoS_{N_n,N_p} \end{cases} \quad (11) $$

(6) Evaluate the reliability index and the initiation probability due to internal erosion on the basis of the probabilistic distribution of the FoS obtained by the MCS. With MCS, the number of simulations must be sufficiently large to obtain the convergence of the variables of interest. The initiation probability can be defined directly as the ratio between the number of BEP initiations and the total number of simulations, which requires a very high number of simulations if the initiation probability is low. In order to limit the number of simulations, this study considered as variables of interest the mean and standard deviation of the FoS whose estimation converges more quickly. By adopting a probability distribution of the FoS with the mean and standard deviation values obtained from the MCS, the initiation probability $P_f$ can be evaluated directly from the FoS distribution as the probability that the FoS is less than 1.

Although the hydraulic models used are standard, the stochastic integration—random field generation, batch simulation, and statistical post-processing—was fully customized. The contribution of this study lies in applying established physical models within a novel spatially probabilistic framework.

## Reliability analysis and results

This section presents a probabilistic seepage analysis based on Monte Carlo simulations using the random finite element model. The Terzaghi criterion was used to calculate the critical hydraulic gradient for the seepage failure of the dam. The location of the dangerous node with the lowest factor of safety (FoS), which is more prone to initiate backward erosion, could be identified. Then, the contours of $P_f$ (initiation probability) were displayed to reflect the risk of backward erosion at each node directly, achieving the goal of data visualisation. Furthermore, a sensitivity analysis was performed to evaluate the significance of each stochastic input parameter on the FoS and the $P_f$. The sensitivity analysis of FoS and $P_f$ considered hydraulic parameters (Upstream water level $H$ and Anisotropy coefficient of permeability $\xi$) as well as random field parameters (COV$_{Ks}$ and correlation length $L_h$ & $L_v$).

### Results of the reliability analysis of backward erosion

#### FoS distribution and initiation probability $P_f$

The convergence analysis of the FoS of the node at the embankment toe was conducted through 1000 hydraulic simulations that incorporated the RFs (see Fig. 6). The estimated mean and standard deviation of the factor of safety (FoS) stabilize as the number of realizations exceeds 300, while the initiation probability $P_f$ converges with more than 500 realizations. Given this convergence behavior, and considering that the estimated failure probability $P_f$ is approximately 16.3% (see Figs. 6 and 11), 500 realizations ($N_p = 500$, as referenced in Sect. 3.3)

<figure><img src="image_8.png" /><figcaption>Fig. 6. Convergence of the estimated mean, standard deviation, coefficient of variation of the FoS and initiation probability.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Spatial variability of permeability for the case using $COV_{K_S}=0.75$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. PDF of FoS against backward erosion obtained from the lognormal distribution of the input parameters.</figcaption></figure>

were selected for each probabilistic seepage analysis. At this level, the standard deviation of the Monte Carlo estimate of $P_f$ is approximately 1.65%, resulting in a 95% confidence interval of ±3.3%, which is considered acceptable for engineering applications.

As done in previous studies, the permeability and critical hydraulic gradient are modeled by means of RFs to consider the spatial variability. The soil permeability of the earth dam has been considered a spatially random field that follows a lognormal distribution with a prescribed mean, variance, and spatial correlation structure⁵³. Similarly, the critical hydraulic gradient has been modeled as RFs that follow a lognormal distribution⁵⁴.

Figure 7 shows the spatial distribution of permeability based on the same $COV_{K_S} = 0.75$ but with different correlation lengths used to generate RFs. A comparison between Fig. 7a,b reveals that the RFs with a higher correlation length correspond to a more homogeneous permeability RF. By comparing Fig. 7a,c (or Fig. 7d), it can be observed that the permeability distribution should follow the main direction of the correlation length. For instance, when the horizontal correlation length is significantly larger than the vertical correlation length, the contour of the permeability field will present a horizontal strip-like distribution, and vice versa.

Figure 8 shows the probability density function of FoS against backward erosion for the configuration with $COV_{K_S} = 0.75$ and $H=9.5$ m. The initiation probability $P_f$ can be estimated from the probability density function fitted to the sample of FoS ($P_f = 16.3\%$, in this case). Figure 9 displays the contours of the estimated mean and the standard deviation of the hydraulic gradient and critical hydraulic gradient. The figures for FoS (Fig. 10) are generated by combining the contours of the hydraulic gradient and critical hydraulic gradient using Terzaghi criteria.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Spatial variability of the hydraulic gradient and critical hydraulic gradient from Monte-Carlo simulations.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Spatial variability of the factor of safety FoS from Monte-Carlo simulations.</figcaption></figure>

The contours of $P_f$ (Fig. 11) can be obtained by estimating the initiation probability from the FoS distribution for each node in the finite element model. The dangerous area in red with the highest backward erosion initiation probability $P_f$ = 16.3% (also with the lowest estimated mean FoS = 1.48) is located at the toe of the embankment, as illustrated in Fig. 11. $P_f$ can reach a high value due to the narrow geometry of the dam and the larger dispersion of input data ($\text{COV}_{\text{Ks}}$, $\text{COV}_{\text{Gs}}$, and $\text{COV}_{\text{np}}$) used in this study. Among these factors, the narrow geometry of the dam is mainly reflected in the narrow width of the dam crest, steep slope, and low dam thickness, which can result in a higher hydraulic gradient. It is noteworthy that, although this study identifies a dominant critical zone near the dam toe, other scenarios may involve the alignment or connection of multiple marginally unstable regions. Such interactions could facilitate the formation of continuous erosion paths, particularly in highly heterogeneous soils.

### Global sensitivity analysis

To assess the influence of key soil parameters ($K_s$, $G_s$, $n_p$) quantitatively on the factor of safety (FoS) against backward erosion initiation, a global sensitivity analysis was performed using the Morris method. This method is widely used in engineering applications because it can effectively screen input variables, capture main effects and potential nonlinear interactions, and has relatively low computational costs⁵⁵.

In this study, the Morris method was applied to three key input variables: saturated permeability ($K_s$), specific gravity ($G_s$), and porosity ($n_p$), all of which play a direct role in the calculation of hydraulic gradient and critical
<figure><img src="image_13.png" /><figcaption>Fig. 11. The contour of initiation probability of a dam in probabilistic backward erosion analysis.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Morris global sensitivity analysis results ($K_s$, $G_s$ & $n_p$).</figcaption></figure>

gradient ($i_{cr}$) based on the Terzaghi criterion (see Sect. 2.2). A total of 5,000 Monte Carlo simulations were used to evaluate the model response in terms of FoS.

Two indices were computed for each input parameter: the mean of the absolute elementary effects ($\mu^*$), which indicates the overall importance of the parameter, and the standard deviation ($\sigma$) of these effects, which reflects nonlinearity or interaction effects. The results are shown in Fig. 12a,b.

Figure 12a shows that $K_s$ is the most influential parameter ($\mu^* \approx 17$), followed by $n_p$ ($\mu^* \approx 9$) and $G_s$ ($\mu^* \approx 6$). Figure 12b illustrates the $\sigma-\mu^*$ scatter plot, where all three parameters exhibit very high $\sigma$ values ($> 190$), indicating significant nonlinearity or interaction effects. In particular, $K_s$ not only dominates in terms of direct influence but also contributes to strong interactions with other variables. The results suggest that variations in $K_s$ largely control the variability of FoS in the model, while $G_s$ & $n_p$ still participate in nonlinear behaviors under uncertain conditions, although with less impact.

These findings provide quantitative evidence supporting the prioritization of $K_s$ in both uncertainty characterization and parameter calibration efforts. Accordingly, the following parametric sensitivity analysis in Sect. 4.2 focuses on the spatial variability of permeability ($K_s$), including its coefficient of variation (COV$_{Ks}$) and correlation lengths ($L_h$ & $L_v$), which were shown to have the most significant influence on the probabilistic assessment of backward erosion initiation.

## Parametric sensitivity analysis

Based on the study in Sect. 4.1.1, the dangerous nodes with the lowest FoS and highest $P_f$ are generally located at the downstream toe of the embankment, where they readily initiate backward erosion. A sensitivity analysis was performed in this study to examine the impact of each stochastic input parameter on the FoS and the $P_f$ of backward erosion. The value of FoS and initiation probability are mainly influenced by hydraulic parameters (water level $H$ and anisotropy coefficient of permeability $\xi$) as well as random field parameters (COV$_{Ks}$ and correlation length $L_h$ & $L_v$). The parameter values used in Sect. 4.1.1 ($H=9.5$ m, $\xi=1$, COV$_{Ks}$ = 0.75, $L_h$ = 20 m & $L_v$ = 2 m) are considered as reference values. The sensitivity analysis considers different configurations of $H$, $\xi$, COV$_{Ks}$, $L_h$ & $L_v$. Each configuration changes only one value while keeping the others at their reference values. To discuss the influence of these parameters on the safety of the dam, this section is divided into two parts. The first part investigates the effect of water level $H$ and the coefficient of variation of permeability COV$_{Ks}$ (Sect. 4.2.1), while the second part examines the influence of anisotropy coefficient $\xi$ and correlation lengths $L_h$.
<table><caption>**Table 4.** The water level *H* and COV<sub>Ks</sub> influence on backward erosion.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Upstream water level H (m)</td><td>8, 8.5, 9, 9.5</td></tr><tr><td>COV<sub>Ks</sub></td><td>0.3, 0.45, 0.6, 0.75, 0.9</td></tr><tr><td>Anisotropy coefficient ξ</td><td>1</td></tr><tr><td>Correlation length L<sub>h</sub> &amp; L<sub>ν</sub> (m)</td><td>L<sub>h</sub> = 20, L<sub>ν</sub> = 2</td></tr></tbody></table>

<table><caption>**Table 5.** The anisotropy coefficient ξ and correlation lengths L<sub>h</sub> &amp; L<sub>ν</sub> influence on backward erosion.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Upstream water level H (m)</td><td>9.5</td></tr><tr><td>COV<sub>Ks</sub></td><td>0.75</td></tr><tr><td>Anisotropy coefficient ξ</td><td>1, 2, 3, 4</td></tr><tr><td>Correlation length (m)</td><td>5 ≤ L<sub>h</sub> ≤ 20, 2 ≤ L<sub>ν</sub> ≤ 8</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Influence of *H* and COV<sub>Ks</sub> on μ<sub>FoS</sub>.</figcaption></figure>

& L<sub>ν</sub> (Sect. 4.2.2). Thus, 20 configurations of water level and COV<sub>Ks</sub> are included in Sect. 4.2.1, and presented in Table 4. Altogether, there are 28 independent configurations of anisotropy coefficient ξ and correlation lengths L<sub>h</sub> &amp; L<sub>ν</sub> (20 configurations in Table 4 and 8 configurations in Table 5). For each configuration, 500 groups of random field realisations were conducted, and the 500 corresponding simulations of seepage analysis were performed. It should be noted that the convergence of the backward erosion analysis results considering 500 realisations was validated previously (see Fig. 6).

### Influence of water level *H* and coefficient of variation of permeability COV<sub>Ks</sub>

This section examines the impact of water level *H* and COV<sub>Ks</sub> on the FoS and P<sub>f</sub> of backward erosion assessed at the downstream toe of the dam. All the cases in this section assumed a correlation length of L<sub>h</sub> = 20 m and L<sub>ν</sub> = 2 m, and isotropic permeability (Anisotropy coefficient ξ = 1) for each node. Water level *H* and COV<sub>Ks</sub> vary from case to case, as shown in Table 4. Figures 13, 14, 15 illustrate the influence of *H* and COV<sub>Ks</sub> on the estimated mean FoS (μ<sub>FoS</sub>), standard deviation of FoS (σ<sub>FoS</sub>), and initiation probability P<sub>f</sub> respectively, which were obtained from the probabilistic seepage analysis. Figure 16 assists in explaining the results shown in Figs. 13, 14, 15 illustrating the influence of *H* and COV<sub>Ks</sub> on COV<sub>FoS</sub> and the influence of FoS and COV<sub>Ks</sub> on P<sub>f</sub> respectively.

Figure 13a displays a significant decrease in the estimated mean FoS (μ<sub>FoS</sub>) with an increase in reservoir water level. The rate of decrease slows down as the water level increases from 8 m to 9.5 m, as shown in the curve graph. Figure 13b shows a slight increase in μ<sub>FoS</sub> with the COV<sub>Ks</sub> at a relative water level of 8 m, but no significant increase for the water level range of 8.5 m to 9.5 m.

Figure 14a,b show that the estimated standard deviation of FoS (σ<sub>FoS</sub>) decreases overall as the water level increases, assuming constant COV<sub>Ks</sub>. However, a slight increase in σ<sub>FoS</sub> is observed when the water level rises from 9 m to 9.5 m. One possible explanation may be related to the position of the point considered in the evaluation of backward erosion (as discussed in Sect. 5.1).

Figure 15a illustrates that the rate of initiation probability P<sub>f</sub> increment accelerates as the water level *H* increases. Figure 15b shows that the initiation probability P<sub>f</sub> increases uniformly with COV<sub>Ks</sub> (with a slightly lower rate of increase for the case with a water level of 8.5 m). In contrast, Fig. 16 shows that the initiation
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Influence of H and $COV_{Ks}$ on $\sigma_{FoS}$.</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Influence of H and $COV_{Ks}$ on $P_f$.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 16.** Influence of FoS and $COV_{Ks}$ on $P_f$.</figcaption></figure>

probability $P_f$ decreases, and the rate of decrease slows down gradually with the increase in the factor of safety (FoS).

*Influence of anisotropy coefficient $\xi$ and correlation length $L_h$ & $L_v$*

For this analysis, the earth dam was assumed to be under the same load condition (Upstream water level is 9.5 m), and the $COV_{Ks}$ of the permeability RFs were kept constant at 0.75 for all subsequent cases. This section examines the anisotropy coefficient $\xi$ and correlation lengths $L_h$ & $L_v$ of the permeability RFs to determine their
<figure><img src="image_21.png" /><figcaption>Fig. 17. Influence of $\xi$, $L_h$ and $L_v$ on $\mu_{FoS}$</figcaption></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 18. Influence of $\xi$, $L_h$ and $L_v$ on $\sigma_{FoS}$</figcaption></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 19. Influence of $\xi$, $L_h$ and $L_v$ on $P_f$</figcaption></figure>

<figure><img src="image_26.png" /></figure>

impact on the estimated mean FoS ($\mu_{FoS}$), standard deviation of FoS ($\sigma_{FoS}$), and initiation probability $P_f$ which are obtained from the probabilistic seepage analysis (see Fig. 17, 18, 19). Table 5 shows the different values used for the anisotropy coefficient $\xi$ and correlation lengths $L_h$ & $L_v$ in this sensitivity analysis.

Figure 17a shows that the estimated mean of the FoS ($\mu_{FoS}$) decreases significantly with an increase in the anisotropy coefficient $\xi$. Additionally, it shows that the rate of decrease of $\mu_{FoS}$ slows down as the anisotropy
coefficient increases from 1 to 4. This is because increasing anisotropy leads to a higher horizontal permeability, which facilitates downstream flow and reduces FoS.

Figure 17b illustrates that $\mu_{FoS}$ may vary slightly with the correlation length, and the curves show the fluctuation of $\mu_{FoS}$. The first four points correspond to an increase in the horizontal correlation length $L_h$, and the last four points correspond to an increase in the vertical correlation length $L_v$. Figure 17a,b show that the anisotropy coefficient has a stronger impact than the correlation length on the mean of FoS.

Figure 18a displays the standard deviation of FoS ($\sigma_{FoS}$) with the variation in the anisotropy coefficient, and this is for different horizontal and vertical correlation length configurations. This figure shows that $\sigma_{FoS}$ is not highly sensitive to the anisotropy coefficient $\xi$. The most significant fluctuations of the standard deviation $\sigma_{FoS}$ as a function of the anisotropy coefficient are obtained for the configuration with a lower correlation length (e.g., the curve colored in black, $L_h = 5$ m, $L_v = 2$ m). Figure 18b illustrates that $\sigma_{FoS}$ decreases as the correlation length (and more particularly the horizontal correlation length $L_h$) increases. This is because a higher correlation length can reduce the dispersion of the permeability random field (as shown in Fig. 7), which in turn decreases the standard deviation of permeability, hydraulic gradient, and FoS step by step. Conversely, using a configuration with a shorter correlation length will result in a greater dispersion of permeability. This, in turn, will increase the dispersion of the FoS at the base of the earth dam.

Figure 19a shows that the initiation probability $P_f$ increases with the anisotropy coefficient, but the rate of increase can be slowed down with the anisotropy coefficient. As shown previously, the anisotropy coefficient has a minor effect on the standard deviation of the FoS (Fig. 18a), but a significant impact on the mean of the FoS (Fig. 17a). Therefore, an increase in the anisotropy coefficient results in a decrease in the mean FoS, which in turn increases the initiation probability due to backward erosion. Figure 19b shows that the initiation probability decreases overall with the correlation length. As observed previously, the correlation length has a minor effect on the mean of the FoS (Fig. 17b), but a significant impact on the standard deviation of the FoS (Fig. 18b). Therefore, an increase in the correlation length results in a decrease in the standard deviation of the FoS, which in turn reduces the initiation probability due to backward erosion.

## Discussion

This section presents a discussion on the of backward erosion initiation position and the assessment of the critical gradient using a local or global approach.

### Discussion on the backward erosion initiation location

Based on the results presented in Sect. 4, the potential initiation position for backward erosion (the node with the lowest FoS in the dam) varies across realisations in probabilistic seepage analyses based on Terzaghi criteria. In practical engineering, it is important to consider not only the initiation probability of backward erosion but also the initiation position. The lowest FoS value is only one of the necessary conditions for backward erosion. Another factor is the position of the node with the lowest FoS in the dam. Among all the realisations, the cases with the lowest FoS near the ground surface are more likely to initiate backward erosion. Therefore, this section highlights two configurations with different upstream water levels (e.g. $H=9$ m & 9.5 m) to comprehensively evaluate the safety of the earth dam. What is more, the potential initiation position of backward erosion can be identified by locating the node with the lowest FoS in the dam.

The results of the 500 numerical simulations conducted in the probabilistic seepage analysis show that the lowest FoS is typically observed in the vicinity of the downstream toe (e.g. points A and B in Fig. 20) or near the surface of the reservoir level (e.g., points C in Fig. 20). However, the initiation of erosion is more critical around the downstream toe (exit gradient). To give an example, the lowest FoS was evaluated at point A for 166 cases out of the 500 realisations (33.2%) when the upstream water level is 9 m (see Fig. 20a).

Figure 20a,b illustrate the position of the points with the lowest FoS for each of the 500 realisations of the numerical simulations (considering the upstream water level at 9 m and 9.5 m in Fig. 20a,b, respectively). By observing Fig. 20a,b, it can be found that the potential initiation position (the point with the lowest FoS) is located more frequently at the downstream toe of the embankment, but the number of realisations in the configuration with a 9 m water level (166 realisations, 33.2%) is higher than that in the configuration with a 9.5 m water level (150 realisations, 25.0%). For the configuration with a water level of 9.5 m, there may be a lower dispersion of the FoS at the embankment toe, resulting in a lower value of $\sigma_{FoS}$. This could explain the slight increase in $\sigma_{FoS}$ when the water level increases from 9 m to 9.5 m (Fig. 14a).

Figure 20 shows that the point with the lowest FoS value in some cases can be located close to the surface of the reservoir level (e.g. point C in Fig. 20): in 18.2% of the cases when the water level is at 9.0 m and in 19.2% of the cases when the water level is at 9.5 m. However, the values of FoS at this node (point C) are higher than 1 in most of the realisations, which makes it difficult to initiate backward erosion from this point. Overall, when using the local approach for backward erosion analysis, it is necessary to calculate the FoS and $P_f$ for each node of the dam and generate the corresponding distribution map as done in Sect. 4.1.1 (Fig. 11) in order to conduct a comprehensive safety assessment that includes the value of $P_f$ and the initiation position.

### Comparison of the local approach and global approach

In addition to a local approach based on the Terzaghi model to evaluate the local critical hydraulic gradient at each point of the dam, practical engineering also considers simplified empirical models based on the evaluation of the global critical hydraulic gradient (such as the Sellmeijer and Hoffmans models). The aim of this section is to compare the results obtained using a local approach and a global approach to evaluate the safety of the dam against backward erosion (factor of safety FoS & initiation probability $P_f$).
<figure><img src="image_27.png" /><figcaption>(a) The configuration with a 9 m water level load</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) The configuration with a 9.5 m water level load.</figcaption></figure>

<table><caption>**Fig. 20.** The distribution of the number of realisations (or probability) for the node with the lowest FoS (500 realisations in total).</caption><thead><tr><td>Result parameters</td><td>Hoffmans method</td><td>Sellmeijer method</td><td>Terzaghi method</td></tr></thead><tbody><tr><td>μ<sub>FoS</sub></td><td>13.664</td><td>4.091</td><td>1.485</td></tr><tr><td>σ<sub>FoS</sub></td><td>41.712</td><td>5.666</td><td>0.535</td></tr><tr><td>COV<sub>FoS</sub></td><td>3.0527</td><td>1.385</td><td>0.360</td></tr><tr><td>P<sub>f</sub></td><td>21.4%</td><td>22.2%</td><td>16.8%</td></tr></tbody></table>

<caption>**Table 6.** The comparison of reliability analysis results using different backward erosion criteria.</caption>

The criteria for backward erosion and the statistical properties of the soil parameters used in the Terzaghi, Sellmeijer and Hoffmans methods are given in Sect. 2.2, Tables 2 and 3. These elements are then used to compute the FoS and estimate the initiation probability $P_f$ using Monte Carlo simulations (see Sect. 3.2).

The results of the FoS distribution (mean and standard deviation) and the initiation probability obtained for the Terzaghi method are presented in Sect. 4.1.1.

For the Sellmeijer and Hoffmans methods, the Monte Carlo simulations consisted of 500 random samples of the input data according to the distribution presented in Table 3 in order to evaluate the global critical hydraulic gradient of each sample and thus obtain a distribution of the global FoS. The results of the probabilistic backward erosion analysis ($μ_{FoS}$, $σ_{FoS}$, & $P_f$) using the Terzaghi, Sellmeijer and Hoffmans methods are summarised in Table 6.

The results presented in Table 6 show that the global approach (Sellmeijer method and Hoffmans method) may overestimate the initiation probability of backward erosion when compared to the local approach (Terzaghi method). The initiation probability for the Sellmeijer and Hoffmans methods is 21.4% and 22.2%, respectively, while for the Terzaghi method it is 16.8%. The Sellmeijer and Hoffmans methods give a high mean FoS, but with a high dispersion (standard deviation and COV), resulting in a higher initiation probability. As presented in
Sect. 2.2, the Sellmeijer and Hoffmans models are based on a global critical hydraulic gradient approach, which takes into account several parameters associated with large uncertainties. This leads to a significant dispersion (e.g., COV) of the global critical hydraulic gradient and the global FoS, resulting in a higher estimated initiation probability in the Hoffmans and Sellmeijer models.

Although this study is based on a simplified dam geometry, it is important to note that backward erosion piping is a scale-sensitive process. Parameters such as hydraulic conductivity, particle size, and layer thickness may behave differently at full scale. Moreover, internal erosion is often a mechanism triggered locally, making it difficult to assess across the entire dam. These limitations are acknowledged, and future work will apply the method to real dam geometries and site-specific data to better evaluate scale effects.

## Conclusion

This study investigated a method for probabilistic backward erosion analysis that takes into account the spatial variability of seepage parameters in earth dams. The probabilistic seepage analysis is based on the random finite element method, which integrates random fields of soil hydraulic properties to obtain a random field of the hydraulic gradient. The method proposed allows obtaining, at any point of the dam, the distribution of the FoS and the initiation probability $P_f$ of backward erosion. The stochastic seepage analysis of the dam against backward erosion included a sensitivity analysis to assess the influence of uncertainty on hydraulic parameters (e.g., reservoir water level, permeability anisotropy coefficient) as well as random field parameters (e.g., coefficient of variation of permeability and correlation lengths). In addition, two different criteria based on local hydraulic gradient (Terzaghi method) and global hydraulic gradient (Sellmeijer and Hoffmans methods) were used to calculate the initiation probability of backward erosion to compare the results between these approaches. Several conclusions can be drawn from this study:

(1) The random finite element method proposed in this study can be used to assess the risk of backward erosion initiation in an earth dam, taking into account the spatial variability of soil properties. The method proposed makes it possible to obtain the FoS distribution (mean and standard deviation) and the initiation probability $P_f$ of internal erosion initiation at each point of the dam. The resulting FoS and $P_f$ distribution maps provide a clear and intuitive visualization of the most vulnerable zones, enabling the direct identification of potential backward erosion initiation areas within the dam body or foundation.

(2) Parametric sensitivity analysis showed that the water level $H$ and permeability anisotropy $\xi$ have a strong influence on the mean factor of safety (FoS), while the coefficient of variation of permeability $COV_{Ks}$ and the correlation lengths primarily affect its variability. Of all the parameters, $H$, $\xi$, and $COV_{Ks}$ significantly impacted the backward erosion initiation probability $P_f$ which increased with increasing values of these parameters, whereas the effect of correlation lengths was comparatively minor. This suggests that in practical applications, while careful characterization of anisotropy, permeability variability, and water level is essential, the selection of correlation lengths for random field generation can be approached with less stringency.

(3) The initiation probability $P_f$ of backward erosion was examined in this study using the local hydraulic gradient approach (Terzaghi method) and compared with results from global methods (Hoffmans and Sellmeijer). In the case study, the Sellmeijer and Hoffmans methods yielded higher mean FoS values, but also showed greater dispersion (standard deviation and COV), resulting in a higher $P_f$. This suggests that, in practical applications, global methods may overestimate the probability of initiation due to compounded uncertainties, and local criteria may offer more stable and spatially resolved risk assessments when supported by detailed input data.

BEP is a progressive, time-dependent process involving detachment, transport, and enlargement phases. This study only focused on the initiation of backward erosion without considering the other phases of the erosion process relating to the progression of erosion until dam failure. Future work may incorporate erosion kinetics, transient flow, and coupled fluid-soil interaction models (e.g., CFD-DEM) to simulate the full backward erosion process.

Despite these limitations, the current framework provides practical insights by identifying zones where the probability of internal erosion initiation is higher. This information can aid engineers not only in prioritizing monitoring efforts, but also in evaluating the effectiveness of design measures—such as drainage systems, slope geometry, or filter configurations—to reduce the probability of backward erosion initiation in real dam projects.

## Data availability

All associated data have been presented in the manuscript which is available from the corresponding author on reasonable request.

Received: 27 July 2024; Accepted: 23 June 2025

Published online: 04 July 2025

## References

1. Yang, M. et al. Uncertainty analysis on flood routing of embankment dam breach due to overtopping failure. *Sci. Rep.* **13**, 20151. https://doi.org/10.1038/s41598-023-47542-6 (2023).

2. Lizárraga, J. J. & Buscarnera, G. Probabilistic modeling of shallow landslide initiation using regional scale random fields. *Landslides* **17** (8), 1979–1988 (2020).

3. Siacara, A. T., Napa-García, G. F., Beck, A. T. & Futai, M. M. Reliability analysis of Earth dams using direct coupling. *J. Rock Mech. Geotech. Eng.* **12** (2), 366–380 (2020).

4. Johari, A. & Heydari, A. Reliability analysis of seepage using an applicable procedure based on stochastic scaled boundary finite element method. *Eng. Anal. Boundary Elem.* **94**, 44–59 (2018).
5. Ji, J., Zhang, C., Gui, Y., Lü, Q. & Kodikara, J. New observations on the application of LS-SVM in slope system reliability analysis. *J. Comput. Civil Eng.* **31** (2), 06016002 (2017).

6. Johari, A., Heydari, A. & Talebi, A. Prediction of discharge flow rate beneath sheet piles using scaled boundary finite element modeling database. *Scientia Iranica*. **28** (2), 645–655 (2021).

7. Zhang, H. et al. Probabilistic analysis of water-sealed performance in underground oil storage considering Spatial variability of hydraulic conductivity. *Sci. Rep.* **12** (1), 13782 (2022).

8. Agbaje, S., Zhang, X., Patelli, E., Ward, D. & Dhimitri, L. Random field failure and post-failure analyses of vertical slopes in soft clays. *Comput. Geotech.* **166**, 106037 (2024).

9. Al-Juboori, M. & Datta, B. Optimum design of hydraulic water retaining structures incorporating uncertainty in estimating heterogeneous hydraulic conductivity utilizing stochastic ensemble surrogate models within a multi-objective multi-realisation optimisation model. *J. Comput. Des. Eng.* **6** (3), 296–315 (2019).

10. Foldager, F. F. et al. Modeling soil aggregate fracture using the discrete element method. *Soil Tillage. Res.* **218**, 105295 (2022).

11. Goharzay, M., Noorzad, A., Ardakani, A. M. & Jalal, M. Computer-aided SPT-based reliability model for probability of liquefaction using hybrid PSO and GA. *J. Comput. Des. Eng.* **7** (1), 107–127 (2020).

12. Terzaghi, K. & Peck, R. B. Soil Mechanics. Engineering Practice (John Wiley and Sons, Inc., 1948).

13. Sellmeijer, J. B. On the Mechanism of Piping Under Impervious Structures (Technische Universiteit Delft, 1988).

14. Schmertmann, J. H. The no-filter factor of safety against piping through sands. In *Judgment and innovation: the heritage and future of the geotechnical engineering profession* 65–132. (2000).

15. Hoffmans, G. J. C. M. An overview of piping models. In *Proceedings of the 7th International Conference on Scour and Erosion* (3–18). (CRC Press, 2014).

16. ICOLD. Bulletin 164 Internal erosion of existing dams, levees and dikes, and their foundations. Volume 1: Internal erosion processes and engineering assessment. Paris: International Commission on Large Dams. (Publication of final preprint in English and French issued on 19 February 2015, first preprint issued in English, 22 January 2013). (2015).

17. Fell, R., Wan, C. F., Cyganiewicz, J. & Foster, M. Time for development of internal erosion and piping in embankment dams. *J. Geotech. GeoEnviron. Eng.* **129** (4), 307–314 (2003).

18. Benaissa, K. et al. Predicting initial erosion during the hole erosion test by using turbulent flow CFD simulation. *Appl. Math. Model.* **36** (8), 3359–3370 (2012).

19. El Shamy, U. & Aydin, F. Multiscale modeling of flood-induced piping in river levees. *J. Geotech. GeoEnviron. Eng.* **134** (9), 1385–1398 (2008).

20. Guo, Y. & Yu, X. Analysis of surface erosion of cohesionless soils using a three-dimensional coupled computational fluid dynamics–discrete element method (CFD–DEM) model. *Can. Geotech.* **J. 56** (5), 687–698 (2019).

21. Peng, S. & Rice, J. D. Inverse analysis of laboratory data and observations for evaluation of backward erosion piping process. *J. Rock Mech. Geotech. Eng.* **12** (5), 1080–1092 (2020).

22. Cho, S. E. Probabilistic analysis of seepage that considers the spatial variability of permeability for an embankment on soil foundation. *Eng. Geol.* **133–134** (none), 30–39. (2012).

23. Liu, L. L. et al. Effects of Spatial autocorrelation structure of permeability on seepage through an embankment on a soil foundation. *Comput. Geotech.* **87**, 62–75 (2017).

24. Tan, X., Wang, X., Khoshnevisan, S., Hou, X. & Zha, F. Seepage analysis of Earth dams considering Spatial variability of hydraulic parameters. *Eng. Geol.* **228**, 260–269 (2017).

25. Liang, Y. et al. Numerical simulation of backward erosion piping in heterogeneous fields. *Water Resour. Res.* **53** (4), 3246–3261 (2017).

26. Hekmatzadeh, A. A., Zarei, F., Johari, A. & Haghighi, A. T. Reliability analysis of stability against piping and sliding in diversion dams, considering four cutoff wall configurations. *Comput. Geotech.* **98**, 217–231 (2018).

27. Chi, F., Breul, P., Carvajal, C. & Peyras, L. Stochastic seepage analysis in embankment dams using different types of random fields. *Comput. Geotech.* **162**, 105689 (2023).

28. Foster, M., Fell, R. & Spannagle, M. The statistics of embankment dam failures and accidents. *Can. Geotech. J.* **37** (5), 1000–1024 (2000).

29. Richards, L. A. Capillary conduction of liquids through porous mediums. *Physics* **1** (5), 318–333 (1931).

30. Van Genuchten, M. T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. *Soil Sci. Soc. Am. J.* **44** (5), 892–898 (1980).

31. Terzaghi, K. A. R. L. Soil mechanics–A new chapter in engineering science. The 45th James Forrest lecture. *J. Institution Civil Eng.* **12**, 106–142 (1939).

32. Sellmeijer, H., de la Cruz, J. L., van Beek, V. M. & Knoeff, H. Fine-tuning of the backward erosion piping model through small-scale, medium-scale and IJkdijk experiments. *Eur. J. Environ. Civil Eng.* **15** (8), 1139–1154 (2011).

33. Hoffmans, G. & Van Rijn, L. Hydraulic approach for predicting piping in dikes. *J. Hydraul. Res.* **56** (2), 268–281 (2018).

34. Wang, R., Chen, G. & Liu, Y. Internal piping erosion evaluation of embankment dam considering the Spatial variability of soil properties. *ASCE-ASME J. Risk Uncertain. Eng. Syst. Part. A: Civil Eng.* **10** (2), 04024006 (2024).

35. Der Kiureghian, A. & Ke, J. B. The stochastic finite element method in structural reliability. *Probab. Eng. Mech.* **3** (2), 83–91 (1988).

36. Carsel, R. F. & Parrish, R. S. Developing joint probability distributions of soil water retention characteristics. *Water Resour. Res.* **24** (5), 755–769 (1988).

37. Mainguenaud, F. et al. A probabilistic approach to levee reliability based on sliding, backward erosion and overflowing mechanisms: application to an inspired Canadian case study. *J. Flood Risk Manag.* e12921. (2023).

38. Mouyeaux, A. et al. Probabilistic analysis of pore water pressures of an Earth dam using a random finite element approach based on field data. *Eng. Geol.* **259**, 105190 (2019).

39. ASTM Committee D-18 on Soil and Rock. *Standard Practice for Classification of Soils for Engineering Purposes (unified Soil Classification system)* 1 (ASTM international, 2017).

40. Cho, S. E. Probabilistic assessment of slope stability that considers the Spatial variability of soil properties. *J. Geotech. GeoEnviron. Eng.* **136** (7), 975–984 (2009).

41. Jiang, S. H., Li, D. Q., Cao, Z. J., Zhou, C. B. & Phoon, K. K. Efficient system reliability analysis of slope stability in spatially variable soils using Monte Carlo simulation. *J. Geotech. GeoEnviron. Eng.* **141** (2), 04014096 (2014).

42. Robbins, B. A. & The bedding angle of sand. Theses and dissertations. 4272. https://scholarsjunction.msstate.edu/td/4272 (2014).

43. Fu, Z. & Bauer, E. (eds). Dam engineering: Recent advances in design and analysis. (2021).

44. Feng, C. et al. Efficient slope reliability analysis under soil Spatial variability using maximum entropy distribution with fractional moments. *J. Rock Mech. Geotech. Eng.* **16** (4), 1140–1152 (2024).

45. Pan, Q. & Dias, D. Probabilistic evaluation of tunnel face stability in spatially random soils using sparse polynomial chaos expansion with global sensitivity analysis. *Acta Geotech.* **12** (6), 1415–1429 (2017).

46. Ahmed, A. S. H. R. A. F. & Soubra, A. H. Application of the subset simulation approach to spatially varying soils. *Risk Reliab. Geotech. Eng.* 591–625. (2018).

47. Phoon, K. K. Numerical recipes for reliability analysis—a primer. In: (eds Phoon, K-K. & Ching, J.) *Reliability-based Design in Geotechnical Engineering*. CRC, New York, 545. (2008).

48. Yuan, Y. L., Hu, C. M., Li, L., Xu, J. & Hou, X. H. Efficient slope reliability analysis using a surrogate-assisted normal search particle swarm optimization algorithm. *J. Comput. Des. Eng.* **11** (1), 173–194 (2024).
49. Johari, A. & Talebi, A. Stochastic analysis of rainfall-induced slope instability and steady-state seepage flow using random finite-element method. *Int. J. Geomech.* **19** (8), 04019085 (2019).

50. Ng, C. W. et al. Risk assessment of soil slope failure considering copula-based rotated anisotropy random fields. *Comput. Geotech.* **136**, 104252 (2021).

51. Liao, K., Wu, Y., Miao, F., Pan, Y. & Beer, M. Probabilistic risk assessment of Earth dams with spatially variable soil properties using random adaptive finite element limit analysis. *Eng. Comput.* **39** (5), 3313–3326 (2023).

52. Sun, R. et al. Failure probability analysis of high fill levee considering multiple uncertainties and correlated failure modes. *Sci. Rep.* **14** (1), 9436 (2024).

53. Griffiths, D. V. & Fenton, G. A. Seepage beneath water retaining structures founded on spatially random soil. *Geotechnique* **43** (4), 577–587 (1993).

54. Robbins, B. A., Griffiths, D. V. & Fenton, G. A. Random finite element analysis of backward erosion piping. *Comput. Geotech.* **138**, 104322 (2021).

55. King, D. M. & Perera, B. J. C. Morris method of sensitivity analysis applied to assess the importance of input variables on urban water supply yield–A case study. *J. Hydrol.* **477**, 17–32 (2013).

## Acknowledgements

The first author thanks the China Scholarship Council, China (CSC No. 202006060027) for providing Fendong Chi with a Ph.D. Scholarship for his research work.

## Author contributions

C. and C.C. prepared the original manuscript including the figures and tables. All the authors revised and reviewed the manuscript. The publication, as well as the listed author information has been approved by all the co-authors.

## Declarations

### Competing interests

The authors declare no competing interests.

### Additional information

Correspondence and requests for materials should be addressed to F.C.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025