Modelling the impact of particle removal on granular material behaviour

Presented here is a numerical method able to predict the mechanical behaviour of granular materials subjected to particle removal. The purpose of this study is to improve understanding of the mechanical behaviour of soils subjected to internal erosion in hydraulic works. The approach is based on a homogenisation technique for deriving the stress–strain relationship of a granular assembly from forces and displacements at the particle level. The soil’s local behaviour is assumed to follow a Hertz–Mindlin elastic law and a Mohr–Coulomb plastic law. Sliding resistance on each plane is made to depend on the actual void ratio of the granular assembly. The solid fraction removal is modelled by progressively increasing the granular assembly void ratio, which provokes a decrease of the sliding resistance of each interparticle contact, leading to macroscopic deformations of the soil specimen. At elevated stress levels, large deformations can develop and lead to soil failure. Numerical simulations also demonstrate that a type of failure called diffuse failure can occur in eroded soil masses whenever an increase in pore pressure is generated within the soil. These numerical results appear to be coherent with observations made on embankment dams that have suffered internal erosion and, in particular, with the description of their modes of failure.


INTRODUCTION
Internal erosion occurs when particles are pulled off by seepage forces and transported downstream.Two fundamental types of internal erosion can be distinguished: piping and suffusion.This paper is concerned with the latter, which corresponds to the extraction of fine particles from a solid matrix.The mechanism of suffusion has been widely studied, mainly with the purpose of defining the criteria for soils likely to develop internal erosion (e.g.Kenney & Lau, 1985;Lafleur et al., 1989;Chapuis et al., 1996;Foster & Fell, 2001;Wan & Fell, 2004;Bendahmane et al., 2008).Several studies have also proposed various numerical approaches to model the erosion process (e.g.Vardoulakis et al., 1996;Cividini & Gioda, 2004;Vardoulakis, 2004).Even though failures of embankments by suffusion have been reported (Fry, 1997), the small number of studies on the mechanical behaviour of eroded soils (Sterpi, 2003;Muir Wood et al., 2010) show that this problem has been insufficiently addressed.Fig. 1 illustrates the different modes of failure associated with the process of internal erosion by suffusion.It has been recognised that the probability of failure of a modern dam due to internal erosion is ten times higher than the probability of failure due to sliding (Fry, 1997).With climate change, this type of accident is likely to occur more frequently.
The aim of this study is to use a numerical method to improve understanding of the mechanical behaviour of soils subjected to internal erosion.The phenomenon of real suffusion is not the subject of investigation.No coupling between water and solid particle is introduced.The present project is to analyse the impact of removing particles from the granular assembly.This has been done by progressively reducing a part of the solid fraction, and by analysing the consequences for the behaviour of a granular material.The approach is based on the use of a homogenisation technique deriving the stress-strain relationship of the granular assembly from forces and displacements at the particle level.The concept comes from Taylor and Budianski in their models for polycrystalline materials (e.g.Batdorf & Budianski, 1949).Similar approaches can also be found in models for soils and rocks (e.g.Calladine, 1971; and multilaminate models by Pande & Sharma, 1982), in models for concrete (e.g. the micro-plane model by Bazant et al., 1995), and in models for granular materials (e.g.Rothenburg & Selvadurai, 1981;Jenkins, 1988;Chang & Liao, 1990;Chang & Gao, 1995;Emeriault & Cambou, 1996;Nicot & Darve, 2005).The macroscopic behaviour of granular materials is governed by their properties at interparticle contacts.The basic idea is to view the packing as represented by a set of micro systems that correspond to the contact planes.The overall stress-strain relationship of the packing is obtained from averaging the contact plane behaviours.Along these lines, the author has developed a stress-strain model that considers the interparticle forces and displacements along a set of contact planes (Chang & Hicher, 2005).In this paper, the capability of the model has been extended to incorporate the structural changes due to particle removal.

