# Internal erosion of sandy gravel and occurrence of open-framework gravels in the subsoil of a river dike

Stéphane Bonelli, Laurence Girolami

► To cite this version:

Stéphane Bonelli, Laurence Girolami. Internal erosion of sandy gravel and occurrence of open-framework gravels in the subsoil of a river dike. Geomechanics for Energy and the Environment, 2025, 42, 10.1016/j.gete.2025.100690. hal-04895821

HAL Id: hal-04895821

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

Submitted on 18 Jan 2025

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
# Internal erosion of sandy gravel and occurrence of open-framework gravels in the subsoil of a river dike

Stéphane Bonelli<fnref n="1" />, Laurence Girolami<fnref n="1" />,<fnref n="2" />

Corresponding author: Stéphane Bonelli; E-mail: stephane.bonelli@inrae.fr

## Abstract

When a river dike is built on a sandy gravel paleo-valley, successive floods can induce internal erosion. This is the subject of this work, with a finite element analysis of a river dike system. This type of analysis makes it possible to find artesian and uplift zones in the protected floodplain, which is an element to be integrated into flood hazard mapping. The study area is the River Agly in southern France, where numerous leaks, sand boils and sinkholes have been observed along the dikes. The aim is to better understand the origin of these surface signatures, as well as the cause of the presence of open-framework gravel in the subsoil. A suffusion model for sandy gravel was used to describe internal erosion. Internal erosion effectively transforms the sandy gravel into gravel, revealing open-framework gravel zone in the paleo-valley. Contact erosion in gravel can be triggered by suffusion, showing that new models coupling suffusion and contact erosion are needed to model internal erosion in sandy gravels.

## Keywords

Internal erosion, suffusion, dike, levee, paleo-valley, paleo-channel, open-framework gravel

## 1 Introduction

Internal erosion is one of the main physical processes leading to the degradation of earthen dams and dikes (Bonelli, 2013). Erosion can occur in the body of the embankment or in its foundation, i.e. in the subsoil beneath the embankment. Subsoils beneath river dikes and in the protected floodplain can be particularly affected by internal erosion processes during flood

<footnote n="1">¹ INRAE Aix-Marseille Université UMR RECOVER, Aix-en-Provence, France</footnote>

<footnote n="2">² GéHCO, Campus Grandmont, Université de Tours, Tours, France</footnote>
periods. These phenomena occur in layers generally made up of highly permeable, coarse-grained
fluvial deposits formed by the presence of paleo-valleys and paleo-channels (Kolb, 1975; Wolff,
2002; Semmens and Zhou, 2019; Girolami et al., 2023). The surface manifestation of these
underlying phenomena generally takes the form of leaks, sand boils, and sinkholes (Semmens and
Zhou, 2019; Marchi et al., 2021; Van et al., 2022; Zwanenburg et al., 2017).

Among the various physical processes involved in internal erosion, suffusion is currently the
most studied, but also the least understood at all levels (Bonelli, 2012): laboratory tests, theoretical
modelling and observations on real hydraulic infrastructures. Suffusion is caused by seepage flow
and corresponds to the detachment and transport of fine particles through the pore domain
formed by the coarse particle matrix.

When the soil under consideration is a sandy gravel, suffusion can transform this bimodal soil
into gravel alone, where the initial sand is no longer present. In geomorphology, a sandy gravel
located in a paleo-valley that has lost its initial sand is called open-framework gravel. Gravelly
fluvial deposits contain open-framework gravel in the form of flat, transverse strata of variable
scale, generally overlain by and resting on bimodal sandy gravels (Lunt and Bridge, 2007). The
presence of open-framework gravel in the subsoil is an important factor in triggering internal
erosion and sinkholes (Luo and Huang, 2020; Fell et al., 2014). This permeable layer can also lead
to the appearance of artesian zones and uplift zones in the protected floodplain during a flood, a
factor that needs to be integrated into flood risk analysis and flood hazard mapping (Julinek et al.,
2020).

This paper deals with the internal erosion finite element analysis of the subsoil of a river dike, when it consists of sandy gravel. Although the finite element method has been a widely used method in engineering for a very long time, it is still very little used to model and analyze suffusion on practical engineering problems such as earth dams and dikes. This type of analysis is particularly powerful for understanding surface observations on dikes, such as sinkholes (Yang et al., 2019), as well as for other hydraulic infrastructures (Yang et al., 2022). In addition, the presence of open framework gravel layers in paleo-valleys is something of a mystery that hasn't been fully explained (Lunt and Bridge, 2007). The possibility that internal erosion could be an element of explanation has not yet been explored: that's what we propose to do here.

The study area covers the Agly river dike system in the southern France, where numerous
leaks, sand boils and sinkholes have been observed (Van et al., 2022; Zwanenburg et al., 2017). In
the work by Girolami et al. (2023), subsurface geophysical images were used as a basis for numerical modeling of Darcy flows, showing that river floods can induce internal erosion of sandy gravel by suffusion. This analysis is continued here with a numerical model of Darcy flow with suffusion.

The paper is organized as follows: Section 2 contains the modelling equations for suffusion in sandy gravel. Section 3 illustrates some of the model's responses to imposed hydraulic gradients, showing the main features of the model. Section 4 deals with finite element numerical modeling of internal erosion in the subsoil of a river dike. Section 5 brings together a series of comments and discussions. Conclusions are presented in Section 6.

# 2 Modelling suffusion in sandy gravel

## 2.1 Modelling flow in porous media

The flow equations in porous media are described by the following three equations

$$ \mathrm{div}\mathbf{q} = 0, \quad -\mathbf{grad}p + \rho_w\mathbf{g} = \mathbf{f}_w, \quad \mathbf{f}_w = \frac{\rho^w\nu^w}{\kappa}\mathbf{q} \qquad (1) $$

where $\mathbf{q}$ is the seepage velocity vector (m/s), $p$ is the pore pressure (Pa), $\mathbf{g} = -g\mathbf{e}_z$ is the vertical gravity vector, $\mathbf{f}_w$ is the interaction vector between fluid and solid phases, $\kappa$ is the intrinsic permeability (m$^2$), $\rho^w$ is the water density (kg/m$^3$), $\nu^w$ is the water kinematic viscosity (m$^2$/s), $g$ is the gravitational constant (m/s$^2$).

