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

HAL
open science

# Probabilistic analysis of three different internal erosion mechanisms - Application to a real earth dam

Fengdong Chi, Claudio Carvajal, Laurent Peyras, Pierre Breul

► To cite this version:

Fengdong Chi, Claudio Carvajal, Laurent Peyras, Pierre Breul. Probabilistic analysis of three different internal erosion mechanisms - Application to a real earth dam. Engineering Geology, 2025, 354, pp.108186. 10.1016/j.enggeo.2025.108186. hal-05460398

HAL Id: hal-05460398

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

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_2.png" /></figure>

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

Technical note

Probabilistic analysis of three different internal erosion mechanisms -
Application to a real earth dam

Fengdong Chi<sup>a,b</sup>, Claudio Carvajal<sup>a</sup>, Laurent Peyras<sup>a,*</sup>, Pierre Breul<sup>b</sup>

ARTICLE INFO

Keywords:

Earth dam

Probabilistic seepage

Backward erosion

Suffusion

Contact erosion

A B S T R A C T

Accident statistics for embankment dams show that internal erosion is a major cause of incidents involving these structures. However, internal erosion is a complex process and remains difficult to assess probabilistically.

This technical note presents a probabilistic analysis of internal erosion for an existing embankment dam. The spatial variability of soil properties is modeled using random field theory and incorporated into a finite element seepage analysis. The probabilistic results yield spatial distributions of the means and variances of key hydraulic parameters, such as flow velocity and hydraulic gradient.

By applying initiation criteria for three internal erosion mechanisms—backward erosion, suffusion, and contact erosion—the safety of the dam is assessed through the initiation probability (*P*<sub>f</sub>), derived from the probability distribution of the factor of safety (FoS) for each mechanism.

This approach is used to evaluate the probability of internal erosion initiation across different zones of the
dam. The results can then be used for risk analysis studies for embankment dams by assessing potential failure
scenarios due to internal erosion.

1. Introduction

Internal erosion corresponds to the detachment and transport of soil
particles by the action of seepage. Internal erosion can be initiated
through four mechanisms: backward erosion piping (BEP), suffusion,
contact erosion, and concentrated leak erosion (Zhang et al., 2016).
Statistics show that 38 % of earth-rock dam failures and 46 % of
embankment failures worldwide are caused by internal erosion (Foster
et al., 2000). Additionally, internal erosion plays a significant role in
river landform evolution. The study and description of the internal
erosion process represents a significant area of research activity within
the fields of geotechnical engineering, hydraulic engineering and geo-
morphology. The core of this issue is the capacity of the soil to resist
internal erosion and the mechanisms associated with it, both of which
have been dealt with by a large number of experiments and numerical
simulations described in the literature (Martínez-Moreno et al., 2018;
Zhang et al., 2022; Ma et al., 2024).

Significant studies have been conducted to evaluate the safety of
different internal erosion mechanisms (backward erosion piping, suf-
fusion, contact erosion, and concentrated leak erosion). Hekmatzadeh

et al. (2018) evaluated stability against backward erosion piping in four cutoff wall configurations using the random finite element method. Safety against piping is defined by the ratio of the critical hydraulic gradient to the exit hydraulic gradient, corresponding to Terzaghi criteria. Robbins et al. (2019) conducted a numerical investigation on the influence of spatial variation in permeability on BEP, and showed that the probability of erosion progression is linked to both the variance and correlation length of permeability random fields (RFs). BEP progression becomes less likely as the variance in permeability increases, and the correlation length in the flow direction decreases. Yang et al. (2019) proposed a finite difference model for suffusion based on the self-filtration process. Their probabilistic study analyzed the effects of variance, spatial correlation length, and cross-correlation of randomly distributed initial porosity and fines content on the eroded mass and permeability evolution. Kovacs (1981) proposed a geometry criterion for suffusion involving the ratio between d '85 (d<sub>85</sub> of the finer fraction of soil) and O<sub>50</sub> (effective constriction size of the coarse fraction). Based on Kovacs geometry criterion and the Terzaghi criterion, Li and Fannin (2008) proposed a hydro-mechanical criterion derived from tests on unstable materials. This criterion considered the reduction factor α,

<footnote><sup>a</sup> Aix-Marseille University, INRAE, RECOVER, 13182 Aix-en-Provence, France</footnote>

<footnote><sup>b</sup> Université Clermont Auvergne, CNRS, Clermont Auvergne INP, Institut Pascal, F-63000 Clermont-Ferrand, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: laurent.peyras@inrae.fr (L. Peyras).</footnote>
which was first proposed by Skempton and Brogan (1994) to illustrate the proportion of effective stress sustained by fines. Shire et al. (2014) established a discrete element model to analyze the fabric and effective stress distribution within idealized gap-graded samples. They found that particle-size distribution, percentage finer fraction, and relative density influence stress distribution. Shire et al. (2014) proposed that soil with a fines content of 35 % or more has an α value of 1, making it over-filled and not suffusive (internally stable). Researchers such as Brauns (1985) and Guidoux et al. (2010) have studied contact erosion and proposed expressions for hydraulic conditions that facilitate particle detachment and transport based on experimental results. Jia et al. (2022) used a two-phase seepage coupling model to study fluid seepage, particle migration, porosity, permeability, and displacement in embankment dam slopes affected by concentrated leak erosion.

This technical note presents the case study of an existing embankment dam to carry out a probabilistic analysis of internal erosion initiation mechanisms that takes into account the spatial variability of soil hydraulic properties. The dam studied has been the subject of previous research focusing on slope stability and seepage analysis (Mouyeaux et al., 2018; Mouyeaux et al., 2019; Guo et al., 2019; Guo et al., 2022). This technical note deals with internal erosion mechanisms and carries on from previous research work.

