# Internal erosion of chemically reinforced granular materials: a mathematical modeling approach

Tony Khalil, Nadia Saiyouri, Bogdan Muresan, Pierre-Yves Hicher

► To cite this version:

Tony Khalil, Nadia Saiyouri, Bogdan Muresan, Pierre-Yves Hicher. Internal erosion of chemically reinforced granular materials: a mathematical modeling approach. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37 (5), pp.491-502. 10.1002/nag.1108 . hal-01006991

HAL Id: hal-01006991

https://hal.science/hal-01006991

Submitted on 11 May 2017

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>
Internal erosion of chemically reinforced granular materials: a
mathematical modeling approach

Tony Khalil¹, Nadia Saiyouri¹, Bogdan Muresan¹,² and Pierre-Yves Hicher¹

¹Ecole Centrale de Nantes, GeM-Institut de Recherche en Génie Civil et Mécanique, Nantes, France

²IFSTTAR, Bouguenais, France

Internal erosion (IE) affects the stability of natural and reinforced materials by causing instability within their granular structure. The dislodgement and transport of eroded particles affect both the particulate con- centration of eroding fluid and the pore network of eroded material. In this study, we examined these mod- ifications using a transport model with a finite element code. First, IE tests on chemically reinforced sand columns were performed to obtain information about eroded material loss of mass, particulate concentration of effluent, porosity and permeability modifications, and existing IE stages. Second, based on experimental results, a mathematical one-dimensional model has been formulated to monitor the evolution and spatial dis- tribution of erodable solids, fluidized particles, porosity, permeability, and seepage stresses. The model con- sists of a set of coupled nonlinear differential equations solved in sequence. It provides valuable information about the extent and the dynamics of structural changes, which can be used to estimate an IE time for the hydraulic work to reach failure.

KEY WORDS: internal erosion; chemically reinforced granular material; erosion tests; transport modeling

# 1. INTRODUCTION

In hydraulic works (like embankment dams or levees), the flow of internal fluids might induce the dislodgement and ensuing migration of solids. In turn, these mechanisms can result in hydraulic and mechanical modifications (like porosity and permeability changes, channelization, loss of tightness, etc.) that might eventually lead to the breakdown of the entire hydraulic work. Such erosion forms are depicted as internal erosion (IE). Close to 46% of accidents concerning embankment dams and sea levees are caused by IE [1]. In France, since 1970, more than 70 incidents of failed dams were directly related to IE [2]. Still, the sum of determining factors initiating IE is not fully understood nor described by any of the numerical models. It is thus relevant to monitor the model changes in the physical characteristics of eroded materials at an early stage of erosion.

Internal erosion is a complex process that includes multiple mechanisms. Four types of mechanisms have been identified in this context [3]: (i) development of defects in the primary fabric of large particles that support imposed stresses (like grain removal inducing stacks faults or dislocations); (ii) suffusion, which is the process whereby the fines fraction is entrained internally through the primary fabric of large particles; (iii) elongated cavities and/or capillaries that are eroded toward the reservoir; and (iv) contact erosion, which occurs at the interface between two solid domains of distinct granulometric distributions, provided the coarse layer is not an appropriate filter for finer
layer. Much research has been conducted in the past to study the IE mechanisms. Significant results are presented in references [4–11]. Most of these works focused on some aspects of IE problem by separately considering two different cases: on the one hand, the effects of IE on the mechanico-physico-chemical properties of eroded materials and, on the other hand, the transport of eroded particles within pores. Yet, because of the mutual relationships between IE mechanisms, the use of a coupled and nonlinear approach is necessary [12].

The objective of this article is to propose a mathematical model of IE and then use it to reproduce experimental test results. The experiment provided information about the loss of mass and hydraulic head, particle concentration in the effluents, and changes in the average porosity (n) and permeability (k) of compacted sand columns (CSCs). To reinforce the CSCs, insoluble silica gel (a binding agent that permits low-permeability material reinforcement: k < $10^{-10}$ m²) was used as a binder. Reinforced CSCs are expected to endure seepage stresses, although their primary fabric of particles proves to be unstable, and/or the binder partially flushed out [13]. The experimental mass and hydraulic head losses were used as inputs to a mathematical one-dimensional IE model that simulates the experiments. The model consisted of a set of coupled nonlinear differential equations solved in sequence. It allowed monitoring the evolution and the spatial distribution of *n*, *k*, pore pressure (*P*), erodible solids content (*C*), and concentration of fluidized particles (*χ*). Current models that treat IE focus on the dislodgement and/or transport of particles that occur during seepage into soils or sands [14–21], whereas, to our knowledge, no data are available on the topic of coupled internal flow effects on pore network of chemically reinforced materials. The present work lies within the scope of a larger research effort intended to probe the mechanisms inherent to IE and ultimately minimize its destructive consequences.

## 2. MODEL DESCRIPTION

A mathematical model is presented to simulate the experimental erosion tests. It is based on the continuum description of a multi-phase system, as detailed in reference [22]. The model depends on two erosion parameters, namely, the erosion coefficient (λ) and the damage coefficient (β). These erosion parameters are determined experimentally, and their values depend on the eroded material. They are assumed to be constant. Unlike the experimental data that provided space-averaged information about ongoing IE mechanisms, the IE simulations allowed to examine the evolution and spatial distributions of *n*, *k*, *P*, *C*, and *χ*.

