Antoine Wautier
email: antoine.wautier@irstea.fr

Stéphane Bonelli

François Nicot

DEM investigations of internal erosion: grain transport in the light of micromechanics

Keywords: Suffusion, DEM, Fluid grain interactions, Grain Transport, Force chain, Granular material

Internal erosion by suffusion can change dramatically the constitutive behavior of granular materials by modifying the fabric of granular materials. In this study, the effect of an internal fluid flow on granular materials is investigated at the material point scale using the numerical coupling between a Discrete Element Method (DEM) and a Pore-scale Finite Volume (PFV) coupling scheme. The influence of the stress state and the hydraulic loading (direction and intensity) on the occurrence of grain transport in dense widely graded granular samples is thus investigated and interpreted in terms of micromechanics. In particular, it is shown that, grain transport is increased when the macroscopic flow direction is aligned with the privileged force chain orientation. The stress induced microstructure modifications are shown to influence the transport distances by controlling the number of rattlers.

INTRODUCTION

Internal erosion is the process by which a saturated soil can lose a part of its solid fraction under the action of an internal fluid flow. If this phenomenon occurs for instance in oil reservoir, the resulting sand production may clog the reservoir and reduce oil production [START_REF] Papamichos | Volumetric sand production model and experiment[END_REF]. If internal erosion occurs in earthen dikes or dams or in natural deposits of loose rocks in mountain rivers, the modifications of the internal structure of the soil can generate mechanical instability that will lead to unexpected failures in case of sudden changes in the loading conditions [START_REF] Seed H Bolton | The failure of Teton dam[END_REF][START_REF] Stéphane | Erosion in geomechanics applied to dams and levees[END_REF][START_REF] Jonathan | The 1985 Bairaman landslide dam and resulting debris flow, Papua New Guinea[END_REF]. Among the four types of internal erosion identified today [START_REF] Stéphane | Erosion in geomechanics applied to dams and levees[END_REF], suffusion, that consists in the selective erosion of the finest particles of a soil, affects directly the internal fabric of granular materials. Its modeling at the continuum scale is still a complex issue as it requires the identification of internal state variables able to account for the microstructure state. To overcome the lack of reliable constitutive models, efforts have been historically devoted to anticipate and prevent the occurrence of grain transport. Consequently, many suffusion susceptibility criteria have been proposed to ensure soil self-filtration. These criteria are usually constructed either on the particle size distribution (PSD) or constriction size distribution (CSD) curves. As a result they often lack of microscale justifications. Moreover, they may prove to be too conservative as they totally dismiss the consequences of suffusion on the mechanical behavior of granular materials. The notion of internal stability derived from these criteria is thus completely different from the notion of mechanical stability. Thanks to ever increasing computation power together with recent efforts to develop numerical schemes accounting for fluid/grain interactions locally [START_REF] Zhao | Lattice Boltzmann method for simulating particle-fluid interactions[END_REF][START_REF] Bruno | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF][START_REF] Jidong | Numerical modeling of fluid-particle interaction in granular media[END_REF][START_REF] Yanhui | LBM-DEM modeling of fluid-solid interaction in porous media[END_REF][START_REF] Martin | Fluid-particle flow simulations using two-way-coupled mesoscale SPH-DEM and validation[END_REF], the full three-dimensional modeling of suffusion at the scale of representative elementary volumes (REV) is now accessible. The results of such numerical simulations can be used to analyze locally the elementary mechanisms responsible for suffusion and their consequences on the mechanical behavior of granular materials. Such a micromechanical approach will provide a physical meaning to additional terms introduced in phenomenological models or drive the definition of physically sound micromechanically-based models. In this paper, REV of granular materials are considered through a discrete element modeling approach and fluid/grains interactions are accounted for using a pore network model developed by Chareyre et al. [START_REF] Bruno | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF]. This numerical scheme offers a local scale modeling at a relatively cheap computational cost in order to enable fully resolved simulations in three dimensions. This contrasts with most of the approaches found in the literature that either consider 2D problems and introduce a hydraulic radius to recover the connectivity of the fluid phase [START_REF] Yanhui | LBM-DEM modeling of fluid-solid interaction in porous media[END_REF][START_REF] Nariman | The impact of fluid flow on force chains in granular media[END_REF][START_REF] Jeff | Two-dimensional numerical simulation of chimney fluidization in a granular medium using a combination of discrete element and lattice Boltzmann methods[END_REF] or use in 3D a non resolved CFD approach [START_REF] Zeghal | A continuum-discrete hydromechanical analysis of granular deposit liquefaction[END_REF][START_REF] El Shamy Usama | Coupled continuum-discrete model for saturated granular soils[END_REF][START_REF] Feng | Coupled discrete element and finite volume solution of two classical soil mechanics problems[END_REF][START_REF] Jidong | Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics[END_REF][START_REF] Tong | A coupled CFD-DEM analysis of granular flow impacting on a water reservoir[END_REF][START_REF] Hama | Analysis of mechanical behaviour and internal stability of granular materials using discrete element method[END_REF][START_REF] Kenichi | Coupled DEM-CFD Analysis of the Initiation of Internal Instability in a Gap-Graded Granular Embankment Filter[END_REF]. The objective of this paper is to investigate the microscale impact on an internal fluid flow on the microstructure of widely graded granular materials for different stress states and hydraulic gradients. In particular, the occurrence of grain transport is analyzed with respect to micromechanical features of the granular assemblies. This paper is organized as follows. In a first section, the numerical method used to tackle the 3D problem of suffusion at the material point scale is reviewed. A particular attention is devoted to the definition of the flow boundary conditions and to an erosion criterion. Then, the relative occurrence of transport and clogging is analyzed for different stress states and hydraulic gradients and related to micromechanical features of granular materials in the second section.