The first equation is the mass balance equation of conservation for water in a saturated, non-deformable porous medium, with incompressible flow. The second equation is the momentum equation for laminar flow in the pore domain. The third equation is the constitutive law describing the viscous friction dissipation with Darcy's law.

## 2.2 Description of the sandy gravel

Gravel and sand grains are of density $\rho^p$ (kg/m$^3$). The quantities that can be directly measured on a laboratory sample are the bulk (total) volume $V$, and the sand and gravel masses $M_s$ and $M_g$. The solid volumes of gravel and sand particles are $V_g = M_g/\rho^p$ and $V_s = M_s/\rho^p$. The gravel porosity is $\phi_g = 1 - V_g/V$ and the sand mass fraction is $f_c = M_s/(M_g + M_s)$.
The mixture is a sandy gravel (Figure 1). The total volume V is defined by the gravel matrix, which transmits the force chains. This gravel matrix is here considered as non-deformable and non-erodible. The sand is located in the pore space of the gravel matrix.

The transition value $f_c^*$ defines the gravel-dominated domain (Figure 1), where $f_c^*$ is empirically set at 30% (Prasomsri et al., 2021). A sandy gravel with a fine fraction greater than $f_c^*$ lies in the interaction range between gravel and sand, as far as the transmission of force chains is concerned. In this case, the confining stress has an influence on sand erosion. For a given stress, part of the sand is not erodible: this is the part that transmits the chains of force (Wautier et al., 2019). The influence of stress state on erosion is a complex subject (Kuwano et al., 2021), which will not be addressed here.

The sandy gravel is assumed to be in the gravel-dominated domain: $f_c < f_c^*$. It is therefore assumed that all sand is in a loose state and does not transmit force chains. This sand is considered erodible. Moreover, when $f_c < f_c^*$ erosion occurs without any change in volume.

We take $\phi_s = V_s/V$ to denote the sand volume fraction, which is used here to model suffusion. The porosity of sandy gravel $\phi$ is then

$$ \phi = \phi_g - \phi_s \qquad (2) $$

The relationship between $f_c$ and $\phi_s$ is as follows

$$ \phi_s = (1 - \phi_g) \frac{f_c}{1 - f_c} \qquad (3) $$

We take $\phi_s^*$ to denote the value of $\phi_s$ corresponding to $f_c = f_c^*$.

## 2.3 Modelling the permeability of sand-gravel mixture

The specific permeability $\kappa$ (m²) of a sand-gravel mixture can be quantified using the Kozeny-Carman equation, as follows

$$ \kappa(\phi_s) = \frac{\phi^3(\phi)}{C S^2(\phi_s)} \qquad (4) $$

where $s$ (m²/m³) is the specific surface of the sandy gravel, and $C$ is a dimensionless dimensionless empirical constant that depends on microscopic parameters, such as grain shape and tortuosity. The sum of the surface area of gravel grains $\partial V_g$ and sand grains $\partial V_s$ gives pore domain surface area $S=\partial V_g + \partial V_s$. The specific surface of the sandy gravel $s = S/V$ is

$$ s = \frac{V_g}{V} \frac{\partial V_g}{\partial V_g} + \frac{V_s}{V} \frac{\partial V_s}{\partial V_s} \qquad (5) $$
Considering that gravel and sand are described by spheres of diameter $d_g$ and $d_s$, respectively, the specific surface of the sandy gravel can be defined as follows (Thies-Weesie and Philipse, 1994; van der Hoef et al., 2005)

$$s(\varphi_s) = s_g + \varphi_s \frac{6}{d_s} , \quad s_g = (1 - \varphi_g) \frac{6}{d_g} \qquad (6)$$

Eq. (4) is rewritten as a function of the gravel permeability $\kappa_g$:

$$\kappa(\varphi_s) = \kappa_g \left( \frac{\phi(\varphi_s)}{\phi_g} \right)^3 \left( \frac{s_g}{s(\varphi_s)} \right)^2, \quad \kappa_g = \frac{\phi_g^3}{C s_g^2} \qquad (7)$$

Eq. (7) shows how the sand volume fraction $\varphi_s$ affects the permeability of the sandy gravel $\kappa$ as the porosity $\phi$ increases and the specific surface $s$ decreases when the volume fraction of sand $\varphi_s$ decreases, $\kappa$ increases when $\varphi_s$ decreases. Consequently, the limits of sandy gravel permeability in the gravel-dominated domain are as follows

$$\kappa(\varphi_s^*) \le \kappa(\varphi_s) \le \kappa_g \qquad (8)$$

The hydraulic conductivity $k$ (m/s) is

$$k(\varphi_s) = \frac{g}{v^w} \kappa(\varphi_s) \qquad (9)$$

## 2.4 Modelling suffusion in a sandy gravel

The erosion law chosen is a linear threshold law as a function of velocity (Vardoulakis et al., 1996)

$$\dot{\varphi}_s = -\lambda \varphi_s \max(0, ||\mathbf{q}|| - q_c) \qquad (10)$$

where $q_c$ (m/s) is the critical velocity (erosion threshold) and $\lambda$ (m⁻¹) is the erosion coefficient.

The presence of $\varphi_s$ as a factor ensures consistency: the stock of erodible material is finite. This factor ensures that $\varphi_s$, which can only decrease, remains positive, and that there is no further erosion when $\varphi_s = 0$.

Determining the erosion parameter $\lambda$ remains a totally open question: at this stage it is purely empirical. Critical velocity can be evaluated by correlation with hydraulic conductivity (Côté, 2010):

$$q_c(\varphi_s) = \beta \phi(\varphi_s) \left( \frac{k(\varphi_s)}{k_{ref}} \right)^n \qquad (11)$$

where $\beta$ (m/s) and $n$ (dimensionless) are two empirical dimensionless constants, and $k_{ref} = 1$ m/s. Since $\phi$ and $k$ increase as the sand fraction $\varphi_s$ decreases, $q_c$ increases as $\varphi_s$ decreases.
Gravel has a high hydraulic conductivity and the gravel/sand particle size ratio $d_g/d_s$ is large. For these reasons, and for the sake of simplification, clogging phenomena are unlikely (Chen et al., 2024; Wautier et al., 2017), and the transport and deposition of eroded sand are not modeled.

The hydraulic gradient is $i = \text{grad}H$ where $H = z + p/\gamma_w$ is the hydraulic head (m), with $\gamma_w = \rho_w g$. The critical hydraulic gradient $i_c$ (erosion threshold) is therefore $i_c = k^{-1}q_c$.