### 2.1. Mathematical formulation

The saturated porous medium is considered to consist of an erodible rigid skeleton, that is, the Fontainebleau sand and binder, and a fluid phase, that is, the injected water and fluidized particles. Because of their nano-micrometric size (see below), most of the fluidized particles are thought to be miscible with the fluid medium. This means that once dislodged, the transport of eroded particle within the pore network is extensively driven by the flow of fluid. In other words, the model does not explicitly take into account filtration mechanisms and group effects. To deal with particle dislodgement and transport of eroded particles within granular materials, the following must be considered: the concentration of fluidized particles of the eroding fluid ($\chi = m_F/V_P$, where $m_F$ is the mass of solids eroded from the rigid skeleton and now in the fluidized state in the pore water per unit volume $V_P$ of pores, that is, volumes of pore water + fluidized particles for fully saturated media), the content of erodible solids of the eroded material ($C = m_R/V_{CSC}$, where $m_R$ is the mass of solids remaining after erosion per unit volume $V_{CSC}$ of material, that is, volumes of pores + CSC rigid skeleton), and changes of porosity and pore pressure (and thus permeability). Hereafter, the terms concentration and content are used for fluid and solid compositions, respectively. Details on the definitions of various fractions and quantities are provided in reference [14]. From the literature about contaminant transport [23–26], a constitutive mass balance equation was formulated for the fluidized particles:
$$ \frac{\partial(n\chi)}{\partial t} + \nabla(\chi\mathbf{q}) - \nabla(\mathbf{D}\nabla\chi) = \lambda\|\mathbf{q}\|\mathbf{C} \quad (1) $$

where $\mathbf{q}$ is the Darcy flux (in m s⁻¹; boldface denotes that $\mathbf{q}$ is a vector and $\|\mathbf{q}\|$ is its norm), and $\mathbf{D}$ is the hydrodynamic dispersion tensor (including the molecular diffusion coefficient, i.e., $D_0$, the longitudinal, i.e., $a_L$, and transversal, i.e., $a_T$, dispersion coefficients and the tortuosity, i.e., $\mathbf{T}$, of the porous medium) [27]:

$$ \mathbf{D} = a_T \|\mathbf{q}\| \mathbf{I} + (a_L - a_T) \frac{\mathbf{q} \otimes \mathbf{q}}{\|\mathbf{q}\|} + nD_0 \mathbf{T} \quad (2) $$

Equation 1 consists of a traditional advection-dispersion equation that accounts for mass transfer from the erodible solids of the granular material toward the eroding fluid. The left-hand terms drive the evolution and space distribution of the concentration of fluidized particles within the eroding fluid, whereas the right-hand term governs the inputs of erodible particles from the rigid skeleton. The latter is set proportional to the erosion coefficient and is expressed as a linear function of the Darcy velocity of the fluid and content of erodible solids.

$$ \frac{\partial \mathbf{C}}{\partial t} = -\lambda \|\mathbf{q}\| \mathbf{C} \quad (3) $$

If the porosity of the material and the velocity of the fluid remain relatively unchanged, the material content of erodible solids would decrease exponentially with time. Yet, the IE affects the rigid skeleton in such a way that the dislodgement of erodible solids increases locally the volume of pores. Thus, the evolution of porosity

$$ n = n_0 + \frac{(C_0 - C)}{\rho_m} \quad (4) $$

is formulated as a linear function of the content of erodible solids. The subscript '0' indicates that quantities under consideration are taken at the initial time, $t=0$ s. $\rho_m$ is the grain density of eroded materials (1300 kg m⁻³). Environmental scanning electron microscope observations (data not showed) indicated that the solids from the effluent consist of nano-micrometric silica gel fragments from the binder and the finest fraction of the Fontainebleau sand (< 5% in mass) [28]. The ratio $(C_0-C)/\rho_m$ stands for the volume of eroded particles per unit volume of material. The fluid velocity is related to the Darcy law

$$ \mathbf{q} = -\frac{\mathbf{k}}{\mu_f} (\nabla \mathbf{P} + \rho_f \mathbf{g}) \quad (5) $$

where $\mu_f$ is the viscosity of the fluid phase, $\rho_f$ is the density of the fluid, and $\mathbf{g}$ the acceleration of gravity. Higher $\mu_f$ and $\rho_f$ values tend to decrease the mean velocity of the fluid. Preliminary simulations indicated that $\|\mathbf{q}\|$ was highly dependent on $\mu_f$ and almost insensitive to $\rho_f$. Yet, changes in $\mu_f$ were ~5 ranges of magnitude lower than in $\rho_f$, that is, $10^{-6}\%$ versus $10^{-1}\%$, respectively. During the IE tests, $\mu_f$ and $\rho_f$ increases resulted in 1–5% $\|\mathbf{q}\|$ decrease. Without precise information about the evolution laws of both variables with regard to $\mathbf{P}$ and $\chi$, two linear functions were proposed:

$$ \mu_f = (\mu_0 - \mu_w) \frac{\chi}{\rho_0} + \mu_w \quad (6) $$

$$ \rho_f = \rho_0 [1 + \varepsilon_P (\mathbf{P} - \mathbf{P}_0) + \varepsilon_\chi (\chi - \chi_0)] \quad (7) $$

