# Analysis of volumetric internal erosion in cohesionless soils: Model, experiments and simulations

Rachel Gelet, Alaa Kodieh, Didier Marot, Ngoc-Son Nguyen

► To cite this version:

Rachel Gelet, Alaa Kodieh, Didier Marot, Ngoc-Son Nguyen. Analysis of volumetric internal erosion in cohesionless soils: Model, experiments and simulations. International Journal for Numerical and Analytical Methods in Geomechanics, 2021, 45 (18), pp.2780-2806. 10.1002/nag.3286. hal-04167924

HAL Id: hal-04167924

https://hal.science/hal-04167924

Submitted on 21 Jul 2023

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 Attribution 4.0 International License
Analysis of volumetric internal erosion in cohesionless soils: model, experiments and simulations

Rachel Gelet$^{a,*}$, Alaa Kodieh$^{a,b}$, Didier Marot$^{a}$, Ngoc-Son Nguyen$^{a}$

$^{a}$Institut GeM, UMR CNRS 6183, Université de Nantes, 44600 Saint-Nazaire, France

$^{b}$Lebanese University, Faculty of Engineering, Laboratory of Mathematics and Applications, Beirut, Lebanon

**Abstract**

A hydro-mechanical continuous model that accounts for volumetric internal erosion (otherwise called suffusion) is developed based on experimental observations and data to fill the gap between laboratory testing and field applications. The model is proposed based on the mixture theory applied to a two-phase four-species porous medium. The erodible soil is partitioned in two phases: one solid phase and one fluid phase. The solid phase is composed of non-erodible grains and erodible particles. The fluid phase is composed of water and fluidized particles. The modelling of internal erosion is contributed directly by mass transfer from the solid phase towards the fluid phase. Two mass transfer relationships are considered, the power-based relationship and the energy-based relationship, and a calibration procedure for all material parameters is proposed. Both relationships assume that the power dissipated by the flow controls the kinetics of suffusion and that detachment prevails upon self-filtration. Both relationships predict reasonably experimental data for several experimental tests on a gap-graded cohesionless soil tested under various hydraulic loading paths. The model has been numerically solved with the finite element method and its predictions are next compared with the experimental results of a physical model of dike. The model along with the energy-based relationship reproduces well the final amount of eroded mass. The spatial distributions of the total water head and the final percentage of fines are also smoothly reproduced.

**Keywords:** internal erosion, suffusion, poromechanics, hydro-mechanical coupling, finite element method, energy approach

1. Introduction

Internal erosion can be categorised into four types, based on the physics of the process and on its location within a hydraulic structure: concentrated leak erosion, backward erosion, contact erosion and suffusion [11]. Suffusion is a complex process that under seepage flow, the fine solid particles can be detached, transported, and for some of them blocked by the constrictions of the solid skeleton [26, 45]. It affects preferably cohesionless soils but silty-sands may also be concerned [1]. The finer fraction eroded and leaving the coarse matrix of the soil will further modify the porosity, the hydraulic conductivity and mechanical parameters of the soil [4]. The development of suffusion gives rise to a wide range of dam incidents [10] which trigger heavy maintenance costs.

Initially, experimental researches on suffusion were devoted to its initiation and its susceptibility [17, 16, 22, 7, 26]. In addition, the focus has been raised to the initiation and the development of suffusion under several hydraulic loadings paths, seepage lengths and stress state conditions [32, 48, 6].

Recently, stakeholders have formulated the need to model internal erosion as a mean to contribute to risk assessment analyses. To date, discrete and continuous numerical methods have been developed. Discrete

<footnote>*Principal corresponding author</footnote>

<footnote>Email addresses: rachel.gelet@univ-nantes.fr (Rachel Gelet), alaa.kodieh@etu.univ-nantes.fr (Alaa Kodieh), didier.marot@univ-nantes.fr (Didier Marot), ngocson.nguyen@univ-nantes.fr (Ngoc-Son Nguyen)</footnote>
approaches are key methods to understand better the physics of this complex coupled process [34, 15, 43, 39]; yet they are still restricted by a relatively low number of particles with respect to that of laboratory specimens or geotechnical structures. On the other hand, continuous models use relationships and averaged variables dedicated to the scale of a representative elementary volume [41, 38, 28, 8, 27, 3, 47, 44, 5]. Most of these studies describe internal erosion by a mass transfer law. Some of these laws describe the detachment and transport of finer soil particles from a porous medium towards an open exit such as a borehole or a cavity, i.e. the detached particles are not concerned by a possible self-filtration. In fact, few models or relationships are dedicated to the simultaneous detachment, transport and self-filtration that cannot be experimentally distinguished [47, 35, 45, 5, 19], and fewer are concerned by the experimental measurement of the averaged parameters [35, 5, 19]. In contrast to other works which lack a clear chain of validation from laboratory specimens, physical models and full-sale structures [47, 45], the present work is based on experimental observations and data to fill the gap between laboratory testing and field applications.

The paper presents the development of a suffusion porous-elastic model that addresses observed suffusion behaviours in laboratory tests (Section 2). In particular, two mass transfer laws dedicated to suffusion are presented [35, 19] which both assume that the power dissipated by the flow controls the kinetics of the suffusion process. Several suffusion tests on a susceptible gap-graded soil are described in Section 3; one test is used to calibrate the suffusion relationships and the other ones are used to assess their prediction capabilities. Finally, the comparison between the model predictions and the experimental results of a physical model of dike is given in Section 4.

## 2. A suffusion porous-elastic model

### 2.1. A two-phase four-species mixture

A suffusive cohesionless saturated soil is viewed as a two-phase, four-species, porous medium. The solid phase S is constituted by non-erodable grains gS and erodable particles pS; while the fluid phase F contains the fluidized particles pF and water wF:

$$S = \{gS, pS\}, \quad F = \{pF, wF\}. \qquad (1)$$

The solid phase is sometimes referred to as the solid skeleton. Motivation for the species segregation are two-fold: both non-erodable grains and erodable particles are assumed to contribute to the overall deformation within the solid phase while both the fluidized particles and water species diffuse within the fluid phase itself. Erodible particles and fluidized particles can transfer from one phase to the other. Indeed, fluidized particles might be redeposited or blocked to become non-erodable grains.

Various macroscopic masses and volumes are required to formulate the constitutive equations. Let the initial volume of the porous medium be $V_0$ and let $V = V(t)$ be its current volume. The current volume and mass of the species $kl \in \{gS, pS, pF, wF\}$ are denoted by $V_{kl}$ and $M_{kl}$, respectively. Various additional entities are attached to species, some are intrinsic like the intrinsic density $\rho_{kl}$, while others are partial like the partial density $\rho^{kl}$. Some are related to the current volume like the volume fraction $n^{kl}$, and others related to the initial volume like the volume content $v^{kl}$:

$$\rho_{kl} = \frac{M_{kl}}{V_{kl}}, \quad \rho^{kl} = \frac{M_{kl}}{V}, \quad n^{kl} = \frac{V_{kl}}{V}, \quad v^{kl} = \frac{V_{kl}}{V_0}. \qquad (2)$$

The corresponding entities associated with the phase $l \in \{S, F\}$ are defined by algebraic summation of individual contributions, e.g. the current volume $V_l$ and mass $M_l$, the partial density $\rho^l$, the volume fraction $n^l$ and the volume content $v^l$. Volume fractions satisfy the compatibility relations $\sum_l n^l = 1$ and $\sum_{kl} n^{kl} = 1$. The total density of the mixture is the sum of its partial contributions $\rho = \sum_{kl} \rho^{kl}$.

Other entities live in their phase such as the concentration of fluidized particles in the fluid phase $C_{pF}$ and the concentration of erodable particles in the solid phase $C_{pS}$,

$$C_{pF} = \frac{V_{pF}}{V_F} \quad \text{and} \quad C_{pS} = \frac{V_{pS}}{V_S}. \qquad (3)$$
Water and fluidized particles are assumed to have the same velocities as that of the fluid phase $v_{wF} = v_{pF} = v_F$; and erodible particles along with non-erodible grains are assumed to move with the solid skeleton $v_{gS} = v_{pS} = v_S$. In addition, all solid particles $gS, pS, pF$ will be assumed incompressible along with the water species $wF$. A variation from this simplified case can be found in Kodieh [18].

Following Biot's approach [2], the fluid phase is viewed as flowing through the solid skeleton and the mass flux per unit current area through the solid skeleton of a species $kl$ is defined as:

$$ M_{kl} = \rho^{kl}(v_{kl} - v_S). \qquad (4) $$

Accordingly, the volumetric flux per unit current area through the solid skeleton is defined by $J_{kl} = M_{kl}/\rho_{kl}$.

The solid phase is endowed by its own (infinitesimal) strain tensor $\epsilon$ and volumetric strain $\text{tr} \, \epsilon$ which are defined from the macroscopic displacement vector $\mathbf{u}$,

$$ \epsilon = \frac{1}{2} [\nabla \mathbf{u} + (\nabla \mathbf{u})^{\text{T}}], \quad \frac{d(\text{tr} \, \epsilon)}{dt} = \text{div} \, \mathbf{v}_S. \qquad (5) $$

Internal erosion may lead to changes of both skeleton deformation and porous space, comprising mechanical induced and erosion induced parts. Since this work focuses on the erosion process, the skeleton deformation is assumed to be reversible $\epsilon^r$, i.e. although suffusion may degrade the elastic material properties of the soil, the strain is assumed to be recoverable upon unloading. Plastic deformations will be considered as perspectives [33]. By contrast with standard poro-elastic models, the volume content evolution $dv^F$ is assumed to contain a reversible part due to mechanical loading $dv^{Fr}$ and explicitly an irreversible part due to internal erosion $dv^{Fi}$. This concept of irreversible porosity was also proposed in poro-elasticity with dissolution by Coussy [9]. The above assumptions are summarized by:

$$ \begin{aligned} d\epsilon &= d\epsilon^r \\ dv^F &= dv^{Fr} + dv^{Fi}. \end{aligned} \qquad (6) $$

These equations imply that in the absence of mechanical loading, the solid skeleton does not experience any deformation, but the fluid content may still be modified due to suffusion.

## 2.2. Balance equations