### 3 Suffusion modelling on a homogeneous evolution

Input data values are shown in Table 1. Values of sandy gravel parameter as a function of sand fraction are shown in Table 2. We take $\varphi_s^0$ to denote the initial sand volume fraction corresponding to $f_c^0$ via Eq.(3).

The influence of the sand fraction in the sandy gravel domain on hydraulic conductivity $k$ Eq.(9), critical velocity $q_c$ Eq.(11), and critical hydraulic gradient $i_c$ is illustrated in Figure 2, for several values of gravel diameter $d_g$. Variations in $k$ cover more than four orders of magnitude in Figure 2(a). Figure 2(b) shows that $q_c$ increases as $f_c$ decreases, by more than two orders of magnitude. Figure 2(c) shows that $i_c$ decreases as $f_c$ decreases, by one order of magnitude, from 0.02 to 0.12 when diameter $d_g = 2$ mm. Figure 2(c) shows that when then hydraulic gradient is constant and erosion is triggered, it does not stop and continues until $\varphi_s = 0$, i.e. when all the sand has been eroded.

The situation encountered on dams and dikes corresponds to imposed hydraulic gradients. The results given by the erosion law Eq.(9) for different imposed hydraulic gradients are illustrated in Figures 3 and 4. The loading shown in Figure 3(a) corresponds to five stages of increasing hydraulic gradient, each lasting 30 min (a total of 150 min). This type of loading corresponds to that used in the laboratory. The initial value of the sand content is $f_s^0 = 15\%$ (Table 1): we are well within the range of gravel-dominated. Results are shown in Figure 3 for several values of the erosion coefficient $\lambda$. Initial and final values of $\varphi_s, f_c, \phi, \kappa, k, q_c$ and $i_c$ are shown in Table 3. The critical hydraulic gradient $i_c$ decreases from 0.05 to 0.02 (Figure 3(b)). The critical velocity $q_c$ increases from $4.3 \times 10^{-5}$ to $4.1 \times 10^{-4}$ m/s (Figure 3(d)). The hydraulic conductivity $k$ increases from $7.8 \times 10^{-4}$ to $1.9 \times 10^{-2}$ m/s (Figure 3(e)). These are the orders of magnitude observed during laboratory suffusion tests on a similar gap-graded soil (Gelet and Marot, 2022; Marot et al., 2024). The $q_c$ values in Figure 3(d) and Table 2 are lower than the low correlation
values of Cote (2010) ($v_c = 0.126k^{0.586}$ with $v_c = q_c/\phi$ and $k$ in cm/s), with a coefficient of 0.64. The $\beta$ value (Table 1) was determined in order to obtain critical gradients (here between 0.02 and 0.05) of the same order of magnitude as the horizontal critical gradient of 0.02 mentioned by USBR (2019) on Mississippi levees. The gravely sand corresponding to Table 1 is therefore highly erodable.

The initial mass of sand per unit volume is $\rho^s\phi_s^0 = 302.10$ kg/m³. This is the maximum asymptotic value reached by the cumulative eroded sand mass per unit volume $\mu_{er}(t) = \rho^s(\phi_s^0 - \phi_s(t))$ in Figure 3(f). Compared to what is observed in the laboratory (Deng et al., 2023; Ke and Takahashi, 2014; Gelet and Marot, 2022; Liu et al., 2021; Marot et al., 2024), this result corresponds to a soil that erodes rapidly. The suffusion kinetics is driven by the erosion coefficient lambda $\lambda$. The value $\lambda = 0.5$ m⁻¹ in Table 1 therefore corresponds to rapid kinetics.

In Figure 4(a), the imposed hydraulic gradient is triangular, with four maximum values, over a total duration of 2700 min (45 h). This corresponds to a simplified description of the hydraulics encountered under a dike during a flood. The cumulative eroded sand mass per unit volume $\mu_{er}$ is reported in Figure 4(b). For this highly erodable and rapidly eroding soil, all of the sand is eroded from the first upward ramp of the hydraulic loading. This observation will be useful to better analyze the results of a succession of floods on a dike.

# 4 Modelling suffusion in the subsoil of a river dike

## 4.1 Description of the two-dimensional model

The study area covers the Agly River dike system located in southern France, where numerous leaks, sand boils and sinkholes have been observed (Zwanenburg et al., 2017; Van et al., 2022). In the work of Girolami et al. (2023), core analysis show that the subsoil is composed of a layer of sandy gravel underlain by low permeability sandy silt. Geophysical measurements reveal a non-tabular geometry of the interface between permeable soil and low-permeable soils at depth (Figure 5). Under the dike, the permeable layer is around ten meters thick, which corresponds to the paleo-valleys. Under the protected floodplain, the thickness of the permeable layer is less than one meter: this corresponds to a paleo-channel. On the surface, a layer of low-permeability topsoil covers this permeable layer.
In order to demonstrate the influence of this geometry on the internal erosion of the soil supporting the dike, a finite element analysis is undertaken. This modeling does not aim to quantitatively reproduce all measurements and observations. Rather, it aims to better understand why sand boils and sinkholes have been observed near the dike, or a few meters or a few dozen meters from it. To achieve this, the model integrates a certain number of elementary physical processes considered important, on a simplified two-dimensional geometry derived from geophysical images (Figures 5 and 6).

The geometry extends horizontally for 400 m from the dike toe, and 100 m in depth from the surface. The dike is 2 m high. On the river side, the width between the dikes is approximately 60 m wide: assuming infiltration symmetry, half (30 m) is modeled. The material parameters for sandy gravel are those given in Table 1. The four types of soil and their hydraulic conductivity are shown in Figure 5 and given in Table 3.

The finite element mesh and the boundary conditions are shown in Figure 6. The mesh comprises between 13,797 and 23,473 quadratic elements, depending on the geometric assumption. The surface condition of the dike and the protected floodplain is a Signorini condition:

$$p \le 0, \quad q \cdot n \ge 0, \quad pq \cdot n = 0 \qquad (12)$$

where $p$ is the pore pressure, and $n$ is the external unit normal to the surface. This boundary condition allows to consider water outlet situations ($q \cdot n > 0, p = 0$)and situations without water outlet ($q \cdot n = 0, p < 0$). The determination between one or the other of these situations is a result of the calculation, and not an input data. This is a source of nonlinearity in the system to be solved. It is only valid if the water outflow is not too important: the thickness of the water layer on the surface is assumed to be negligible ($p = 0$). If this is not the case, for example in the case of flooding of the protected floodplain by rising water table, it is then necessary to make a coupling with the shallow water equations.