where $\mu_w$ (1.00 10⁻³ Pa s) is the viscosity of water at 1 atm and 20 °C, $\mu_0$ and $\rho_0$ are the initial viscosity (measured with a rheometer; 1.10 10⁻³ Pa s) and density (1050 kg m⁻³) of fluid phase, respectively. The
$\epsilon_P$ (3.8 $10^{-8}$ Pa$^{-1}$) and $\epsilon_\chi$ (2.9 $10^{-4}$ m$^3$ kg$^{-1}$) coefficients were experimentally determined by gauging the volume (resp. weight) of 12-hour degassed ultrapure water ($\rho > 18.2$ M$\Omega$ cm) at various air pressure (resp. silica gel fragment concentrations). These coefficients introduce the effect of density change of the fluid as a result of changes in P and $\chi$, respectively. They also allowed to formulate the fluid mass balance equation.

$$ \frac{\partial n}{\partial t} + \epsilon_P \frac{\partial (nP)}{\partial t} + \epsilon_\chi \frac{\partial (n\chi)}{\partial t} + \nabla \mathbf{q} = 0 \quad (8) $$

Finally, because of the drain of erodible solids, the intrinsic permeability of the granular material is increased during the passing-through of the fluid. Unlike the porosity that accounts for the volume of pores, intrinsic permeability depends on the geometry and spatial/size distribution of constrictions (the small windows between interconnected pores). Hence, n and k cannot be firmly related: slight porosity increase may result in marked permeability changes (either increases, e.g., because of creation of preferential pathways for fluid, or decreases, e.g., clogging of constriction by fluidized particles). For that reason, we decoupled the evolution of both variables by using the linear function below:

$$ \mathbf{k} = \mathbf{k}_0[1 + \beta(C_0 - C)] \quad (9) $$

Equation 8 depends on the initial permeability ($\mathbf{k}_0$), the initial content of erodible solids, and is proportional to the damage coefficient.

## 2.2. Calibration and limitations of the model

The proposed IE model depends on two experimentally determined parameters: $\lambda$ and $\beta$ (Table I). The calibration procedure is detailed in reference [29]. In summary, $\lambda$ values, which determine the IE kinetics, were approached (at constant fluid phase viscosity and without dispersion) from Equation 10 and the monitored mass losses ($\Delta m / \Delta t$) after saturation, that is, when the mass variations are primarily driven by the dislodgement and transport of particles.

$$ \frac{\Delta m}{\Delta t} = \frac{V_{CSC} \|\mathbf{q}\| C_b}{L} \left(1 - \frac{\rho_w}{\rho_m}\right) [1 - \exp(-\lambda L)] \quad (10) $$

$\lambda$ (in m$^{-1}$) is the only unknown of Equation 10 where $C_b=C(0,t)$ (~$C_0$ at early IE stages) is the boundary condition, L is the length of the column (0.1 m), and $\rho_w$ is the density of injected water.

<table><caption>Table I. Parameters used in the numerical internal erosion simulation of $N_{50}$ test.</caption><thead><tr><th></th><th>Notation</th><th>Value of N<sub>50</sub> test</th></tr></thead><tbody><tr><td>Sample size (diameter, length)</td><td>(Ø, L)</td><td>(0.1 m, 0.1 m)</td></tr><tr><td>Initial content of erodible particles</td><td>C<sub>0</sub></td><td>1927 kg m<sup>-3</sup></td></tr><tr><td>Initial concentrations of fluidized particles in pores</td><td>χ<sub>0</sub></td><td>0 kg m<sup>-3</sup></td></tr><tr><td>Initial porosity</td><td>n<sub>0</sub></td><td>0.242</td></tr><tr><td>Initial intrinsic permeability</td><td>k<sub>0</sub></td><td>1.35 10<sup>-12</sup> m<sup>2</sup></td></tr><tr><td>Gravity</td><td>g</td><td>9.81 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Inflow</td><td>Q</td><td>0.1 L min<sup>-1</sup></td></tr><tr><td>Coefficient of longitudinal dispersion</td><td>a<sub>L</sub></td><td>10<sup>-2</sup> m</td></tr><tr><td>Coefficient of transversal dispersion</td><td>a<sub>T</sub></td><td>10<sup>-3</sup> m</td></tr><tr><td>Coefficient of molecular diffusion</td><td>D<sub>0</sub></td><td>10<sup>-10</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Viscosity of injected water</td><td>μ<sub>w</sub></td><td>1.00 10<sup>-3</sup> Pa s</td></tr><tr><td>Initial viscosity of pore water</td><td>μ<sub>0</sub></td><td>1.10 10<sup>-3</sup> Pa s</td></tr><tr><td>Initial density of pore water</td><td>ρ<sub>0</sub></td><td>1050 kg m<sup>-3</sup></td></tr><tr><td>Pressure coefficient</td><td>ε<sub>P</sub></td><td>3.8 10<sup>-8</sup> Pa<sup>-1</sup></td></tr><tr><td>Concentration coefficient</td><td>ε<sub>χ</sub></td><td>2.9 10<sup>-4</sup> m<sup>3</sup> kg<sup>-1</sup></td></tr><tr><td>Density of eroded materials</td><td>ρ<sub>m</sub></td><td>1300 kg m<sup>-3</sup></td></tr><tr><td>Erosion coefficient</td><td>λ</td><td>2.5 10<sup>-4</sup> m<sup>-1</sup></td></tr><tr><td>Damage coefficient</td><td>β</td><td>9 10<sup>-3</sup> m<sup>3</sup> kg<sup>-1</sup></td></tr></tbody></table>
To calculate $\beta$, we used the experimental slope of the hydraulic head ($\Delta P / \Delta t$) with the integrated Darcy law (over [0, L], Equation 11).