A single balance of momentum is required to describe the mechanical behaviour of the mixture as a whole in which $\sigma$ is the total stress and $g$ is the gravity,

$$ \text{div} \, \sigma + \rho g = 0. \qquad (7) $$

To highlight the mass interactions between the species or the phases, the balance of mass is introduced at the level of the species $kl$, with respect to the solid skeleton,

$$ \frac{d\rho^{kl}}{dt} + \rho^{kl}\text{div}\mathbf{v}_S + \text{div}\mathbf{M}_{kl} = \hat{\rho}^{kl}. \qquad (8) $$

The rate of mass supply $\hat{\rho}^{kl}$ represents the mass supply due to internal mass transfer, commonly due to physical changes between the phases under the action of flow during suffusion. Thus, the total rate of mass supply $\hat{\rho}$ is the sum of the mass supplies due to internal contributions which vanishes to zero for this closed system,

$$ \hat{\rho} = \sum_{kl} \hat{\rho}^{kl} = 0. \qquad (9) $$

Since non-erodible grains and water species do not transfer $\hat{\rho}^{gS} = \hat{\rho}^{wF} = 0$. On the other hand, erodible particles are transferred into fluidized particles so that $\hat{\rho}^{pS} = -\hat{\rho}^{pF}$, which is equivalent to stating that $\hat{\rho}^{S} = -\hat{\rho}^{F}$ at the phase level. Importantly, the computation of the erosion rate per unit volume $\hat{\rho}^{F}$ is physically constrained by the presence of erodible particles, i.e. $C_{pS} > 0$.
Four balance of mass equations are required to describe this two-phase four-species mixture. Since, the fluid pressure $p_F$, the fluid porosity $n^F$, the concentration of fluidized particles in the fluid phase $C_{pF}$ and the concentration of erodible particles in the solid phase $C_{pS}$ are of particular interest; the balance of mass for the fluid phase, for the solid phase, for the fluidized particles and for the erodible particles are expressed by:

$$ \frac{dv^F}{dt} + \mathrm{div}\mathbf{J}_F = \frac{\hat{\rho}^F}{\rho_S}, \qquad (10) $$

$$ \frac{dn^F}{dt} = (1 - n^F)\mathrm{div}\mathbf{v}_S + \frac{\hat{\rho}^F}{\rho_S}, \qquad (11) $$

$$ n^F \frac{dC_{pF}}{dt} + C_{pF} \frac{dv^F}{dt} + \mathrm{div}(C_{pF} \mathbf{J}_F) = \frac{\hat{\rho}^F}{\rho_S}, \qquad (12) $$

$$ \frac{dC_{pS}}{dt} = -\frac{1 - C_{pS}}{1 - n^F} \frac{\hat{\rho}^F}{\rho_S}. \qquad (13) $$

The derivation of these four equations uses the definitions introduced in eqs. (2), (3) and (4), and the balance of mass eq. (8). In particular, the balance of mass for the fluid phase is obtained by summing up the balances of mass of the fluidized particles $pF$ with that of the water $wF$. Similarly, the balance of mass for the solid phase is obtained by summing up the balance of mass of the erodible particles $pS$ with that of the non-erodible grains $gS$. Additional details may be found in Kodieh [18].

## 2.3. Poro-elastic coupling

Considering non-compressible grains, the poro-elastic constitutive equations can be expressed as a function of Lame's constants of the drained solid $\lambda^{DS}$ and $\mu^{DS}$:

$$ \begin{aligned} d\boldsymbol{\sigma} &= \lambda^{DS} d(\mathrm{tr}\boldsymbol{\epsilon}) \mathbf{I} + 2\mu^{DS} d\boldsymbol{\epsilon} - dp_F \mathbf{I} \\ dv^{Fr} &= d(\mathrm{tr}\boldsymbol{\epsilon}) \end{aligned} \qquad (14) $$

The erosion process is herein described by the irreversible volume content $v^{Fi}$ and by the rate of mass supply to the fluid phase $\hat{\rho}^F$ which may be related by comparing the mass balance equation for the whole solid phase:

$$ \frac{d\rho^S}{dt} + \rho^S \mathrm{div}\mathbf{v}_S = -\hat{\rho}^F. \qquad (15) $$

with eq. (14)$_2$. Hence, the erosion induced-part of the volume content variation is directly linked to the rate of mass supplied to the fluidized particles,

$$ dv^{Fi} = \frac{\hat{\rho}^F}{\rho_S}. \qquad (16) $$

## 2.4. Darcy's law

The flow in the porous medium is governed by Darcy's law which states that the flow rate is driven by the gradient of the hydraulic head and is proportional to the ratio of the intrinsic permeability $k_F$ over the fluid dynamic viscosity $\mu_F$:

$$ \mathbf{J}_F = -\frac{k_F}{\mu_F} [\nabla p_F - \rho_F \mathbf{g}]. \qquad (17) $$

Since the fluid phase is a mixture, the fluid density $\rho_F$ depends on the fluidized particles concentration by $\rho_F = C_{pF}\rho_{pF} + (1 - C_{pF})\rho_{wF}$. Yet in view of experimental observations and other numerical studies [1, 24, 32, 45], it is reasonable to consider that $C_{pF} \ll 1$, so that $\rho_F \approx \rho_{wF}$. For the same reason, the fluid
dynamic viscosity will be assumed independent of the fluidized particle concentration with $\mu_F \approx \mu_{wF}$.

With time the intrinsic permeability $k_F$ is considered to be dependent on the porosity $n^F$ via the Kozeny-Carman equation,

$$k_F = k_0 \frac{n^F}{(1-n^F)^2} \qquad (18)$$

where $k_0$ is the Kozeny-Carman coefficient [28].

## 2.5. Eroded mass relationships for suffusion

To describe the suffusive behaviour of a porous medium, the mass exchange term $\hat{\rho}^F$ needs to be formulated. Considerations based on the thermodynamics of irreversible processes can be found in several references [5, 18, 33, 47]. The committed position of this work is to favour relationships that use measurable material parameters and that take into account experimentally-based observations. In this respect, two non-linear relationships are presented: the power-based relationship and the energy-based relationship. These relationships are both based on three ideas: (i) the power dissipated by the flow is thought to drive the suffusion kinetics [30, 34], (ii) the history of the hydraulic loading also influences the suffusion kinetics [32] and (iii) the rate of eroded mass may only be positive or null. This latter concept arise from the hypothesis that the whole specimen is considered to be a Representative Elementary Volume (REV). On this REV and during the current experimental procedures, more particles are exiting the specimen than entering the specimen [22, 7, 32]. Inside this REV, suffusion comprises three mechanisms (detachment of some fine particles, transport and sometimes self-filtration) that cannot be distinguished from each other since they are taking place first successively and rapidly simultaneously. Once particles exit a REV, they are assumed to remain fluidized and the properties of the fluid which enters in a REV are assumed to be constant. To summarize, self-filtration (which is a kind of reversible detachment) is hence accounted for in each REV. By enforcing the rate of eroded mass to be positive or null, we mimic the conditions of the experimental tests. The validity of the above hypothesis, i.e. the whole specimen is considered to be a REV, would deserve a dedicated work [45, 43] that is not in the scope of this study.

### 2.5.1. The power-based relationship

This phenomenological model has been proposed by Sibille et al. [35] to describe the development of the suffusion process. This model relates the erosion rate per unit volume $\hat{\rho}^F$ to the instantaneous volumetric power dissipated by the flow $\bar{P}_{flow}$ by using two material parameters $\alpha_{ref}$ and $b_S$ which characterise the material's erodibility,

$$\hat{\rho}^F = \alpha_{ref} \left[ \frac{1}{\frac{\Delta \bar{E}_{stage}}{\bar{P}_{flow} t^*} + 1} \right] (\bar{P}_{flow})^{b_S}. \qquad (19)$$

$t^*$ is a characteristic time relative to the tested material. The history of the hydraulic loading (i.e. the amplitude and the duration of each loading stage) is accounted for by $\Delta \bar{E}_{stage}$, which is the volumetric flow energy cumulated from the initiation of each hydraulic loading stage characterised by $t_{init}$,

$$\Delta \bar{E}_{stage} = \int_{t_{init}}^{t>t_{init}} \bar{P}_{flow} dt. \qquad (20)$$

Throughout, the overhead bar will refer to a volumetric quantity with respect to the initial volume $V_0$. The volumetric flow power $\bar{P}_{flow}(t)$ is computed from the fluid flow and the flow head gradient [5],

$$\bar{P}_{flow}(t) = -\mathbf{J}_F \cdot [\nabla p_F - \rho_w \mathbf{g}] = \frac{k_F}{\mu_F} \|\nabla p_F - \rho_w \mathbf{g}\|^2. \qquad (21)$$
The parameters $\alpha_{ref}$ and $b_S$ are determined from the erosion rate at the initiation of each hydraulic loading stage [35]. This initiation corresponds to the largest erosion rate experimentally observed on each stage and characterises well suffusion when detachment prevails upon filtration. It is represented by the upper limit envelop of data approximated by the following power law [35],

$$ \hat{\rho}_{\text{upperlimit}}^{F} = \alpha_{\text{ref}} (\bar{P}_{\text{flow}})^{b_S}. \qquad (22) $$

It should be noted that the determination of $\alpha_{ref}$ and $b_S$ requires several eroded mass measurements over each loading stage. This metrology is challenging to achieve without closing and opening the downstream valve that can produce uncontrolled pressure surges. If only one measurement is done for each loading stage, $\alpha_{ref}$ and $b_S$ can be slightly underestimated, i.e. they will characterise averaged erosion rates rather than the largest ones.

### 2.5.2. The energy-based relationship

Experimentally, the development of suffusion can be characterized by the cumulative eroded mass, per unit volume, collected at the outlet of the specimen,

$$ \bar{m}_{\text{cum}}(t) = \int_{t_0}^{t} \bar{m}(t) dt. \qquad (23) $$

so that the erosion rate per unit volume $\hat{\rho}^F$ is the time derivative of the volumetric cumulative eroded mass $\bar{m}_{cum}$. Accordingly, the volumetric flow power may be cumulated to obtain the cumulative volumetric energy expanded by flow,

$$ \bar{E}_{\text{cum}}(t) = \int_{t_0}^{t} \bar{P}_{\text{flow}}(t) dt. \qquad (24) $$