The originality of this study is to compare and evaluate the safety of a real earth dam – the Sillonnière Dam located in France – based on the criteria of three internal erosion mechanisms: BEP, suffusion, and contact erosion. The spatial variability of embankment properties is modeled with RFs based on a geostatistical analysis of construction control data. This study is novel in that it provides a comprehensive probabilistic framework integrating multiple internal erosion mechanisms on a real dam using spatially variable input parameters. Previous studies often focus on single mechanisms or simplified hypothetical dams. Our approach enables a spatially-resolved safety assessment (distribution of the FoS and initiation probability assessed at each point of the dam), identifying the most vulnerable zones within the dam structure for different internal erosion mechanisms.

The findings have important implications in the field of engineering geology, as they improve the understanding of how internal erosion may initiate under the spatial variability of material properties. The approach proposed provides more informed dam safety assessments and enhances decision-making regarding surveillance strategies and rehabilitation measures. Finally, this study aims to provide a probabilistic reliability analysis of the Sillonnière dam, and more generally provide a probabilistic assessment of internal erosion mechanisms contributing to dam engineering and risk assessment.

As seen below, some studies have proposed to evaluate dam safety by estimating the probability of internal erosion initiation, but by focusing on a single internal erosion mechanism or on a theoretical dam case.

Thus, the main objectives of this study are as follows: (1) to consider stationary RFs to represent the spatial variability of hydraulic parameters (dry density, saturated density, and porosity); (2) to estimate the spatial variability of permeability using available data on soil particle size distribution and dry density. (3) to evaluate the safety of a real dam in terms of initiation position and probability based on the criterion of three internal erosion mechanisms.

## 2. Materials and methods

### 2.1. Presentation of the dam

The Sillonnière dam (see Fig. 1), located in the west of France, was constructed in 1998. It is a 23.8 m high earth-filled dam fed by a catchment area of 50 km². Its width is around 140 m, and the length in the third dimension is 170 m for the crest. The dam closes a valley covered with alluvial deposits and can retain a reservoir of about 5 × 10⁶ m³. The normal and maximal water level is 50 m and 51.6 m, respectively. This work focuses on dam stability under steady-state flow conditions, and the water level is fixed at a level of 50 m.

Fig. 2 presents the main cross-section of the dam. It can be seen that the dam is formed by three different zones, including a core consisting of lean clay and two backfill zones consisting of well-graded gravel with silt in the upstream and downstream parts of the core respectively. In this work these three zones are named Core, Upstream shoulder (UPS), and Downstream shoulder (DOS), respectively. The difference in permeability between the core and the shoulders is less significant than with a zoned dam. A vertical and a horizontal drain were installed in the downstream part in order to collect the water flows and to lower the phreatic surface inside the dam. Although the foundation under the core wall was treated with a grout curtain, this part has been simplified in the present model. This choice was made after comparing simulations with and without the grouted curtain, which showed that the location of internal erosion initiation remained similar. As the grouted curtain provides a safer condition, removing it represents a conservative assumption.

The grain size distribution curves of these soil materials used for the UPS shoulder (in red colour), DOS shoulder (in green colour), core (in blue colour), and foundation (in purple colour) can be distinguished in Fig. 3.

### 2.2. Seepage modeling

This study adopts Darcy's law, according to which the flow velocity is proportional and linearly related to the hydraulic gradient:

$$i = \frac{v}{K_s} \qquad (1)$$

<figure><img src="image_4.png" /><figcaption>Fig. 1. Aerial view of the Sillonnière dam.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 2. The main cross-section of the dam and the location of pore water pressure cells.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. The grain size distribution curve of the zones in the Sillonnière dam.</figcaption></figure>

where: $i$ = hydraulic gradient, $v$ = flow velocity, $K_s$ = saturated permeability.

The hydraulic model in this research uses the finite element method to solve governing equations, primarily Richards' eq. (Richards, 1931), characterizing variably saturated flow through porous media:

$$ \frac{\partial}{\partial x} \left[ K_x(\theta) \frac{\partial h}{\partial x} \right] + \frac{\partial}{\partial y} \left[ K_y(\theta) \frac{\partial h}{\partial y} + 1 \right] = 0 \quad (2) $$

where $h$ = water pressure head, $\theta$ = volumetric water content and $K_x$ and $K_y$ are the permeability in the horizontal and vertical directions, respectively.

The unsaturated condition in earth dams typically requires unsaturated seepage analysis, as adopted in most studies (Tan et al., 2017), by incorporating the soil-water retention curve (SWRC) through the van Genuchten model (van Genuchten, 1980). In this study, the van Genuchten model is also adopted to represent the unsaturated behavior of soils in different dam zones.

$$ \theta_e(\psi) = \frac{\theta - \theta_r}{\theta_s - \theta_r} = \frac{1}{[1 + (\psi/a)^n]^m} \quad (3) $$

$$ K = K_s \theta_e^{1/2} \left[ 1 - \left( 1 - \theta_e^{1/m} \right)^m \right]^2 \quad (4) $$

where $\theta$ is the volumetric water content of the soil, $\psi$ is the matric suction, $\theta_s$ is the saturated volumetric water content, $\theta_r$ is the residual volumetric water content, and parameters $a$, $n$, and $m$ are the fitting parameters of the volumetric water content function, $\theta_e(\psi)$ is a function of the volumetric moisture content of the soil.

## 2.3. Internal erosion mechanisms

This section presents three common internal erosion mechanisms - backward erosion piping (BEP), suffusion, and contact erosion - and their initiation criteria. Concentrated leakage erosion is excluded as it is irrelevant to the Sillonnière dam, which has no cracks or other discontinuities.

### 2.3.1. Backward erosion piping (BEP)

BEP typically occurs in a dam's foundation but can also develop within the embankment. It starts at the downstream free surface and erodes upstream, forming a pipe. BEP exposes the foundation's sand layer, providing a free outlet at the dam toe. Terzaghi (1922) proposed a critical exit gradient for piping initiation in internally stable soils with good grain size distributions, see Eq. (5).

