Author

PhD Jie Yang
email: doc.jie.yang@gmail.com

Zhen-Yu Yin
email: zhenyu.yin@polyu.edu.hk

PhD Farid Laouafa
email: farid.laouafa@ineris.fr

Emeritus Pierre-Yves Hicher
email: pierre-yves.hicher@ec-nantes.fr

Hydro-mechanical modeling of granular soils considering internal erosion

Keywords: Granular soils, internal erosion, hydro-mechanical coupling, critical state, fines content

This paper attempts to formulate a coupled practical model in the framework of continuum mechanics to evaluate the phenomenon of internal erosion and its consequences on the mechanical behavior of soils. For this purpose, a four-constituent numerical approach has been developed to describe the internal erosion process. The detachment and transport of the fine particles have been described by a mass exchange formulation between the solid and fluid phases. The stress-strain relationship of the soil is represented by a non-linear incremental model. Based on experimental data, this constitutive model has been enhanced by the introduction of a fines content dependent critical state, which allows accounting for the influence of fines on soil deformation and strength. The applicability of the practical approach to capture the main features of the internal erosion process and its impact on the mechanical behavior of the eroded soil has been validated by comparing numerical and experimental results of internal erosion tests on HK-CDG mixtures, which demonstrated that the practical model was able to reproduce, with a reasonable success, the experimental data.

Furthermore, the influence of the stress state, the initial soil density, and the initial fraction of fines have been analyzed through numerical simulations using the proposed model.

Introduction

Internal erosion occurs when fine particles are plucked off by hydraulic forces and transported through the coarse matrix. The known causes are either a concentration of leak erosion, backward erosion, soil contact erosion, or suffusion [START_REF] Bonelli | On the modelling of internal soil erosion[END_REF][START_REF] Fell | Time for development of internal erosion and piping in embankment dams[END_REF][START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF]). This paper focuses on suffusion, which corresponds to the detachment and migration of fine particles within the voids of coarse particles by seepage flow. Subsequently, the induced modification of the soil microstructure may lead to deformations at the macroscopic scale and, therefore, influence significantly the mechanical behavior of the soil. Many damages and failures of embankment dams can be associated with internal erosion [START_REF] Fell | Time for development of internal erosion and piping in embankment dams[END_REF][START_REF] Xu | Breaching parameters for earth and rockfill dams[END_REF][START_REF] Zhang | Seepage failure mechanism of the Gouhou rockfill dam during reservoir water infiltration[END_REF][START_REF] Zhang | Analysis of earth dam failures: A database approach[END_REF][START_REF] Rönnqvist | On the use of empirical methods for assessment of filters in embankment dams[END_REF]. [START_REF] Foster | The statistics of embankment dam failures and accidents[END_REF] showed that nearly 46% of the damages of 128 embankment dams could be attributed to internal erosion.

Besides, sinkholes and cavities triggered by internal erosion are frequently observed within dams and dikes [START_REF] Muir Wood | The magic of sands-the 20th Bjerrum Lecture presented in Oslo, 25 November 2005[END_REF][START_REF] Sterpi | Effects of the erosion and transport of fine particles due to seepage flow[END_REF]. Finally, internal erosion can also trigger slope failures and generate significant landslides [START_REF] Crosta | On slope instability induced by seepage erosion[END_REF][START_REF] Hu | Seismic precursor to instability induced by internal erosion in loose granular slopes[END_REF].

To evaluate the internal instability of filters in embankment dam, engineers have traditionally employed empirical methods [START_REF] Burenkova | Assessment of suffusion in non-cohesive and graded soils[END_REF][START_REF] Foster | Assessing embankment dam filters that do not satisfy design criteria[END_REF]Kenney andLau 1985, 1986;[START_REF] Kézdi | Soil physics-selected topics-developments in geotechnical engineering-25[END_REF][START_REF] Li | Comparison of two criteria for internal stability of granular soil[END_REF][START_REF] Sherard | Sinkholes in dams of coarse, broadly graded soils. 13th[END_REF][START_REF] Wan | Assessing the potential of internal instability and suffusion in embankment dams and their foundations[END_REF]. [START_REF] Rönnqvist | On the use of empirical methods for assessment of filters in embankment dams[END_REF] compiled a database of 80 dams of which 23 dams suffered from documented internal erosion. Their gradation analysis demonstrated that empirical methods are capable of distinguishing dams with a probable occurrence of erosion from dams with no up-to-date indication of internal erosion.