For some suffusion tests, the suffusion development is observed to reach a stable state of erodibility characterized by a constant permeability and a decreasing erosion rate [21]. This stable state is characterised by the maximum cumulative eroded mass and the maximum cumulative expended energy,

$$ \bar{m}_{\max} = \int_{t_0}^{t_{\text{end}}} \bar{m}(t) dt, \quad \bar{E}_{\max} = \int_{t_0}^{t_{\text{end}}} \bar{P}_{\text{flow}}(t) dt. \qquad (25) $$

When this stable state is reached, experimental evidences have shown that the maximum cumulative eroded mass $\bar{m}_{max}$ and the maximum cumulative expended energy $\bar{E}_{max}$ were proportional in a repetitive manner [31]. This proportionality introduces the suffusion resistance index $I_\alpha$,

$$ \bar{m}_{\max} = 10^{-I_{\alpha}} \bar{E}_{\max}. \qquad (26) $$

The energy-based relationship [19] relates the cumulative eroded mass to the cumulative expended energy up to the stable state of suffusion by a non-linear power function,

$$ \frac{\bar{m}_{\text{cum}}(t) - \bar{m}_{\text{sat}}}{\bar{m}_{\max} - \bar{m}_{\text{sat}}} = \left( \frac{\bar{E}_{\text{cum}}(t)}{\bar{E}_{\max}} \right)^{b(t)} \qquad (27) $$

where $\bar{m}_{sat}$ is the mass lost during the saturation phase and $b(t)$ is a term that controls the kinetics of suffusion. Small values of $b(t)$ indicate a fast kinetics and conversely. Aside from $\bar{m}_{sat}$, which is initially measured, and from $b(t)$ which will be discussed below, the energy-based relationship uses two parameters: the suffusion resistance index $I_\alpha$ and the maximum cumulative expended energy $\bar{E}_{max}$.

Since the erosion rate is observed to be large at the beginning of each loading step and to decrease during the step, $b(t)$ is expected to vary with time. From eq. (27), the definition of $b(t)$ is constrained by two restrictions: (i) $b(t)$ should be a dimensionless variable and (ii) the cumulative eroded mass $\bar{m}_{cum}(t)$
should never decrease to remain physically admissible.

(i) In the spirit of the energy-based approach, the kinetic term $b(t)$ is defined as a function of the power dissipated by the flow $\bar{P}_{\text{flow}}$ and a smoothed value of this volumetric power determined based on a moving average method $\bar{P}_{\text{smoothed}}(t)$,

$$b(t) = \frac{\bar{P}_{\text{smoothed}}(t)}{\bar{P}_{\text{flow}}(t)}. \qquad (28)$$

The above equation is inspired from signal processing techniques and implies that each sudden increase in $\bar{P}_{\text{flow}}(t)$ is associated with a decrease in $b(t)$. The new smoothed series of values is obtained by averaging the original series over a number of instantaneous power values $n_b$ within a chosen smoothing time $t_{\text{smoothed}}$,

$$\bar{P}_{\text{smoothed}}(t) = [\bar{P}_{\text{flow}}(t) + \bar{P}_{\text{flow}}(t - \Delta t) + \bar{P}_{\text{flow}}(t - 2\Delta t) + \dots + \bar{P}_{\text{flow}}(t - n_b\Delta t)]/n_b. \qquad (29)$$

$\Delta t$ begin the time interval between two power values. Importantly, the knowledge of $\bar{P}_{\text{flow}}(t)$ values ahead of the time $t$ are not required to compute $\bar{P}_{\text{smoothed}}(t)$. In fact, the smoothing time $t_{\text{smoothed}}$ is the material parameter that controls the $b(t)$ term and hence the kinetics. A larger smoothing time indicates a faster kinetics and conversely.

(ii) To prevent the cumulative eroded mass from decreasing, the following rule is enforced to further constrain $b(t)$ between times $t_n$ and $t_{n+1}$: $\bar{m}_{\text{cum}}(t_{n+1}) \ge \bar{m}_{\text{cum}}(t_n)$.

In other words, the following rule is applied:

$$b(t_{n+1}) = \begin{cases} \frac{\bar{P}_{\text{smoothed}}(t)}{\bar{P}_{\text{flow}}(t)}, & \text{if } b(t_{n+1}) \ge b(t_n) \frac{\log\left(\frac{\bar{E}_{\text{cum}}(t_n)}{\bar{E}_{\text{max}}}\right)}{\log\left(\frac{\bar{E}_{\text{cum}}(t_{n+1})}{\bar{E}_{\text{max}}}\right)} \\ b(t_n) \frac{\log\left(\frac{\bar{E}_{\text{cum}}(t_n)}{\bar{E}_{\text{max}}}\right)}{\log\left(\frac{\bar{E}_{\text{cum}}(t_{n+1})}{\bar{E}_{\text{max}}}\right)}, & \text{otherwise.} \end{cases} \qquad (30)$$

## 2.6. Governing equations

The following governing equations were obtained by combining the balance and constitutive equations to describe the coupled suffusion poro-elastic process:

- Balance of momentum for the mixture:

$$\mu^{\text{DS}} \nabla^2 \mathbf{u} + (\lambda^{\text{DS}} + \mu^{\text{DS}}) \nabla(\text{div}\mathbf{u}) - \nabla p_F + \rho \mathbf{g} = 0 \qquad (31)$$

- Balance of mass for the fluid phase:

$$\text{div} \left( -\frac{k_F}{\mu_F} (\nabla p_F - \rho_F \mathbf{g}) \right) + \frac{d(\text{tr}\epsilon)}{dt} = 0 \qquad (32)$$

- Balance of mass for the fluidized particles:

$$n^F \frac{dC_{pF}}{dt} + C_{pF} \frac{d(\text{tr}\epsilon)}{dt} + \text{div} \left( -C_{pF} \frac{k_F}{\mu_F} (\nabla p_F - \rho_F \mathbf{g}) \right) = (1 - C_{pF}) \frac{\hat{\rho}^F}{\rho_S} \qquad (33)$$

- Balance of mass for the solid phase:

$$\frac{dn^F}{dt} = (1 - n^F) \frac{d(\text{tr}\epsilon)}{dt} + \frac{\hat{\rho}^F}{\rho_S} \qquad (34)$$
- Balance of mass for the erodible particles:

$$ \frac{dC_{pS}}{dt} = -\frac{1 - C_{pS}}{1 - n^F} \frac{\hat{\rho}^F}{\rho_S} \qquad (35) $$

to which should be added the Kozeny-Carman equation (18) and one of the two mass transfer relationships (19) or (27). Details of the derivations of eqs. (31) to (35) can be found in Kodieh [18]. In the above system, one may observe that the fluidized particles concentration $C_{pF}$ does not appear in the balance of momentum, nor in the balance of mass equations for the fluid phase, the solid phase and the erodible particles. This is due to the fact that $\hat{\rho}^F$ does not depend on $C_{pF}$, and that $\rho_F$ and $\mu_F$ are assumed independent of $C_{pF}$. In this context, the fluidized particles concentration is not needed to obtain the unknowns: $\mathbf{u}$, $p_F$, $n^F$ and $C_{pS}$; and the system to be solved reduces to four equations:

$$ \mu^{DS} \nabla^2 \mathbf{u} + (\lambda^{DS} + \mu^{DS}) \nabla(\mathrm{div} \mathbf{u}) - \nabla p_F + \rho \mathbf{g} = 0 \qquad (36) $$

$$ \mathrm{div} \left( -\frac{k_F}{\mu_F} (\nabla p_F - \rho_F \mathbf{g}) \right) + \frac{d(\mathrm{tr} \mathbf{\epsilon})}{dt} = 0 \qquad (37) $$

$$ \frac{dn^F}{dt} = (1 - n^F) \frac{d(\mathrm{tr} \mathbf{\epsilon})}{dt} + \frac{\hat{\rho}^F}{\rho_S} \qquad (38) $$

$$ \frac{dC_{pS}}{dt} = -\frac{1 - C_{pS}}{1 - n^F} \frac{\hat{\rho}^F}{\rho_S} \qquad (39) $$

The system of partial differential equations involves 6 scalar fields in 3D problems and 5 scalar fields in 2D problems. They are highly non linear and not easy to solve. Yet, no boundary conditions for the porosity $n^F$ and the erodible particles concentration $C_{pS}$ are required. For this reason, these variables, along with the permeability $k_F$, will be considered as state variables rather than primary unknowns. In addition, the differential equations describing the fluid porosity and the erodible particles concentration will be determined explicitly. Hence, the numerical integration is organised in two steps: (i) an implicit step to solve the displacement and the pore fluid pressure and (ii) an explicit step to update the porosity, the permeability and the erodible particles concentration. The numerical integration of the implicit step is a simplified version of the work detailed in Gelet et al. [13] which was verified against analytical models in Gelet [12]. The explicit step has been developed as part of this work in the FE code written in FORTRAN [18].

### 3. Suffusion tests on cohesionless specimens

Suffusion has been investigated in the past [32, 21, 48] on various soils, with specimens of several sizes, subjected to various hydraulic loading paths. However, the conclusions obtained at the scale of a specimen cannot directly be applied to a structure. As an intermediate step, suffusion tests were performed on a physical model of dike for one selected soil, named soil B [31]. For this reason, the capabilities of both constitutive models for suffusion, eqs. (19) and (27), are now investigated specifically for soil B. Other soils are considered in Kodieh et al. [19]. The specimen-scale is addressed herein while the scale of the physical model is postponed to Section 4.

This soil B, created by a mixture of sand and gravel from Sablière Palvadeau [31], is a cohesionless soil with a gap-graded distribution (Figure 1). Although contradictory, susceptibility criteria are still a matter of interest for dam engineers. According to Kenney and Lau's criterion [17, 16], this soil is potentially unstable (Table 3). As the gap ratio $G_r$ is less than 3, this soil is classified as internally stable by Chang and Zhang [7]. Finally, according to the energy approach classification [31], the specimens of soil B suffused up to a stable state of suffusion are identified as being moderately erodible, i.e. $3 < I_\alpha < 4$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Grain size distribution for the selected gap-graded soil B.</figcaption></figure>