$$ i_{cr} = (1 - n_p) \frac{\rho_s - \rho_w}{\rho_w} = \frac{\gamma'}{\gamma_w} \quad (5) $$
where $n_p$ is the porosity; $\rho_s$ and $\rho_w$ are the densities of the soil solids and water, respectively; $\gamma'$ is the effective unit weight of soil, and $\gamma_w$ is the unit weight of water. In most cases, this equation yields values of $i_{cr}$ around 1.0–1.1.

The factor of safety (FoS) is calculated by comparing the local hydraulic gradient $i$ to the critical hydraulic gradient $i_{cr}$:

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

### 2.3.2. Suffusion

Suffusion involves the transport of fine particles through coarser material pores due to seepage. Unlike BEP, suffusion initiates with the migration of fine particles through the pore spaces between larger particles under seepage flow. As seepage velocity increases, these preferential flow paths may evolve into continuous channels, potentially leading to internal erosion (Yang et al., 2019). Skempton and Brogan (1994) noted that piping in unstable sandy gravels can occur at lower hydraulic gradients than Terzaghi's theoretical critical gradient, introducing a stress-reduction factor to the critical exit gradient:

$$ i_{cr} = \alpha \frac{\gamma'}{\gamma_w} \qquad (7) $$

The reduction stress factor $\alpha$ is related to the mobility of the finer particles and is dependent on Kovacs's geometric criterion $d$ '85/$O_{50}$ (Kovacs, 1981), where $d$ '85 is the $d_{85}$ of the finer fraction of soil and $O_{50}$ is the effective constriction size of the coarse fraction. Li and Fannin (2008) proposed a set of empirical formulas:

$$ \alpha = 3.85 \frac{d_{85}}{O_{50}} - 0.616 \qquad (8) $$

$$ O_{50} = 4 * \frac{n_c}{1 - n_c} * D_h^c \qquad (9) $$

$$ n_c = n_p + F_n(1 - n_p) \qquad (10) $$

$$ D_h^c = \frac{1}{\sum \frac{\Delta F_c}{D_i^c}} \qquad (11) $$

where $\alpha_D$ = shape coefficient (6 for rounded particle, 7–9 for angular particle); $n_c$ = porosity of skeleton; $n_p$ = porosity of soil; $F_n$ = mass passing (%) at diameter $D_n$; $D_n$ is the separating diameter (mm) which separates the finer fraction and coarser fraction. $D c h$ = the effective diameter of the coarse skeleton; $D c i$ = average diameter in the $i$-th interval of the particle size distribution curve of coarse fraction; $\Delta F c i$ = the weight of grains in the $i$-th interval of the particle size distribution curve of the coarse fraction. Kézdi geometry criteria were used to divide the soil into two components which are fine fraction and coarse fraction (Kézdi, 1979).

### 2.3.3. Contact erosion

Contact erosion occurs when fine particles are eroded through a coarse soil layer due to flow velocity differences. It requires geometric and hydraulic conditions: particles must pass through the coarse layer, and flow velocity must be sufficient (Guidoux et al., 2010). This process, similar to concentrated leak erosion, involves fine soil eroding into gravel and requires specific assessment methods. Brauns (1985) derived a formula to estimate the critical velocity based on layer characteristics, which influences hydraulic conditions and potential erosion:

$$ U_{crit} = 0.65 * n_D \sqrt{\left(\frac{\rho_s - \rho_w}{\rho_w}\right)gd_{50}} \qquad (12) $$

with $n_D$ being the porosity of the gravel layer, $\rho_s$ (kg/m³) the density of the sand particles, $\rho_w$ (kg/m³) the water density, and $d_{50}$ (m) the median

diameter of the sand grading curve.

### 2.4. Data for probabilistic analysis

This section presents the data available for the probabilistic modeling of the soil properties involved in the assessment of internal erosion initiation mechanisms. Parameters for which spatially distributed data are available from field investigations (e.g., dry density, saturated density, and porosity) are modeled as spatial random fields to account for their heterogeneity within the dam body. Permeability is also treated as a spatial random field, as it is estimated using the empirical relationship proposed by Chapuis and Aubertin (2003), which depends on density and porosity. In the suffusion mechanism, a reduction factor derived from an empirical function of porosity (Li and Fannin, 2008) is likewise considered as a spatial random field.

For the contact erosion mechanism, the grain diameter $d_{so}$ is modeled as a random variable, since spatially distributed measurements are not available. Its statistical parameters are derived from our previous study (Mouyeaux et al., 2019), which analyzed 27 to 30 grain size distribution curves to characterize material variability across dam zones.

Table 1 provides a summary of the parameters used in the probabilistic analysis and indicates whether they are modeled as random fields (RFs) or random variables (RVs).

#### 2.4.1. Data for probabilistic modeling of permeability

The probabilistic modeling of permeability is based on the method proposed by (Mouyeaux et al., 2019), based on the expressions developed by Chapuis and Aubertin (2003) (see Eq. (13)).

$$ K_s = \sqrt{10} * \left( \frac{e^3}{(1 + e)G_s^2 S_s^2} \right) \qquad (13) $$

$$ e = \frac{n_p}{1 - n_p} \qquad (14) $$

where $K_s$ is the saturated permeability, $e$ is the void ratio, $n_p$ is the porosity, $G_s$ is specific gravity, and $S_s$ is the specific surface. The specific surface $S_s$ can be estimated using the method proposed by Chapuis and Legare (1992), which links the grain size distribution curve and specific gravity. For the Sillonnière dam, the core is made of lean clay whose void ratio is 0.54. The Upstream, Downstream, and Foundation are made of well-graded gravel with silt whose void ratio is 0.44. The estimated permeability values and the corresponding parameters for the van Genuchten unsaturated model are presented in Table 2.

<table><caption>**Table 1**
Summary of parameters used in the probabilistic internal erosion analysis.</caption><thead><tr><th>Types of parameters</th><th>Parameters</th><th>Modeling types</th><th>Remarks</th></tr></thead><tbody><tr><td rowspan="5">Input parameters in seepage analysis</td><td>Dry density ρ<sub>d</sub></td><td>RFs</td><td>Based on field investigation</td></tr><tr><td>Saturated density ρ<sub>sat</sub></td><td>RFs</td><td>Based on field investigation</td></tr><tr><td>Porosity n<sub>p</sub></td><td>RFs</td><td>Based on field investigation</td></tr><tr><td>Permeability K<sub>s</sub></td><td>RFs</td><td>Estimated by empirical function (Chapuis and Aubertin, 2003)</td></tr><tr><td>Reduction factor α (suffusion mechanism)</td><td>RFs</td><td>Estimated by empirical function (Li and Fannin, 2008)</td></tr><tr><td rowspan="2">Parameters for internal erosion criteria</td><td>Median diameter of the soil grading curve d<sub>so</sub></td><td rowspan="2">RVs</td><td rowspan="2">Based 27–30 grain size distribution curves (Mouyeaux et al., 2019)</td></tr><tr><td>(contact erosion mechanism)</td></tr></tbody></table>
<table><caption>**Table 2**
Soil parameters for the hydraulic model.</caption><thead><tr><th>Zone</th><th>Ss (m²/kg)</th><th>G<sub>s</sub></th><th>K<sub>x</sub>: K<sub>y</sub></th><th>K<sub>x</sub> (m/s)</th><th>K<sub>y</sub> (m/s)</th><th>θ<sub>f</sub></th><th>θ<sub>s</sub></th><th>α (kPa)</th><th>n</th></tr></thead><tbody><tr><td>Upstream</td><td>72.0</td><td>2.7</td><td>2</td><td>4.4E-6</td><td>2.2E-6</td><td>1E-4</td><td>0.25</td><td>50</td><td>1.441</td></tr><tr><td>Downstream</td><td>45.3</td><td>2.7</td><td>2</td><td>1.6E-5</td><td>7.8E-6</td><td>1E-4</td><td>0.24</td><td>50</td><td>1.441</td></tr><tr><td>Core</td><td>930.5</td><td>2.7</td><td>5</td><td>5.0E-8</td><td>1.0E-8</td><td>1E-4</td><td>0.31</td><td>40</td><td>1.332</td></tr><tr><td>Foundation</td><td>49.3</td><td>2.7</td><td>1</td><td>1.0E-5</td><td>1.0E-5</td><td>1E-4</td><td>0.23</td><td>50</td><td>1.441</td></tr></tbody></table>

This study uses the Karhunen-Loève (K-L) method to generate random fields of soil parameters as it has been acknowledged as a highly efficient method. It uses the single exponential auto-correlation function (SNX) to simulate soil property variability.

### 2.4.2. Data for probabilistic modeling of the BEP mechanism

As mentioned in Section 2.4.1, estimated permeability is influenced by porosity $n_p$, while the Terzaghi critical gradient (see Section 2.3.1) depends on sand particle density $\rho_s$ and porosity $n_p$. The experiment shows $\rho_s = 2.7$ g/cm³ for both materials (lean clay and well-graded gravel with silt). Therefore, the spatial variability of dry density $\rho_d$, saturated density $\rho_{sat}$, and porosity $n_p$ is considered (see Eq. (15) and Eq. (16)) to assess the spatial variability of the critical gradient (see Eq. (5)).

Log-normal random fields of dry density and saturated density were generated for each zone of the dam (UPS, Core and DOS) on the basis of means, standard deviations and correlation lengths calculated from the distributions of compaction control measures. The correlation lengths in the vertical and horizontal directions are determined using the auto-correlation function in the geostatistical analysis of field data, shown in Table 3.

$$n_p = 1 - \frac{\rho_d}{\rho_s} \qquad (15)$$

$$\rho_{sat} = \rho_s + \left(1 - \frac{\rho_s}{\rho_w}\right) \bullet \rho_w \bullet n_p \qquad (16)$$

Fig. 4 presents the histogram and probability density function (PDF) of the generated values for porosity $n_p$ in the UPS, DOS, and Core zone.

### 2.4.3. Data for probabilistic modeling of the suffusion mechanism

The criterion for suffusion initiation (see Section 2.3.2) employed the reduction factor $\alpha$, which is relative to the soil particle size. The reduction factor $\alpha$, discussed here, is a key parameter calculated from porosity and grain size distribution using Skempton and Brogan's (1994) method.

For the well-graded soil materials (e.g. UPS, DOS and Foundation), the RFs for $\alpha$ can be derived from RFs for porosity $n_p$. Based on the filtering rule proposed by Kézdi (1979), the soil can be separated into two components with different sizes of soil particles and the suitable separating diameter $D_n$ can be identified. In view of the above, the initial reduction factor is $\alpha_{l\_UPS} = 0.604$, $\alpha_{l\_DOS} = 1.052$, $\alpha_{l\_Foundation} = 1.271$ (see Table 4). Because the maximum value of $\alpha$ is 1, the final modified value is $\alpha_{l\_UPS} = 0.604$, $\alpha_{l\_DOS} = 1$, $\alpha_{l\_Foundation} = 1$.

For the gap-graded soil material (e.g. Core), the gradation is

considered internally unstable by the Kézdi criteria, but it appears to be a little conservative in evaluating the potential for internal instability. Using the Discrete Element modeling (DEM) of spheres, Shire et al. (2014) estimated by that soils with fines content of 35 % or more have a value of 1.0 and are over-filled and not suffusive (internally stable). Therefore, the reduction factor $\alpha = 1$ is for the Core zone.

To sum up, the UPS zone is the only part made of unstable material with a reduction factor smaller than 1. Compared with the BEP mechanism, the UPS zone should taken into account to study the suffusion mechanism.

As mentioned above, the probabilistic modeling of the reduction factor $\alpha$ is done directly by linking the reduction factor with the porosity. Fig. 5 illustrates the histogram and fitted CDF curve for the reduction factor $\alpha$ in the UPS zone obtained in Fig. 4 which illustrates the histogram for porosity $n_p$.

### 2.4.4. Data for probabilistic modeling of the contact erosion mechanism

The contact erosion initiation criterion considers that the porosity and soil particle size and the FoS of contact erosion are linked to critical velocity (see Section 2.3.3). According to the experiment data, the density of the sand particles is equal to 2.7 g/cm³. The distribution of critical velocity is therefore influenced by the porosity $n_D$ (coarser material: Drain zone) and the median diameter $d_{50}$ of the sand grading curve (finer material: DOS or Foundation zone).

The spatial variability of median diameter $d_{50}$ can be considered in this study, e.g., $d_{50}$ for Foundation material can be considered as a random variable following normal distribution with the mean value at 17 mm, and COV with a value of 1.5. Therefore, the probabilistic modeling of critical velocity is done directly by linking the critical velocity with median diameter and porosity, shown in Eq. (12).

Contact erosion occurs between coarse and fine soils due to differences in flow velocity, causing fine soils to deteriorate and pass through the coarse layer. In this case study, the interfaces between the horizontal Drain and DOS, as well as between the Drain and the Foundation (see Fig. 6), are examined for the following reasons:

1) Contact erosion is more likely between coarse and fine soils, with a higher initiation probability when there are significant differences in particle size. The boundary of the Drain is studied due to the large material differences with DOS and the Foundation.

2) J. Brauns' research on horizontal flow erosion (Brauns, 1985) provides a criterion for contact erosion, making the horizontal drain boundary more suitable for analysis than the vertical boundary.