Laboratory tests further investigated the mechanics of internal erosion and the susceptibility of soil to seepage induced internal instability by using primarily the conventional rigid-walled permeameter under self-weight conditions (e.g. [START_REF] Kenney | Internal stability of granular filters[END_REF][START_REF] Skempton | Experiments on piping in sandy gravels[END_REF][START_REF] Sterpi | Effects of the erosion and transport of fine particles due to seepage flow[END_REF], Wan andFell 2008, Sibille et al. 2015a;[START_REF] Rochim | Effects of Hydraulic Loading History on Suffusion Susceptibility of Cohesionless Soils[END_REF]. The eroded mass, the hydraulic gradient and the flow rate were monitored to determine the susceptibility of seepage induced instability by incorporating the post-erosion analysis of the particle size distribution. The initial fines content and the history of the hydraulic loading appeared to play an important role in the development of internal erosion (Sibille et al. 2015a;[START_REF] Rochim | Effects of Hydraulic Loading History on Suffusion Susceptibility of Cohesionless Soils[END_REF]. It is important to add that the soils constituting embankment dams or levees are usually subjected to complex stress states [START_REF] Zhang | Effects of abutment slopes on the performance of high rockfill dams[END_REF]. To address this phenomenon, Moffat andFannin (2006, 2011) developed a rigid-wall permeameter with an axial loading system so that the onset of internal instability of cohesionless soils under stress state could be investigated. More recently, advanced flexible wall 0 K permeameters capable of controlling the stress on the tested specimen have been developed so that the development of internal erosion under different stress states and the mechanical behavior of soils during and after erosion could be investigated (e.g., [START_REF] Bendahmane | Experimental parametric study of suffusion and backward erosion[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF][START_REF] Marot | Assessing the susceptibility of gap-graded soils to internal erosion: proposition of a new experimental methodology[END_REF]Ke and Takahashi 2014a;[START_REF] Slangen | A Flexible Wall Permeameter for Investigating Suffusion and Suffosion[END_REF]. Most of the studies were conducted under an isotropic stress state. The results showed that the larger effective confining pressure led to a less extent of internal erosion [START_REF] Bendahmane | Experimental parametric study of suffusion and backward erosion[END_REF][START_REF] Mehdizadeh | Progressive internal erosion in a gap-graded internally unstable soil: mechanical and geometrical effects[END_REF], even if the deformation of the specimen caused by erosion was limited. [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF] showed that, by performing internal erosion tests under different anisotropic stress states, the maximum erosion rate, the variations in soil permeability, and the erosion induced deformation of the specimen will all increase with the increase of the deviatoric stress. For the purpose of investigating the influence of the loss of fine particles on the mechanical behavior of eroded soils, they performed triaxial tests on eroded specimen under different stress paths. Drained triaxial tests on eroded specimens indicated that the initially dilative stress-strain behavior changed to a contractive behavior and the soil strength decreased under drained shearing [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF]Ke and Takahashi 2014a, 2014b, 2015). The undrained behavior of the original specimen changed from a strain hardening behavior to a flow-type behavior with limited deformation after internal erosion [START_REF] Ke | Soil deformation due to suffusion and its consequences on undrained behavior under various confining pressures[END_REF][START_REF] Mehdizadeh | Progressive internal erosion in a gap-graded internally unstable soil: mechanical and geometrical effects[END_REF].

Parallel to experimental studies, numerical approaches have also been proposed to describe the erosion process, either within the framework of continuum mechanics by using finite element or finite difference methods [START_REF] Cividini | Finite-element approach to the erosion and transport of fine particles in granular soils[END_REF][START_REF] Papamichos | Volumetric sand production model and experiment[END_REF][START_REF] Stavropoulou | Coupled wellbore erosion and stability analysis[END_REF][START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF][START_REF] Zhang | A thermodynamicsbased model on the internal erosion of earth structures[END_REF], or within a discrete framework by associating the discrete element method to describe the dynamics of the granular solid phase to the computational fluid dynamics to solve the interstitial fluid flow [START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: application to piping erosion[END_REF][START_REF] Mansouri | Numerical simulation of the quicksand phenomenon by a 3D coupled Discrete Element -Lattice Boltzmann hydromechanical model[END_REF][START_REF] Reboul | Transport de particules dans les milieux granulaires: Application à l'érosion interne[END_REF][START_REF] Sari | Investigation of internal erosion processes using a coupled dem-fluid method[END_REF]Sibille et al. 2015b;[START_REF] Zhang | Three-dimensional DEM modeling of the stress-strain behavior for the gap-graded soils subjected to internal erosion[END_REF][START_REF] Zhao | Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics[END_REF], mainly focused on the hydraulic process, i.e. the evolution of the eroded mass and permeability. The mechanical part was simplified by using an elastic model with induced damage parameters. Some attempts have also been made to model the mechanical response to internal erosion by removing particles in granular materials. [START_REF] Muir Wood | Modelling mechanical consequences of erosion[END_REF] modeled the mechanical consequences of internal erosion by a two-dimensional discrete element analysis and found that internal erosion changes the density state of the material from dense to loose as they compared the current soil state to the critical state line. [START_REF] Scholtès | Multiscale approaches to describe mechanical responses induced by particle removal in granular materials[END_REF] developed a three-dimensional discrete element model and an analytical micromechanical model to describe the mechanical responses induced by particle removal in granular materials. They found that the mechanical behavior of the soil changes from dilative to contractive with the removal of soil particles. However, this removal was based on the size of the particles and their degree of interlocking, and the fluid phase was not considered. [START_REF] Hosn | A discrete numerical model involving partial fluid-solid coupling to describe suffusion effects in soils[END_REF] studied the macroscopic response during internal erosion and the post-internal erosion properties of granular materials by particle removal, considering a one-way coupling between the fluid and the solid phases via the discrete element and the pore scale finite volume methods. The erosion process was simplified into the definition of the detached and transportable particles by calculating the unbalance fluid force and by defining a controlling constriction size. [START_REF] Hicher | Modelling the impact of particle removal on granular material behaviour[END_REF] presented a numerical method to predict the mechanical behavior of granular materials subjected to particle removal and concluded that the removal of soil particles may cause diffuse failure in eroded soil masses. However, the process of real erosion was not considered. Therefore, it is of interest to develop a numerical method that will provide the modeling of the hydro-mechanical response of soils subjected to internal erosion.

Despite the relevance of the discrete element method to shed light onto the constitutive properties of granular materials, its computational cost remains a restrictive aspect which could become an obstacle wherever real soil structures such as embankment dams are concerned. Hence, it is necessary to develop the approaches that can take better into account the engineering boundary value problems. The objective of this paper is to formulate a numerical model within the framework of continuum mechanics to investigate the coupled hydro-mechanical behavior of soils subjected to internal erosion. First, an enhanced four-constituent model has been developed to describe the detachment and transport of fine particles induced by the fluid flow with more realistic consideration of the effects of the fines content and the shear stress level. This process has been coupled with an advanced critical state-based mechanical model, which has been enhanced to consider the effects of the evolution of the fines content. By assuming a homogenous sample, a practical simulation tool for triaxial erosion tests under mechanical loadings has also been proposed. Then, the predictive ability of the approach has been examined by simulating internal erosion tests on HK-CDG mixtures under mechanical loadings. Finally, the influence of the stress state, the initial density, and the initial fines content of the soil have been investigated by the developed numerical model.

Coupled hydro-mechanical model

According to [START_REF] Schaufler | Infiltration processes in cohesionless soils[END_REF], it is possible to consider the saturated porous medium as a material system composed of 4 constituents in 2 phases: the stable fabric of the solid skeleton, the erodible fines, the fluidized particles and the pure fluid. The fines can behave either as a fluid-like (described as fluidized particles) or as a solid-like (described as erodible fines) material.

It is assumed that the seepage velocities is given as the relative motion of the fluid phase with w l v respect to the solid phase . The relative motion of the pore fluid and the fluidized particles is s v neglected. The total discharge of the liquid is then given as:

w q (1)         , , , with ,     w l s x t x t x t x t q w w v v
where is the porosity. Furthermore, in the framework of small transformations (small strains   , x t  and small displacements), the strain of the solid skeleton is given by the symmetric part of the displacement gradient.

Mass balances and transported particles

The balance equations are derived based on the porous medium theory (de Boer 2000; [START_REF] Steeb | A thermodynamic-consistent model describing growth and remodeling phenomena[END_REF][START_REF] Uzuoka | Dynamics of unsaturated poroelastic solids at finite strain[END_REF]. The local form of the volume balance of a constituent is deduced from the partial balance of mass [START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF], given as:

(2)

    div i i i n n n t     v
where and denote, respectively, the volume fraction and the velocity of the corresponding

i n i v constituent.
The liquid-solid phase transition process has been accounted for by introducing a volumetric mass exchange term in the corresponding mass balances for the 4 constituents. More n details can be found in the Appendix A1.

A system of non-linear partial differential equations is deduced from the mass balance of the mixture system:

(3)

    div div 0 w s   q v (4)     div div s s n t         v v (5)         div div c c c s c s f f f f n t t           v v (6)       div s c c c n t t           w v q
where , and are the porosity, the amount of erodible fines and the  

, x t    , c f x t   , c x t
concentration of the fluidized particles, respectively; is the source term describing the exchange n between the erodible fines and the fluidized particles; is the velocity of the soil skeleton;

s v   , w x t q
is the total discharge of the pore fluid.

Balance of linear momentum

The model has been developed for geotechnical applications and it is not the aim of the present paper to present the thermodynamical setting of the erosion model. Therefore, the linear momentum balance equation of a constituent is given as follows under isothermal and quasi-static conditions: [START_REF] Borja | On the mechanical energy and effective stress in saturated and unsaturated porous continua[END_REF]   ˆdiv 0, 0

i i i i     σ p

Porous fluid flow

In this study, it is assumed that the flow in the porous medium is governed by Darcy's law which states that the flow rate is driven by the gradient of the pore fluid pressure :

  , w p x t (10)       , c w k k f p c      w q grad
where denotes the intrinsic permeability of the medium; is the kinematic viscosity of the

  , c k f  k  fluid;
is the pore fluid pressure; and is the density of the mixture defined as:

w p   c  (11)   1 s f c c      
with the density of the solid and the density of the fluid. For a material made of a mixture of s  f  coarse and fine grains, the intrinsic permeability of the porous medium depends on the   , k x t current porosity and the fines content as [START_REF] Revil | Permeability of shaly sands[END_REF]:

  , x t    , c f x t 10 / 35 (12)   3 0 1 1 k m c k k f        
where is a power index which varies with the pore geometry [START_REF] Revil | Permeability of shaly sands[END_REF]. k m

Mass exchange term

The variable in Eqs. ( 4)-( 6) is the volumetric mass exchange at any point in time, given by the n following relation (Uzuoka et al. 2012):

  

ˆ1 e c c n f f        w q
where is the ultimate fines content after a long seepage period; is a material parameter. The is assumed to decrease when the hydraulic gradient increases as follows (Cividini et

  c c f f   c f  al. 2009): (14)       3 2 1 1 0 1 exp 10 exp exp c c p p f f M M                             w q
where is the initial fines content, , and are material parameters. A mechanical

0 c f 1  2  3 
component " " defined later on is introduced into the definition of , in order to take into account  c f  the influence of the stress state on the erosion rate.

Note that the mass exchange term is defined as a function of the porosity, the initial fines content, n the flow rate and the stress state. There isn't any explicit expression to define the critical hydraulic gradient required to initiate the erosion or the soil skeleton deformation. All the unknown variables are coupled implicitly and solved by the system of the partial differential equations.

Critical state based soil model considering fines content effect

A non-linear incremental constitutive model including the critical state concept has been adopted (Yin et al. 2018a). The non-linear frictional behavior of the granular materials can be described by an asymptotic relation between the stress ratio and the shear strain as follows:

 d  (15) 1 d p G M p M e              
where is the maximum strength ratio; is the initial stiffness; is the effective

p M 3 G G p    p pressure;
is the shear modulus. The model has been enhanced by introducing a non-linear G elasticity, a non-linear stress dilatancy, and the critical state concept [START_REF] Jin | Estimation of critical state-related formula in advanced constitutive modeling of granular material[END_REF][START_REF] Jin | Selection of sand models and identification of parameters using an enhanced genetic algorithm[END_REF]Jin and Yin 2018). More details of the model can be found in Yin et al. (2018a).

The constitutive equations under axisymmetric triaxial condition are briefly summarized as follows:

(16)

 

1 exp                            v d pt d q p K A M d p (17) 3 1 d p q q q G p p M p                  with (18)         2 2 0 0 2.97 2.97 , 1 1                     n n at at e e p p G G K K e p e p (19) 6sin 6sin , 3 sin 3 sin p pt p pt p pt M M        
where and are determined by

p  pt  (20) tan tan , tan tan                    p d n n c c p c pt c e e e e
with the critical state void ratio given by c e (21)

0 c cr at p e e p           
where the stress ratio under axisymmetric triaxial condition is defined by the deviatoric = ' q p  stress (

) and the mean effective stress ( ), with and the

q a r q     ' p ' 2 3 ( ) a r p       a  r 
axial and radial stress, respectively; the deviatoric strain increment is a

  2 3 d a r     
function of the axial and radial strain increments; is the volumetric strain 2 v a r

  

  increment; and are the bulk and shear moduli with the reference bulk and shear moduli and

K G 0 K ; Poisson's ratio
; the non-linear elastic constant (typically 0.5-

0 G     0 0 0 0 3 2 6 2     K G K G n 0.7);
is the atmospheric pressure chosen as the reference stress; and are the 101.325 Since the fines content evolves during internal erosion, the constitutive model needs to consider such internal soil mass evolution. According to [START_REF] Yin | A micromechanics-based model for sand-silt mixtures[END_REF]2016), is a function of the evolutive 0 cr e fines content , in order to unify the mechanical behavior of a sand-silt mixture for different fines

c f contents: (22)             0 , 0 , 0 1 tanh 1 tanh 1 1 2 2 c th c th c cr hc cr c c hf cr c m d f f f f f e e f af e f R                        
where and are the reference critical void ratios for the pure sand and the pure silt, obtainable at low experimental cost, as suggested by [START_REF] Yin | A micromechanics-based model for sand-silt mixtures[END_REF]2016). It should be pointed out that Eq. ( 22) is mathematically structured based on the hyperbolic tangent function to ensure a continuous and smooth transition domain; more validation and discussion can be found in [START_REF] Yin | A micromechanics-based model for sand-silt mixtures[END_REF]2016).