MATERIAL POINT SCALE SIMULATIONS OF SUFFUSION

Down to the microscale, soil can be seen as a collection of individual particles interacting together at contact points. However, at the engineering scale, this discrete nature cannot be accounted for and is hidden behind a macroscopic continuum description. Between the macroscopic and microscopic scales, an intermediate scale exist for which a continuum media can be defined by averaging discrete quantities. This scale is known as the Representative Elementary Volume (REV) scale in the sense of Hill [START_REF] Rodney | Elastic properties of reinforced solids: some theoretical principles[END_REF] and a quantitative definition with respect to suffusion can be found in Wautier et al. [START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF]. In the following sections, the modeling of suffusion at this scale is detailed.

DEM/PFV numerical modeling

At the REV scale, individual soil particles are accounted for using a discrete element method (DEM). In DEM, only two ingredients are needed to account for the evolution of a mechanical system: the geometry of the particles and the interaction forces between them. In this study only non-cohesive granular materials are considered in the form of poly-disperse assemblies of spheres. Particle interactions are modeled with the classical elasto-frictional contact law proposed by Cundall and Strack [START_REF] Peter | A discrete numerical model for granular assemblies[END_REF] and illustrated in Figure 1 . In this law, a contact is created as two particles overlap. The two contacting bodies are then assumed to be linked by two fictitious linear springs and a slider. Based on the direction of the vector joining the particle centers, a normal force is defined as proportional to the overlapping distance between the two particles where the contact normal stiffness is proportional to the material's Young modulus and to the harmonic average of the two particle radii and ′ ( = 2 ′ + ′ ). In addition to the normal force, a tangential force is introduced. This tangential force is proportional to the relative tangential displacement between the two particles where the horizontal contact stiffness is a fraction of its normal counterpart ( = ). Compared to the normal displacement, the tangential displacement cannot be defined directly on the current geometry. It is indeed defined in an incremental form with Δ being the relative tangential displacement of the contact point between two successive geometrical configurations. The linear behavior of the tangential force is limited by a Coulomb friction condition based on an internal friction angle , which defines the largest accessible ratio ∕ (see Figure 1 ). The input parameters used in this elasto-frictional contact law are reported in Table 1 . They are chosen equal to those of Hadda et al. [START_REF] Nejib | Micromechanical analysis of second order work in granular media[END_REF] and Wautier et al. [START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF]. After computing all inter-particle contact forces, the induced particle displacements are integrated based on Newton's second law of motion using the DEM open source code YADE (23). In order to provide a full microscopic modeling of suffusion, the local fluid impact has to be simulated also at the local scale. From a conceptual point of view, accounting for fluid/grain interaction consists in the simulation of a laminar flow in a pore space which geometry is explicitly known with no-slip boundary conditions on the fluid/grain interfaces. Under the additional assumption of fluid incompressibility, the local fluid problem consists in solving Stokes and mass conservation equations

∇ 2 -∇ + = ∇ ⋅ = 0 , (1)
where is the fluid velocity, the fluid pressure field, the fluid dynamic viscosity and a body force taken here equal to zero as gravity forces are dismissed.

Given the fluid velocity field, the action of the fluid on the different particles of the material is deduced by integrating the shear and normal stresses on grain interfaces. Eventually, these additional forces can be used in DEM to update the grain positions and the pore space geometry. Different numerical schemes exist in the literature in order to solve for the fluid phase. The four main approaches are summarized below:

-A first approach is to simulate the fluid flow with resolved Computational Fluid Dynamics (CFD) techniques. In this approach, the pore space is triangulated in order to solve Equation (1) using Finite Element Methods (FEM). As tens to hundreds FEM elements per pore are required to get accurate results, the fluid solving step is extremely expensive and published studies are extremely scarce [START_REF] Bruno | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF]. Moreover, frequent remeshing is needed to account for microstructure evolutions.