STRESS-STRAIN MODEL BASED ON MICROMECHANICAL APPROACH
The microstructural model developed by Chang & Hicher (2005) views a soil as a collection of non-cohesive particles.The deformation of a representative volume of the material is generated by mobilising particle contacts in various orientations.On each contact plane, an auxiliary local coordinate can be established by means of three orthogonal unit vectors n, s and t.The vector n is outward normal to the contact plane.Vectors s and t are on the contact plane.
Interparticle behaviour Elastic stiffness.The contact stiffness of contact plane of orientation AE includes normal stiffness k AE n and shear stiffness k AE r : The elastic stiffness tensor is defined by which can be related to the contact normal and shear stiffness by The value of the stiffness for two elastic spheres can be estimated from the Hertz-Mindlin formulation (Mindlin & Deresiewicz, 1953).For sand grains a revised form was adopted (Chang et al., 1989), given by where G g is the elastic modulus for the grains; f n is the contact force in the normal direction; l is the branch length between two particles; and k n0 , k r0 and n are material constants.
Plastic yield function.The yield function is assumed to be of the Mohr-Coulomb type, defined in a contact-force space (e.g.f n , f s , f t ) where k(˜P) is a hardening/softening parameter.The shear force T and the rate of plastic sliding ˜p are defined as The hardening function is defined by a hyperbolic curve in the k-˜p plane, which involves two material constants: ö p and k p0 : Plastic flow rule.Plastic sliding often occurs along the tangential direction of the contact plane with an upward or downward movement: thus shear-induced dilation or contraction takes place.The dilatancy effect can be described by where the material constant ö 0 can be considered in most cases to be equal to the interparticle friction angle ö ì : On the yield surface, under a loading condition, the shear plastic flow is determined by a normality rule applied to the yield function.However, the plastic flow in the direction normal to the contact plane is governed by the stressdilatancy equation in equation ( 7).Thus the flow rule is non-associated.

Influence of void ratio on mobilised friction angle
Resistance against sliding on a contact plane depends on the degree of interlocking by neighbouring particles.The resistance can be related to the packing void ratio e by where m is a material constant.e c corresponds to the critical void ratio for a given state of stress.For dense packing e c /e is greater than 1, and therefore the apparent interparticle friction angle ö p is greater than the internal friction angle ö ì : When the packing structure dilates, the degree of interlocking and the apparent friction angle are reduced, which results in a strain-softening phenomenon.For loose packing, the apparent friction angle ö p is smaller than the internal friction angle ö ì , and increases during the material contraction.
The critical void ratio e c is a function of the mean stress applied to the overall assembly, and can be written as where ˆand º are two material constants, p9 is the mean effective stress of the packing, and (e ref , p ref ) is a reference point on the critical state line.

Local elasto-plastic relationship
With the elements discussed above, the final incremental stress-strain relation of the material can be derived that includes both elastic and plastic behaviours, and is given by A detailed expression for the elasto-plastic stiffness tensor is given in Chang & Hicher (2005).

Micro-macro relationship
The stress-strain relationship for an assembly can be determined by integrating the behaviour of interparticle contacts in all orientations.In the integration process, a micro-macro relationship is required.Using the static hypothesis, the relation between the global strain and interparticle displacement is obtained as where the branch vector l AE k is defined as the vector joining the centres of two particles, and the fabric tensor is defined as The mean force on the contact plane of each orientation is The stress increment can be obtained by the contact forces and branch vectors for all contacts (Christofferson et al., 1981;Rothenburg & Selvadurai, 1981), as Stress-strain relationship From equations ( 11)-( 14), the relationship between stress increment and strain increment can be obtained as where When the contact number N is sufficiently large in an isotropic packing, the summation of the flexibility tensor in equation ( 15) and the summation of the fabric tensor in equation ( 12) can be written in integral form, given by and The integration of equations ( 16) and ( 17) in spherical coordinates can be carried out numerically by using Gauss integration points over the surface of the sphere (Chang & Hicher, 2005).