<table><caption>**Table 3**
The statistical properties of soil (dry density, porosity, and saturated density).</caption><thead><tr><th rowspan="2">Zone</th><th colspan="3">Dry density ρ<sub>d</sub> (g/cm<sup>3</sup>)</th><th colspan="3">Porosity n<sub>p</sub></th><th colspan="3">Saturated density ρ<sub>sat</sub> (g/cm<sup>3</sup>)</th><th rowspan="2">Correlation length (m)</th></tr><tr><th>Mean</th><th>Variation</th><th>COV</th><th>Mean</th><th>Variation</th><th>COV</th><th>Mean</th><th>Variation</th><th>COV</th></tr></thead><tbody><tr><td>UPS</td><td>1.99</td><td>0.004</td><td>3.18 %</td><td>0.263</td><td>0.00055</td><td>8.92 %</td><td>2.25</td><td>0.00159</td><td>1.77 %</td><td>L<sub>h</sub> = 78.1<br>L<sub>v</sub> = 7.8<br>L<sub>h</sub> = 13.0<br>L<sub>v</sub> = 1.5<br>L<sub>h</sub> = 4.9<br>L<sub>v</sub> = 1.9</td></tr><tr><td>Core</td><td>1.83</td><td>0.004</td><td>3.46 %</td><td>0.322</td><td>0.00055</td><td>7.28 %</td><td>2.11</td><td>0.00159</td><td>1.89 %</td><td></td></tr><tr><td>DOS</td><td>2.05</td><td>0.005</td><td>3.45 %</td><td>0.241</td><td>0.00069</td><td>10.9 %</td><td>2.29</td><td>0.00199</td><td>1.5 %</td><td></td></tr><tr><td>Foundation</td><td>2.03</td><td>-</td><td>-</td><td>0.25</td><td>-</td><td>-</td><td>2.28</td><td></td><td></td><td></td></tr><tr><td>Drain</td><td>2.16</td><td>-</td><td>-</td><td>0.2</td><td>-</td><td>-</td><td>2.36</td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Fig. 4. Histogram and PDF of the porosity in the UPS, DOS, and Core zone.</figcaption></figure>