-As a result, much of the CFD approaches found in the litterature are indeed unresolved CFD approaches [START_REF] Zeghal | A continuum-discrete hydromechanical analysis of granular deposit liquefaction[END_REF][START_REF] El Shamy Usama | Coupled continuum-discrete model for saturated granular soils[END_REF][START_REF] Feng | Coupled discrete element and finite volume solution of two classical soil mechanics problems[END_REF][START_REF] Jidong | Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics[END_REF][START_REF] Tong | A coupled CFD-DEM analysis of granular flow impacting on a water reservoir[END_REF][START_REF] Hama | Analysis of mechanical behaviour and internal stability of granular materials using discrete element method[END_REF][START_REF] Kenichi | Coupled DEM-CFD Analysis of the Initiation of Internal Instability in a Gap-Graded Granular Embankment Filter[END_REF]. Unresolved CFD is a hybrid approach as the fluid problem is not solved at the pore scale but on fluid cells much larger than the particles. As a result, Equation (1) has to be changed to account for the local porosity of the granular media. Compared to resolved CFD, unresolved CFD is computationally less expensive and requires no remeshing but in the meantime the microscale description of the problem is lost and ad hoc macroscopic assumptions are needed to account for the fluid flow (Darcy's law and Kozeny-Carman expression for the permeability for instance).

-A third option to account for the fluid at the local scale is to use Lattice Boltzmann Methods (LBM) [START_REF] Pierre | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF][START_REF] Mansouri Mouloud | A 3D DEM-LBM approach for the assessment of the quick condition for sands[END_REF][START_REF] Zhao | Lattice Boltzmann method for simulating particle-fluid interactions[END_REF][START_REF] Pablo | Micromechanical analysis of the surface erosion of a cohesive soil by means of a coupled LBM-DEM model[END_REF][START_REF] Jaeger Robin | Mechanism behind erosive bursts in porous media[END_REF][START_REF] Jeff | Two-dimensional numerical simulation of chimney fluidization in a granular medium using a combination of discrete element and lattice Boltzmann methods[END_REF]. These methods rely on a statistical physics approach and Equation (1) is replaced by Boltzmann discrete equation. The fluid behavior is described through a collision-propagation scheme [START_REF] Pierre | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]. As for resolved CFD approaches, this numerical scheme requires a discretization of the pore space but LBM relies on a fixed lattice that is not changed as particles move.

Even if this numerical technique can be easily parallelized, its implementation in 3D is still very expensive. In addition, it should be noted that, incompressible flow represents a high difficulty for the conventional LBM, which is entirely based on density fluctuations.

-A fourth class of approach lies in the use of pore network models [START_REF] Steven | Network model evaluation of permeability and spatial correlation in a real random sphere packing[END_REF][START_REF] Thompson Karsten | Modeling flow in disordered packed beds from pore-scale fluid mechanics[END_REF][START_REF] Stig | 3D pore-scale modelling of sandstones and flow simulations in the pore networks[END_REF][START_REF] Markus | Calibration of a pore-network model by a pore-morphological analysis[END_REF][START_REF] Tarek | Network model for hydraulic conductivity of sand-bentonite mixtures[END_REF]. Indeed, these methods (initially developed to predict the permeability of porous media) are based on a discretization of the void space into connected pores. The complexity of the fluid flow is then reduced to fluid exchange laws between adjacent pores. As a result, the accuracy of these approaches strongly relies on adequate definitions of how fluids are exchanged between pores in terms of the local pore geometry. A second difficulty consists in computing the fluid forces on the surrounding particles as the fluid state is defined as piecewise constant within each pore. However, provided that these difficulties can be overcome, these approaches provide a microscale fully coupled fluid/grain modeling at a relatively cheap computational cost.

In this paper, the efficient Pore-scale Finite Volume method introduced by Chareyre et al. ( 6) is used to account for internal fluid flows in granular materials at the local pore scale and in three dimensions.

In the PFV method, a regular Delaunay triangulation is used to simplify the geometry of the void space by defining a connected network of pores. Based on this partition of the void space, the mass conservation equation gives for an incompressible fluid (and for no-slip boundary conditions on the fluid/grain interfaces) a relation between the time derivative of the fluid volume contained in pore Ω and fluid fluxes from pores to as

d d = - 4 ∑ = 1 . ( 2
)
As Stokes equation ( 1) predicts a linear relationship between the fluid velocity and the pressure gradient, the fluid exchange between two adjacent pores is sought in a Darcy (or Poiseuille) form. By assuming the fluid pressure to be constant in each pore, fluid fluxes are computed as

= - (3)
where is the hydraulic conductivity of constriction linking pore to pore and is the inter-pore distance. In the PFV model, Equation ( 1) is thus reduced to the two Equations ( 2) and ( 3) and the unknowns of the problem are the pore pressures and the fluid fluxes where ( , ) ∈ [[1, pores ]], pores being the the number of pores. As a result, the fluid flow problem is reduced to a linear algebra problem and basically to the inversion of a sparse matrix for which several efficient algorithms exist. In Equation (3), the definitions of the geometrical parameters and are not given. The accuracy of the PFV model relies indeed on astute definitions of these two parameters based on the local geometry of the constrictions joining two adjacent pores. A detailed definition of and as well as a validation of the PFV scheme against fully resolved CFD simulations can be found in Chareyre et al. [START_REF] Bruno | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF].