MATERIAL PARAMETERS
One can summarise the material parameters as Except for the critical state parameters, all the parameters are linked to the interparticle relationship.The interparticle elastic constant k n0 is assumed to be equal to 61 000 N/mm.The value of k t0 /k n0 is commonly about 0 .4, corresponding to a Poisson's ratio of the soil í ¼ 0 . 2 and the exponent n ¼ 0 .5. (Hicher, 1996;Hicher & Chang, 2006).
Standard values for k p0 and ö 0 are k p0 ¼ k n0 and ö 0 ¼ ö ì (Chang & Hicher, 2005).Therefore, for dry or saturated samples, only five parameters have to be obtained from experimental results, and these can all be determined from the stress-strain curves obtained from triaxial tests.
The number of contacts per unit volume, Nl 3 /V, changes during the deformation.Using the experimental data by Oda (1977) for three mixtures of spheres, the total number of contact per unit volume related to the void ratio can be approximated through the expression This equation is used to describe the evolution of the contact number per unit volume.The initial contact number per unit volume can be obtained by matching the predicted and experimentally measured elastic moduli for specimens with different void ratios (Hicher & Chang, 2005, 2006).
The soil selected for this study was determined by the conclusions made by several studies on the erodability potential of soils.Soils that are prone to suffusion are said to be internally unstable.This instability is caused by the coexistence of coarse and fine grains inside the soil structure, corresponding to a soil mass that does not meet selffiltering conditions.These soils are clay-silt-sand-gravel or silt-sand-gravel mixtures.Different criteria have been proposed in order to relate the risk of suffusion to the grainsize distribution of the soil (see the synthesis made by Wan & Fell, 2004).Fig. 2 shows two grain-size distributions that are typical of soils prone to developing internal erosion: a coarse, widely graded soil and a gap-graded soil.The grainsize distribution can take into account the geometrical conditions that lead to potential erosion, which are related mainly to the capability of the fine particles to move freely within the pores of the coarse matrix.This requires the dimensions of the fines to be smaller than the conscription -that is, the opening connecting two pores -and also insufficient fine particles to fill up the voids of the coarse matrix completely.Beyond these geometrical conditions, the velocity of the flow has to be high enough in order to drag the fine particles through the conscriptions.
For this study, a soil made of non-cohesive particles was selected, typically a silt-sand-gravel mixture, with a grainsize distribution as one of those presented in Fig. 2. The material parameters were selected by using the correlations proposed by Biarez & Hicher (1994) relating the physical properties of a granular assembly to its mechanical properties.The maximum and minimum void ratios of a widely graded granular material are typically e max ¼ 0 .6 and e min ¼ 0 .25.The values can then be derived of the two parameters corresponding to the position of the critical state in the e-p9 plane: º ¼ 0 .05 and p ref ¼ 0 .01 MPa for e ref ¼ e max ¼ 0 .6.The friction angle at critical state, ö ì , is considered equal to 308.In equation ( 8), the value of m is taken equal to 0 .5, which is a typical value for this type of soil.The set of parameters for the selected soil is presented in Table 1.
Two different specimens have been selected for the study: a well-compacted material with an initial void ratio e 0 ¼ 0 .3, corresponding to a relative density D r ¼ 85%; and a lesscompacted material with an initial void ratio e 0 ¼ 0 .45, corresponding to a relative density D r ¼ 43%.Fig. 3 shows the results of the simulation of a triaxial test on each specimen subjected to an initial isotropic confining stress equal to 300 kPa.These results demonstrate the ability of the model to reproduce the main features of granular materials: contractive or dilative behaviour, depending on the initial material density; an increase in the maximum strength when the density decreases; and softening behaviour of a dense, dilative specimen.