Outside flood periods, the water table is assumed to be at a depth of 4.5 m in the floodplain, corresponding to a depth of 50 cm in the minor bed. The Agly is indeed a perched river, i.e. the ground level in the protected floodplain decreases as one moves away from the dike, and it is assumed that the river is not fed by the water table. During flood periods, the river level reaches a maximum of +2m above ground level in the protected floodplain, and +6m above ground level in the minor bed: the water is at the level of the crest of the dike. Analysis of the largest historical floods shows that the peak flow is reached in around 22.5 hours. The flood hydrograph is modeled
in a simplified way, with a triangular evolution over a duration of 45 h (Figure 6), and a maximum
height at the level of the dike crest. The rate of water level rise is 15 cm/h.

Four geometries are studied (Figure 8). Configuration G1 is close to a tabular geometry, usually considered in works on backward erosion and sand boils. Configurations G2, G3 and G4 incorporate a reduction in the thickness of the permeable layer, from 10 m to 2 m, at x-coordinate $x_1 = x_0$, $x_0 + 12$ m, and $x_0 + 24$ m respectively, where $x_0$ is the x-coordinate of the dike toe on the protected zone side.

**4.2 Numerical modeling results**

The pore pressure field at the flood peak ($H_{up} = H_{max}$) corresponding to G3 geometry before internal erosion is shown in Figure 9. The unsaturated zones are considered here in a simplified manner: the hydraulic conductivity is constant and equal to the saturated hydraulic conductivity, and the transient term driven by the capillary capacity is not taken into account.

The volume fraction field of eroded sand $1 - \phi_s(t)/\phi_s^0$ is shown in Figure 10 for the four geometries, at different times $t$, where $\phi_s^0$ is the initial volume fraction of the sand.

The comments are as follows:

- Erosion appears after the first flood at three locations:

i) at the dike toe on the river side (water inlet) for all configurations;

ii) at the dike toe on the protected area side (water outlet) for all configurations;

iii) at location $x_1$ for G3 and G4, inside the permeable layer, at the location of the decrease in the thickness of the permeable zone.

- Erosion develops in two locations, essentially in a progressive manner (from the river to the protected floodplain), and to a lesser extent, in a regressive manner:

i) under the dike for all configurations;

ii) at position $x_1$ for G2, G3 and G4, at the location of the decrease in the thickness of the permeable zone.

Figure 10 illustrates a situation where the low permeability topsoil (thickness $e$) covers the permeable layer. During a flood, river water infiltrates into the permeable zone and excess pore pressure can act on the base of the topsoil layer.

We take $p_e$ to denote the pore pressure and $H_e$ to denote the piezometric level below the topsoil, where $H_e = p_e/\gamma_w - e$ is defined at depth $z = -e$. The artesian condition occurs when the
piezometric level is above ground level. It is defined by $p_e > \gamma_w e$. Since the ground level is here at $z = 0$, the artesian condition becomes $H_e > 0$.

The total vertical stress acting on the base of the topsoil is $\sigma_z = -\gamma e$ where $\gamma$ is the soil density. The uplift condition occurs when the pore pressure exceeds the submerged unit weight of the soil, $p_e > |\sigma_z|$, where the uplift pressure is $p_e = \gamma_w (H_e + e)$. In a simplified analysis, the order of magnitude of the saturated soil density is $\gamma \approx 2\gamma_w$. The uplift condition becomes $H_e > e$.

Figure 11 shows values of $H_e$ for the four geometries studied, plotted as a function of the horizontal position $x$ at the initial time (no erosion) and after several floods. The observations are as follows:

- The geometric hypothesis (G1-G4) has a limited influence on the decreasing shape of the piezometric level profiles in the protected floodplain.

- All piezometric levels increase with the number of floods.

- The geometric hypothesis (G1-G4) has an influence on the extent of the artesian and uplift zones in the protected floodplain.

Figure 11 shows the increase in the piezometric level $\Delta H_e$ under the topsoil layer for the four geometries studied. It is plotted as a function of the horizontal position $x$ after several floods, where $\Delta H_e(t) = (p(t) - p(0))/\gamma_w$ is defined at depth $z = -e$. The observations are as follows:

- G4 is the only geometry for which the increase $\Delta H_e$ is significant after the first flood.

- G4 is the only geometry for which internal erosion leads to a decrease in pore pressure under the dike ($\Delta H_e < 0$), after the first flood.

- G4 is the geometry for which the increase is the greatest (of the order of 1 m after 20 floods at $x = 35$ m approximately, or 23 m from the dike toe).

- The maximum increase is of the same order of magnitude for G1, G2 and G3 (approximately 60 cm), but it is not located in the same place. For G2 (resp. G3), it is closer (resp. further away) from the dike than G1.

Figure 13 shows the piezometric level $H_e$ under the topsoil as a function of the number of floods, at several positions $x_0$ from the dike toe. The observations are as follows:

- At the dike toe ($x_0=0$), the piezometric levels are initially artesian and above uplift threshold ($H_e > 1$ m). They increase with the number of floods, then reach a maximum which is of the same order of magnitude for the four geometries (of the order of 1.60 m). The increase is slowest for G3 geometry.
- At $x_0=12$ m from the dike toe, the piezometric levels are initially artesian ($H_e > 0$) but below the uplift threshold ($H_e < 1$ m). They increase with the number of floods (more slowly than at $x_0=0$) to reach values greater than the uplift threshold ($H_e > 1$ m). A date therefore appears from which there is uplift at $x_0=12$ m. This date depends on the geometry, and is between 4 and 9 floods. The piezometric levels stabilize at a maximum which is of the same order of magnitude for the four geometries (around 1.30 m, above uplift threshold).

- At $x_0=24$ and 36 m from the dike toe, the piezometric levels are artesian ($H_e > 0$) but below the uplift threshold ($H_e < 1$ m). Over the duration studied (20 floods), the maximum stabilized levels are not reached and the situation remains evolving.

Figure 14 shows the evolution of the artesian zone size $l_a$ and the evolution of the uplift zone size $l_u$ in the protected area for the four geometries studied. It is plotted as a function of number of floods. The observations are as follows:

- The artesian zone size $l_a$ (defined by $H_e > 0$) can extends over several tens of meters in the protected area.

- The uplift zone size $l_u$ (defined by $H_e > e$) initially extends over a few meters in the protected area (less than 10 m); it increases very significantly (by approximately 10 to 15 m) with the number of floods.

- The largest artesian zone is given by G1 (above 60 m).

- The largest uplift zone is given by G3 up to 3 floods, then by G4 beyond that (up to 23 m).

- The geometric hypothesis (G1-G4) has a significant influence on $l_a$ and $l_u$. This influence is complex, and varies according to the number of floods.

# 5 Comments on the news results obtained and open questions

## *Open-framework gravels created by internal erosion in the paleo-valley*

The presence of open-framework gravel layers in paleovaleys is a mystery that has not been completely elucidated (Lunt and Bridge, 2007). The present results show that internal erosion transforms sandy gravels into gravels, revealing open-framework gravel zones in the paleo-valley. These are the red areas in Figure 10. This is a new finding. Until now, the causes put forward were based on surface hydraulics and sediment transport. Supplementing these hypotheses with arguments based on internal erosion opens up new perspectives on the subject. This result requires confirmation by more in-depth analyses.
## Regarding leaks occurrence in the protected floodplain

Due to the high permeability of the sandy gravel layer, pressure losses are low and pore pressures are high. In the presence of low-permeability topsoil, high pore pressures result, from the first flood, in the appearance of artesian zones and uplift zones in the protected area.

In the artesian zones, the slightest defect in the topsoil can lead to leaks (concentrated outflow of seepage water, also known as boils) in the protected area: for example, a fluidized sand vein (heave) or a pre-existing vertical crack. In the uplift zones, the cohesive topsoil is subjected to hydraulic failure (cracking).

## Regarding sand boils occurrence scenarios

When a topsoil defect causes leakage, this can initiate erosion of the sand contained in the sandy gravel (Figure 15(d)). The sand is fluidized and transported by the outgoing flow, and is deposited on the soil surface. This erosion creates a pipe through the gravel matrix (not a conventional water-filled pipe). This is a backward suffusion, which is also a normal contact erosion at the sand/flow interface (Figure 18). This is a different process than that commonly considered to explain sand boils and their consequences. The usual assumption is that the permeable layer is sand, and that backward erosion creates a pipe under the topsoil (backward erosion piping, van Beek et al., 2011).

This result raises new open questions. Fluidization and sedimentation are intimately linked (Amin et al., 2021; Girolami et al., 2024). In this configuration, little is known about sedimentation (Chen et al., 2024). What are the conditions for fluidization of sand within the gravel matrix?

How does such a pipe progress through the gravel matrix? What are the conditions required for this pipe to connect the river and the protected floodplain? Can this process lead to a levee failure and flooding of the protected floodplain?

## On erosion within the permeable layer and sinkhole occurrence

The initiation and development of erosion in the subsoil beneath the dike is an expected result: this is where the hydraulic gradient is highest. On the other hand, the initiation and development of erosion within the permeable layer, far from the boundary conditions, and far from the dike, is a novel result. This is a crucial result, driven by the non-tabular geometry of the interface between the permeable layer and the low-permeability layer. This result needs to be confirmed by further analysis.
Erosion occurring at the entrance to the paleo-channel, at the level of the reduction in the thickness of the sandy gravel layer (Figure 15(e)) can induce localized settlement, which will result in a surface sinkhole, in the protected floodplain. It is important to note that these sinkholes can occur without the prior appearance of sand boils, contrary to what is often assumed (Zwanenburg et al., 2017; Van et al., 2022).

### On the appearance of surface signatures on the long term (leaks, sand boils, sinkholes)

The presence of a low-permeability top soil is the cause of the uplift zone. A new uplift zone can appear in the protected floodplain after several floods, even though this zone did not initially exist (Figure 13). This is a novel result. This important result shows that even if a dike has withstood several floods, and no disorder has been observed, it should not be extrapolated to future times that this dike will resist other floods, and that no disorder will appear. If no surface signatures (leaks, sand boils, sinkholes) have been observed to date at a given location, this doesn't mean that they will not occur after a future flood.

### Regarding erosion kinetics at field scale

Suffusion is considered to be a very slow internal erosion process, evolving over several years or even decades (Bonelli 2012; Bonelli, 2013). This suggests the need for very long-duration laboratory tests to better study suffusion, or the need to find scale effects for suffusion. The present work shows that the question can be considered differently in the present case.

The results in Figure 4 show that, if we consider erodible and rapidly eroding sandy gravel, all the sand is eroded from the first flood. This is a local observation, based on a homogeneous trend. This type of analysis corresponds to a laboratory suffusion test on a small odometer or triaxial cell.

The results in Figure 10 show the evolution kinetics of the eroded zone (the open-framework gravel zone in red) over 20 floods. The results obtained in Figures 13 and 14 show that the evolution kinetics do not stabilize after 20 floods. These are large-scale results.

These results indicate that the kinetics of internal erosion on site, on large time scales, are not driven by the erosion law in the present case (and in particular not by the erosion coefficient $\lambda$), if we consider erodible and rapidly eroding sandy gravel. The site kinetics are of a different nature. Understanding the parameters influencing the site kinetics, and its formalization, is a work that remains to be done. This is an important result, but it remains to be confirmed.
On the modelling of real historical situations

Although these models take up characteristics of the real dikes, the present numerical models
are too simplified to reproduce exactly the real historical situations. The most important
improvements are:

- The geometry of the subsoil layers, and their main physical characteristics (subsoil stratigraphy). For this, the classic geotechnical techniques (surveys, core sampling, CPT) must be supplemented by geophysical images, among which the equally classic electrical methods are effective (EMI, ERT) (Girolami et al., 2023). This is one of the important points of this work, which shows the influence of stratigraphy on the initiation and development of internal erosion.

- The embankment and its drainage system. They must be modelled in a manner close to the existing system.

- Hydraulic boundary conditions. They must be derived from measurements: actual water levels in the river and in the protected flood plain (in depth).

- Calibration and validation of the model. The presence of an on-site measurement system would make it possible to have measurements during a flood. Piezometric measurements, for example, would make it possible to calibrate the model and then validate it. Unfortunately, river dikes are rarely systematically equipped with a flood monitoring system (unlike dams, this is not mandatory).