Representative elementary volume definition and stress boundary conditions

In this study, non periodic cubic volumes of granular materials are considered. Such REVs are composed of spherical particles following the widely graded particle size distribution (PSD) illustrated in Figure 2 . This PSD is such that the ratio between the largest and smallest particles is 10 and such that the coefficient of uniformity Cu = 60 ∕ 20 = 1.98 where 20 and 60 are the particle diameters corresponding to respectively 20 % and 60 % of the mass passing fractions. After generating a cloud of 5,000 non-overlapping spheres surrounded by six bounding planes defining a cube, the particles are inflated and allowed to rearrange according to the radius expansion technique. This process is stopped when the confining pressure applied to the bounding planes reaches 20 kPa and the normalized unbalanced force ( unb ) of the system decreases below 1.5 10 -2 . By definition unb is equal to the mean summary force on the particles divided by the mean contact force magnitude on the contacts:

unb = 1 ∑ =1 || || 1 ∑ =1 || || (4)
where = ∑ is the sum of all contact forces acting on particle . unb has no units and provides an intrinsic measure of how close to equilibrium a mechanical system is. During this process, the inter-particle friction angle is maintained to 35 • in order to prepare a dense granular material with a void index = 0.55. The resulting sample shown in Figure 2 is then rescaled such that mean = 2.10 -4 m. It should be highlighted here that despite the large aspect ratio between the largest and smallest particles, the used PSD is not identified as internally unstable according to Kezdi's filter rule [START_REF] Kézdi | Soil physics (selected topics) developments in geotechnical engineering[END_REF] or Kenney and Lau criterion [START_REF] Kenney | Internal stability of granular filters[END_REF]. At the REV scale, a Cauchy stress tensor can be defined from the forces applied on the bounding walls as long as the system remains quasi-static. According to soil mechanics convention, compressions are counted positive. The stress ratio is then introduced as the ratio between the deviatoric stress and the mean pressure

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ = = √ 3 2 dev ∶ dev dev = - = 1 3 Tr( ) . ( 5
)
where stands for the identity tensor and ":" stands for the double dot contraction product.

From the bounding walls displacements, a macroscopic homogeneous strain tensor is also defined. The volumetric strain is then obtained as

= Tr( ) (6)
where Tr is the trace operator.

From the obtained equilibrium state, a drained triaxial loading is imposed to the sample in the form of a double-step procedure. First, the confining pressure ( 0 ) is increased from 20 kPa to 100 kPa by allowing the bounding walls to move. Once a new equilibrium state is reached ( unb < 10 -4 , dimensionless), a vertical compression strain rate ̇ = 0.01 s -1 is applied up to 3 % of deformation while keeping the same lateral confining pressure 0 . This strain rate is chosen similar to the one used in previous numerical studies [START_REF] Nejib | Micromechanical analysis of second order work in granular media[END_REF][START_REF] François | Failure in rate-independent granular materials as a bifurcation toward a dynamic regime[END_REF] and is supposed to be sufficiently small so that the loading can be considered as quasi-static.

In Figure 3 , the stress ratio and volumetric strain responses are shown. A typical contractive/dilative behavior of a dense soil is observed. While the stress ratio increases, two samples are saved for ∈ {0.7, 0.9} (corresponding to diamonds in Figure 3 ). As a result, two granular material REVs with slightly different stress and microstructure states are obtained. This ends the sample preparation procedure.

Flow boundary value problem and bounce back erosion criterion

By adjusting the displacements of the six bounding planes around the considered granular samples, the stress or strain state of the material is controlled at the REV scale. Concerning the fluid boundary conditions, hydraulic gradients are imposed by controlling the water pressure or flux on the six bounding planes. In the following, two flow directions are considered:

-parallel to the principal loading direction (along axis) and referred to as vertical flow case hereafter;

-perpendicular to the principal loading direction (along axis) and referred to as horizontal flow case hereafter.