NUMERICAL SIMULATIONS OF SOIL SUBJECTED TO SUFFUSION Modelling particle removal
What could be the impact of removing a part of the solid fraction on the behaviour of granular materials?Considering the homogenisation technique presented in the previous section, two major influences can be noted.First, the number of interparticle contacts, N, can be reduced, leading to a decrease of the normalised number of contacts per unit volume, Nl 3 /V.This decrease in the number of contacts will concentrate the contact forces on the remaining contacts, and therefore create more interparticle displacements.By integrating these micro-displacements all over the elementary volume, additional deformation of the granular assembly can be obtained.The second consequence is linked to what was previously referred to as the interlocking influence.The sliding resistance at each contact has two origins.One corresponds to the friction between two particles, defined by the internal friction angle ö ì ; the other is linked to a geometrical effect of the neighbouring particles, which will resist the sliding displacement of two particles in contact, and increasingly so if the packing is dense.Therefore, the pulling-off of a fraction of the particles will release the interlocking influence and allow the contacting particles to slide together.In the model, this interlocking influence is taken into account by equations ( 8) and ( 9), which relate the sliding resistance, controlled by the apparent friction angle ö p , to the void ratio of the assembly.So, when a part of the solid fraction is removed from the material, the void ratio increases, and as a consequence the sliding resistance decreases, which can create additional interparticle displacements, and a greater deformation of the granular assembly.
The model cannot take into account the size of individual particles; only the mean size value is considered as a model parameter.It is therefore not possible to model exactly the suffusion process that corresponds to the removal of the finest fraction.This study simply considers the removal of a given fraction of the solid particles.For this purpose the eroded fraction f e is defined as where W f is the weight of the eroded particles, and W s0 is the initial total solid weight per unit volume.If one can assume that the particle density is the same for any particle size, equation ( 19) can be written in terms of volume as where V f is the eroded solid volume, and V s0 is the initial solid volume per unit volume of soil.One can then relate f e to the removed solid fraction according to where V s is the solid volume after extraction.When no deformation takes place during the extraction of solid particles, the total volume remains constant, and the volume occupied by the extracted particles is replaced by the same volume of voids.Therefore it is possible to write Combining equations ( 21) and ( 22),  where e 0 is the initial soil void ratio before the extraction process begins, and (˜e) er is the change in void ratio due to the extraction process.
When the extracted fraction f e increases progressively, it creates a change in the void ratio.This change is taken into account in the model by using equations ( 8) and ( 9).If the material is subjected to a constant state of external stresses, the evolution of the sliding resistance introduced by equation ( 8) when the void ratio is changed creates a disequilibrium at each contact point, leading to local sliding.All the local displacements are then integrated to produce the macroscopic deformation of the soil specimen.This macroscopic deformation induces a volumetric change å v , and therefore a change in the void ratio, which is added to the void ratio change caused by particle removal The void ratio is calculated at the end of each step, which consists, first, in imposing a change in void ratio corresponding to an incremental increase of the eroded fraction f e , and then in calculating the induced deformations.
The change of void ratio at the end of each step can be higher or lower than that due to the sole removed solid fraction, depending on the sign of the volume change due to the induced deformations, as will be shown below.A contraction during this process will contribute to delay its impact, whereas a dilation will accelerate the deformation progress.
The change of grading could influence the position of the critical state line in the e-log p9 plane.Strictly speaking, this aspect should have been incorporated in the model, which would have been done by linking the value of e c in equation ( 8) to the grading change due to the increase of the eroded fraction.However, in the case of widely graded materials such as those prone to suffusion (Fig. 2), it was assumed that a limited removal of the fine fraction would not influence the position of the critical state.

Mechanical behaviour of eroded material during triaxial loading
Several numerical tests were performed on initially dense and loose samples.Each specimen was initially subjected to triaxial loading up to a given state of stress, and then to progressive particle removal while the external stresses were kept constant.Figs 4 and 5 present the results obtained for the dense and loose samples respectively.All the specimens are isotropically consolidated up to p9 ¼ 300 kPa and then sheared at different stress levels.It can be seen that the erosion-induced deformations are larger when the stress ratio is higher.For low stress ratios, the strain amplitude remains limited to values lower than 1% when the eroded fraction increases up to 12%.Under these conditions the soil remains stable, and the damage induced within the earth structure is limited.For higher stress ratios deformation increases much faster, and large deformations can occur for eroded fractions measuring more than 10% for the studied cases.If similar conditions are experienced by soil masses within the earth structure, heavy damage and even a general collapse can be the consequence of the suffusion process.
The critical stress level depends on the initial density of the soil.By comparing the results in Figs 4 and 5, it can be concluded that, at any given deviatoric stress, a denser sample is more resistant and less inclined to deform when particles are progressively removed.If the results are compared at the same ratio between the applied deviatoric stress and the maximum soil resistance, the difference in behaviour is slight; apparently the two samples develop a similar strain evolution.