<table><caption>**Table 4**

Calculation of the stress reduction factor $\alpha$ with separating diameters.</caption><thead><tr><th>Parameters</th><th>UPS</th><th>DOS</th><th>Foundation</th></tr></thead><tbody><tr><td>The median diameter of the soil grading curve:<br>d<sub>50</sub> (mm)</td><td>8</td><td>12</td><td>17</td></tr><tr><td>Shape coefficient: a<sub>D</sub></td><td>8</td><td>8</td><td>8</td></tr><tr><td>d '85 — d<sub>85</sub> of the finer fraction of soil: d<sub>2</sub> (mm)</td><td>1.81</td><td>2.74</td><td>3.295</td></tr><tr><td>d '15 — d<sub>15</sub> of the finer fraction of soil: d<sub>1</sub> (mm)</td><td>0.0015</td><td>0.0018</td><td>0.006</td></tr><tr><td>d (mm)</td><td>5.744</td><td>8.972</td><td>10.5</td></tr><tr><td>D c h — Effective Diameter (mm)</td><td>11.30</td><td>14.604</td><td>18.216</td></tr><tr><td>Effective constriction size of the coarse fraction<br>O<sub>50</sub> (mm)</td><td>5.712</td><td>6.326</td><td>6.7233</td></tr><tr><td>Diameter D<sub>n</sub>, which separates finer and coarser<br>soil (mm)</td><td>3</td><td>4</td><td>5</td></tr><tr><td>F<sub>n</sub> — Mass passing (%) at diameter D<sub>n</sub></td><td>0.328</td><td>0.295</td><td>0.243</td></tr><tr><td>The porosity of coarser fraction n<sub>c</sub></td><td>0.503</td><td>0.4642</td><td>0.4247</td></tr><tr><td>The initial reduction factor α<sub>i</sub></td><td>0.604</td><td>1.0515</td><td>1.271</td></tr><tr><td>The modified reduction factor α</td><td>0.604</td><td>1.000</td><td>1.000</td></tr></tbody></table>

3) Tha factor of safety (FoS) for contact erosion depends on seepage flow velocity. Contact erosion occurs when the flow velocity exceeds the critical flow velocity (FoS < 1). Therefore, since the horizontal drain boundary has the highest flow velocity (see Fig. 9), it is considered the most vulnerable area.

## 3. Numerical calculations and results

### 3.1. Deterministic analysis results

The deterministic seepage analysis was performed using the finite element software OPTUM G2. This software enables the integration of user-developed procedures, which is highly beneficial for probabilistic analysis. Before implementing the probabilistic model, it was important to verify that the hydraulic model, based on the finite element method, produced reliable results during a deterministic seepage calculation with constant permeability (see Table 2).