$$ \frac{\Delta P}{\Delta t} = \frac{\mu_r \|\mathbf{q}\|^2 C_b \beta}{\|\mathbf{k}_0\|} [1 - \exp(-\lambda L)] \quad (11) $$

If $\lambda$ is determined first, $\beta$ (in $m^3 kg^{-1}$) is the only unknown of Equation 11. Overall, whereas $\lambda$ controls the amount of eroded particles, $\beta$ modulates the passing through of the eroding fluid.

Actually, as a result of granular filtration mechanisms, a fraction of dislodged particles does not make it through the sample. Sizes of constrictions within the pore network are too small and, therefore, prevent the coarse fluidized particles from being flushed out by seepage. Yet, since the lambda and beta parameters have been experimentally determined, they take into account (at least partly) the effects of granular filtration.

## 2.3. The numerical algorithm

The algorithm used to numerically solve the set of nonlinear and coupled equations is detailed in references [29,30]. The coupled problem is solved by a Gauss-Seidel type method for nonlinear partial differential equation systems as explained in reference [31]. In brief, the equations of flow and transport with viscosity and density effects are implemented in finite elements code. The problem is discretized in time by finite differences using a classical trapezoidal rule (or $\theta$-rule) and in space by finite element method. The discretized system is solved in sequences with an outer iteration. It consists of solving the initial system within an iterative procedure for each state variable separately, leading to a decoupled solution. The pressure equation is first solved with respect to $P$. The pressure equation is still nonlinear and is solved by successive substitutions. Then, the transport equation is solved with respect to $\chi$, using the most recently computed $p$ values. The attractive feature of the Gauss-Seidel iteration approach in nonlinear problem is that we only need to solve standard scalar partial differential equations. For instance, if $i$ denotes the iterative index, Equation 5 is first solved with respect to $P^{i,j+1}$ (where $j$ represents the time level) considering the other state variables at the previous iteration, $i-1$. Once $P^{i,j+1}$ is known, the transport and the erosion equations are solved with respect of $C^{i,j+1}$ and $\chi^{i,j+1}$, respectively. Iterations stop when a convergence criterion is fulfilled: that is, the erodible solid content of the material stabilizes ($|\partial C/\partial t| < 10^{-7} kg m^{-3} s^{-1}$).

# 3. INTERNAL EROSION EXPERIMENTS

Details about IE experiments and quality insurance are given in reference [28]. In brief, three IE tests were performed at a constant inflow rate (0.1 L min$^{-1}$, fluid velocities around $10^{-3} m s^{-1}$) flowing from the bottom to the top of a transparent tube of 0.1 m in length and 0.1 m in diameter filled with reinforced compacted Fontainebleau sand (Figure 1). The tests are different according to the reinforcement of the CSCs, that is, amounts of binding agent. During every test, the weight of the CSC, the injection pressure ($P_{inj}$, and hence average intrinsic permeability), the concentration and the granulometric distribution (from 0.01 to 100 µm) of effluent eroded particles (and hence average porosity) were measured at different time intervals.

The Fontainebleau sand is a reference granular material that consists of sub-rounded grains with intermediate sphericity and a grain size range of 50–400 µm. Its grain and bulk densities are ~2600 and 1500–1600 kg m$^{-3}$, respectively, with an internal friction angle of 40–45° and cohesion of 70 Pa [32]. A colloidal dispersion of high molecular polysilicic acids was used to reinforce the sand matrix. Indeed, preliminary IE tests done on unreinforced Fontainebleau sand columns highlighted the prompt destabilization of their granular structure. Because of its colloidal nature and its high water content (~70%), the dispersion allows the chemical reinforcement of low-permeability materials (< $10^{-10} m^2$). The colloid size lies between 5 and 75 nm, and the transformation of the colloidal silica dispersion into insoluble silica gel (the binder; ~1300 kg m$^{-3}$ grain density), upon addition of a saline solution, is irreversible. The gelation time was set to ~1 hour. After impregnation
<figure><img src="image_2.png" /><figcaption>Figure 1. Experimental setup for internal erosion tests.</figcaption></figure>

with the dispersion, and prior to gelation, sand columns were compacted using a dynamic impact method (sand piles: 0.01 m, dropped load of 2.5 N; drop height of 0.1 m). Finally, prepared samples were kept at 20 °C and 50% relative humidity for 1 week to gain strength and harden fully until IE tests. In the experiments, the amounts of colloidal silica dispersion were set for the silica gel to occupy 33, 50, and 66% of the porosity of the compacted Fontainebleau sand (n=0.38). We thus referred to the prepared CSCs as N₃₃, N₅₀, and N₆₆, respectively. The obtained (n;k) (as determined from the added amount of binder and water permeability measurements) for N₃₃, N₅₀, and N₆₆ equal to (0.25;10⁻¹³ m²), (0.19;5 10⁻¹⁴ m²), and (0.13;5 10⁻¹⁵ m²), respectively. The bulk density of the sand matrix was kept constant (close to 1600 kg m⁻³).