Figure 4 summarizes the boundary value problem solved for the two flow cases. On upstream faces, the fluid pressure is set to 0 while on downstream faces null fluid pressure and null displacement are imposed. On the remaining faces of the sample, zero fluxes conditions are imposed ( = = 0 for the vertical case and = = 0 for the horizontal one). The upstream pressure 0 is prescribed such that the hydraulic gradient = 0 is set either to 0.1 or 1 with = 1000 kg.m -1 the fluid density, = 9.81 m.s -2 the earth gravity and ∈ { , } the sample length in the flow direction (since gravity forces are not considered, is simply a dimensionless rescaling of 0 ). It should be emphasized that the additional fluid pressure 0 applied only on the upstream wall generates an effective stress gradient between the upstream and downstream walls of the sample. As a result, in the direction of the pressure gradient, the downstream effective stress is 0 higher than the upstream one. To ensure the controllability of the system in YADE, zero displacement conditions = 0 or = 0 are imposed to the downstream walls. Since the effective stress is now varying at the scale of the sample, the problem of interest is no longer a homogenization problem and the concept of REV vanishes. However, as the fluid pressure drop 0 is small compared to the effective stress 0 and if we neglect the flow induced microstructure heterogeneities at the scale of the sample, the concept of REV is restored approximatively. Under this hypothesis, the fluid can be seen as an external perturbation that makes the microstructure evolve without any change in the stress state. In total, six DEM/PFV simulations are defined by considering different stress/microstructure states and different hydraulic loading conditions. They are summarized in Table 2 . In this paper, non-periodic boundary conditions are used in order to avoid collective motion of the particles under the action of the fluid. As a result fluid forces can be applied to all particles. However, non periodic boundary conditions generate strong boundary effects when grain transport is considered. Indeed, transported particles are either i) extracted out from the core of the FIGURE 4 Boundary value problems solved for two flow directions. In the vertical case (a), zero flux boundary conditions are applied in directions and while on the top face a pressure 0 is imposed while keeping = 0 on the bottom face. In the horizontal case (b), zero flux boundary conditions are applied in directions and while on the left face a pressure 0 is imposed while keeping = 0 on the right face. sample as they reach the downstream wall or ii) filtrated by the primary fabric as they are blocked by a too narrow constriction. In the first case, the downstream wall induces artificial clogging as illustrated in Figure 5 . To avoid such boundary effects, it is thus important to detect and delete particles reaching the downstream wall. A naive and computationally expensive approach is to look for particles with a single contact with the downstream wall. The main problem with this criterion is that it has to be run at every time step as extracted particles often bounce several times and the contact with the downstream wall never lasts for long. As a result, a prediction/correction algorithm was preferred to this first approach. This algorithm is illustrated in Figure 5 and consists in the following two steps:

-At time , rattlers are identified. By assuming that their velocity ( ) remains constant over the interval [ , + Δ ], rattlers that will virtually go through the downstream wall during this interval are identified as candidates to particle removal. This is the prediction step.

-At time + Δ , the scalar product between the updated velocity ( + Δ ) and the outward normal to the downstream wall is computed for rattlers candidates to particle removal. ( + Δ ) ⋅ < 0 means that the considered particle bounced back on the downstream wall during the time interval [ , + Δ ] and will be removed from the simulation. This is the correction step.

This approach reduces the computational cost needed to detect extracted particles as loops over the number of particles are required only from time to time and not at every DEM time step. In practice, Δ is set to be equal to 1,000 DEM time steps.

To improve computational efficiency, when extracted particles are identified, they are not deleted immediately. Indeed, as the PFV scheme relies on a particle based tessellation of the space (see Section 2.1), remeshing is needed each time a particle is deleted. As a result, extracted particles are kept in memory and are only deleted after several runs of the prediction/correction algorithm (every 5,000 iterations in practice which is small enough to avoid that particle that bounced back collide with an other particle).

INTERNAL FLUID FLOW IMPACT ON GRANULAR MICROSTRUCTURES

For cohesionless granular materials, internal fluid flows can modify the constitutive mechanical behavior at the continuum scale through two kinds of mechanisms:

-a first change lies in the transport of rattlers, i.e. particles with no contacts with their neighbors in the absence of gravity.

Under the action of an internal fluid flow, these particles move through the pore space and are either eroded or clogged depending on their size with respect to those of pores and constrictions.

-a second flow impact consists in introducing additional forces on the primary skeleton of the granular material [START_REF] Nariman | The impact of fluid flow on force chains in granular media[END_REF][START_REF] Rémi | Pore-scale flow measurements at the interface between a sandy layer and a model porous medium: Application to statistical modeling of contact erosion[END_REF][START_REF] Antoine | Flow impact on granular force chains and induced instability[END_REF]. This results in stress redistribution and possibly to grain detachment if the existing contact network is forced to rearrange.

In the following, the first aspect is emphasized while the second one has been isolated by artificially removing rattlers and was discussed in a previous paper [START_REF] Antoine | Flow impact on granular force chains and induced instability[END_REF].

Numerical results for grain transport

By tracking the positions of all particles during the fluid flow simulations, travel distances can be estimated for all the particles. The effective travel distances may be computed by integrating the norm of the particle velocities over time, but if transient evolutions are dismissed, the shortest displacements of all particles can be defined by simply comparing the initial and final positions of all particles. For a given particle , moving an initial position ( 0 ) to a final one ( ) in DEM/PFV simulations, the associated shortest travel distance reads