The hydraulic model geometry is illustrated in Fig. 7, where the structure and foundation were meshed into 8510 triangular elements, each with three nodes. We consider steady-state analysis under normal operating conditions. The downstream water level was set at 31.35 m,

<figure><img src="image_8.png" /><figcaption>Fig. 5. Histogram and fitted CDF for the reduction factor $\alpha$ in the UPS zone.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 6. Studied contact erosion interface.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. The Finite element model in OPTUM G2.</figcaption></figure>

based on the average PZ4 piezometer readings, while the normal reservoir level of 50 m was used as the upstream boundary condition.

The results of the deterministic analysis (Figs. 8 to 11) provide key insights for subsequent probabilistic analysis.

Fig. 8 illustrates the distribution of pore water pressure and the saturation line, which are used to delineate the saturated zone. Fig. 9 shows maximum flow velocity around the horizontal drain, which is useful for identifying the potential initiation position for contact erosion. Fig. 10 reveals a higher hydraulic gradient in the core zone, indicating potential initiation points for BEP and suffusion, crucial for identifying weak spots in the dam for probabilistic analysis.

The critical hydraulic gradient was calculated using Eq. (5), and the Factor of Safety (FoS) was then derived for each internal erosion mechanism accordingly. Figs. 11(a) and 11(b) show the deterministic FoS contours for the BEP and suffusion mechanisms, respectively. These results are based on the mean values of the input parameters and provide an initial identification of potential initiation zones under nominal, steady-state conditions. The upstream toe and the core are highlighted as the most vulnerable areas for BEP and suffusion, respectively.

In addition, special attention should be given to zones where the FoS is close to 1. Although these areas are not classified as unstable in a deterministic sense, they represent regions of marginal stability and are

more susceptible to failure under slight parameter fluctuations. As discussed in Section 3.2, these zones often correspond to areas with higher local initiation probability ($P_f$), and thus warrant detailed probabilistic evaluation and, where necessary, targeted risk mitigation measures.

## 3.2. Probabilistic analysis

This study employed the Monte Carlo Simulation (MCS) method using MATLAB to assess the reliability of a finite element model developed in OPTUM G2. Random fields for key parameters such as dry density, porosity, and saturated density were generated and incorporated into the model to simulate spatial variability. We can obtain at each mesh point an evaluation of the various variables of interest for analyzing internal erosion (hydraulic gradient, flow velocity, etc.) by taking into account the spatial variability of the soil properties for each numerical simulation.

The Factor of Safety (FoS) for internal erosion was evaluated using 1000 runs to ensure result stability. Convergence was observed after approximately 100 runs, with stable statistical metrics achieved by 500. Thus, 1000 simulations provide robust and reliable estimates for probabilistic analysis (see Fig. 12).

<figure><img src="image_11.png" /><figcaption>Fig. 8. The contour of pore water pressure.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 9. The contour of flow velocity.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. The contour of the hydraulic gradient.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. The contours of the deterministic internal erosion analysis results.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Convergence of the estimated mean and standard deviation of the FoS.</figcaption></figure>

### 3.2.1. Probabilistic analysis for BEP

The contour of the mean value of probabilistic BEP analysis results – hydraulic gradient and FoS against BEP - are illustrated in Figs. 13 and 14. Based on the FoS from 1000 runs for each node, the FoS was found to

follow a log-normal distribution, as determined using MATLAB's Distribution Fitter app. The initiation probability at each node was then calculated, and the contour of BEP initiation probability is shown in Fig. 15.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Mean of the hydraulic gradient.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Mean of FoS (BEP).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Initiation probability $P_f$ (BEP).</figcaption></figure>

For the BEP mechanism, we can see that the higher $P_f$ is obtained at the toe of the core where the FoS is close to 1. This highlights the importance of considering both the mean safety level and the associated variability when identifying critical zones for internal erosion risk (see detailed discussion in Section 4.1).

### 3.2.2. Probabilistic analysis for suffusion

The contour of the mean value of probabilistic BEP analysis results - FoS against suffusion - is illustrated in Fig. 16. The suffusion initiation probability contour is presented in Fig. 17.

For suffusion mechanism, We can see that the higher $P_f$ is also obtained at the toe of the core where the FoS is close to 1. A secondary risk

<figure><img src="image_19.png" /><figcaption>Fig. 16. Mean of the FoS (Suffusion).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 17. Initiation probability $P_f$ (Suffusion).</figcaption></figure>

zone is also identified at the upstream toe, though its associated initiation probability is substantially lower. Nevertheless, this area may still be of interest in a comprehensive risk assessment (see detailed discussion in Section 4.2).

### 3.2.3. Probabilistic analysis for contact erosion

Contact erosion was analyzed at the boundary of the horizontal drain, where the maximum mean flow velocity is located (see Fig. 18).

The interface boundary between the horizontal drain and the foundation comprises 106 nodes (see Fig. 6). For each node, the flow velocity, critical flow velocity, mean FoS, standard deviation of FoS, and initiation probability $P_f$ were calculated. Fig. 19 shows the probability density function (PDF) of FoS for the node with the lowest FoS, indicating a higher likelihood of contact erosion initiation. Fig. 20 presents a similar analysis for the interface boundary between the horizontal drain and the downstream section (see Fig. 6).

## 4. Discussion

### 4.1. Discussion of the results for BEP

The contours of the Factor of Safety (FoS) and initiation probability ($P_f$) for the backward erosion piping (BEP) mechanism (Figs. 14 and 15) indicate that the most vulnerable area is located at the toe of the Core zone, where both a lower FoS and a higher $P_f$ are observed. Although the Core material is classified as a stable soil (as mentioned in Section 2.4.3), resulting in a reduction factor $\alpha$ equal to 1, the hydraulic gradient in this zone remains high due to boundary conditions and geometric constraints.