# 4. RESULTS AND DISCUSSION

## 4.1. *Experimental features*

Figures 2A, B, and C show the evolution of mass of N₃₃, N₅₀, and N₆₆ CSCs during IE tests, respectively. *t*=0 second denotes the outflowing of the effluent. Mass evolution curves developed up to three zones: an early mass intake zone, which increases with the amount of binder; a later zone of plateau; and a final mass loss zone. This observation is consistent with the initiation (loss of impervious function), continuation (creation of pathways for eroding fluid), and progression

<figure><img src="image_3.png" /><figcaption>Figure 2. Evolution of the mass of A) N₃₃, B) N₅₀, and C) N₆₆ compacted sand columns (CSCs) during internal erosion (IE) tests (time = 0 seconds denotes the outflowing of the effluent).</figcaption></figure>
(destabilization of the granular structure) phases for embankment dam filters [33]. The mass intake, which is observable only for $N_{50}$ and $N_{66}$ [average ± SD; $(5 ± 2) \times 10^{-7}$ and $(29 ± 18) \times 10^{-7}$ kg s$^{-1}$, respectively], proves the lasting infiltration of the injected water into the pore network. The plateau zone is developed particularly in $N_{50}$ [(85–130) $\times 10^3$ seconds] and to a lesser extent in $N_{66}$ [(30–50) then $\times 10^3$ seconds] and $N_{33}$ [(0–1) $\times 10^3$ seconds]. It supports the existence of an IE stage during which both mass intake via water infiltration and flushing of eroded particles by seepage were low or balanced each other. In the last zone, the mass of the reinforced CSCs decreased slower with increasing amounts of binder: $(8.3 ± 1.1) \times 10^{-7}$, $(0.8 ± 0.1) \times 10^{-7}$ and $(0.4 ± 0.2) \times 10^{-7}$ kg s$^{-1}$ for $N_{33}$, $N_{50}$, and $N_{66}$, respectively.

Figures 3A and B display the evolution of eroded volumes of 1- to 100-µm particles (as determined from the granulometry of effluent particles) and n, respectively. n of $N_{33}$, $N_{50}$, and $N_{66}$ CSCs exhibited analogous trends (Figure 3B): an initial period of rapid linear increase [$(12.3 ± 1.5, 5.6 ± 0.4$, and $45 ± 7) \times 10^{-7}$ s$^{-1}$, respectively], which, except for $N_{33}$, corresponds to the mass intake zone, followed by a slower increase in the plateau and mass loss zones [$(6.0 ± 0.7, 0.8 ± 0.1$, and $0.4 ± 0.3) \times 10^{-7}$ s$^{-1}$, respectively]. The rapid initial increase of n accounts for the filling of accessible pores with injected water and, as the volumes of 1- to 100-µm eroded particles indicate [$(32.2 ± 0.3, 1.4 ± 0.1$, and $0.8 ± 0.1) \times 10^{-10}$ m$^3$ s$^{-1}$, respectively] (Figure 3A), the flushing out of most movable particles (like loose and/or small particles that can easily exit the material). By comparison, in the plateau and mass loss zones, the slower increases of n likely account for slower filling of the less accessible pores and lower effluent loading with 1- to 100-µm particles [$(1.4 ± 0.1, 0.3 ± 0.1$, and $0.2 ± 0.1) \times 10^{-10}$ m$^3$ s$^{-1}$, respectively] and rigid skeleton constituents.

During all the IE tests, head losses between the bottom and top of the CSCs markedly decreased ($ΔP$ of 10, 107, and 506 kPa, respectively) (Figure 3C). This resulted in significant k increases (k:k$_0$ ratio values reached up to 40, 28, and 9, respectively) (Figure 3D). Assuming that infiltration of injected water and removal of silica gel lead to local drain of the finest sand grains and further destabilization of CSCs structures, we used exponential regression models to fit $P_{inj}$ and k evolution: $P_{inj}(Pa) = 3 \times 10^7 e^{-3n}$ ($r^2$ is

<figure><img src="image_4.png" /><figcaption>Figure 3. Evolution of the following: A) eroded volume of 1- to 100-µm size particles, B) average porosity, C) injection pressure, and D) average intrinsic permeability of $N_{33}$, $N_{50}$, and $N_{66}$ eroded CSCs, respectively.</figcaption></figure>
the correlation coefficient = 0.87) and k(m²) = 6 × 10⁻¹⁶ e³⁰n (r² = 0.81) (p-value = 0.05). Because intrinsic permeability depends on the geometry of constrictions and interconnections between pores, the obtained relationships suggest that the filling of pores and subsequent drain of fluidized particles could have initiated (in terms of size or quantities) pathways and/or capillaries for water flowing throughout the material. This also is consistent with the previous assumption stating that slight porosity increase may result in marked permeability changes (see Section 2.1).

## 4.2. Internal erosion simulations