<table><caption>Table 1: Properties of the tested gradation soil B (Figure 1). P = percentage of particles smaller than 0.063mm; $G_r$ = $d_{max}/d_{min}$ ($d_{max}$ and $d_{min}$: maximal and minimal particle sizes characterizing the gap in the grading curve); $C_u$ = uniformity coefficient; F and H are the mass percentages of the grains with a size lower than a given particle diameter d and between d and 4d, respectively; $D(H/F)_{min}$ is the corresponding diameter for the minimum value of ratio H/F.</caption><thead><tr><th>P (%)</th><th>G<sub>r</sub></th><th>C<sub>u</sub></th><th>(H/F)<sub>min</sub></th><th>D(H/F)<sub>min</sub> (mm)</th><th>Kenney and Lau's criterion</th><th>Chang and Zhang's criterion</th><th>Energy approach</th></tr></thead><tbody><tr><td>1.533</td><td>2.14</td><td>19.52</td><td>0.035</td><td>0.400</td><td>Unstable</td><td>Stable</td><td>Moderately erodible</td></tr></tbody></table>

### 3.1. *Experimental program and procedures*

For soil B, six specimens were prepared using a single layer semi-static compaction technique in a triaxial erodimeter device (used in oedometric conditions) and are compressed until the initial targeted dry unit weight is reached. All specimens are 50 mm diameter and length cylinders. Downstream of each specimen, a sieve is placed with a pore opening size selected according to the tested soil. The pore opening of the downstream sieve must be selected to avoid the migration of coarse grains, but also to permit the migration of all fine grains. To achieve that, fine and coarse fractions have to be distinguished. Skempton and Brogan [36] suggested that this distinction is best represented by the minimum value of the Kenney and Lau's ratio H/F (see Fig. 1). However, Marot et al. [23] showed that the use of a downstream filter with a pore opening equal to the maximum diameter of the fine fraction can restrict the suffusion development. It can be assumed that this limitation is due to a geometric filtration by arching effect, which increases with the number of detached particles reaching the filter at the same time. This effect depends on the hydraulic loading path and can lead to an overestimation of the soil resistance towards suffusion. In this context, the authors suggest to identify the smallest value of the pore sieve opening by the maximum value of the ratio H/F that follows the minimal value of H/F [40]. For the tested soil, Figure 1 shows that this method permits to select a pore opening of the downstream sieve greater or equal to 1.2 mm. The pore opening of the downstream sieve used for these tests was equal to 4 mm.