Interestingly, for both the BEP and suffusion mechanisms, the lowest mean FoS values are found at the same location near the toe of the Core zone. The most critical node, located at coordinates (6.52, 29.42), exhibits a mean FoS of 1.282 and an associated initiation probability $P_f$ of 1.59 %. This highlights a key insight of the probabilistic framework: the presence of uncertainty can lead to a non-negligible initiation probability even when the mean FoS is greater than 1.

<figure><img src="image_21.png" /><figcaption>Fig. 19. PDF of FoS (Foundation - Drain).</figcaption></figure>

Therefore, relying solely on the mean FoS may be misleading in assessing the risk of internal erosion. A more comprehensive evaluation should incorporate the initiation probability, which integrates both the central tendency and the variability of input parameters. This reinforces the value of the probabilistic approach adopted in this study.

Fig. 21 shows the histogram and fitted CDF curve for the FoS of this node. It illustrates that the FoS results for the most dangerous node have higher relative dispersion over 1000 runs, ranging from 0.8 to 2.6. The areas of relative higher variability of FoS (i.e., with the higher standard deviations) are mostly located in the Core zone. This is due to the gradient of permeability in the Core zone which can be high or not, depending on the permeability values obtained during the simulations.

<figure><img src="image_22.png" /><figcaption>Fig. 18. Mean of flow velocity.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 20. PDF of FoS (DOS - Drain).</figcaption></figure>

## 4.2. Discussion of the results for suffusion mechanisms

In steady-state conditions, the risk profile for the suffusion mechanism differs from that of BEP in the Upstream (UPS) zone, where the material is classified as unstable based on the filtering rule. Despite this classification, the node with the highest initiation probability for suffusion coincides with that of the BEP mechanism—located at the toe of the Core zone, as shown in Figs. 15 and 17. This outcome suggests that the hydraulic conditions and geometric configuration of the Core zone are more unfavorable, although the UPS zone contains more erosion-prone material.

The second most critical node for suffusion is located at the toe of the UPS zone (coordinates: -20.99, 30.49), with a mean FoS of 1.278 and a corresponding initiation probability $P_f$ of 0.44 %. Again, this illustrates the importance of complementing deterministic FoS values with a probabilistic perspective. Although the mean FoS is greater than 1, the associated non-zero $P_f$ highlights potential susceptibility to suffusion due to variability in hydraulic and geotechnical parameters.

These findings further support the notion that risk assessment based solely on FoS can be misleading, particularly in marginally stable zones. By incorporating initiation probability, the probabilistic framework offers a more complete and nuanced understanding of internal erosion

risks across different dam zones.

For transient conditions, such as a rapid drawdown of the upstream water level, the risk of suffusion initiation in the UPS zone may increase due to the reversal of seepage gradients. Although not addressed in this study, transient flow modeling could offer deeper insights into upstream infiltration failure. Future work may explore this scenario using unsteady-state seepage analyses.

Fig. 22 presents the histogram and fitted CDF curve for the FoS of the second most critical node. The FoS values for this node range from 0.9 to 1.7 across 1000 runs, exhibiting relatively low dispersion. The Upstream zone shows less variability in FoS, with lower standard deviations, due to the low permeability gradient in this area. Water flow upstream of the dam is relatively stable compared to the core, and this stability, even when considering the spatial variability of permeability, leads to less fluctuation in the FoS.

## 4.3. Discussion of the results for contact erosion

Fig. 13 and Fig. 14 highlight the highest-risk nodes along the upper and lower boundaries of the horizontal drain. The node with the lowest FoS at the interface between the Drain and Foundation was further analyzed. The factor of safety (FoS) for contact erosion, specifically the node with the lowest FoS on the interface between the Drain and Foundation, is calculated using the following formula:

$$ 
\begin{aligned}
\text{FoS} = \frac{V_{er}}{V} &= \frac{0.65 \cdot n_D \sqrt{\left(\frac{\rho_s - \rho_w}{\rho_w}\right) g d_{50}}}{V_{max}} = \frac{0.65 \cdot 0.2^* \sqrt{1.7 \cdot 9.81 \text{ m/s}^2 \cdot 0.017 \text{ m}}}{5.42 \cdot 10^{-5} \text{ m/s}} \\
&= 1277
\end{aligned}
$$

These results indicate that contact erosion is unlikely to occur in the Sillonnière dam for the following reasons:

1) With the protection function of the Core zone, the flow velocity in the horizontal Drain is far lower than the critical flow velocity, even if the flow velocity near the horizontal Drain is higher than in other areas of the dam. Based on the above formula, the critical flow velocity required to initiate contact erosion is 0.069 m/s. The maximum flow velocity obtained in the numerical analysis is only $5.42 \cdot 10^{-5}$ m/s (corresponding to 4.68 m/day in Fig. 18), far below the threshold for contact erosion.

<figure><img src="image_24.png" /><figcaption>Fig. 21. Histogram and fitted CDF for the FoS of the most critical node (BEP & suffusion).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 22. Histogram and fitted CDF for the FoS of the 2nd most critical node (suffusion).</figcaption></figure>

2) The flow velocity in the DOS zone and horizontal drain is significantly reduced due to the Core zone's protective function. Additionally, the particle size distribution of the materials in the Core, DOS, and Foundation zones (see Fig. 3 and Table 4) further enhances stability. Overall, given the current configuration, contact erosion is unlikely to occur in the earth dam studied.

## 5. Conclusions

This paper presented a probabilistic analysis of three internal erosion initiation mechanisms - backward erosion piping (BEP), suffusion, and contact erosion - while accounting for the spatial variability of soil hydraulic properties in an existing dam. The study used a finite element hydraulic model developed with OPTUM G2, incorporating random fields derived from geostatistical analysis of available dam data. Monte Carlo simulations were employed to assess flow variability, and the probabilistic seepage analysis was used to evaluate the probability of internal erosion initiation.

The method presented in this technical note was applied to an existing dam. The key conclusions are as follows:

1) The method accounts for spatial variability in soil properties, providing maps showing spatial probabilistic distributions for hydraulic results, such as flow velocities, pore pressures, and hydraulic gradients at every point in the dam and its foundation.