Here, we present IE simulations for N₅₀ test, that is, a CSC that displayed delayed occurrence of significant amounts of eroded particles into the effluent (Figure 3A). Simulations were restricted to the mass loss stage starting 130 10³ seconds after the outflowing of the effluent (Figure 2B). During this IE stage, the less accessible pores are thought to be filled with fluid, and most of the movable particles flushed out (as the experimental data indicated). The numerical data show that C (resp. n and k) increases (resp. decrease) slower with distance from injection point (located at x=0 m) and decreases (resp. increase) linearly with increasing IE duration (as observed from changes in the shape and height of curves in Figures 4A, C, and D). This indicates that particle dislodgement and ensuing formation of pores and opening of constrictions first occurred near the input of the CSC. Unlike for C, the spatial distribution of fluidized particles exhibits a marked nonlinear behaviour: χ increases quasi-exponential with the distance from the injection point (Figure 4B). This indicates that the fluidized particles are readily transferred toward the output of the CSC and concentrate as the fluid flows throughout the material. Thus, depending on the material (resp. on λ and β values),

<figure><img src="image_5.png" /><figcaption>Figure 4. Simulations at different times of A) content of erodable solids, B) concentration of fluidized particles at t=(10, 500, and 1000) 10³ seconds, C) porosity, and D) intrinsic permeability of N₅₀ CSC against the distance from the injection point. Because internal erosion tests did not exceed 2.5 10⁵ seconds, curves at t=5 10⁵ and 10⁶ seconds could not be compared with experimental results (time=0 second denotes the beginning of the mass loss stage).</figcaption></figure>
the shape of $\chi$ profile may notably differ from that of C. As a matter of fact, we calculated that > 99% of the fluidized particles concentrate in the last 0.02 m of the CSC, and this proportion increases with time. The numerical simulations agree reasonably well with ongoing IE mechanisms via the formation of pathways for fluid (as *n* and *k* profiles indicated) that propagate forward (as can be observed from changes in the shape of C and $\chi$ profiles).

Simulated data for $N_{50}$ mass losses, injection pressure, *n*, and *k* have been compared with experimental results. The resulting curves are presented in Figures 5A, B, C, and D, respectively. Simulations exhibit curves close to experimental data ($r^2 > 0.6, p=0.05$), thereby indicating that the model overall reproduces well the IE effects. Yet, discrepancies (up to 30 %) for mass loss, *n*, and *k* usually occur at late IE stages when the CSCs become less stable. Indeed, because the filtration of fluidized particles has not explicitly been taken into account (see Section 2.2), the model tends to overestimate the experimental values. Furthermore, local shifts (like at $t=33 \times 10^3$ seconds), which correspond to swift discharges of fine Fontainebleau sand grains [28], were neglected and thus consist another limit to the model. Hence, the proposed model only provides estimates of the extent and dynamics of local structural changes (like *n*, *p* increase) that occur during IE. Such information is still important. By numerically evaluating $\chi$ of effluent waters and/or *n* and *k* values of the superficial layer of close to failing internally eroded materials (as tested in the laboratory), obtained values can be compared with *in situ* measurements to estimate an IE time for the hydraulic work to reach failure.

### 4.3. Exploiting the outcomes of the model

To assess the spatial distribution of structural changes and seepage stresses, we considered a unit volume of eroded material with many parallel cylindrical capillaries of mean hydraulic diameter $D_H$.

<figure><img src="image_6.png" /><figcaption>Figure 5. Comparisons between the evolution of modeled and experimental A) mass loss, B) injection pressure, C) average porosity, and D) average permeability of $N_{50}$ CSC (time = 0 second denotes the beginning of the mass loss stage).</figcaption></figure>
crossing it. Using the Hagen-Poiseuille equation with the Darcy's law (the flow in capillaries is laminar, Re = 0.013 ± 0.001) as well as modeled n and k values, D<sub>H</sub> is calculated as follows:

$$D_H^2 = \frac{32k}{n} \qquad (12)$$

The obtained D<sub>H</sub> values (average of 13.7 ± 0.2 µm) are consistent with sizes of eroded silica gel fragments (average diameter 11 ± 5 µm) and smallest constriction based on the granulometry of the effluent particles (around 16 µm) [28]. D<sub>H</sub> slowly decreases with the distance from the injection point, and when averaged over all the sample length (D̄<sub>H</sub> = L<sup>-1</sup> ∫<sub>x=0</sub><sup>x=L</sup> D<sub>H</sub>(x)dx), it increases quasi-linearly with IE duration [r<sup>2</sup> = 0.98, p = 0.05; D̄<sub>H</sub>(µm) = 10<sup>-5</sup> t (s) + 13.4] (Figure 6A). This result confirms that fluidized particles typically consist of silica gel fragments, which mobility (in terms of D<sub>H</sub>) decreases toward the output of the sample, whereas it increases with IE duration. Based on D<sub>H</sub> data, the Hagen-Poiseuille equation for laminar flow was used to estimate the mean flow velocity, n, of the fluid within a cylindrical capillary:

$$n = \frac{1}{2} n_{\max} = \frac{D_H^2}{32\mu_f} \nabla P \qquad (13)$$

∇P is the hydraulic gradient. Calculated [(9.31–10.05) 10<sup>-4</sup> m s<sup>-1</sup>] and experimental [(9.37–9.83) 10<sup>-4</sup> m s<sup>-1</sup>] as determined from porosity measurements and the flow of discharging effluent] ||n|| values are similar. This supports the assumption that cylinders fit the overall shape of capillaries. n is essentially driven by ∇P: ||n|| decreases slower with distance from the injection point and quasi-linearly with IE duration when averaged over all the sample length [p = 0.05, r<sup>2</sup> = 0.99; ||n|| (m s<sup>-1</sup>) = -1.26 10<sup>-10</sup> t (s) + 9.72 10<sup>-4</sup>] (Figure 6B). From the mean flow velocity, the mean shear stress (t) incurred by capillary boundary was calculated using the formulation below:

$$\mathbf{t} = \mu_f \frac{\partial n}{\partial x} \qquad (14)$$

It results that ||t|| and ||t|| display quasi-linear decreases with distance from the injection point and IE duration [p = 0.05, r<sup>2</sup> = 0.99; ||t|| (Pa) = -2.24 10<sup>-12</sup> t (s) + 1.17 10<sup>-6</sup>], respectively (Figure 6C). This implies that IE efforts and particle dislodgement (see Section 4.2) primarily concentrate in the vicinity of the water input. From an engineering point of view, this zone would eventually constitute an unstable portion of the tested material. It is noteworthy that because of the filtration mechanisms and CSC heterogeneity, obtained ||t|| values likely underestimate local shear stresses that actually developed in N<sub>50</sub>. Failure of similarly prepared then eroded CSCs usually occurred at > 1.5 10<sup>5</sup> seconds from the beginning of the mass loss stage. According to the model calculations, this corresponds to > 30% decrease (resp. > 5% and > 30%

<figure><img src="image_7.png" /><figcaption>Figure 6. Simulations at different times of A) mean capillary hydraulic diameter, B) mean flow velocity in capillaries, and C) mean shear stress incurred by capillary boundary in N<sub>50</sub> CSC versus the distance from the injection point (time = 0 second denotes the beginning of the mass loss stage).</figcaption></figure>
increases) of $||\bar{\tau}||$ (resp. of $n$ and $k$) (Figures 6C, 4C, and 4D). It is therefore assumed that the decrease in the amplitude of seepage stresses is not fast enough to prevent further destabilization of $N_{50}$ CSC, that is, because of the removal of binder and fine sand grains that maintain the fragile coarse sand grains matrix still. In other words, the kinetic of silica gel removal and ensuing swift releases of fines prevailed over the diminution of the amplitude of seepage stresses. Hence, we question the fact that similarly reinforced soils or man-made infrastructures can undergo long-term development of pathways for passing-through water.

## 5. CONCLUDING REMARKS

The aim of this paper is to improve our knowledge on the transport of particle-loaded fluids in chemically reinforced granular materials affected by IE. It discusses the evolution of key parameters for chemically reinforced materials during IE. As a matter of fact, significant variations in the pore pressure, concentration of fluidized particles, and content of erodible solids were observed. Thus, it demonstrates that IE mechanisms modify both the structure of the rigid skeleton and pore network, for example, through porosity and permeability increases.

At a constant inflow rate, IE tests performed on reinforced CSCs highlighted the existence of up to three IE stages: (i) an early CSCs mass intake stage, which corresponds to filling of accessible pores with injected water and flushing of the most movable particles; (ii) a subsequent stage of low binder erodability coupled with and lower effluent loading with fines; and (iii) a later CSCs mass loss stage of silica gel removal and ensuing swift discharges of fines. Overall, the observed correlation between $n$ and $k$ indicated that the filling of pores and the subsequent drain of fluidized particles could have initiated pathways for water flowing throughout the material.

The mathematical model reproduces well the effects of IE. It allowed following the evolution and spatial distribution of structural changes and seepage stresses. The simulations indicated that, during the mass loss stage, dislodgement of particles from the rigid skeleton and channelization occurred mostly near the water injection point where the shear stresses were high. Furthermore, the mobility of fluidized particle was shown to decrease toward the output of the sample. This is consistent with cavities and/or capillaries that are eroded toward the output of the CSCs. One limit of the model is not explicitly taking into account granular filtration mechanisms. It thus eventually overestimates the actual mass losses that are the sum of dynamic releases and filtering of silica gel fragments and sand grains. Hence, the model only provides estimates of the extent and dynamics of the spatial distribution of structural changes occurring during IE. Based upon IE tests, such values can be used to estimate an IE time for the hydraulic work to reach failure. Ongoing research are led in our laboratory to take into consideration granular filtration mechanisms.

## ACKNOWLEDGEMENTS

The research work described herein was led in the framework of the POLESUR project (Federal multi-disciplinary research on the pollution of soil particles and water in underground aqueous medium). It was financially supported by the 'Pays de la Loire (PDL)' Region and the French National Research Agency (ANR) through the Erinoh Project.

## REFERENCES

1. Foster M, Fell R, Spannagle M. Statistics of embankment dam failures and accidents. *Canadian Geotechnical Journal* 2000; **37**:1000-1024.

2. Fry J-J. Internal Erosion: Typology, Detection, Repair. *Barrages and Reservoirs* N° 6. Comité Francais des Grands Barrages, Le Bourget-du-lac Cedex, 1997.

3. Bonelli S, Brivois O. The scaling law in the hole erosion test with a constant pressure drop. *International Journal for Numerical and Analytical Methods in Geomechanics* 2008; **32**:1573-1595.

4. US Army Corps of Engineers. Filter experiments and design criteria. *Technical memorandum of the U.S. Army. Vicksburg, 1953; 3-360*.

5. Loebotsjkov EA. The calculation of suffusion properties of nocohesive soils when using non-suffosion analog (in Russia). *Publication of the technical university of Brno*, Svazek, Czechoslovakia 1969; 135-148.