short num ( ) = || ( ) -( 0 )||. ( 7
)
Because of the finite dimension of the sample, short num ( ) does not always correspond to a maximum travel distance. Indeed, transported rattlers are either blocked by a too narrow constriction or reach the downstream side of the sample before being removed from the sample. In the latter case, short num ( ) corresponds to the initial distance between the particle and the downstream wall. If the sample would have been larger, these particles would have traveled farther away. In Figures 6 a and6 b short num ( ) are plotted for all particles for the sample corresponding to = 0.9 and subjected to a hydraulic gradient = 0.1 applied in vertical or horizontal directions, respectively. In Figures 6 a and6 b the clouds of points corresponding to extracted (red diamonds) and non extracted (blue crosses) particles superimpose, which shows that the samples are sufficiently large to compute mean travel distances representative of material properties and not influenced by sample boundaries. Compared to the sample dimensions ( , , ) = (43.7 mean , 43.7 mean , 43.5 mean ) some particles are transported over quite large distances up to 20 times the mean radius, but the average travel distances ̄ short num are limited ( ̄ short num < 2 mean ). Indeed, too few pores and constrictions large enough to allow for transport exist, and clogging dominates over erosion. This is consistent with the fact that the used PSD in Figure 2 is not identified as internally unstable according to Kezdi's filter rule [START_REF] Kézdi | Soil physics (selected topics) developments in geotechnical engineering[END_REF] or Kenney and Lau criterion [START_REF] Kenney | Internal stability of granular filters[END_REF]. The mean transport distances ̄ short num are shown in blue solid lines. In black solid lines the numerical predictions corresponding the transport model of Section 3.2 are shown (pressure map computed using the PFV scheme). Dashed and dot-dashed lines correspond to the predictions computed from the linear approximation of the macroscopic pressure gradient respectively on the pore network computed with all particles and on the pore network computed with only non free particles.

Pore network model predictions

In this subsection, the mean travel distances ̄ short num computed numerically from DEM/PFV simulations are compared with pore network model predictions ̄ short . To this respect, three predictions of ̄ short based on the three weighted pore networks illustrated in Figure 7 are considered.

-The first pore network relies on a triangulation of the pore space based on all particles and a pressure map computed using the PFV scheme (Figure 7 (a)).

-The second pore network relies on a triangulation of the pore space based on all particles and a pressure map computed using the linear approximation of the macroscopic pressure gradient (Figure 7 (b)).

-The third pore network relies on a triangulation of the pore space based only on non free particles and a pressure map computed using the linear approximation of the macroscopic pressure gradient (Figure 7 (c)).

Mean transport distances are then estimated on these weighted pore network without running any numerical simulations using the propagation criterion introduced in (20). This criterion is based on the definition of a flow intensity indicator in every constriction linking adjacent pores. Based on the Hagen-Poiseuille flow velocity profile for a circular tube and the expression of the drag force of a uniform flow acting on a sphere for low Reynolds numbers (Figure 8 ), the drag force acting on a particle of radius in a constriction of radius and subjected to a pressure gradient is indeed assumed to scale as Given a pore connected to four neighboring pores { 1 , 2 , 3 , 4 }, a particle of radius in the pore will move to a neighboring pore ∈ { 1 , 2 , 3 , 4 } if and only if

∼ 2 . ( 8
)
> || || 2 ≤ || || 2 ∀ ∈ { 1 , 2 , 3 , 4 } , (9)
where is the constriction radius of the connection between pores and . In other words, the particle propagates to the next accessible pore in the direction of maximum drag force. This propagation criterion initially introduced in (20) enables to compute the expected shortest travel distance short ( ) using the propagation criterion constructed from Equations ( 9) and ( 8) for any radius and pore . The expected mean travel distance ̄ short ( ) is then computed by averaging short ( ) over the number of pores. Mean expected shortest travel distances ̄ short computed for the three pore networks of Figure 7 are plotted in Figures 6 a and 6 b. While ̄ short computed on the pore network build only on non-rattlers overestimates the real travel distances, ̄ short computed on the pore network build on all particles underestimates the real travel distances. Mean expected shortest travel distances ̄ short give indeed lower and upper bounds of the transport distances estimated using fully coupled DEM/PFV simulations:

-When the pore network is computed with all particles, pore and constriction sizes are underestimated. For instance, for a free particle in a tetrahedron composed of four particles, four small pores are identified in the corresponding pore network whereas the free particle moves in a single large pore (see Figure 7 ). As a result, mean travel distances computed from this pore network provide lower bounds for the effective mean travel distances.

-When the pore network is computed with only non free particles, pore and constriction sizes are more relevant but collective motions of the free particles are not accounted for (in particular collective clogging). As a result, mean travel distances computed from this pore network provide upper bounds for the effective mean travel distances.

In Figures 6 a and 6 b, the comparison between the predictions obtained with the pore networks of Figure 7 (a) (pressure map computed using PFV scheme) and Figure 7 (b) (pressure map computed using the linear approximation of the macroscopic pressure gradient) shows that local pore pressure fluctuations tend to increase the travel distances. This result justifies the use of a resolved approach to model the fluid flow among the modeling options listed in Section 2.1.