A new understanding of the zones of internal erosion in the subsoil of a river dike

From Figure 10, four main erosion zones can be identified (Figure 15):

(A) At the interface between the gravel formed by internal erosion and the sandy gravel;

(B) At the interface between the gravel formed by internal erosion and the base of the sandy silt of the topsoil;

(C) In a defect crossing the topsoil (a hole);

(D) At the interface between the gravel formed by internal erosion and the roof of the sandy silt at depth.

These are new results concerning the internal erosion of the subsoil of a river dike, when this
dike is located on a paleo-valley filled with sandy gravels: the state of knowledge does not reflect
this (Bonelli 2012; Bonelli 2013; Van et al., 2022; Zwanenburg et al., 2017).
*Coupling suffusion and contact erosion and new erosion processes?*

Figure 15 illustrates four erosion processes, detailed in Figures 16, 17 and 18:

(b) Flow in the gravel erodes the sand of the gravely sand (Figure 16);

(b) Water flowing from the sand into the gravel erodes the sand of the gravely sand (Figure 16);

(c) Flow in the gravel erodes the silt of the topsoil and the silt at depth (Figure 17);

(d) Water flowing from the sand into the defect fluidizes and transports the sand of the gravely sand (Figure 18). A pipe through the gravel matrix is formed by suffusion.

Starting from suffusion, we finally arrive at identifying four erosion processes which are not explicitly identified by the state of knowledge. In particular, it appears that contact erosion can be triggered by suffusion. This is a new result, which must be analyzed in depth by complementary work that requires new models coupling suffusion and contact erosion.

*On the modelling assumptions and the erosion law*

The modeling hypotheses can be eliminated one after the other, making the model increasingly complex. Here are some examples that can be made from what is already present in the literature: the transport of eroded soil particles through the pore domain and their possible deposition; unsaturated areas; a nonlinear law for flow in gravel (e.g. Ergün-type law); a three-dimensional geometry to better consider the geometry of the paleo-valley and model leaks and sandboils; a soil stress-strain behaviour law to model the consequences of internal erosion (e.g. sinkholes).

However, there are hypotheses that correspond to totally open questions. There are convincing models that reproduce experimental observations well, such as the suffusion of sandy gravels (Deng et al, 2023), but no model has reached consensus. There is no consensus on the very notion of critical hydraulic gradient or critical velocity, and there are several approaches (Gelet and Marot, 2022). The question of the erosion law for suffusion is totally open. The driving force to be associated with the erosion rate by a behavior law is still unknown! We still don't know whether it's pressure gradient, Darcy velocity or pore velocity, hydraulic energy or power, or something else.

The suffusion law Eq.(10) is deliberately simple, as this paper is primarily a theoretical and qualitative study of a real-life situation. This model has two main drawbacks. For a constant hydraulic gradient, the erosion process will not cease until the complete loss of all fine particles. Observations show that the quantity of particles eroded depends on the hydraulic gradient and the
confinement stress. This can be taken into account in the suffusion law (Bonelli and Marot, 2011; Deng et al., 2023). In addition, analysis of the results obtained in this work shows that it is necessary to develop new models coupling suffusion and contact erosion. Sampling sandy gravel soils is extremely difficult, if not impossible. Internal erosion tests on intact samples in the laboratory are out of the question. For this reason, the parameters of internal models must be quantified physically, without laboratory testing, or with large-scale modelling, by comparing results with observations and field measurements.

# 6 Conclusion

This numerical study concerns the finite element modelling of a dike located on a sandy gravel paleo-valley, and subjected to successive floods. The internal erosion of sand by suffusion leads to the creation and development of gravel zones in the subsoil, similar to open-framework gravel. The results obtained lead to a new understanding of the internal erosion zones in the subsoil of a river dike. In particular, it appears that suffusion can trigger contact erosion. Regarding sand boils occurrence scenarios, the observed process is a different process than that commonly referred as to backward erosion piping. The main limitation of this work is the erosion law when the eroded material is highly permeable, which opens up new questions for internal erosion models. Coupling suffusion with contact erosion is essential to model flow-induced erosion in open-framework gravel. Coupling with fluidization of sand within the gravel matrix is necessary to model sand-boil initiation.

# 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.

# Fundings

This research was supported by the Water departement of INRAE (Aqua department, Pari Scientifique), by Aix-Marseille University (ECCOREV/ITEM 2023 PALERO), by the Directorate General for Risk Prevention (DGPR), by the Water and Environment INRAE Carnot Institute (AquaGeoPhi project), and by the Region SUD in France (AquaGeoPhi project).
References

Amin A., Girolami L., Risso F. (2021) On the fluidization/sedimentation velocity of a homogeneous suspension in a low-inertia fluid, *Powder Technology* **391**:1-10.

van Beek V. M., Knoeff H., Sellmeijer H. (2011) Observations on the process of backward erosion piping in small-, medium-and full-scale experiments, *European Journal of Environmental and Civil Engineering* **15**(8):1115-113.

Bonelli S. (2012) *Erosion of Geomaterials*, (Ed), John Wiley & Sons, 372 p.

Bonelli S. (2013) *Erosion in Geomechanics Applied to Dams and Levees*, (Ed), Wiley-ISTE, 388 p.

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

Chen F., Wautier A., Philippe P., Benahmed N., Nicot F. (2024). Micromechanics of fine-grain infiltration in coarse grain sands, *Acta Geotechnica*.

Côté S. (2010) Internal stability criteria for low plasticity materials subjected to water flow (in French), PhD Thesis, Université Laval, 2010.

Deng Z., Wang G., Jin W., Tang N., Ren H., Chen X. (2023) Characteristics and quantification of fine particle loss in internally unstable sandy gravels induced by seepage flow, *Engineering Geology*, **321**:107150.

Fell R., MacGregor P., Stapleton D., Bell G., Foster M., 2014. *Geotechnical Engineering of Dams*, 2nd Edition. 419 ed. CRC Press.

Gelet R., Marot D. (2022). Internal erosion by suffusion on cohesionless gap-graded soils: Model and sensibility analysis, *Geomechanics for Energy and the Environment* **31**:100313.

Girolami L., Bonelli S., Valois R., Chaouch N., Burgat J. (2023) On Internal Erosion of the Pervious Foundation of Flood Protection Dikes, *Water* **15**(21):3747.