6. Arulanandan K, Gillogley E, Tully R. Development of a quantitative method to predict critical shear stress and rate of erosion of natural undisturbed cohesive soils. *Technical report of the U.S. Army*. Washington DC, 1980; GL-80-5,
7. Arulanandan K, Perry EB. Erosion in relation to filter design criteria in earth dams. *Journal of Geotechnical Engineering* 1983; **109**:682–698.

8. Sterpi D. Effects of the erosion and transport of fine particles due to seepage flow. *The International Journal of Geomechanics* 2003; **3**:111–122.

9. Wan CF, Fell R. Investigation of internal erosion and piping of soils in embankment dams by the slot erosion test and the hole erosion test. *UNICIV Technical report of The University of New South Wales*, Sydney, Australia 2002; No. R-412.

10. Wan CF, Fell R. Investigation of rate of erosion of soils in embankment dams. *Journal of Geotechnical & Geoenvironmental Engineering* 2004; **130**:373–380.

11. Bendahmane F, Marot D, Alexis A. Experimental parametric study of suffusion and backward erosion. *Journal of Geotechnical & Geoenvironmental Engineering* 2008; **134**:57–67.

12. Bonelli S, Marot D. On the modelling of internal soil erosion. *The 12 International conference of International Association for Computer Methods and Advances in Geomechanics*, Goa, India, 2008; 2544–2550.

13. Tam JG. Grout expectations. Description of and applications for grouts used in the civil engineering field. *Unpublished report*. (Available from: http://www.ocf.berkeley.edu/~tamjo/grout/Grout-rep.html, 01/03/2011) 1–22.

14. Vardoulakis I, Stavropoulou M, Papanastasiou P. Hydromechanical aspects of sand production problem. *Transport in Porous Media* 1996; **22**:225–244.

15. Vardoulakis I, Papanastasiou P, Stavropoulou M. Sand erosion in axial flow conditions. *Transport in Porous Media* 2001; **45**:267–281.

16. 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**:789–808.

17. Wan RG, Wang J. Modelling sand production within a continuum mechanics framework. *Journal of Canadian Petroleum Technology* 2002; **41**:46–52.

18. Papamichos E, Vardoulakis I. Sand erosion with a porosity diffusion law. *Computers and Geotechnics* 2005; **32**:47–58.

19. Muir-Wood D, Maeda K, Nukudani E. Discrete element modelling of soil erosion. *Proceedings of the 4th International Conference on Scour and Erosion*, Tokyo, Japan, vol. **29**, 2008; 47–65.

20. Fujisawa K, Murakami A, Nishimura S. Numerical analysis of the erosion and the transport of fine particles within soils leading to the piping phenomenon. *International Symposium on Prediction and Simulation Methods for Mitigation of Geohazard*, Kyoto, Japan, vol. **50**, 2010; 471–482.

21. Sibille L, Scholtès L. Effects of internal erosion on mechanical behavior of soils: a DEM approach. In *Proceedings of the International Conference on Particle-Based Methods*, Onate E, Owen DRJ (eds). CIMNE Publisher: Barcelona, Spain, 2009; 167–170.

22. Bear J, Bachmat Y. *Introduction to modelling phenomena of transport in porous media*. Kluwer Academic Publishers: Dordrecht, 1991; 554.

23. Haverkamp R, Vauclin M, Touma J, Wierenga PJ, Vachaud G. A comparison of numerical simulation models for one dimensional infiltration. *Soil Science Society of America Journal* 1977; **41**:285–294.

24. Hassanizadeh M, Gray WG. General conservation equations for multi-phase systems: 1. averaging procedure. *Advances in Water Resources* 1979; **2**:131–144.

25. Haggerty R, Gorelick SM. Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity. *Water Resources Research* 1995; **31**:2383–2400.

26. Baeumer B, Meerschaert MM, Benson DA, Wheatcraft SW. Subordinated advection-dispersion equation for contaminant transport. *Water Resources Research* 2001; **37**:1543–1550.

27. Bear J. *Dynamics of fluids in porous media*. Elsevier Publications: McGraw-Hill, NY, 1972; 768.

28. Muresan N, Saiyouri N, Guefrech A, Hicher P-Y. Internal erosion of chemically reinforced granular materials: a granulometric approach. *Journal of Hydrology* 2011, DOI:10.1016/j.jhydrol.2011.09.009

29. Chupin O, Saiyouri N, Hicher P-Y. The effects of filtration on the injection of cement-based grouts in sand columns. *Transport in Porous Media* 2007; **72**:227–240.

30. Kolditz O, Ratke R, Diersch H-JG, Zielke W. Coupled groundwater flow and transport: 1. verification of variable density flow and transport models. *Advances in Water Resources* 1998; **21**:27–46.

31. Langtangen HP. *Computational Partial Differential Equations: Numerical Methods and Diffpack Programming*. Springer-Verlag: Berlin, 1999; 855.

32. van Mechelen JLM. Strength of moist sand controlled by surface tension for tectonic analogue modelling. *Tectono-physics* 2004; **384**:275–284.

33. Foster MA, Fell R. Assessing Embankement Dam Filters that do not satisfy Design Criteria. *Journal of Geotechnical and Geoenvironmental Engineering* 2001; **127**:398–407.