Micromechanical insight on particle transport

When comparing Figures 6 a and 6 b, transport and erosion seem to be more important for the vertical than for the horizontal flow direction. This can be quantified in Figure 9 by checking the time evolution of the particle mass fraction removed from the simulation for the six DEM/PFV simulations of Table 2 . In order to present the results in a dimensionless form, a characteristic time scale is defined as = √ 2 0 with being the length of the sample in the flow direction ( ∈ { , } depending on the flow case), being the particle density and 0 the fluid pressure drop. corresponds indeed to the time needed for a particle initially at rest to move from a distance under the action of a pressure gradient 0 ∕ . Figure 9 , shows that only 1.4-1.9 % of the total mass of the particles is removed from the samples in the end of the simulations. This corresponds to 10-14% of the initial mass of rattlers as for ∈ {0.7, 0.9} the initial rattlers population represents respectively 13.41 % and 13.96 % of the total mass of particles. This result is consistent with the fact that the used particle size distribution is self-filtrating and that clogging dominates over erosion. A larger eroded mass is observed i) for the vertical flow case, ii) for the larger hydraulic gradient = 1 and iii) for the largest stress ratio = For 3.5 , all particles have the downstream size the sample which shows that the local tortuosity slows down particle extraction by a factor of 3.5. These observations indeed be interpreted through micromechanical analysis. It is now well that stresses in granular materials concentrate on a limited number of organized in mesostructures called force chains [START_REF] Drescher | Photoelastic verification of a mechanical model for the flow of a granular material[END_REF][START_REF] Liu | Force fluctuations in bead packs[END_REF][START_REF] Farhang | Bimodal character of stress transmission in granular packings[END_REF][START_REF] Bardenhagen | Numerical study of stress distribution in sheared granular material in two dimensions[END_REF][START_REF] Bernard | Micromechanics of granular materials[END_REF]. It has also been shown the fraction of particles involved in force chains chained is closely related to the macroscopic stress state. From drained triaxial tests, it was shown that, when the axial loading starts, force chains are in the loading directions [START_REF] Antoine | Scale separation between grain detachment and grain transport in granular media subjected to an internal flow[END_REF]. As long as the material contracts, more and more are to new force chains. For dense materials characterized by a contractive/dilative behavior, the number of force chains decreases after the characteristic point ( is minimum) up to the peak stress point ( is maximum). Simultaneously to the decrease in the number of chained particles, more and more rattlers (i.e. particles with no contacts) are released in the pore space. These trends are illustrated in Figure 10 in which the chained particles are identified using the force chain definition proposed by Peters et al. [START_REF] Peters | Characterization of force chains in granular material[END_REF]. As illustrated in Figure 10 , for the two samples considered in this study, force chains are aligned with the vertical direction (the major compression direction). This load induced anisotropy explains why more particles are extracted from the sample in the vertical flow case. The flow is indeed parallel to the force chains which facilitates grain transport.

The two samples considered in this paper are saved in the dilation regime (see Figure 3 ). As a result, the number of chained particles decreases when increases from 0.7 to 0.9 and conversely the number of rattlers for = 0.7 is smaller than for = 0.9 (see Figure 10 ). In addition, the macroscopic dilation imposes locally an increase in pore and constriction sizes. All these observations justify why erosion is a little more pronounced in the sample corresponding to = 0.9 than in = 0.7 despite the fact that it is more stressed (see Figure 9 ). It should also be underlined that the very large proportion of rattlers found in the sample (64% of the particles) is linked to the fact that i) the used PSD is very wide ( max ∕ min = 10) and ii) the volume fraction of the small particles is quite large (Cu = 1.98). As a result the primary fabric of the considered material relies on few particles and many smaller particles fill the interstitial voids. As illustrated by the small travel distances observed in Figures 6 a and 6 b the rattlers remain trapped because of too small constrictions and collective clogging. In Figure 9 , when the hydraulic gradient is increased, the fluid forces applied on the particles are larger which explains why eroded tot increases faster for = 1 than for = 0.1. If particles are locally transported always in the direction of the maximum flow intensity (as assumed previously in the pore network transport models), the linearity of the PFV scheme (see Section 2.1) would result in similar transport and erosion whatever the intensity of the hydraulic gradient. In Figure 9 , this is obviously not the case as more erosion occurs when hydraulic gradients are more intense. This can result from the fact that fluid forces induce stress redistribution in the primary network of granular materials. The primary fabric of the material deforms which slightly change the constriction sizes. If the existing contact network is not resilient enough, fluid forces can even lead to force chain collapse and production of new rattlers. As no initially jammed particles (with coordination number > 3) become rattlers at some point in the six DEM/PFV simulations, the flow induced deformation of the primary fabric provides a rational basis to interpret the results in Figure 9 .