COMPARISON WITH A DEM APPROACH
At present the author is not aware of the existence of well-documented experimental tests studying the impact of internal erosion on the mechanical behaviour of soil samples.To gain more confidence in the numerical simulations, it was decided to compare the results to other numerical results obtained by the discrete-element method (DEM) where the process of internal erosion was modelled by a progressive extraction of the finer particles from a granular assembly composed of spheres while the sample was maintained under constant stresses (Scholte `s et al., 2010).A similar procedure was used by Muir Wood et al. (2010) on assemblies of discs in a two-dimensional condition.The most erosion-prone particles were progressively removed from the assembly, leading to a rearrangement of the microstructure.The removal of a given particle left some unequilibrated interparticle forces in the assembly.Under constant external stresses, a deformation of the granular assembly, due to these unequilibrated forces, was observed.Upon reaching a new state of equilibrium, this process of particle removal was repeated.Two different extraction processes have been defined.
The first is based on particle sizes only, and consists in removing the smallest particle contained in the assembly, as also performed by Muir Wood et al. (2010).The second is based on both the particle's sizes and the degree of interlocking, estimated through the assessment of the particle's mean internal moment M p , which involves the removal of the less loaded of the smallest particles.
Two criteria were defined for terminating the particle removal: when the axial strain exceeded 25%, or when the eroded fraction became equal to 5% of the total solid mass.The specimen deformations during the erosion process are plotted in Fig. 6 for different constant stress levels.Even if this particle removal procedure does not correspond to any ε 1 : % f e : % q 400 kPa ϭ q 450 kPa ϭ q 500 kPa ϭ Fig. 5. Evolution of soil deformation with eroded fraction for different states of stress; initial void ratio e 0 0 .45 realistic erosion process, the description of its consequences for the granular assembly stability appears strikingly similar to what has been described in the previous section concerning the internal mechanism of erosion considered by the micromechanical model.It would therefore be interesting to see whether comparative results can be obtained by simulating the numerical tests performed by DEM.In this study, the normal stiffness and the ratio between tangential and normal stiffness are taken equal to the values considered by the DEM approach.The parameters of the plastic part of the contact law are given in Table 2. Parameters determining the critical state line of the DEM simulations are, as well as the value of m, calibrated by curve-fitting using the test results shown in Fig. 6.
Figure 7 presents the simulations done by the micromechanical model.Several erosion tests were simulated at various stress levels.The results demonstrate that, for elevated deviatoric stresses, the model gives similar results as the ones obtained by DEM.Progressive deformations develop within the specimens during the erosion process.For high stress levels, the axial strains increase rapidly and lead to specimen failure (Fig. 8).The strain path is linked to the stress ratio.The test at a stress ratio equal to ç cs ¼ 0 .72, corresponding to the value of q/p at critical state obtained from the drained triaxial test, deforms roughly at constant volume.At higher stress ratios the volume increases during erosion, whereas at smaller stress ratios it decreases.This behaviour is in complete agreement with the DEM simulations, as can be seen in Fig. 6.
For smaller deviatoric stresses, however, it seems that the DEM simulations produce larger deformations during erosion than those produced by the microstructural model.Indeed, rather large values of the eroded fraction f e are required in order to obtain a significant straining of the specimens (Fig. 8).
As in the DEM approach, in accordance with the increase of porosity induced by particle removals, degraded specimens behave like a loose material, with significant contraction and monotonic increase of the stress ratio.In particular, the internal friction angle decreases from 24 .18 for the intact medium to 19 .38 for the degraded ones, which corroborates well the shear strength reduction given by DEM simulations (Fig. 9).
The reduction in shear strength and the contractive behaviour of degraded samples can lead to an unstable response to external loading, as will be shown in the following section.

MATERIAL INSTABILITY AND DIFFUSE FAILURE
An important aspect of soil behaviour that has to be studied in connection with earth structure safety is instability.This phenomenon falls into two categories: material instability (also known as intrinsic/constitutive instability) and geometrical instability (e.g.Goddard, 2003).Strain localisation in a finite-size specimen can be studied as a boundary value problem for which initial inhomogeneities, as well as boundary conditions, play an important role.However, as has been shown by Rudnicki & Rice (1975), the occurrence of localisation can be predicted on a constitutive level.More recently, several studies have shown that other modes of instability can occur within granular materials.One of these, called diffuse failure, has been studied by Darve and co-authors (Darve & Roguiez, 1998;Laouafa & Darve, 2002;Darve et al., 2004).Nova (1994) has developed similar concepts.Experimental evidence supporting these theoretical approaches -as, for example, in loose sand under undrained conditions -shows that an unstable condition can be obtained at a low shear stress level and, subsequently, that the strength is reduced to almost zero, corresponding to a material state known as static liquefaction.
This particular mode of instability will be examined for the case of an embankment subjected to internal erosion.

Constant-q loading
This type of test consists of shearing the specimen to a prescribed stress ratio along a drained compression triaxial path, and then decreasing the mean effective stress while keeping the deviatoric stress constant.This stress path can simulate the loading condition of a soil element within a slope when pore pressure increases progressively.
Several investigations have demonstrated that instability can occur during a constant-q stress path (Sasitharan et al., 1993;Nova & Imposimato, 1997;Gajo et al., 2000;Lade, 2002;Chu & Leong, 2003;Darve et al., 2007).These studies observed that sand specimens experienced a sudden collapse for stress states located well below the critical state failure line when a decrease of the mean effective stress p9 was applied while q was kept constant.At a given point of the test, the axial strain rate started increasing very rapidly, and the deviatoric stress could no longer be kept constant.The test lost its controllability (as defined by Nova, 1994), in the sense that the imposed loading programme could no longer be maintained.
This phenomenon can be linked to Hill's sufficient condition of stability (Hill, 1958), which states that a material, progressing from one stress state to another, is stable if the second-order work is strictly positive Thus, according to Hill's condition, whether a material is stable or not depends not only on the current stress state but also on the direction of the stress increment.For loading conditions, such as those considered in this study, the expression for the second-order work can be rearranged to give For constant-q tests (dq ¼ 0), according to the above equation the second-order work is reduced to d 2 W ¼ dp9då v : Since the mean stress is being progressively decreased (i.e.dp9 , 0), the second-order work becomes negative if and only if då v .0 (i.e. the volume contracts).Thus the onset of instability corresponds to the minimum of the p9-å v curve.Indeed, experimental results on loose sand specimens have shown that the volumetric strain passes through a minimum (då v ¼ 0), which corresponds to the vanishing of the secondorder work.At that point, any small perturbation imposed on the specimen could provoke its sudden collapse.

Numerical simulations of constant-q tests on intact soil
The material parameters in Table 1 were used to predict the results of constant-q tests.A typical result is presented in Fig. 10.The test was performed on a loose specimen (e 0 ¼ 0 .6) consolidated up to an isotropic stress p9 equal to 300 kPa, and then loaded up to a deviatoric stress q equal to 300 kPa.The deviatoric stress q was then kept constant 9. Mohr-Coulomb failure envelopes before and after particle extraction while the effective mean stress p9 was progressively decreased.Only the constant-q part of the test is presented.The initial part of the p9-å v curve shows that, as the mean stress p9 decreases, the volume increases.This trend continues up to a certain point, where the volume starts to decrease.
A number of similar numerical tests were performed with different values of q and different initial void ratios e 0 : All the stress states corresponding to the onset of instability are plotted in the p9-q plane (Fig. 11).One can see that, for any given initial void ratio, the instability criterion corresponds to a straight line starting from the origin of the axes.
The slope of the instability line increases when the void ratio decreases.For void ratios smaller than 0 .4, the instability condition cannot be reached ahead of the plastic limit condition.The position of the plastic limit depends on the material void ratio.For dense specimens this is located above the critical state line, which corresponds, for this material, to a stress ratio M ¼ q/p9 ¼ 1 .2 (i.e. a friction angle ö ¼ 308 at critical state).The results of constant-q tests on dense specimens (e 0 ¼ 0 .3) are shown in Fig. 12.One can see that the decrease of the mean effective stress is accompanied by an increase of the volumetric strain, which corresponds to a positive value of the second-order work.This condition persists all along the constant-q stress path.At the end of the test, the axial strain starts to increase rapidly, leading to the failure of the soil specimen.
This failure condition is obtained when the stress state reaches the plastic limit, as shown in Fig. 13.In this case, a classic mode of failure in dense granular materials is obtained, above the critical state line.This type of result agrees with the experimental results obtained from tests on sand.

Numerical simulations of constant-q tests on eroded specimens
Now consider a specimen that is unconditionally stable under the considered loading condition, which means that its void ratio is initially low enough to prevent an unstable state from developing along constant-q tests, as shown in Fig. 13.The specimen is then subjected to suffusion, leading to an increase in its void ratio.One can see in Fig. 11 that this increase in the void ratio, without a change in the stress state, can cause a condition of stress that could lead to instability.To illustrate this, a specimen was selected with an initial void ratio e 0 ¼ 0 .3, and was placed under conditions similar to those defined previously: a mechanical loading that corresponds to a triaxial loading with a confining stress equal to 300 kPa, and a deviatoric stress q ¼ 300 kPa.It was then subjected to particle removal, which led to an increase in its void ratio.Two cases were examined: the first with an increase up to e 1 ¼ 0 .45 ( f e1 ¼ 10 .5%), and the second with an increase up to e 2 ¼ 0 .55 ( f e2 ¼ 16 .5%).For dense materials, large deformations will develop during the process of removing part of the solid phase only at much higher deviatoric stress levels (see Fig. 4).Therefore the material remains stable so long as the effective stress state does not evolve.Now assume an increase in pore pressure due to changes in hydraulic conditions (e.g. rise of the water table, or heavy rain).The two eroded specimens are therefore subjected to a decrease of the mean effective stress p9 at constant deviatoric stress q.The results are presented in Fig. 14, showing the volumetric change with the mean effective stress p9.One can see that an instability condition will result in both cases.For the less eroded material, a pore pressure increase of 120 kPa is needed, whereas for the more eroded one this increase is only 65 kPa.When instability sets in (zero slope, i.e. ˜åv ¼ 0), the axial strain will start to increase rapidly, indicating that the material is leading to a diffuse failure state (see Fig. 14).These results can be compared with those obtained on intact samples subjected to the same constant-q Fig. 12. Axial and volumetric strain evolution during constant-q tests on dense specimens (e 0 0 .3) Fig. 13.Failure line for constant-q tests on dense specimens (e 0 0 .3) loading (Figs 12 and 13).One can see that large deformations in eroded samples tend to develop for a much smaller decrease of the mean effective stress p9.This difference in behaviour corresponds to the difference in the slope of the plastic limit of the intact samples and the slopes of the instability lines of looser samples (Fig. 11).
A second example can be considered, in which the deviatoric stress q has been changed from 300 kPa to 400 kPa.Under this condition, the less eroded specimen reaches an unstable state for a pore pressure increase of 53 kPa, which is considerably lower than in the first case (Fig. 15).The more eroded specimen is a very interesting example, because the erosion process causes it to act inside a domain where it is unstable in the loading direction corresponding to a constant-q test and decreasing mean effective stress p9, as shown in Fig. 10.Under these conditions, the axial strain will develop very rapidly as soon as the mean stress p9 starts to decrease.Therefore any perturbation in the hydraulic conditions, leading to the slightest increase in pore pressure, will cause the soil to be unstable, which in turn can lead to a collapse of the earth structure.
This study has focused on constant-q stress paths, but it can be extended to other loading conditions, such as, for example, undrained loading.As demonstrated by Darve et al. (2004), the condition of instability for constant-q tests coincides with the condition for undrained tests, so the stress states in the p9-q plane corresponding to the vanishing of the volumetric strain increment during constant-q tests are the same as the stress states corresponding to the peaks of the stress-strain curves for undrained tests.Chang et al. (2009) showed that the instability line in the p9-q plane obtained by the micromechanical model for a given material at a given density is the same for both constant-q and undrained loadings, in agreement with experimental testing on the same material.Therefore, in the second example described above, if a slight perturbation is applied to the material in undrained condition, it will fail.This was the case in the Dycho ´w Dam accident in 1997.The downstream part of the dam lost its stability.The analyses provided by various experts concluded that the accident was caused mainly by the liquefaction of a silty sand layer due to the vibrations produced by a heavy vehicle's passing on the road at the top of the embankment (Dłuz ˙ewski & Hrabowski, 1998;Wolski et al., 1999).This silty sand layer was subjected to suffusion from the initial operation of the dam in 1940.The erosion created a loose material structure, as demonstrated by in situ penetrometer testing performed 2 years prior to the accident, which suggested that the relative density of the soil was around 15-20%.This soil layer was therefore in a domain where it was already unstable, as discussed previously, and the disturbance caused by the vehicle was sufficient to trigger its liquefaction, which produced the collapse of the embankment.

CONCLUSION
Based on a microstructural approach, this study has examined the influence on specific soil properties of removing a part of the solid fraction.The material was selected for its typical grain-size distribution in the family of soils likely to develop internal erosion; the parameter values typical for this kind of material were assumed, using correlations between the micro and macro properties of granular materials.
To obtain an elasto-plastic stress-strain relationship for granular materials, a homogenisation technique was used, which consisted in mobilising particle contacts in various orientations to obtain the deformation of a representative volume of the material.Thus a stress-strain relationship could be derived as an average of the behaviour of these local contact planes.The local behaviour is assumed to follow a Hertz-Mindlin elastic law and a Mohr-Coulomb plastic law.The friction angle on each plane is made dependent on the actual void ratio of the granular assembly.On the whole, the model requires only a limited number of parameters, which can easily be determined from conventional triaxial testing.
The progressive removal of the solid fraction leads to a decrease of the sliding resistance of each interparticle contact, which creates a disequilibrium between the external applied loading and the internal contact forces.As a consequence, local slidings occur, which lead to macroscopic deformations of the soil specimen.The amplitude of the induced deformations depends on the quantity of particles removed, as well as on the applied stress level.At elevated stress levels, large deformations can develop when the removed fraction increases.A comparison with numerical results obtained by DEM showed that similar trends of behaviour could occur when progressive removal of particles developed within the granular material.
A mechanism called diffuse failure can occur in eroded soil masses.Granular assemblies can develop instability at a shear stress level much lower than the critical state failure line.The appearance of an unstable state can be analysed by taking into consideration Hill's sufficient condition of stability along a given stress path.Numerical simulations of constant-q loading, corresponding to a pore pressure increase within the soil mass, show that the condition of instability can be reached when the removed solid fraction is high enough.These results demonstrate that changes in the hydraulic conditions within the embankment dam can lead to a sudden collapse of the earth structure due to the existence of unstable states within the eroded soil mass.
This study is an attempt to understand the behaviour of Fig. 15.Axial and volumetric strain evolution during constant-q tests on eroded specimens (q 400 kPa); arrows show sign changes for volumetric strain, corresponding to vanishing of second-order work soils subjected to internal erosion.Even if the solid fraction removal process used here does not reproduce the exact conditions experienced by the soil during suffusion, the numerical results obtained appear to agree with observations made on embankment dams that have suffered internal erosion, and in particular with the description of their modes of failure.More experimental studies are necessary, however, to validate the hypotheses underlying the modelling developed in this study.In particular, laboratory experiments at the specimen size would be highly useful in order to confirm the conclusions of this work not only qualitatively but also quantitatively, and would allow the model to be calibrated.

ACKNOWLEDGEMENTS
The financial support for this study came from the French Ministry of Ecology within the ERINOH and LEVEES projects, and from the Region Pays de la Loire within the project EMERMOD.

Fig. 2 .
Fig. 2. Grain-size distribution of soils prone to suffusion (after Wan & Fell, 2004): (a) soil with continuous gradation and a given percentage of fines; (b) gap-graded soil Fig. 3. Triaxial tests on selected soil at two different initial void ratios (ó 3 300 kPa)

Fig. 4 .
Fig. 4. Evolution of soil deformation with eroded fraction for different states of stress; initial void ratio e 0 0 .3

Fig. 6 .
Fig. 6.Summary of simulated loadings applied to DEM model (Scholte `s et al., 2010).Thick line: response of intact specimen to triaxial compression.Thin lines: responses of assembly when subjected to particle extraction for different levels of q and to subsequent triaxial compressions when eroded samples have stabilised.Symbols indicate transition between extraction and compression processes