Girolami L., Risso F., Amin A., Rousseau L., Bondessan A., Bonelli S. (2024) Sedimentation of short-lived fluid-solid suspensions, *Physics of Fluids* **36**:113308.

van der Hoef M.A., Beetstra R., Kuipers J.A.M. (2005) Lattice-Boltzmann simulations of low-Reynolds-number flow past mono- and bidisperse arrays of spheres: results for the permeability and drag force, *Journal of Fluid Mechanics* **528**:233-254.

Julinek T., Duchan D., Riha J. (2020). Mapping of uplift hazard due to rising groundwater level during floods, *Journal of Flood Risk Management* **13**:e12601.
Kolb C.R. (1975) Geologic Control of Sand Boils Along Mississippi River levees. USACE Vicksburg, Report S-75-22, 35 p.

Ke L., Takahashi A. (2014). Experimental Investigations on Suffusion Characteristics and Its Mechanical Consequences on Saturated Cohesionless Soil, *Soils and Foundations*, **54**(4):713-730.

Kuwano R., Santa Spitia L. F., Bedja M., Otsubo M. (2021). Change in Mechanical Behaviour of Gap-Graded Soil Subjected to Internal Erosion Observed in Triaxial Compression and Torsional Shear, *Geomechanics for Energy and the Environment*, **27**:100197.

Liu K., Qiu R., Su Q., Ni P., Liu B., Gao J., Wang T. (2021). Suffusion response of well graded gravels in roadbed of non-ballasted high speed railway, *Construction and Building Materials*, **284**:122848.

Lunt I.A., Bridge J.S. (2007) Formation and preservation of openframework gravel strata in unidirectional flows, *Sedimentology* **54**:71-87.

Luo Y., Huang Y. (2020) Effect of open-framework gravel on suffusion in sandy gravel alluvium, *Acta Geotechnica* **15**:2649-2664.

Marchi M., Martinez M.F.G., Gottardi G., Tonni G.L. (2021) Field measurements on a large natural sand boil along the river Po, Italy. *Quarterly Journal of Engineering Geology and Hydrogeology* 54(4):qjegh2020-097.

Marot D., Oli B., Bendahmane F., Gelet R., Leroy P. (2024) A hydraulic gradient and stress-controlled erosion apparatus for assessing soil internal erosion, *Geotechnical Testing Journal* **47** (6).

Prasomsri J., Shire T., Takahashi A. (2021) Effect of fines content on onset of internal instability and suffusion of sand mixtures, *Géotechnique Letters* **11**(3):209-214.

Semmens S.N., Zhou W. (2019) Evaluation of environmental predictors for sand boil formation: Rhine-Meuse Delta, Netherlands, *Env. Earth Sci.* **78**:1-11.

Thies-Weesie D.M.E., Philipse A.P. (1994) Liquid permeation of bidisperse colloidal hard-sphere packings and the Kozeny-Carman scaling relation, *Journal of Colloid and Interface Science* **162**(2):470-480.

USBR, Best Practices in Dam and Levee Safety Risk Analysis, Bureau of Reclamation & US Army Corps of Engineers. Dam Safety Office. Security, Safety and Law Enforcement Office, 2019, 1337 p.

Van M.A., Rosenbrand E., Tourment R., Smith P., Zwanenburg C. (2022) Failure paths for levees, *International Society of Soil Mechanics and Geotechnical Engineering (ISSMGE)*, TC201.
Vardoulakis, I., Stavropoulou, M. & Papanastasiou, P. (1996) Hydro-mechanical aspects of the sand production problem. *Transport in porous media* **22**(2):225-244.

Wautier A., Bonelli S., Nicot F. (2017) Scale separation between grain detachment and grain transport in granular media subjected to an internal flow, *Granular Matter* **19**(22).

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

Wolff T. F. (2002) Performance of levee underseepage controls: A critical review. Rep. No. ERDC/GSL TR-02-19. Washington, DC: USACE.

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

Yang, J., Yin, Z. -Y., Laouafa, F., & Hicher, P. -Y. (2022). Numerical analysis of internal erosion impact on underground structures: Application to tunnel leakage. *Geomechanics for Energy and the Environment, 31*:100378.

Zwanenburg C., López-Acosta N.P., Tourment R., Tarantino A., Pozzato A., Pinto A. (2017) Lessons Learned from Dike Failures in Recent Decades. *International Journal of Geoengineering Case Histories, 4*(2), 203-229.
<table><caption>Table 1. Material parameters for sandy gravel</caption><thead><tr><td>Material parameter</td><td>Notation</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Soil particle density</td><td>ρ<sup>s</sup></td><td>2,650</td><td>kg/m<sup>3</sup></td></tr><tr><td>Water density</td><td>ρ<sup>w</sup></td><td>1,000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Water kinematik viscosity</td><td>η<sup>w</sup></td><td>1.15×10<sup>-6</sup></td><td>m<sup>2</sup>/s</td></tr><tr><td>Gravitational constant</td><td>g</td><td>9.81</td><td>m/s<sup>2</sup></td></tr><tr><td>Gravel porosity</td><td>φ<sub>g</sub></td><td>0.35</td><td>-</td></tr><tr><td>Diameter of gravel particles</td><td>d<sub>g</sub></td><td>2×10<sup>-3</sup></td><td>m</td></tr><tr><td>Diameter of sand particles</td><td>d<sub>s</sub></td><td>2×10<sup>-4</sup></td><td>m</td></tr><tr><td>Critical velocity coefficient</td><td>β</td><td>1.2×10<sup>-2</sup></td><td>m/s</td></tr><tr><td>Critical velocity exponent</td><td>n</td><td>0.586</td><td>-</td></tr><tr><td>Erosion coefficient</td><td>λ</td><td>0.5</td><td>m<sup>-1</sup></td></tr><tr><td>Kozeny-Carman coefficient</td><td>C</td><td>5</td><td>-</td></tr><tr><td>Initial value of the sand mass fraction</td><td>f<sub>c</sub><sup>0</sup></td><td>0.15</td><td>-</td></tr></tbody></table>