Hydro-mechanical responses of soils during internal erosion

Two main influences can be noted concerning the mechanical impact of the decrease of the fines content induced by internal erosion. First, according to Eq.( 18), the increase of the void ratio can reduce the elastic moduli and , leading to increasing deformation of the granular assembly. The G K second consequence is linked to the interlocking features. The sliding resistance between particles has two origins: one is the friction between two particles, the other is related to the geometrical (topological) effect of the neighboring particles which will resist the sliding displacement of two particles in contact, i.e. the interlocking effect. From a macroscopic standpoint, the first aspect is characterized by the internal friction angle , corresponding to the critical state; the second aspect is c  related to the density state, which will in turn influence the peak strength and the dilative/contractive behavior of the soil. In the model, this interlocking influence is taken into account by Eqs. ( 20)-( 22 (2014,2016), in the case of a fine grain skeleton, the loss of fine particles will decrease the critical void ratio . From a microscopic standpoint, the loss of fine particles makes it easier for the force c e chains which control the soil stability to buckle. As observed by [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF], the dilative behavior of the non-eroded specimen changed to a contractive behavior after erosion. In the case of a coarse grain skeleton, the loss of fine particles will make the critical state line move upwards in the e-logp' plane. The interlocking effect could, therefore, be enhanced when c e e increases. This might explain the observed increase of the initial peak strength of the eroded specimen in some cases (e.g. [START_REF] Mehdizadeh | Progressive internal erosion in a gap-graded internally unstable soil: mechanical and geometrical effects[END_REF]. Thus, when a part of the fine particles is detached from the soil skeleton, the void ratio increases, and the critical void ratio evolves e c e according to the fines content , creating a disequilibrium of internal and external forces and, c f consequently, a deformation of the granular assembly.

During an internal erosion experimental test, the soil specimen is subjected to constant boundary stresses, whereas the amount of eroded fine particles increases progressively. Therefore, the stress ratio is kept constant while the peak strength ratio and the initial stiffness evolve as a  p function of the density state of the soil. The shear strain can be written in an incremental form from Eq. ( 15): ( 23)

                1 1 1 1 ln 1 ln 1 i i i i p p d i i i i p p M M G M G M                              
with the superscript (i) representing the current loading step and the previous loading step.

 

1 i 
The volumetric strain increment can then be computed as:

(24)

  2 ' 1 exp                      v d pt d p q K A M d K p
The incremental modification of the void ratio induced by the body deformation is added to the void ratio change caused by erosion: (25)

      total er 1 v e e e      
Under triaxial condition, the axial and radial strain increments ( and ) are functions of the a  r  shear and volumetric strain increments:

(26) 1 1 1 , 3 3 2           a d v r v d
Thus, and can be used alternatively instead of and .

a  r   d  v
The influence of the effective stress on the erosion of fine particles is considered in two parts. First, the volume change induced by the effective stress ( ) has been used in the mass

    div s v t     v
balance equations (3)-(6). Thus, the impact of the effective stress is implied in the volume change, which will further induce changes in the permeability as well as in the flow rate and the erosion rate, according to Eqs. ( 10)-( 14). Secondly, the mass exchange term for the mass balance equations (3)-n ( 6) is formulated to be dependent directly on the shear stress ratio , according to Eqs. ( 13

Simulations of tests

Description of laboratory tests

Chang and Zhang (2011) performed several hydraulic-gradient controlled downward erosion tests on gap-graded HK-CDG mixtures with a newly developed stress-controlled erosion apparatus. These experiments were performed in order to investigate the initiation and the development of internal erosion in soils subjected to multi-step fluid flow under complex stress states, as well as the stressstrain behavior of soils subjected to internal erosion. The experimental set-up consisted of a triaxial apparatus, a pressurized water supply system, a particle collection system, and a water collection system. A hollow base pedestal was designed for the eroded soil and the seepage water to pass through in the downward direction. The tested soil was obtained by mixing the commercial washed

Leighton Buzzard sand (fraction E), and the completely decomposed granite (CDG) extracted from a construction site on Beacon Hill, Hong Kong, following the grain size distribution shown in Figure 1.

Each triaxial specimen, 100 mm in diameter and 100 mm in high, was prepared by moist tamping to prevent soil segregation. The erosion tests were carried out under initial isotropic and anisotropic stresses with the major principal stress being parallel to the seepage direction. After a saturation stage, the consolidation of the specimen was realized by increasing gradually the confining stress.

For the tests under an anisotropic stress condition, the vertical stress was gradually increased up to the target value after the specimen had been isotropically consolidated to a given stress state. After the stress state being applied to the specimen, the pressurized water supply and the water collection system were connected to the triaxial system. The hydraulic gradient, , was increased in stages to i the final value (i.e., 0.15 per 10 min for , 0.25 per 10 min for , and 0.50 per 10

1.0 i  1.0 2.0 i  

min for

). For 10 min, the hydraulic gradient was increased within the first 2 min and then 2.0 i  kept constant for the following 8 min. Once the soil erosion was initiated, the hydraulic gradient was kept constant until no further soil loss was observed (i.e. < 25.0g/m 2 per 10 min). Then, the hydraulic gradient was increased to another level. A typical increase of the hydraulic gradient is illustrated in Figure 2. During the internal erosion testing, the applied hydraulic gradient, the vertical displacement, the outflow rate, and the eroded soil mass were continuously measured. The radial deformation was deduced from planar deformation observed on digital photographs. Conventional drained triaxial compression tests were then carried out to study the stress-strain behavior of the soil having been subjected to internal erosion. Note that the samples of gap-graded or broadly graded granular soils are normally more or less heterogeneous with non-uniform distribution of initial void ratio and fines content in spatial space. During the erosion process, fluidized fine particles may be blocked when passing the interface between the local "loose" area and the local "dense" area, which results in heterogeneous particle size distributions. However, not enough detail (i.e. the evolution of pore pressure along the height of the specimen) was given by [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF] to be able to consider the initial heterogeneity of the soil specimen; as a consequence, the specimens in the following simulations were assumed to be initially homogeneous.

A series of tests on the gap-graded HK-CDG mixtures samples [START_REF] Chang | Internal erosion and overtopping erosion of earth dams and landslide dams[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF] were selected for simulations. BM-C-1, BM-C-2, and BM-C-3 were conventional drained triaxial tests on isotropically consolidated samples without erosion, under confining stresses of 50, 100, and 200 kPa, respectively. GS-C-1, GS-C-2, and GS-C-3 were erosion tests under isotropic stress states with confining stress of 50, 100 and 200 kPa, respectively. GS-C-4, GS-C-5 and GS-C-6 were erosion tests under anisotropic stress states with confining stress of 50 kPa and deviatoric stress of 50, 100 kPa and 150 kPa, respectively. For the eroded samples, drained triaxial tests were carried out starting at the stress state applied during the erosion tests. More details about the laboratory test soil can be found in [START_REF] Chang | Internal erosion and overtopping erosion of earth dams and landslide dams[END_REF].

Simulation schemes

The numerical modeling scheme can be summarized in Figure 3 for granular media under constant stresses considering internal erosion.

represents the current step and the previous step. The

  i   1 i 
set of partial differential equations were solved explicitly in time and by the finite difference method in space, the primary unknowns being the soil skeleton displacement (

), the pore pressure ( x Based on the above numerical scheme, the simulation could be conducted with the real dimensions of the tests by using the finite element method to investigate the temporal and spatial variations of the fines content and hydraulic conductivity during internal erosion. However, under the framework of continuum mechanics, the erosion process is uniform under the assumption of an initially uniform homogeneous specimen. In this study, the example of a well-controlled erosion coupled with the mechanical loading triaxial tests on reconstructed specimens can be seen as a case of a representative elementary volume with the hydraulic gradient across the specimen as input boundary condition.

Hence, a very interesting simplified erosion and mechanics coupled simulation model can be formulated, which should be useful in evaluating the eroded mass as well as the induced deformation of the soil element. In conventional drained triaxial test simulations before and after erosion, the confining stress was maintained constant. Thus, the ratio between the incremental deviatoric stress and the incremental mean effective stress was equal to 3. From Eqs. ( 16) and ( 17), the deviatoric stress and the volumetric strain could be obtained by imposing an incremental deviatoric strain as follows:

(27) 3 1 1 3 d p q q q G p M p                      (28) 3 q p     (29)   ' 1 exp v d pt d p q A M d K p                    
Note that the above conventional drained triaxial test simulations were performed on a representative elementary volume, thus the time and the space were not considered in the simulations. t

x

Calibration of parameters

The physical properties of the soil mixtures are summarized in Table 1. The initial void ratio 0 e before the external stress being applied is 0.53. The mechanical parameters were calibrated as follows (see Figure 4): (a) the elastic parameters ( , and ) were determined from an isotropic 2 andTable 3 were used to predict the hydro-mechanical behavior during the erosion tests and the following drained triaxial tests after erosion. Note that the parameters can also be identified by using optimization methods (Jin et al. 2018;[START_REF] Yin | An efficient optimization method for identifying parameters of soft structured clay by an enhanced genetic algorithm and elasticviscoplastic model[END_REF]Yin et al. 2018b).

Figure 5 shows the cumulative eroded soil mass and the variations in hydraulic conductivity during the internal erosion tests under different stress states. The experimental results show that the stress state during internal erosion affects significantly the hydraulic behavior of the soil. Specimens GS-C-1, GS-C-5, and GS-C-6 were tested under the same confining pressure (50 kPa), but with different deviatoric stresses: 0 kPa, 100 kPa and 150 kPa, respectively. Figure 5(a) shows the cumulative weights of the eroded soil in the three tests. It was found that the average erosion rate within the testing period increased with the increase of the deviatoric stress. [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF] attributed this observation to the fact that the primary structure formed by the coarse particles became more unstable with the increase of the applied stress ratio. Combined with the effect of the loss of fine particles that provided lateral "support" for the primary soil skeleton structure, the erosion rate in the highest stress ratio case (specimen GS-C-6) was much larger than that in the isotropic stress case (specimen GS-C-1). Moreover, due to the loss of fine particles, the permeability increased significantly in the three tests as shown in Figure 5(b). The more the fine particles were lost, the larger the permeability became. Finally, the permeability reached a relatively constant value.

The proposed model was able to reproduce the general trend of the erosion evolution until a stable step was reached. However, some differences between experimental and numerical results could be found, especially for the evolution of the hydraulic conductivity. After the loss of a significant amount of fine particles, the soil specimen reached a new equilibrium state, where the increment of the eroded mass was generally limited. During this period, the soil permeability slightly decreased because part of the transported fine particles was assumed to be captured within the soil skeleton during the erosion process, which induced a clogging of the constrictions among the coarse particles, i.e. a self-filtration phenomenon. The decreased diameters of the effective pore throats resulted in a reduction of permeability. Only few data are available in the literature concerning the self-filtration phenomenon during an erosion test. Ke andTakahashi (2014a, 2014b) attributed the deviation of the hydraulic conductivity to the difference in homogeneity along the reconstituted soil specimens , induced probably by the segregation during preparation and placement of the specimen. Another aspect which was not taken into account was the unknown influence of the saturation stage, which might also lead to the heterogeneity of the soil sample before erosion. Not enough detail on this particular aspect was given by [START_REF] Chang | Internal erosion and overtopping erosion of earth dams and landslide dams[END_REF] to calibrate the filtration term in the mass exchange equations.

Figure 6 presents the applied increasing hydraulic gradients for the three internal erosion simulations under different stress states. When the hydraulic gradient is increased to 0.6, i.e. the initiating hydraulic gradient , the fine particles will start to be eroded, as shown in Figure 5. According to start i the above cited experimental tests of [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF], once the rate of the eroded weight became smaller than 0.02g/min, the hydraulic gradient was increased to another level. Sudden increases in the erosion rate, hydraulic conductivity and deformations of the specimen were observed at a gradient of 1.5, i.e. the soil skeleton deformation hydraulic gradient , as shown in Figures 5 sd i and 7. In the numerical simulations, the final equilibrium states were reached at a gradient around 2.2 at 200 min for GS-C-1 and at 280 min for GS-C-5 and GS-C-6. Once the self-filtration mechanism is added into the simulations, the critical hydraulic gradients could be expected to increase, and the appearance of the final equilibrium state to be delayed. With the increase of the hydraulic gradient, the pore channels previously clogged could be unblocked. This aspect will be discussed in future work based on better documented experimental data that will consider the spatial variation of the fine particles during internal erosion.

Figure 7 shows the axial and radial strains of the soil specimens caused by internal erosion. Here, the negative values mean that the sample diameter increased. Larger axial and radial strains were associated with higher stress ratios during the erosion process. The model was able to reproduce the initial small deformation stage, the skeleton deformation stage and the final equilibrium stage. A higher stress ratio was associated with larger axial and radial strains during the erosion process.

Moreover, from a macroscopic point of view, the jumps in the evolution of the specimen deformation should be attributed to the changes in the hydraulic gradient. The results of the simulations under the applied increasing hydraulic gradients during internal erosion shown in Figures 6 and7 demonstrate that the numerical approach is also able to reproduce the temporary equilibrium state reached under the applied stresses with a relatively low hydraulic gradient.

Mechanical responses after erosion

Figure 8 illustrates a comparison between the experimental and the numerical results of drained triaxial tests on samples before and after erosion. Both experimental and numerical results show that the loss of fine particles during internal erosion affects significantly the mechanical behavior of the soil. When the soil loses a significant number of fine particles, its stress-strain behavior changes from dilative to contractive. When a part of the solid fraction is eroded, the density state ( ) c e e decreases and, as a consequence, the shear resistance decreases, creating additional deformation of the granular assembly.

The proposed numerical approach was able to reproduce the general trend of the soil behavior up to the peak strength which corresponded in the experiments to the development of strain localization within the specimens. In some cases, discrepancies between experimental and numerical results could be seen, possibly due to the difficulty of measuring accurately certain material properties, such as the elastic stiffness or the friction angle, as a function of the fines content. Another aspect which was not taken into account in the simulations was that heterogeneity within the specimen could be produced by the combined effects of transport and filtration of the fine particles. Not enough detail of this particular aspect was given by Chang (2012).

Discussion

Several numerical tests were performed on the HK-CDG samples to investigate systematically the soil behavior which lost fines content during the erosion process, as well as the influences of soil's initial density and initial fines content. After applying the external stresses, each specimen was subjected to a hydraulic gradient at constant external stresses, as shown in Figure 2.

Influence of the soil's initial state on the initiation of erosion

According to Figures 5 and6, the proposed model is able to distinguish the initiating hydraulic gradient and the soil skeleton deformation hydraulic gradient of the erosion process. ' 100 p  levels. Afterwards, each specimen was subjected to a hydraulic gradient while the external stresses were kept constant. Figure 10 presents the axial and the volumetric strains versus the eroded fraction of fine particles. For small stress ratios ( ), a slight contraction occurred with the increase

1 q p    
of the eroded fraction and the samples reached a final equilibrium state with an axial strain smaller than 5% after the entire erosion process. For intermediate values ( =1.0; 1.25), two stages could   be identified. Slight contractions occur firstly, and the axial strain developed slowly. After 20% and 16.5% of the total mass being eroded for the samples under =1.0 and 1.25, the deformation  increased much faster and the samples exhibited a large contraction. This result indicates that, at this stage, important grain rearrangements occurred in the microstructure so as to reach a new stable mechanical state. However, for , no new equilibrium state could be reached; more details can 1.5   be found in the next section for similar numerical tests on denser samples. Obviously, given (or c f  the fines loss in terms of percentage of the total fines ), values can be identified as the 0 100%

c c f f  
triggering factor for the development of large macroscopic deformations. Furthermore, the higher the stress ratio is during erosion, the smaller an eroded fraction will be needed in order to lead the sample to failure. The predicted axial strains of GS-C-1 ( ), GS-C-5 ( ) and GS-C-6 0

  1.2   (
) with large fines content loss are compared in Figure 11 (the solid lines correspond to loss 1.5   of fines content observed in experiments, and the dashed lines predict the erosion evolution in case the loss of fines content continues). It is shown that for the specimen GS-C-1 ( ), even if all the 0   fines are eroded, the deformation remains below 1%; for specimens GS-C-5 (

) and GS-C-6 1.2  

(

), a temporary equilibrium state is reached at , and then remains stable up to 15% 1.5

  5% c f  
(for GS-C-6) and 18% (for GS-C-5) of eroded fines. The axial strain then develops very rapidly with a small loss of fines.

As presented in Figure 12(c) and Figure 12(f), the erosion process tends to increase the sample porosity without producing significant macroscopic deformation, especially for low and medium  values. As shown by the DEM simulations in [START_REF] Scholtès | Multiscale approaches to describe mechanical responses induced by particle removal in granular materials[END_REF]Wood et al. 2008;Wood and Maeda 2008), a certain amount of fine particles can be removed without significant consequence at the macroscopic scale due to the presence of floating grains not involved in the force chain network.

Furthermore, the creation of a more open microstructure due to the loss of fines content dominates the effect of sample compression induced by local destabilizations.

In order to bring all the specimens to failure, triaxial compression was performed on the samples which stabilized after a certain fraction of the eroded mass. Figure 12 ). Even if the soil became denser and, therefore, 1.68 cs   somewhat stronger, large deformations might occur, leading to heavy damages to the earth structures.

As stated in [START_REF] Hicher | Modelling the impact of particle removal on granular material behaviour[END_REF], granular assemblies can develop instability at a shear stress level much lower than the critical state failure line; the so-called diffuse failure can occur in highly eroded soils.

Besides, the results demonstrate that the proposed model was able to take into account the evolution of the critical state which was different for the degraded specimens in comparison to the intact ones.

Conclusion

This study provides a novel contribution to the numerical approach of quantifying the impact of internal erosion on the mechanical behavior of granular soils. The detachment and transport of fine particles induced by internal erosion has been modeled by a mass exchange between the solid and the fluid phases. The governing equations were formulated based on the mass balance of four assumed constituents: the stable fabric of the solid skeleton, the erodible fines, the fluidized particles, and the fluid. A non-linear incremental model was adopted in order to obtain a non-linear stress-strain relationship for the soil mixture. The model was enhanced by the introduction of the critical state concept and by formulating the influence of the fines content on the critical state line in the e-p' plane based on experimental results in order to take into account the impact of the amount of eroded fines on soil deformability and strength.

An erosion and mechanics coupled model was formulated which is practical useful to evaluate the eroded mass as well as deformation of soil specimen during internal erosion. This practical model was used to simulate the behavior of HK-CDG mixtures before, during, and after erosion through triaxial tests. After a significant amount of fine particle loss, the stress-strain behavior of the soil changed from dilative to contractive with the shear resistance decreasing. The comparison between experimental results and simulations demonstrated that the model can reproduce with a reasonable success the initiation of internal erosion, the sudden increases in erosion rate, the change in hydraulic conductivity and the deformations of the specimen. The results confirmed that the deformation was linked to the stress ratio under which the erosion process was active. Note that it is in our interest to understand better the development of internal erosion. The experimental analysis on erosion induced deformations under controlled stress states is still in its infancy, not to mention the systematic analysis of the influence of the soil's initial states (initial relative density and fines content) under different stress states. Therefore, the influence of the stress state and the soil's initial states (relative density and fines content) were investigated from a numerical standpoint. Interestingly, besides the stress ratio, the amount of loss of fines indicating the yielding from a stable to an unstable mechanical response appeared to be related to the initial density as well as to the initial fines content of the soil mixture. These results should provide some inspiration for the design of future experiments in order to better understand the physics of internal erosion and its impact on the soil behavior.

It is worth mentioning that this study does contain some limitations which could be remediated in future work. The consideration of the self-filtration process and the initial heterogeneity of the specimen might influence the critical hydraulic gradient and delay the final equilibrium state. More advanced laboratory testing in a triaxial permeameter with well-documented data on the spatial variations during internal erosion will be needed for the development, calibration, validation and verification of the model. Tables Table 1 Physical

Figure captions

  not only on the total discharge of fluid , but also on the residual fines content w q .

  constants controlling the non-linearity of the critical state line; is 

  controls the evolution rate of the transition zone ( after[START_REF] Yin | A micromechanics-based model for sand-silt mixtures[END_REF] 2016),  20   is the threshold fines content from a coarse grain skeleton to a fine grain skeleton, varying th f usually from 0.2 to 0.35 on the basis of existing experimental data; is the ratio of the mean size of d R the coarse grains D 50 to the mean size of the fine grains d 50 . The material constants , and can

  ), which relate the sliding resistance, controlled by the mobilized peak friction angle to the density p  state ( ) of the assembly. Therefore, the loss of fine particles affects the interlocking, i.e. the c e e density state, by changing simultaneously the void ratio and the critical void ratio . When a part e c e of the fine particles is detached from the soil skeleton, the void ratio increases, the critical void e ratio evolves according to the fines content evaluated by Eqs. (21)-(22). According to Yin et al.

  , as shown by experimental results (e.g.[START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF]. However, the impact of the confining pressure on the erodibility of fine particles still remains unclear since the results in the literature[START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF][START_REF] Bendahmane | Experimental parametric study of suffusion and backward erosion[END_REF][START_REF] Luo | Hydro-mechanical experiments on suffusion under long-term large hydraulic heads[END_REF], and Ke and Takahashi 2014b[START_REF] Moffat | Spatial and temporal progression of internal erosion in cohesionless soil[END_REF] present contradictions. More advanced experimental tests with visualization of particles are needed to calibrate more precisely the relationship between the erosion rate and the stress state.

  -(12). For clarity, time and space variables are omitted in the equations. t

  stress paths of the triaxial tests; (c) the CSL related parameters ( , ) were   measured from drained triaxial tests. It was difficult to determine these two parameters because strain localizations appeared systematically at the end of the tests (Chang 2012); due to a lack of experimental data, the assumption of parallel critical state lines was made, and the values suggested by Yin et al. (2014) were adopted. (d) due to the lack of data concerning Chang's tests, ); (e) finally, the dilatancy parameters ( and ) were determined by fitting the stress-strain d d A relations obtained in the triaxial tests BM-C-1, BM-C-2, and BM-C-3. The erosion parameters were calibrated by fitting the time evolution of the hydraulic conductivity and the cumulative loss of dry mass in tests GS-C-1 and GS-C-5. All determined parameters summarized in Table

  performed on the HK-CDG samples to further study the influence of the soil's initial state on the initiation of erosion. Figure9presents the cumulative weight of the eroded soil as a function of hydraulic gradient for (a) different initial fines contents and (b) different initial void ratios. It is obvious that the increase of the initial fines content or the decrease of the initial void ratio increases the hydraulic gradient required to cause erosion, which appears reasonable in view of the experimental phenomenon.4.2. Effect of internal erosion on the hydro-mechanical responseNumerical tests were performed on the HK-CDG samples, were isotropically consolidated up to kPa and then sheared at different stress

  presents the simulations of the erosion for different stress ratios and the following triaxial compression of the samples = q porosity produced a change in the specimen behavior from dilative to contractive. As shown also in Figure10, the strain path was linked to the stress ratio. It is remarkable that the limiting threshold of under which the sample failed during erosion was generally smaller than the  value of at critical state ( ) obtained from the drained triaxial tests and dependent on the q p cs  amount of the loss of fines. For a high degree of erosion (i.e.as shown in Figure12d-f), 0.25 c f   the specimen developed large compaction during erosion, even though the stress ratio was far =1.0  from its value at the critical state (

4. 3 .

 3 Influence of the soil's initial density on the hydro-mechanical responseNumerical tests were performed on dense samples of HK-CDG mixture, with and 0 0.35 e  , in order to evaluate the influence of the initial density of the soil on the mechanical erosion. By comparing Figure13with Figure10, it can be concluded that, at a given stress ratio, a denser sample is more resistant and less inclined to deform when fine particles are progressively eroded. A larger value is, therefore, necessary to trigger the development of c f  large macroscopic deformations. For high stress ratios ( ), the specimens dilated with the loss 1.5   of fines. Let us consider the specimen with as an example. The volumetric and axial strains 1.6   of the specimen increased slowly to 0.2% and 1.4%, respectively, until . specimen increased rapidly, indicating that the mechanical state of the soil changed from stable to unstable, leading to the transition of the development of the axial strain from a slow to a large increase. At , a sudden decrease of the volume indicated the failure of the 18the specimen was linked to the collapse of the initial microstructure, leading to the rapid decrease of the specimen volume.

Figure

  Figure 14 presents the simulations of erosion for different stress ratios and the following = q p 

Figure 1 .Figure 7 .

 17 Figure 1. Grain size distribution of the tested HK-CDG mixture

Figure 8 .

 8 Figure 8. Comparison between experimental results and simulations for HK-CDG mixture before

Figure 14 .

 14 Figure 14. Simulations of erosion on initially dense HK-CDG mixture with 0 0.35 e  for different

Figure 2 .Figure 5 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 12 .Figure 13 .Figure 14 .Figure 15 .

 257891012131415 Figure 2. Illustration of increasing hydraulic gradient applied during internal erosion test

Table 2 .

 2 Mechanical constants of HK-CDG mixture

		Elastic		Dilatancy		CSL-related
	parameters K 0 /MPa G 0 /MPa	n	parameters d A d	parameters e hc,cr0 e fc,cr0 		 c
	50	16.7	0.95	18 0.35	0.8	0.86 0.45 0.03 40.5

Table 3 .

 3 Erosion constants of HK-CDG mixture

		Fines		Erosion
	a	parameters m	f th	parameters  e  1  2  3 m k
	0.79 0.46 0.28	1 0.8 2.7 0.05 4

Acknowledgement

The financial support provided by both the National Natural Science Foundation of China (51579179) and the National Institute for Industrial Environment and Risks of France (INERIS) is gratefully acknowledged.

Appendix A1. Four-constituent continuum model

In a given Representative Elementary Volume (REV), dV consisting of the four constituents, the volume fraction of a single constituent x is expressed as:

with  

, , , i ss se ff fp  denoting the 4 constituents: the stable fabric of the solid skeleton, the erodible fines, the fluidized particles and the pure fluid phase;

x dV denoting the volume of the corresponding constituent.

At a material point level, the mass balance for the i phase is given, neglecting the hydro-mechanical dispersion tensor, by [START_REF] Schaufler | Infiltration processes in cohesionless soils[END_REF]:

where , ex i  and i v denote, respectively, the mass exchange term and the velocity of the corresponding constituent. The mass balances for the four constituents can be reduced to the corresponding volume balance:

n is the term of the volume of mass exchange.

The porosity field  