Fig. 7 .
Fig. 7. Microstructural model simulations of particle extraction for different mobilised friction ç q/p following triaxial compressions realised on initial dense granular specimen (thick black line).Solid thin lines: extraction processes leading to stable configurations after entire 5% mass extraction (located with symbols).Dashed thin lines: extraction processes leading to sample failure

Fig. 8 .
Fig. 8. Evolution of vertical deformation during erosion at different stress levels: (a) DEM and (b) micromechanical model simulations

Fig. 14 .
Fig.14.Axial and volumetric strain evolution during constant-q tests on eroded specimens (q 300 kPa); arrows show sign changes for volumetric strain, corresponding to vanishing of second-order work

Table 1 .
Model parameters for selected soil

Table 2 .
Model parameters for DEM assembly of spheres

NOTATION
A ik fabric tensor C ijmp flexibility tensor D r relative density e packing void ratio e c critical void ratio e max maximum void ratio e min minimum void ratio e ref reference critical void ratio e 0 initial void ratio ˜e change in void ratio (˜e) er change in void ratio due to extraction process f e eroded fraction f AE contact force in the contact plane in the orientation AE f n contact force in normal direction f s contact force in direction s in sliding plane f t contact force in direction t in sliding plane G g elastic modulus for the grains k AE component of the tangential plastic displacement in direction t å d deviatoric strain å ij component of the strain tensor å v volumetric strain å 1 major principal strain ç stress ratio ç cs critical stress ratio k hardening function k(˜p) hardening/softening parameter º slope of the critical state line í Poisson's ratio ó stress ó ij component of the stress tensor _ ó ij stress increment ô shear stress ö p apparent interparticle friction angle ö 0 material constant ö ì interparticle friction angle