Next, the specimens were saturated by injection of carbon dioxide and subsequently demineralized water, at a low rate. The dry unit weight of each specimen *after* the saturation phase is indicated in Table 2. Denser specimens are referred to as B-97 (97 corresponds to the percentage of the dry density in comparison with the Proctor's optimum), while others are named B-90. Afterwards, each specimen is subjected to a downward flow driven by a specific hydraulic loading path illustrated in Figure 2. Each hydraulic loading path is named
by a letter. These letters are then attached to the name of the specimens. For instance, specimens B-90a and B-97a are both subjected to the hydraulic loading 'a'. It is worth highlighting that the hydraulic gradients are either single-staged, such as loading 'c', or multi-staged with different stage durations or ranges. The hydraulic loading 'a' consists of increasing the hydraulic gradient by steps of 0.1 until 2, then by steps of 0.5 between 2 and 4, and by steps of 1 beyond. Steps are directly equal to 1 for the hydraulic loading 'b'. For both loadings 'a' and 'b', stage durations were equal to 10 minutes. The hydraulic loading 'h' consists of increasing the hydraulic gradient by steps of 0.2 until 1, then by steps of 1 beyond, for 60 minutes each; whereas the hydraulic loading 'k' only consists of hydraulic gradients 0.5 and 1 which are kept constant during 12 hours.

<table><caption>Table 2: Properties of tested specimens and summary of testing program.</caption><thead><tr><th>Specimen</th><th>Initial dry unit weight (kN/m³)</th><th>Applied hydraulic gradient i (m/m)</th><th>Initial hydraulic conductivity 10<sup>-5</sup>(m/s)</th><th>Test duration (min)</th></tr></thead><tbody><tr><td>B-90a</td><td>17.38</td><td>From 0.10 to 6.00</td><td>0.97</td><td>180</td></tr><tr><td>B-90c</td><td>17.32</td><td>4</td><td>2.82</td><td>300</td></tr><tr><td>B-90h</td><td>17.15</td><td>From 0.20 to 10.0</td><td>6.04</td><td>720</td></tr><tr><td>B-90k</td><td>17.14</td><td>From 0.50 to 1.00</td><td>5.19</td><td>1440</td></tr><tr><td>B-97a</td><td>18.74</td><td>From 0.10 to 12.0</td><td>1.39</td><td>240</td></tr><tr><td>B-97b</td><td>18.74</td><td>From 1.00 to 9.00</td><td>3.16</td><td>90</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Multi-stage and single-stage hydraulic gradients applied to the specimens of soil B.</figcaption></figure>

The suffusion behaviour of each specimen is best represented by the hydraulic conductivity (Figure 3, left), computed from Darcy's law, and the erosion rate (Figure 3, right), which is the mean erosion rate over each hydraulic stage. According to Figure 3, some tests display a decrease in the hydraulic conductivity either at the beginning of the test (B-90a, B-97a B-97b and B-90k) or during the whole test (B-90c). This decrease is attributed to some fine particles detached and transported under the imposed water seepage and then filtered within the soil itself. This filtration induces a local physical clogging thus decreasing the hydraulic conductivity of the soil. From the erosion rate point of view, this local clogging aspect can be observed as a decrease in the erosion rate despite the increase in hydraulic gradients. The subsequent or late
<figure><img src="image_4.png" /><figcaption>Figure 3: Variation with time of the hydraulic conductivity (left) and the erosion rate (right), for six specimens of soil B [31].</figcaption></figure>

increase in the hydraulic conductivity and erosion rate (B-90a, B-97a, B-97b and B-90h) is attributed to the washing out of these fine particles caused by larger hydraulic gradients, overcoming the clogging aspect. Especially for some specimens (B-90a, B-97a and B-97b), the hydraulic conductivity and the erosion rate both increase to a limit level, where the former reaches a steady state while the latter starts to decrease. By following the energy approach [26], these conditions characterize the stable state of erodibility.

The cumulative eroded mass per unit volume of each specimen is displayed in Figure 4. Importantly, the eroded mass is currently collected separately for each hydraulic stage, so that its acquisition frequency is much lower than that used to compute the hydraulic conductivity. Additional information relative to the experimental device and procedure as well as the testing program and results can be found in Rochim [31] and Rochim et al. [32].

## 3.2. Calibration procedure

In Section 2.5, two constitutive laws for suffusion were presented. Both of them are built upon the experimental observation that the power dissipated by the flow is controlling the evolution of suffusion, up to the end of each loading stage for the power-based approach and up to a stable state of erodibility for the
<figure><img src="image_5.png" /><figcaption>Figure 4: Experimental measurements of the cumulative eroded mass per unit volume, for six specimens of soil B [31].</figcaption></figure>

energy-based approach. As stated beforehand, this stable state of erodibility is characterized by a constant permeability and a decreasing erosion rate [26]. However, this stable state of erodibility may not be reached for each test, e.g. B-90c and B-90k. Recall that the energy-based relation uses three material parameters: the suffusion resistance index $I_α$, the maximum volumetric cumulative expended energy $\bar{E}_{max}$ and the smoothing time $t_{smoothed}$. Since $I_α$ and $\bar{E}_{max}$ characterise the stable state of erodibility, the calibration procedure should use a specimen that reaches the stable state of erodibility. Among all specimens, three fulfil this condition (B-90a, B-97a and B-97b) upon which B-90a is arbitrarily chosen to calibrate the material parameters.

Next, the calibration procedure should account, if possible, for the effect of the dry unit weight on the material parameters. Yet, the available experimental results [21, 32] are not numerous enough to safely infer a possible correlation between the maximum volumetric cumulative eroded mass $\bar{m}_{max}$ and the dry unit weight. In addition, for the initial amount of fines at hand, i.e. 25%, fine particles-coarse particles contacts are created and contribute to carrying the stress, yet contacts between coarse particles primarily carry the stress so that the fine particles are slightly stressed [39]. Consequently, the maximum cumulative eroded mass per unit volume is assumed to be constant for the range of $\gamma_d$ at hand:

$$ \bar{m}_{max} B = \bar{m}_{max} B-90a. \qquad (40) $$

Loose specimens are more susceptible to suffusion than denser ones [32, 21]. Adapting the work of Le et al. [20, 21] to soil B, a linear correlation is assumed between the suffusion resistance index $I_α$ and the dry unit weight $\gamma_d$:

$$ I_{\alpha} B = \frac{\gamma_d}{\gamma_d B-90a} \times I_{\alpha} B-90a. \qquad (41) $$

The maximum volumetric cumulative expended energy for soil B is then expressed by replacing eqs. (40) and (41) into the stable state relationship (26),

$$ \bar{E}_{max} B = 10 \left[ I_{\alpha} B-90a \left( \frac{\gamma_d}{\gamma_d B-90a} - 1 \right) \right] \bar{E}_{max} B-90a. \qquad (42) $$

Finally, based on the available experimental results, no obvious correlation can yet be made between the development of suffusion and the dry unit weight (Figures 4 and 5). Hence, the smoothing time is also assumed independent of the dry unit weight,
<figure><img src="image_6.png" /><figcaption>Figure 5: Variation of the cumulative eroded mass per unit volume with the cumulative expended energy per unit volume, both measured by Rochim [31].</figcaption></figure>

$$t_{\text{smoothed B}} = t_{\text{smoothed B-90a}}. \qquad (43)$$

Regarding the power-based approach, the determination procedure for $\alpha_{\text{ref}}$ and $b_S$ was initially proposed by Sibille et al. [35]. Again, the calibration procedure should account, when possible, for the effect of the dry unit weight on the material parameters introduced in the power-based approach: $\alpha_{\text{ref}}$ and $b_S$. Yet, due to the lack of experimental data, the determination procedure initially proposed by Sibille et al. [35] was adopted. This procedure equally considers all the specimens of one grain size distribution at once, and material parameters are assumed to be constant in the range of variation of the tested dry density. This procedure relies on the maximum erosion rate measured for each hydraulic loading stage and uses the upper limit envelop of the data approximated by a power law (Figure 6),

$$\hat{\rho}_{\text{upperlimit}}^{F} = \alpha_{\text{ref}} (\bar{P}_{\text{flow}})^{b_S}. \qquad (44)$$

Here, the parameters $\alpha_{\text{ref}}$ and $b_S$ are measured by using the data of all tests realised on soil B. The potential influence of the dry unit weight on these parameters is not accounted for, in line with the work of Sibille et al. [35]. The current experimental apparatus gives access to one erosion rate per hydraulic stage, unlike the data initially used [35] for which several erosion rates were available for each hydraulic stage. This difference may affect the precision of $\alpha_{\text{ref}}$ and $b_S$, especially $\alpha_{\text{ref}}$ which might be underestimated.

Next, the characteristic time $t^*$ for the reference specimen B-90a is best fitted. All other parameters being constant, the increase in $t^*$ is found to increase the amount of eroded mass which suggests that $t^*$ should be inversely correlated with the dry unit weight $\gamma_d$. Following the calibration approach developed for $I_\alpha$, the characteristic time $t^*$ of soil B is calibrated with,

$$t_B^* = \frac{\gamma_d B - 90a}{\gamma_d} \times t_{B-90a}^*. \qquad (45)$$

Based on the above calibration procedure, i.e. eqs. (40) to (45), the material parameters $I_\alpha$, $\bar{E}_{\text{max}}$, $t_{\text{smoothed}}$, $\alpha_{\text{ref}}$, $b_S$ and $t^*$ are calibrated for each specimen (Table 3) by taking the parameters of specimen B-90a as reference. To summarise, three parameters have the same value for all specimens ($t_{\text{smoothed}}$, $\alpha_{\text{ref}}$ and $b_S$) while the other three vary slightly with the initial dry unit weight of each specimen ($I_\alpha$, $\bar{E}_{\text{max}}$ and $t^*$). The parameters obtained from the calibration procedure and displayed in Table 3, are the only one being used to obtain the numerical predictions presented in Section 3.3. In other words, these material parameters are not being changed with the applied hydraulic loading path.
<figure><img src="image_7.png" /><figcaption>Figure 6: Identification of the maximum erosion rate per unit volume as a function of the flow power per unit volume for the specimens of soil B.</figcaption></figure>

<table><caption>Table 3: Suffusion parameters for each specimen of soil B upon calibration. The reference specimen and values used for the calibration are highlighted in bold.</caption><thead><tr><th>Parameters</th><th>Units</th><th>B-90k</th><th>B-90h</th><th>B-90c</th><th><b>B-90a</b></th><th>B-97a</th><th>B-97b</th></tr></thead><tbody><tr><td>γ<sub>d</sub></td><td>kN/m<sup>3</sup></td><td>17.14</td><td>17.15</td><td>17.32</td><td><b>17.38</b></td><td>18.74</td><td>18.74</td></tr><tr><td>I<sub>α</sub></td><td>SI units</td><td>3.18</td><td>3.18</td><td>3.21</td><td><b>3.22</b></td><td>3.47</td><td>3.47</td></tr><tr><td>E&#x0302;<sub>max</sub></td><td>10<sup>4</sup>(J/m<sup>3</sup>)</td><td>6.45</td><td>6.48</td><td>6.97</td><td><b>7.15</b></td><td>12.78</td><td>12.78</td></tr><tr><td>t<sub>smoothed</sub></td><td>(s)</td><td>60</td><td>60</td><td>60</td><td><b>60</b></td><td>60</td><td>60</td></tr><tr><td>α<sub>ref</sub></td><td>10<sup>-3</sup> SI units</td><td>7.8</td><td>7.8</td><td>7.8</td><td><b>7.8</b></td><td>7.8</td><td>7.8</td></tr><tr><td>b<sub>S</sub></td><td>10<sup>-2</sup>(-)</td><td>44.46</td><td>44.46</td><td>44.46</td><td><b>44.46</b></td><td>44.46</td><td>44.46</td></tr><tr><td>t*</td><td>(s)</td><td>121.6</td><td>121.6</td><td>120.4</td><td><b>120.0</b></td><td>111.3</td><td>111.3</td></tr></tbody></table>

### 3.3. Experimental and numerical results

The aim of this section is to illustrate the abilities of both models to qualitatively and quantitatively reproduce the cumulative eroded mass of a suffusive soil during a suffusion test, at the scale of a material point. Numerical integrations of both constitutive models, eqs. (19) and (27), are compared with the collected mass measured during laboratory tests [31]. It is worth highlighting that the measured flow rate and applied hydraulic gradient (Figure 2) have been directly used to compute the power dissipated by the flow and the hydraulic conductivity.

The correspondence between the cumulative eroded mass predicted using the power-based relationship and the experimental results is qualitatively good for specimens with a multi-stage hydraulic loading: B-90a and B-97b (Figure 7). This can be explained by the fact that these two specimens are controlling the upper limit envelop (Figure 6). Moreover, this correspondence is fairly good for the specimens with a single-stage or two-stage hydraulic loading: B-90c and B-90k (Figures 8 and 9). The power-based approach is able to reproduce a step-like evolution of the cumulative eroded mass for some tests (B-90c and B-90k) which is attributed to the parameter $\Delta\overline{E}_{stage}$ that integrates the hydraulic loading history over each loading stage. On the other hand, this approach may severely overestimates the maximum cumulative eroded mass towards the end of some tests (B-97a and B-90h in Figures 8 and 9), since no cut-off value exists to limit
the maximum cumulative eroded mass. When considering the specimen B-90h alone in Figure 6, the upper limit envelop should have a decreasing slope. Hence, the spirit of the power-based approach does not allow a good representation of this test.

<figure><img src="image_8.png" /><figcaption>Figure 7: Comparison of the cumulative eroded mass predicted using the energy-based relationship and the power-based relationship with the collected mass measured for tests B-90a and B-97b.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison of the cumulative eroded mass predicted using the energy-based relationship and the power-based relationship with the collected mass measured for tests B-97a and B-90k.</figcaption></figure>

For the energy-based relationship, the predicted cumulative eroded mass captures fairly well the step-like evolution of specimens B-90a, B-90k and B-97a and a delayed response for B-97b. Nevertheless, the behaviour of specimens B-90c and B-90h could not be well captured. Those discrepancies may be attributed to two issues. (i) For the specimen B-90h, the parameter $\bar{E}_{max}$ that partly characterize the stable state of erodibility is calibrated to $6.48 \times 10^4$ J/m$^3$, as detailed in Section 3.2. Yet, upon this calibration, this parameter might not be characterized with enough precision (Figure 5), thus causing the mentioned discrepancies.
<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison of the cumulative eroded mass predicted using the energy-based relationship and the power-based relationship with the collected mass measured for tests B-90c and B-90h.</figcaption></figure>

Hence, the dependency of the maximum cumulative eroded mass per unit volume $\bar{m}_{max}$ with the dry unit weight should be further studied. (ii) The smoothing time that partly controls the kinetics of our prediction is assumed equal for all tests. All other parameters being constant, the decrease in the smoothing time slows the kinetics down and conversely. In fact, this smoothing time might be related to the characteristic time of hydraulic diffusion. In view of Figure 3, one can see that the permeabilities of specimens B-90c, B-90k, and B-90h are larger in average than that of specimen B-90a; so that, the corresponding average characteristic time of diffusion of these specimens is smaller compared with B-90a. This might suggest that the smoothing time should also be smaller for these specimens (B-90c, B-90h and B-90k) than for B-90a. Hence, further investigations should be pursued to study the effect of the smoothing time on the prediction of the cumulative eroded mass. Importantly, the predicted cumulative eroded mass, using the energy-based relationship, for these specimens does not grow infinitely, due to the presence of a cut-off value, as presented in eq. (40).

To enhance the comparison between the energy-based and the power-based capabilities, the absolute errors between the predicted and the measured values are calculated for (i) the mass, at each measured instant, and (ii) for the time, at each measured mass. In fact, these absolute errors measure the vertical and the horizontal shifting of the prediction. The average of these absolute errors is next determined separately for mass and time, for each prediction (Table 4). Both values are useful to categorise the quality of the prediction due to the small amount of experimental measures. Low average errors on the mass (< 3 kg/m³) and on the time (< 40 min) indicate a relatively good accuracy. A low average error on the mass only corresponds to a good prediction of the eroded mass range with an early or a delayed kinetics. A low average error on the time only reveals a good prediction of the overall kinetics with an overestimation or an underestimation on the eroded mass. Finally, large average errors on both mass and time point towards less accurate predictions. Overall, both approaches are able to predict quite well the cumulative eroded mass for various hydraulic loading paths. Based on the aforementioned simulations, the power based law seems to be the most relevant relation to model suffusion induced by a large constant hydraulic gradient that mainly involves a decrease in the hydraulic conductivity (B90-c). On the other hand, the energy based law seems to be the most relevant relation to model suffusion induced by multi-step hydraulic gradient paths (B97-a). This latter asset is due to the ability of the energy based law to cap the total cumulative eroded mass by using the cut-off value $\bar{m}_{max}$. Additional simulations should be realized to strengthen these conclusions.
<table><caption>Table 4: Comparison of the average absolute errors (AE) of mass and time for each prediction.</caption><thead><tr><th>Specimen</th><th>Relationship</th><th>Average AE of mass (kg/m³)</th><th>Average AE of time (min)</th><th>Comments</th></tr></thead><tbody><tr><td>B-90a</td><td>Energy-based</td><td>1.7</td><td>17</td><td rowspan="3">Good accuracy on the mass and the kinetics</td></tr><tr><td>B-97a</td><td>Energy-based</td><td>1.6</td><td>15</td></tr><tr><td>B-90a</td><td>Power-based</td><td>2.8</td><td>32</td></tr><tr><td>B-90k</td><td>Energy-based</td><td>2.9</td><td>661</td><td rowspan="3">Good eroded mass range with an early or delayed kinetics</td></tr><tr><td>B-90c</td><td>Power-based</td><td>1.9</td><td>240</td></tr><tr><td>B-90k</td><td>Power-based</td><td>1.8</td><td>240</td></tr><tr><td>B-97b</td><td>Energy-based</td><td>11.6</td><td>20</td><td rowspan="2">Good kinetics with over or underestimation of the eroded mass</td></tr><tr><td>B-97b</td><td>Power-based</td><td>5.0</td><td>7</td></tr><tr><td>B-90c</td><td>Energy-based</td><td>10.1</td><td>120</td><td rowspan="4">Low accuracy</td></tr><tr><td>B-90h</td><td>Energy-based</td><td>5.0</td><td>115</td></tr><tr><td>B-90h</td><td>Power-based</td><td>9.3</td><td>136</td></tr><tr><td>B-97a</td><td>Power-based</td><td>10.6</td><td>53</td></tr></tbody></table>

## 4. Suffusion test and simulation of a physical model of dike

Embankment structures are usually characterized by a simple geometry, in particular when focusing on the cross section. Hence when modelling suffusion in such structures, the main challenge is to achieve a good representation of the physical process involved. Our long term objective is to validate the ability of our numerical model to represent the suffusive behaviour of real scale earth structures [46]. As an intermediate step, our current objective is to reproduce numerically an experiment realized on a physical model of dike [31] by using the poro-elastic model extended to suffusion developed in Section 2 and both relationships previously studied.

### 4.1. *Experimental equipment and procedures*

The testing apparatus of the half-dike physical model consists of a transparent box, a water supply system, a soil collection system, a water collection system and an acquisition system (Figure 10). A series of pressure ports is installed on the back side of the box to measure local hydraulic heads. The box has inlet and outlet holes on both sides: a water supply system to the right and a funnel-shaped drainage system to the left. On the vertical sidewall of the water supply system, a fiber mesh is installed which allows the water to seep into the dike, yet avoids the fine fraction of the soil from migrating into the upstream retaining reservoir.

The dimensions of the dike are 100 cm in length, 45 cm in height and 15 cm in width (Figure 10). The crest has a length of 22 cm and the downstream slope is inclined with a 30° angle. The seepage test is divided into two phases: (i) the saturation phase at which the dike box is filled at a very small rate by injecting water from both the upstream and the downstream side of the soil specimen to saturate it over 12 hours, (ii) the seepage phase that starts by the de-watering of the downstream side. The water seeping through the body of the dike eventually passes the funnel-shaped drainage outlet, where the discharged eroded particles and water flow are delivered into an effluent tank. The eroded particles are being collected in 8 beakers and the water flow is continuously weighed by a balance to determine the water flow rate. The bottom side of the box is equipped with 7 pressure ports connected to piezometers. Two digital cameras are placed at the front and at the back of the box to estimate the displacement field of the dike and local water heads from image processing. Additional details concerning the dike testing apparatus can be found in Rochim [31].
<figure><img src="image_11.png" /><figcaption>Figure 10: Schematic diagram of the dike apparatus and dimensions (in cm) of the half-dike physical model [31].</figcaption></figure>

## 4.2. Finite element mesh

The finite element mesh (Figure 11) is characterized with 136 elements that have 5 cm lengths and heights in average, i.e. alike the seepage length experimentally used to calibrate the material parameters. Hence, the “model length / mesh size” ratio lies between 9 and 20 which is similar to that used by Zhang et al. [47] and greater in average than that used in Yang et al. [44]. The relatively coarse characteristic of this mesh is currently constrained by the probable dependency of some material parameters on the seepage length. In fact, previous studies [18, 19] have indicated that the maximum cumulative expended energy per unit volume $E'_{max}$ tends to increase with the seepage length. This is due to the experimental measurement of the whole suffusion process which involves three mechanisms: detachment, transport and self-filtration. While the volumetric energy used for the detachment should not be influenced by the flow length, it is expected that more volumetric energy will be required for the transportation of fine particles through the solid skeleton when the seepage length increases. In the future, further studies will be carried out to improve our physical understanding of the dependency of $E'_{max}$ upon the seepage length. To date, the relation linking the maximum cumulative expended energy per unit volume $E'_{max}$ to the seepage length remains unknown. To obtain this relation, additional experimental tests are required for different seepage lengths. By implementing this relation in the FEM code, the current constraint on the mesh size would be alleviated. This particular point will be addressed in future works.

In fact, the asset of our work lies in the careful calibration of the material parameters and the full awareness of the current limitations of our knowledge regarding the potential dependency of some material parameters on the seepage length. For example, the critical hydraulic gradient was found to be decreasing with the seepage length [48, 25]; and both the coefficient of erosion and the critical shear stress were found to be characterized by a noticeable too large dispersion of their relative error with the flow length [19].
<figure><img src="image_12.png" /><figcaption>Figure 11: Finite element mesh and boundary conditions used for the numerical simulations of the half-dike physical model.</figcaption></figure>

### 4.3. Initial conditions

Numerically speaking, the initial time corresponds to the end of the saturation process. Upon saturation, a single water level is imposed to the dike that corresponds to a total head of 43 cm. Accordingly, the initial pore pressures are supposed to be hydrostatic everywhere. Above the water level, the soil is assumed to be saturated by capillarity and pore pressures are characterized with negative pore water pressures. Although, the displacements and strains are identically null everywhere, initial total stresses are assigned to geostatic stresses.

### 4.4. Boundary conditions

The seepage test is conducted by maintaining the upstream water level constant and decreasing the downstream water level to the base level. The hydraulic (H) and mechanical (M) boundary conditions (BC) are explained in turn for each dike boundary starting by the bottom side and following a counterclockwise order (Figure 11).

(H) Along the bottom side, a no flow BC is imposed. Along the upstream side, the pore pressures corresponding to the experimentally applied total head are imposed (43 cm up to 120 minutes and 44 cm after). Along the crest of the dike, initial negative pore pressures are maintained corresponding to hydrostatic BC. Finally, along the slope of the dike, transient pore pressures are imposed during the first 8 minutes. For each node, the pore pressure corresponds to the experimentally applied water height. Upon emptying the downstream side, null pore pressures are reached along the slope.

(M) Both horizontal and vertical displacements are fixed along the bottom side. Along the upstream side, horizontal displacements are fixed and zero vertical stresses are imposed. Along the crest of the dike, zero total stresses are imposed. Finally, along the slope of the dike, transient total stresses are imposed during the first 8 minutes. For each node, the compressive total vertical stress corresponds to the experimentally applied water height. Again, upon de-watering the downstream side, null total stresses are reached all along the slope.

### 4.5. Material parameters

Table 5 summarizes the properties of the tested half-dike (B-90s-rep [31]) which was carefully compacted in five layers to the dry unit weight of 17.39 kN/m³. The suffusion process is modelled by using both the energy-based and the power-based approach. The related parameters are determined using the calibration technique detailed in Section 3. The intrinsic permeability $k_F$ was obtained by fitting the downstream volumetric flow rate (Figure 12). Finally, due to a lack of data, the volumetric amount of eroded particles detached during the saturation phase $\bar{m}_{sat}$ is assumed equal to 0.

### 4.6. Simulations versus experimental results

Throughout, numerical results are systematically compared against experimental data. In addition, numerical results are provided for both the energy-based approach and the power-based approach.
<table><caption>Table 5: Material parameters for the physical model of dike. † Fitted from Figure 12. ★ Calibrated from B-90a.</caption><thead><tr><th>Material parameter</th><th>Value</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>Drained Young's modulus <i>E</i><sup>DS</sup></td><td>7.777</td><td>MPa</td><td>[31]</td></tr><tr><td>Drained Poisson's ratio <i>ν</i></td><td>0.3</td><td>-</td><td>[31]</td></tr><tr><td>Dry unit weight <i>γ</i><sub>d</sub></td><td>17.39</td><td>kN/m<sup>3</sup></td><td>[31]</td></tr><tr><td>Porosity <i>n</i><sup>F</sup></td><td>0.33</td><td>-</td><td>[31]</td></tr><tr><td>Solid density <i>ρ</i><sub>s</sub></td><td>2650</td><td>kg/m<sup>3</sup></td><td>-</td></tr><tr><td>Initial concentration of erodible particles in the solid phase <i>C</i><sub>pS</sub>(<i>t</i> = 0)</td><td>0.25</td><td>-</td><td>Figure 1</td></tr><tr><td>Permeability <i>k</i><sub>F</sub></td><td>2.22 × 10<sup>-10</sup></td><td>m<sup>2</sup></td><td>†</td></tr><tr><td>Fluid dynamic viscosity <i>μ</i><sub>F</sub></td><td>1.0005 × 10<sup>-3</sup></td><td>Pa.s</td><td>-</td></tr><tr><td>Fluid density <i>ρ</i><sub>F</sub></td><td>998.2</td><td>kg/m<sup>3</sup></td><td>-</td></tr><tr><td>Suffusion resistance index <i>I</i><sub>α</sub></td><td>3.23</td><td>SI units</td><td>★</td></tr><tr><td>Maximum expended energy per unit volume <i>E</i><sub>max</sub></td><td>7.204 × 10<sup>4</sup></td><td>J/m<sup>3</sup></td><td>★</td></tr><tr><td>Smoothing time <i>t</i><sub>smoothed</sub></td><td>60</td><td>s</td><td>★</td></tr><tr><td>Power-based relationship parameter <i>α</i><sub>ref</sub></td><td>7.8 × 10<sup>-3</sup></td><td>SI units</td><td>★</td></tr><tr><td>Power-based relationship parameter <i>b</i><sub>S</sub></td><td>44.46 × 10<sup>-2</sup></td><td>-</td><td>★</td></tr><tr><td>Characteristic time <i>t</i><sup>*</sup></td><td>120</td><td>s</td><td>★</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 12: The intrinsic permeability $k_F = 2.22 \times 10^{-10}$ m² was obtained by fitting the predicted downstream volumetric flow rates with the experimental one.</figcaption></figure>

Figure 13 shows the variation of the cumulative eroded mass with time for the whole physical model of dike. Major fines erosion takes place in the early stage of the seepage test (until around 20 minutes). In this period, the discharge flow rate is also relatively large. After that, the erosion rate gets smaller in agreement with the discharge rate which gradually decreases with time to 1-1.5 l/min. Both numerical methods cannot reproduce the experimental kinetics. On one hand, the energy-based approach is able to predict the total amount of cumulative eroded mass, while the power based approach underestimates this total amount. On the other hand, the power-based approach is able to predict a two-step evolution, which is similar to the
real experimental evolution. This ability can be attributed to the $\Delta\bar{E}_{stage}$ parameter that takes into account the hydraulic loading history during each hydraulic step, see eq. (20). It should be noted that the great increase in the eroded mass observed at the beginning of the test might result from the sliding and the surface erosion on the slope during the de-watering at the downstream side.

<figure><img src="image_14.png" /><figcaption>Figure 13: Variation of the cumulative eroded mass with time: comparison between the experimental test [31] (symbols) and simulated data (discontinuous lines) using the energy-based relationship and the power-based relationship.</figcaption></figure>

Figure 14 shows the spatial distribution of the variation of the percentage of fines inside the dike after suffusion. This variation is defined with respect to the initial value of 25%. A positive variation indicates a decrease in fine content and conversely. Overall variations of the percentage of fines predicted by both approaches are in good agreement with that measured experimentally: the lowest variations are close to the upstream side, the strongest variations are located along the downstream slope and intermediate values are obtained in the core of the dike. While the experiment displays a nearly-smooth contour with some particular points, the numerical predictions display a smoother version of this contour. In terms of magnitude, the experimental variations are usually stronger than that obtained with both approaches. As expected from Figure 13, the variations predicted by the power-based approach tend to be smaller in magnitude than the ones obtained with the energy-based approach. The main differences between the numerical predictions and the experimental data can be discussed. (i) No negative variations are predicted along the upstream side since the rate of erosion $\dot{\rho}^F$ cannot become negative, i.e. both approaches assume that overall detachment prevails upon self-filtration. (ii) The small variation of the percentage of fines at point X9 is not predicted by both numerical approaches (Figure 14). This may be attributed to a local filtration induced by the loss of fines upstream of this point. (iii) Along the downstream slope, a large erosion takes place that is not reproduced by both approaches. This may be attributed to the sliding observed at the beginning of the experimental test (described below), which probably induced a loose zone that localised the flow and favoured larger variations of the percentage of fines. (iv) The time evolution of the permeability is obtained by using the semi-empirical relationship of Kozeny-Carman (18) which might not be accurate for soils with an evolving fine content [14]. In particular, the effect of the evolution of the constriction size distribution on the permeability is not accounted for [29].

Figure 15 presents the comparison between the experimental data (left) and numerical results obtained with the energy-based approach (right) for the spatial distribution of the total water head, at three instants. Recall that the water level was initially equal on both sides of the dike to 0.43 m. Next, the downstream level was lowered within the first 8 minutes while the upstream level was controlled. Generally speaking, the total water head is influenced by both the hydraulic boundary conditions and the permeability variation.
<figure><img src="image_15.png" /><figcaption>Figure 14: Spatial distribution of the variation of the percentage of fines after suffusion, with respect to the initial value of 25 %: (top) Experimental test [31], (middle) simulation with the energy-based relationship and (bottom) simulation with the power-based relationship. [To view this figure in color, please visit the online version of this Journal.]</figcaption></figure>

Experimentally, most of the total water head evolution took place at the very beginning of the test (before $t = 17$ minutes) which indicates that the permeability variation after this time is reduced. The overall tendency is well reproduced by the FEM predictions, where both approaches give the same results. As expected for this homogeneous dike, a smoother version of the total head spatial distribution is numerically obtained in comparison with the experimental data. Differences could be attributed to the lack of experimental data. Indeed, the experimental contours of Figure 15 have been obtained from data measured at ten pressure ports.

Figure 16 presents the spatial distributions of the displacement at 2.25 mins and 6.75 mins. Note that, both the energy-based and the power-based simulations predict the same results so that only the energy-based predictions are presented. The comparison between results simulated with a poro-elastoplastic model that does not account for suffusion (Plaxis) and our poro-elastic model that accounts for suffusion informs on the predominant mechanism that triggers the experimentally observed deformations. Since the Plaxis predictions are qualitatively better than those computed with our model, one may conclude that the experimentally observed displacements can be primarily identified as plastic induced deformations rather
<figure><img src="image_16.png" /></figure>

Distribution of total water head [m]
at time = 17 min (Experimental)

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

Distribution of total water head [m]
at time = 360 min (Experimental)

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

Distribution of total water head [m]
at time = 3.75 min (FEM)

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

Distribution of total water head [m]
at time = 17 min (FEM)

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

Distribution of total water head [m]
at time = 360 min (FEM)

<figure><img src="image_21.png" /><figcaption>Figure 15: Spatial distribution of total water head at different times: (left) experimental data [31] and (right) simulated data using the energy-based approach. [To view this figure in color, please visit the online version of this Journal.]</figcaption></figure>

than suffusion induced elastic deformations. Experimentally, a progressive sliding along the downstream slope is observed within the first seven minutes. The magnitude of this displacement is quite large (up to 4 centimeters) and cannot be easily reproduced by a conventional poro-elastoplastic model which assumes small strains and small deformations. According to Rochim [31], this sliding was attributed to an over-interstitial pressure generated by the water seepage and/or a loss of fine fraction near the downstream.

With Plaxis, the predicted displacements tend to decrease with time. In fact, the early induced displacements are mainly controlled by the mechanical dis-equilibrium that is triggered by the decrease in the applied mechanical pressure along the downstream slope. This dis-equilibrium tends to induce a large displacement which is next mitigated by the decrease in the pore water pressure induced by the flow (see Figure 15).

As expected, during this sliding, the magnitude of the displacement predicted by our poro-elastic model is also much lower compared with the experimental data. By focusing on several points inside the dike
<figure><img src="image_22.png" /><figcaption>Figure 16: Spatial distribution of the displacement, at instants 2.25 mins (left) and 6.75 mins (right). (top) Experimental data [31], (middle) poro-plastic simulation with no erosion [31] and (bottom) simulation data with the energy-based approach. [To view this figure in color, please visit the online version of this Journal.]</figcaption></figure>

(located in Figure 15), our poro-elastic model predicts larger vertical downward displacements along the downstream slope, points B3 and E4, in comparison with those located near the upstream side, points G5, G4, D3 and F3 (Figure 17). In addition, opposite horizontal displacements are predicted along the slope of the dike, see B3 and E4, which mimic a sliding process (Figure 18). In all cases, vertical displacements are approximately ten times larger than horizontal ones. Hence, the poro-elastic simulations suggest an overall displacement towards the toe of the dike, similar to the experimental response to within a lower magnitude.

The effects of an orthotropic permeability tensor and of a loose zone were discussed in Kodieh [18] (not shown here). By comparing the numerical data with the experimental ones, better results were found for an isotropic permeability tensor. In addition, the numerical results displayed an improved behaviour with respect to the kinetics of the cumulative eroded mass when artificially enforcing a loose zone, designed to mimic the effect of the experimentally observed early sliding.
<figure><img src="image_23.png" /><figcaption>Figure 17: History of the vertical displacement, for points B3, D3 and F3 (left) ; and E4, G4 and G5 (right). [To view this figure in color, please visit the online version of this Journal.]</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 18: History of the horizontal displacement, for points B3, D3 and F3 (left) ; and E4, G4 and G5 (right). Negative horizontal displacements are oriented towards the downstream slope of the dike. [To view this figure in color, please visit the online version of this Journal.]</figcaption></figure>

Based on the FEM simulations of this boundary value problem, the pros and cons of each relationship
can be drafted. The power based law has the advantage of reproducing a two-step cumulative eroded mass
with time. On the other hand, the energy based law better tackles the total amount of eroded mass obtained
towards the end of the test. From a computational point of view, the power based approach is put at a
disadvantage by the dependence of the Δ*E*<sub>stage</sub> parameter on the time evolution of the hydraulic loading
path. Finally, in terms of material parameters, the energy based law presents the downside of a dependence
of the *Ē*<sub>max</sub> parameter to the seepage length, which implies that at least two experimental tests should be
realized to obtain this relation and to free the element size from any related constraint.

The above results and discussions raise the need for the following improvements that will be tackled
in future works: (i) perform additional experimental tests on a physical model of dike with a softer slope
[31, 14], (ii) account for the effect of the fine content on the permeability [37, 14] and (iii) provide an additional suffusion law dedicated to the case where self-filtration prevails upon detachment [29, 42].

## 5. Conclusions

This study has been an attempt to provide a novel contribution to the numerical modelling of soil suffusion based on experimentally measurable parameters. During suffusion experimental tests, three mechanisms occur simultaneously that cannot be easily distinguished: the detachment of fine particles under the action of flow, the transport of these particles by the fluid and the self-filtration of some but not all of these fluidized particles. The work consisted in modelling this complex volumetric internal erosion through the mass exchange between the solid and the fluid phases. The mass exchange term accounted for the summation of both positive and negative contributions, i.e. detachment and self-filtration, into a unique term. In addition, this mass exchange term was restricted to positive or null values implying that overall detachment prevailed upon self-filtration. Each phase was composed of two species: the solid phase comprised non-erodible grains and erodible particles, while the fluid phase gathered the water and fluidized particles. Governing differential equations were formulated based on the mass balance for each species and on the momentum balance for the mixture as a whole. The hydro-mechanical coupling was restricted to a classic poro-elastic constitutive model while the main focus was placed on two suffusion behaviour models: the power-based relationship and the energy-based relationship.

The presented experiments have shown that the applied hydraulic loading path influenced the kinetics of suffusion. In particular multi-stages hydraulic paths favour the evolution of the suffusion process up to a stable state, which is physically interpreted as the end of the suffusion process. One specimen, which was suffused up to this stable state and which displayed an averaged dry unit weight, was chosen as the reference specimen for calibration. For the energy-based relationship, two material parameters, namely the erosion resistance index and the volumetric maximum energy were measured on this reference specimen while the smoothing time was fitted. For the power-based relationship, two material parameters were measured by capping the six specimens responses while the intrinsic time was fitted on the reference specimen. A calibration procedure was also proposed to account for the effect of the initial dry unit weight on some material parameters. The comparison of models results with experimental data showed that for realistic suffusion predictions non-linear behaviour laws are required. Although both approaches were able to predict qualitatively well the cumulative eroded mass for various hydraulic loading paths, the energy-based relation displayed the ability to control the total cumulative eroded mass.

The numerically predicted suffusion of a physical model of dike showed a reasonable agreement with the experimental data; especially with respect to the total amount of eroded mass and the ultimate spatial distribution of fines, when using the energy-based relationship. Experimentally, the suffusion test on the physical model of dike led to a sharp initial peak of the flow rate followed by a smoother increase, to a three stage evolution of the cumulative eroded mass and to a heterogeneous ultimate spatial distribution of the percentage of fines. Indeed more fines are eroded near the downstream slope, while self-filtration was observed near the upstream side. In addition, a significant early sliding was experimentally observed which probably favoured a loose preferential flow path. The modelling of the physical model of dike demonstrated the applicability of the energy-based relationship to simulate the evolution of suffusion with time for the case where detachment prevails upon self-filtration.

## Acknowledgments

This work is supported by two PhD fellowships provided by the Lebanese University (UL) and the National Council for Scientific Research of Lebanon (CNRS-L). This work also benefited from the support of a PHC CEDRE grant. The authors also wish to thank the sponsors of EDF-CIH for supporting this research.
**Data Availability Statement**

The data that support the findings of this study are available from the corresponding author upon
reasonable request.

**References**

[1] Bendahmane F, Marot D, Alexis A. Experimental parametric study of suffusion and backward erosion, Journal of Geotechnical and Geoenvironmental Engineering, 2008,134(1):57-67.

[2] Biot MA. General theory of three-dimensional consolidation, Journal of Applied Physics, 1941,12:155-164.

[3] Bonelli S, Marot D. Micromechanical modeling of internal erosion, European Journal of Environmental and Civil Engineering, 2011,15(8):1207-1224.

[4] Bridle R, Fell R. Internal erosion of existing dams, levees and dikes, and their foundations, ICOLD (International Commission on Large Dams), Bulletin 164, volume 1: internal erosion processes and engineering assessment,2013.

[5] Bui TA, Gelet R, Marot D. Modeling of internal erosion based on mixture theory: General framework and a case study of soil suffusion, International Journal for Numerical and Analytical Methods in Geomechanics, 2019,43(15):2407-2430.

[6] Chang DS, Zhang LM. Critical hydraulic gradients of internal erosion under complex stress states, Journal of Geotechnical and Geoenvironmental Engineering, 2012,139(9):1454-1467.

[7] Chang D, Zhang LM. Extended internal stability criteria for soils under seepage, Soils and Foundations, 2013,53(4):569-583.

[8] Cividini A, Gioda G. Finite-element approach to the erosion and transport of fine particles in granular soils, International Journal of Geomechanics, 2004,4(3):191-198.

[9] Coussy O. Poromechanics, John Wiley & Sons, 2004.

[10] Foster M, Fell R, Spannagle M. The statistics of embankment dam failures and accidents, Canadian Geotechnical Journal, 2000,37(5):1000-1024.

[11] Fry J-J. Introduction to the Process of Internal Erosion in Hydraulic Structures: Embankment Dams and Dikes, Wiley, In Erosion of Geomaterials, Bonelli S, 2013.

[12] Gelet R. Thermo-hydro-mechanical study of deformable porous media with double porosity in local thermal non-equilibrium, Université Grenoble Alpes and the University of New South Wales, 2011.

[13] Gelet R, Loret B, Khalili N. Thermal recovery from a fractured medium in local thermal non-equilibrium, International Journal for Numerical and Analytical Methods in Geomechanics, 2013,37(15):2471-2501.

[14] Horikoshi K, Takahashi A. Suffusion-induced change in spatial distribution of fine fractions in embankment subjected to seepage flow, Soils and Foundations, 2015,55(5):1293-1304.

[15] Hosn RA, Sibille L, Benahmed N, Chareyre B. A discrete numerical model involving partial fluid-solid coupling to describe suffusion effects in soils, Computers and Geotechnics, 2018,95:30-39.

[16] Kenney TC, Lau D. Internal stability of granular filters: Reply, Canadian Geotechnical Journal, 1986,23(3):420-423.

[17] Kenney TC, Lau D. Internal stability of granular filters, Canadian Geotechnical Journal, 1985,22(2):215-225.

[18] Kodieh A. Poro-elastic numerical modeling of internal erosion: development of an energy-based approach, Université de Nantes, 2020.

[19] Kodieh A, Gelet R, Marot D, Fino A. A study of suffusion kinetics inspired from experimental data: comparison of three different approaches, Acta Geotechnica, 2020,https://doi.org/10.1007/s11440-020-01016-5

[20] Le VT. Development of a new device and statistical analysis for characterizing soil sensibility face suffusion process, Université de Nantes, 2017.

[21] Le VT, Marot D, Rochim A, Bendahmane F, Nguyen H-H. Suffusion susceptibility investigation by energy-based method and statistical analysis, Canadian Geotechnical Journal, 2018,55(1):57-68.

[22] Li M, Fannin RJ. Comparison of two criteria for internal stability of granular soil, Canadian Geotechnical Journal, 2008,45(9):1303-1309.

[23] Marot D, Bendahmane F, Rosquoet F, Alexis A. Internal flow effects on isotropic confined sand-clay mixtures, Soil & sediment contamination, 2009, 18(3):294-306.

[24] Marot D, Regazzoni P-L, Wahl T. Energy-based method for providing soil surface erodibility rankings, Journal of Geotechnical and Geoenvironmental Engineering, 2011,137(12):1290-1293.

[25] Marot D, Le V-D, Garnier J, Thorel L, Audrain P. Study of scale effect in an internal erosion mechanism: centrifuge model and energy analysis, European Journal of Environmental and Civil Engineering, 2012,16(1):1-19.

[26] Marot D, Rochim A, Nguyen H-H, Bendahmane F, Sibille L. Assessing the susceptibility of gap-graded soils to internal erosion: proposition of a new experimental methodology, Natural Hazards, 2016,83(1):365-388.

[27] Papamichos E, Vardoulakis I. Sand erosion with a porosity diffusion law, Computers and Geotechnics, 2005,32(1):47-58.

[28] Papamichos E, Vardoulakis I, Tronvoll J, Skjaerstein A. Volumetric sand production model and experiment, International Journal for Numerical and Analytical Methods in Geomechanics, 2001,25(8):789-808.

[29] Reboul N, Vincens E, Cambou B. A computational procedure to assess the distribution of constriction sizes for an assembly of spheres, Computers and Geotechnics, 2010,37(1-2):195-206.

[30] Reddi LN, Lee IM, Bonala MVS. Comparison of internal and surface erosion using flow pump tests on a sand-kaolinite mixture, Geotechnical Testing Journal, 2000,23(1):116-122.

[31] Rochim A. Characterization of suffusion susceptibility of granular soils, Université de Nantes, 2015.
[32] Rochim A, Marot D, Sibille L, Le VT. Effects of Hydraulic Loading History on Suffusion Susceptibility of Cohesionless Soils, Journal of Geotechnical and Geoenvironmental Engineering, 2017,143(7):04017025.

[33] Rousseau Q, Sciarra G, Gelet R, Marot D. Modelling the poroelastoplastic behaviour of soils subjected to internal erosion by suffusion, International Journal for Numerical and Analytical Methods in Geomechanics, 2020,44(1):117-136.

[34] Sibille L, Lominé F, Poullain P, Sail Y, Marot D. Internal erosion in granular media: direct numerical simulations and energy interpretation, Hydrological Processes, 2015,29(9):2149-2163.

[35] Sibille L, Marot D, Sail Y. A description of internal erosion by suffusion and induced settlements on cohesionless granular matter, Acta Geotechnica, 2015,10(6):735-748.

[36] Skempton A, Brogan J. Experiments on piping in sandy gravels, Geotechnique, 1994, 44(3):449-460.

[37] Smith M, Konrad J-M. Assessing hydraulic conductivities of a compacted dam core using geostatistical analysis of construction control data, Canadian Geotechnical Journal, 2011,48(9):1314-1327.

[38] Stavropoulou M, Papanastasiou P, Vardoulakis I. Coupled wellbore erosion and stability analysis, International journal for numerical and analytical methods in geomechanics, year=1998,22(2):749-769.

[39] Taha H, Nguyen NS, Marot D, Hijazi A, Abou-Saleh K. Micro-scale investigation of the role of finer grains in the behavior of bidisperse granular materials, Granular Matter, 2019,21(2):1-28,

[40] Tran D-M. Suffusion susceptibility characterization by taking into consideration the context of engineering practice, 2020, Université de Nantes

[41] Vardoulakis I, Stavropoulou M, Papanastasiou P. Hydro-mechanical aspects of the sand production problem, Transport in porous media, 1996,22(2):225-244

[42] Vincens E, Witt KJ, Homberg U. Approaches to determine the constriction size distribution for understanding filtration phenomena in granular materials, Acta Geotechnica, 2015,10(3):291-303.

[43] Wautier A, Bonelli S, Nicot F. DEM investigations of internal erosion: Grain transport in the light of micromechanics, International Journal for Numerical and Analytical Methods in Geomechanics, 2019,43(1):339-352.

[44] Yang J, Yin Z-Y, Laouafa F, Hicher P-Y. Internal erosion in dike-on-foundation modeled by a coupled hydromechanical approach, International Journal for Numerical and Analytical Methods in Geomechanics, 2019,43(3):663-683.

[45] Yang J, Yin Z-Y, Laouafa F, Hicher P-Y. Modeling coupled erosion and filtration of fine particles in granular media, Acta Geotechnica, 2019,14(6):1615-1627.

[46] Zhang L, Gelet R, Marot D, Smith M, Konrad J-M. A method to assess the suffusion susceptibility of low permeability core soils in compacted dams based on construction data, European Journal of Environmental and Civil Engineering, 2019,23(5):626-644.

[47] Zhang XS, Wong H, Leo CJ, Bui TA, Wang JX, Sun WH, Huang ZQ. A thermodynamics-based model on the internal erosion of earth structures, Geotechnical and Geological Engineering, 2013,31(2):479-492.

[48] Zhong C, Le VT, Bendahmane F, Marot D, Yin Z-Y. Investigation of spatial scale effects on suffusion susceptibility, Journal of Geotechnical and Geoenvironmental Engineering, 2018,144(9):04018067.