CONCLUSION AND OUTLOOK

Thanks to fully coupled DEM/PFV simulations, the occurrence of grain transport has been investigated at the REV scale in widely graded granular materials. By considering several hydraulic loadings (two directions and two intensities) and two stress states, the influence micromechanical features on the travel distances has been highlighted. The stress induced microstructure modifications are shown to influence grain transport as -for anisotropic stress states, force chains are built in privileged orientations. As a result, transport distances and erosion increase if the macroscopic flow direction is aligned with the principal force chains orientation;

-while subjected to drained triaxial tests, dense granular materials exhibit a contractive/dilative behavior. As they dilate, contact networks weaken and more rattlers are released in the pore space. This feature combined with larger pore sizes results in larger transport distances in dense granular materials subjected to large deviatoric stress.

By comparing the results obtained from simple transport models to fully coupled DEM/PFV simulations, the influence of collective grain motion and local pressure fluctuations has been highlighted. By over or under estimating these collective effects, the proposed models provide lower and upper bounds for grain transport. In addition, by modifying the local pressure maps used in the transport models, local pressure fluctuations around macroscopic pressure gradients have been shown to increase transport distances. This impact underlines the importance of using fully resolved fluid grain coupling schemes while modeling suffusion at the representative elementary volume scale. In this study, mechanical loading is applied prior to hydraulic loading and no post-erosion mechanical loading is considered. While considering the problem of suffusion in dams or dikes, it will be of primary interest indeed to analyze the fluid flow effect on the mechanical stability of the structure undergoing internal erosion. A proper analysis can be achieved thanks to the combined use of the second-order work criterion introduced by Hill (44) and directional analyses [START_REF] Gudehus | A comparison of some constitutive laws for soils under radially symmetric loading and unloading[END_REF]. Such a work is currently in progress.

FIGURE 1

 1 FIGURE 1Elasto-frictional contact law used in DEM simulations.

FIGURE 2

 2 FIGURE 2 Cumulative distribution of the mass passing fraction with respect to the particle sizes. Particle radii are normalized by the mean radius value mean = 2.10 -4 m. The generated sample in the isotropic state is shown.

FIGURE 5

 5 FIGURE 5 Illustration of the non-physical clogging induced by the downstream wall.

FIGURE 6

 6 FIGURE 6Shortest travel distances of all the particles of the initially unstable sample corresponding to = 0.9 subjected to a hydraulic gradient = 0.1 in the vertical (a) or horizontal (b) directions. Extracted particles are identified by red diamonds. The mean transport distances ̄ short num are shown in blue solid lines. In black solid lines the numerical predictions corresponding the transport model of Section 3.2 are shown (pressure map computed using the PFV scheme). Dashed and dot-dashed lines correspond to the predictions computed from the linear approximation of the macroscopic pressure gradient respectively on the pore network computed with all particles and on the pore network computed with only non free particles.

FIGURE 7

 7 FIGURE 7 Three different pore networks used to estimate mean travel distances. The network can be constructed either from all particles ((a) and (b)) or only from non free particles (c). The pore pressure map can be derived from the use of the PFV scheme (a) or from the linear approximation of the macroscopic pressure gradient ((b) and (c)).

FIGURE 8

 8 FIGURE 8 Hagen-Poiseuille flow velocity profile in a circular tube (a) and Stokes drag force (b). In both cases viscous flow conditions are assumed.

7 , 1 FIGURE 9

 719 FIGURE 9Dimensionless time evolution of the mass (corresponding to the mass of deleted particles) for the six DEM/PFV simulations considered.

FIGURE 10

 10 FIGURE 10 Strain evolution of the proportions of rattlers (bottom left) and chained particles (top right) during the first increments of the drained triaxial loading. The two states corresponding to ∈ {0.7, 0.9} are shown with diamond points. For these two states, visualization of the force chains are provided (top right) and the force chain definition of Peters et al. (43) is illustrated (top left).

TABLE 1

 1 Mechanical parameters used in the elasto-frictional contact law implemented in YADE.

	Parameters	Value
	Density ( )	3,000 kg.m -3
	Young Modulus ( )	356 MPa
	Stiffness ratio ( )	0.42
	Inter-particle friction angle ( ) 35 •
	Particle-wall friction angle	0 •
	Number of particles	5,000

TABLE 2

 2 Drained triaxial test response. The diamond points correspond to the mechanical states considered for the stability analysis. The six DEM/PFV simulations carried out.

					0.005
	0.8				0.004
					0.003
	0.6				0.002
	0.4				0.001
					0.000
	0.2				0.001
	0.00 0.0	0.01	0.02	0.03	0.00 0.002	0.01	0.02	0.03
	FIGURE 3 Stress	Hydraulic gradient	Flow direction
	state		= 0.1	= 1	Vertical ( ) Horizontal ( )
	= 0.7	x		x	x
	= 0.9	x	x	x	x