2) Using hydraulic parameter results and internal erosion criteria from the literature, the method enables mapping of the factor of safety and initiation probabilities for different internal erosion mechanisms.

3) It highlights the dam's most vulnerable areas with the highest initiation probabilities and can be integrated into fault tree analyses typically used in dam risk assessments.

The study also suggested areas for improvement, including incorporating additional criteria for internal erosion assessment, considering the filtration phase, and refining the numerical modeling of erosion mechanisms. Furthermore, since this work focused on the initiation phase of internal erosion, the identification of high-risk zones provided valuable insights for improving monitoring strategies. More specifically, the outcomes of the probabilistic analysis can guide the placement of pore water pressure (PWP) sensors and flowrate monitoring devices in critical areas. Integrating such risk-informed monitoring measures could enhance early detection capabilities and aid more effective dam safety management.

## CRediT authorship contribution statement

**Fengdong Chi:** Writing – original draft, Visualization, Validation, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Claudio Carvajal:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Laurent Peyras:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Pierre Breul:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgement

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

## Data availability

Data will be made available on request.

## References

Brauns, J., 1985. Erosionsverhalten geschichteten bodens bei horizontaler durchströmung. Wasserwirtschaft 75 (10), 448–453.

Chapuis, R.P., Aubertin, M., 2003. On the use of the Kozeny Carman equation to predict the hydraulic conductivity of soils. Can. Geotech. J. 40 (3), 616–628.

Chapuis, R.P., Legare, P.P., 1992. A Simple Method for Determining the Surface Area of Fine Aggregates and Fillers in Bituminous Mixtures.

Foster, M., Fell, R., Spannagle, M., 2000. The statistics of embankment dam failures and accidents. Can. Geotech. J. 37 (5), 1000–1024.
Guidoux, C., Faure, Y.H., Beguin, R., Ho, C.C., 2010. Contact erosion at the interface between granular coarse soil and various base soils under tangential flow condition. *J. Geotech. Geoenviron. Eng.* 136, 741–750.

Guo, X., Dias, D., Carvajal, C., Peyras, L., Breul, P., 2019. A comparative study of different reliability methods for high dimensional stochastic problems related to earth dam stability analyses. *Eng. Struct.* 188, 591–602.

Guo, X., Dias, D., Carvajal, C., Peyras, L., Breul, P., 2022. Modelling and comparison of different types of random fields: case of a real earth dam. *Eng. Comput.* 38 (Suppl. 5), 4529–4543.

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

Jia, Y., Ding, Y., Wang, X., Zhang, J., Chen, X., 2022. A numerical analysis of the leakage characteristics of an embankment dam slope with internal erosion. *Front. Earth Sci.* 10, 866238.

Kézdi, A., 1979. Soil Physics-Selected Topics-Developments in Geotechnical Engineering-25 (No. Monograph).

Kovacs, G., 1981. Seepage Hydraulics. Elsevier Scientific Publishing Company, Amsterdam, 730p.

Li, M., Fannin, R.J., 2008. Comparison of two criteria for internal stability of granular soil. Can. Geotech. J. 45 (9), 1303–1309.

Ma, C., Peng, M., Zhang, L., Shi, Z., Zhou, J., Chen, H., Li, Z., 2024. Erosion, deposition and breach evolution of landslide dams composed of various dam material types based on flume tests. *Eng. Geol.* 337, 107598.

Martínez-Moreno, F.J., Delgado-Ramos, F., Galindo-Zaldívar, J., Martín-Rosales, W., López-Chicano, M., González-Castillo, L., 2018. Identification of leakage and potential areas for internal erosion combining ERT and IP techniques at the Negratín Dam left abutment (Granada, southern Spain). *Eng. Geol.* 240, 74–80.

Mouyeaux, A., Carvajal, C., Bressolette, P., Peyras, L., Breul, P., Bacconnet, C., 2018. Probabilistic stability analysis of an earth dam by Stochastic Finite Element Method based on field data. *Comput. Geotech.* 101, 34–47.

Mouyeaux, A., Carvajal, C., Bressolette, P., Peyras, L., Breul, P., Bacconnet, C., 2019. Probabilistic analysis of pore water pressures of an earth dam using a random finite element approach based on field data. *Eng. Geol.* 259, 105190.

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

Robbins, B.A., Griffiths, D.V., Fenton, G.A., 2019. Influence of spatially variable soil permeability on backward erosion piping. In: 7th International Symposium on Geotechnical Safety and Risk (ISGSR 2019) Taipei, Taiwan.

Shire, T., O'Sullivan, C., Hanley, K.J., Fannin, R.J., 2014. Fabric and effective stress distribution in internally unstable soils. *J. Geotech. Geoenviron. Eng.* 140 (12), 04014072.

Skempton, A.W., Brogan, J.M., 1994. Experiments on piping in sandy gravels. *Geotechnique* 44 (3), 449–460.

Tan, X., Wang, X., Khoshnevisan, S., Hou, X., Zha, F., 2017. Seepage analysis of earth dams considering spatial variability of hydraulic parameters. *Eng. Geol.* 228, 260–269.

Terzaghi, K., 1922. Der grundbruch an stauwerken und seine verhutung. *Wasserkraft* 17, 445–449.

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

Yang, J., Yin, Z.Y., Laouafa, F., Hicher, P.Y., 2019. Analysis of suffusion in cohesionless soils with randomly distributed porosity and fines content. *Comput. Geotech.* 111, 157–171.

Zhang, L., Peng, M., Chang, D., Xu, Y., 2016. Internal Erosion in Dams and their Foundations. *Dam Failure Mechanisms and Risk Assessment. John Wiley & Sons* Singapore Pte. Ltd, Singapore, pp. 67–126.

Zhang, D., Zhang, X., Du, W., 2022. DEM-FEM based numerical analysis on mechanical responses of sandy soil and pipeline to seepage erosion. *Eng. Geol.* 310, 106868.