<table><caption>Table 2. Sandy gravel parameter values as a function of sand fraction</caption><thead><tr><td>Material parameter</td><td>Notation</td><td>Initial value (at t = 0)</td><td>Value at φ<sub>s</sub> =0</td><td>Unit</td></tr></thead><tbody><tr><td>Sand mass fraction</td><td>f<sub>c</sub></td><td>0.150</td><td>0</td><td>-</td></tr><tr><td>Sand volume fraction</td><td>φ<sub>s</sub></td><td>0.114</td><td>0</td><td>-</td></tr><tr><td>Sandy gravel porosity</td><td>φ(φ<sub>s</sub>)</td><td>0.236</td><td>0.350</td><td>-</td></tr><tr><td>Sandy gravel specific permeability</td><td>κ(φ<sub>s</sub>)</td><td>9.17×10<sup>-11</sup></td><td>2.26×10<sup>-9</sup></td><td>m<sup>2</sup></td></tr><tr><td>Sandy gravel hydraulic conductivity</td><td>k(φ<sub>s</sub>)</td><td>7.82×10<sup>-4</sup></td><td>1.92×10<sup>-2</sup></td><td>m/s</td></tr><tr><td>Critical velocity</td><td>q<sub>c</sub>(φ<sub>s</sub>)</td><td>4.29×10<sup>-5</sup></td><td>4.15×10<sup>-4</sup></td><td>m/s</td></tr><tr><td>Critical gradient i<sub>c</sub> = q<sub>c</sub>/k</td><td>i<sub>c</sub>(φ<sub>s</sub>)</td><td>5.48×10<sup>-2</sup></td><td>2.16×10<sup>-2</sup></td><td>-</td></tr></tbody></table>

<table><caption>Table 3. Hydraulic conductivity of dike soil and subsoil</caption><thead><tr><td>Soil</td><td>Value (m/s)</td></tr></thead><tbody><tr><td>Dike soil</td><td>10<sup>-5</sup></td></tr><tr><td>Topsoil</td><td>10<sup>-6</sup></td></tr><tr><td>Deep sandy silt</td><td>10<sup>-6</sup></td></tr><tr><td>Sandy gravel</td><td>4×10<sup>-5</sup> ≤ k ≤ 2×10<sup>-2</sup></td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1. Diagram of a sand gravel mixture. The transitional fine fractions are determined empirically as $f_c^* = 30\%$ and $f_c^{max} = 35\%$ (Prasomsri et al., 2021).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2. Influence of the sand fraction $f_c$ in the gravely sand domain on the hydraulic conductivity $k$, the critical velocity $q_c$, the critical hydraulic gradient $i_c$, with the parameters of Table 1, and four gravel diameter values $d_g = 1$ mm (C1), 2 mm (C2), 4 mm (C3), 8 mm (C4).</figcaption></figure>
<figure><img src="image_4.png" /></figure>

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

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

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

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

<figure><img src="image_9.png" /><figcaption>Figure 3. Results of a homogeneous loading in hydraulic gradient (a) with the parameters of Table 1, and four erosion coefficient values $\lambda$ = 1 m⁻¹ (C1), 0.5 m⁻¹ (C2), 0.25 m⁻¹ (C3), 0.125 m⁻¹ (C4).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 4. Homogeneous loading in hydraulic gradient (a) of the "flood" type (duration 2700 s = 45 h) with four peak values, and result in terms of cumulated eroded mass per unit volume (b). Calculations carried out with the parameters of Table 1.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 5. Cross-section of the dike and subsoil, with the four types of soil and their hydraulic conductivity. The topsoil is assumed to have the same hydraulic conductivity as the sandy silt at depth.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 6. Hydrographs of four historical floods of River Agly, and simplified model of the flood hydrograph. The rate of water level rise is 15 cm/h. The considered flood duration if 45 h.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7. Finite element mesh and boundary conditions for the numerical simulations.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8. The four geometries of paleo-valley and paleo-channel studied.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9. Water pore pressure at the peak flood ($H_{up} = H_{max}$) corresponding to the geometry G3 before internal erosion.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 10. Volume fraction of eroded sand $1 - \varphi_s/\varphi_s^0$ for the four geometries, at different times. Red areas correspond to open-framework gravel created by suffusion (gravel, no sand, $\varphi_s = 0$). Blue areas correspond to non-eroded areas.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 11. Piezometric level $H_e$ under the topsoil layer as a function of the horizontal position $x$ at the initial time (no erosion) and after 1, 5, 10 and 20 floods, where $H_e = p/\gamma_w - e$ is defined at depth $z = -e$ ($e = 1$ m).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 12. Increase in the piezometric level $\Delta H_e$ under the topsoil layer as a function of the horizontal position $x$ after 1, 5, 10 and 20 floods, where $\Delta H_e(t) = (p(t) - p(0))/\gamma_w$ is defined at depth $z = -e$ ($e = 1$ m).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 13. Piezometric level $H_e$ under the topsoil layer as a function of the number of floods (and according to time in years if there is one flood per year), at positions $x_0 = 0, 12, 24$ and $36$ m from the dike toe located at $x=12$ m, where $H_e = p/\gamma_w - e$ is defined at depth $z = -e$ ($e = 1$ m).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 14. Artesian zone size $l_a$ and uplift zone size $l_u$ in the protected area as a function of the number of floods (and according to time in years if there is one flood per year), from the dike toe located at $x=12$ m.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

Figure 15. Four main erosion zones in the subsoil of a river dike. Three erosion zones at interfaces between the gravel formed by internal erosion and: (A) the sandy gravel, (B) the base of the topsoil, (C) the roof of the sandy silt at depth. An erosion zone at a fault crossing the topsoil (D). Four processes of subsoil erosion. Water flowing from sand to gravel erodes the sand (a). Flow in gravel erodes sand from gravely sand (b), silt from topsoil (c) and silt at depth (e). If there is a defect crossing the topsoil, water flowing from sand to the defect fluidizes and transports the sand (d).

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

Figure 16. Schematic of the occurrence of an open-framework gravel by internal erosion. (a) Flow in the porous domain of gravel erodes sand. This is selective tangential contact erosion at the sand/flow interface. (b) Water flowing out of the sand erodes sand. This is backward suffusion, which also corresponds to normal contact erosion at the sand/flow interface.
<figure><img src="image_23.png" /><figcaption>Figure 17. Diagram of contact erosion of silt by flow in gravel. (c) Flow in the porous domain of gravel erodes silt. This is a classical tangential contact erosion.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 18. Schematic of the occurrence of sand boil with internal erosion. (a) Water leaks out of the sandy gravel through the defect. (b) Fluidization of sand located at the surface in the defect. (c) Flow transports sand to the surface. A pipe forms through the gravel matrix.</figcaption></figure>