Micromechanical assessment of an internal stability criterion

The internal stability of a soil is a measure of its susceptibility to suffusion and suffosion, two forms of internal erosion. The internal stability of granular filters must be carefully considered when designing new embankment dams and assessing the risk associated with existing embankment dams. Current guidelines for assessing the internal stability of such filters were empirically derived from macroscale observations and consider the shape of the particle-size distribution curve. These guidelines lack a fundamental, scientific micromechanical basis. The initiation and propagation of internal erosion is clearly a particle-scale phenomenon, and this paper applies particulate mechanics to provide a micromechanical justification for one currently used stability criterion. The study used discrete element simulations of idealised virtual soil samples that had various degrees of internal stability when assessed using the criterion proposed by Kézdi [10]. The internal topologies of stable and unstable samples were analysed by considering the distributions of inter-particle contact forces, the number of particle–particle contacts and the average particle stresses. Clear correlations are observed between the filter stability criterion and the average number of contacts per particle and the probability that a given particle participates in stress transmission. The phenomenon of a critical fines content, at which the existing guidelines are no longer considered to be valid, is also considered.


Introduction
Internal erosion presents a significant threat to the stability of embankment dams. There are several distinct mechanisms of internal erosion, which include backward erosion, concentrated leaks or suffusion [4]. Internal stability refers to the geometric capacity of a granular material to prevent its fine fraction migrating through the void spaces between the larger particles under the action of seepage flow. ''Suffusion'' and ''suffosion'' are two similar, but distinct, mechanisms of internal erosion that are encountered in internally unstable soils [15]. When suffusion occurs, seepage flow carries the fine particles out of the soil without volumetric strain or settlement. There may be increases in local permeability at the point where the fine material is eroded and possible reductions in permeability downstream where the fines come to rest. When there is suffosion, the loss of fine particles under seepage load causes a change in the overall soil volume. Suffosion leads to local changes in permeability, and it can also lead to settlement with consequent loss of freeboard, the appearance of sinkholes and potential dam failure. Both mechanisms lead to a change in soil grading, which will have consequences for material behaviour, for example, by altering the critical state line [16] or the available peak strength [22].
Soils that are susceptible to either suffusion or suffosion are described as ''internally unstable''. When a soil is identified as being internally unstable, it does not always follow that either suffusion or suffosion will actually take place. A critical hydraulic gradient must also be applied to initiate the mechanism. Skempton and Brogan [25] showed that the critical hydraulic gradient for suffusion failure in internally unstable soils is considerably lower than the value of 1 typically considered necessary for surface heave in internally stable soils.
The initiation and propagation of suffusion and suffosion are particle-scale mechanisms. The current guidelines that use the particle-size distribution (PSD) to identify whether a soil is stable or unstable are empirical and based upon macroscale observations. This paper aims to identify, in a scientific manner, the differences in the internal structures of materials with different levels of stability. Idealised virtual samples of soil with differing PSDs were created to allow consideration of various levels of internal stability as assessed according to existing empirical guidelines. Following a general discussion of the criteria used to assess filter stability, the paper describes the discrete element method (DEM) simulations used and the analysis approaches adopted. The key differences in the material structure are then outlined.

Background
Current empirical guidelines to assess the internal stability of a filter consider the shape of the PSD by mass curve over a discrete interval [11]. No account is taken of the material fabric or structure, the particle shape or the void ratio. Two of the most commonly used methods are those proposed by Kenney and Lau [8,9] and by Kézdi [10], subsequently advocated by Sherard [24]. Both criteria were developed for the assessment of soils with non-plastic fines. Consequently, the effect of cohesion on internal stability is not considered here.
In this paper, the internal stability of soils is quantified using the (D 0 15 =d 0 85 ) max ratio, as proposed by Kézdi [10]. Using this approach, the PSD is divided into coarse and fine fractions. The D 0 15 of the coarse fraction is compared with the d 0 85 of the fine fraction (as illustrated in Fig. 1).
The ''fine'' and ''coarse'' fractions are selected systematically splitting the PSD at different points along the curve (or at the gap in the case of a gap-graded soil). A soil is considered to be internally stable if (D 0 15 =d 0 85 ) max \ 4. This approach for assessing stability was developed by extending Terzaghi's criterion for protection of a base soil by a filter to a single soil that could be considered to be comprised of coarse and fine parts. It is important to note that the ratio (D 0 15 =d 0 85 ) max = 4 was empirically chosen to provide safe filters for new dams and does not represent a clear transition between stable and unstable states. For example, in a series of permeameter tests, Fannin and Moffat [3] showed that the fines lost under seepage increased progressively with (D 0 15 =d 0 85 ) max , but that no clear transition between stable and unstable states could be identified. Interestingly, comparisons of characteristic grain size for coarse and fine fractions are used elsewhere in the soil mechanics of mixed soils, for example, to define the inter-grain contact parameter b for the assessment of fines contribution to shear strength [18].
Kézdi [10] provided an analytical micromechanical justification for this rule, by considering uniform spheres with an ideal packing. This is a very idealised system that differs significantly from a real soil. A valid micromechanical analysis would need to consider a material with a range of particle sizes and a heterogeneous, random packing.
Several researchers have considered, hypothetically, the internal stability of soils from a particle-scale perspective. For example, Kenney and Lau [8] proposed that the following three conditions apply to an internally unstable soil: 1. A primary fabric of particles, which transfers stresses. 2. Loose particles, present in the voids between the primary fabric, which do not carry effective stress and can be moved by seepage. A soil in which the voids between the primary fabric are not filled with particles is described as ''under-filled'' [30]. 3. Constrictions within the primary fabric are large enough to allow the loose particles to be transported from void to void by seepage.
These conditions were defined based on an intuitive understanding of the system, rather than particle-scale observations. Kenney and Lau [8] proposed a critical fines content, above which the fine particles completely fill the voids between the coarse particles. In this case, the soil will be predominantly supported by the fine particles, which will participate in load transfer and will not be susceptible to erosion. Kenney and Lau state that even for very loose soils, when the fines content exceeds 30 %, the soil will be inherently stable. Skempton and Brogan [25] stated that for fines contents above 35 %, the coarse particles will not be in contact with one another, and the soil properties will be  a soil is considered to be internally stable dominated by the fine fraction. Vaughan [30] described this as an ''over-filled'' fabric.
Prior discussion on the particle-based justification for the grading criteria has been largely hypothetical. The detailed measurements of the structure needed to give a rational explanation as to why the gradings are stable or not have not been made. DEM is well suited to providing the necessary data to analyse the influence of grading on the internal topology of a particulate material. The application of DEM in geomechanics was originally proposed by Cundall and Strack [2], and its use in basic soil mechanics research is now well established (e.g. [1,17,27]).

Simulation approach
A series of three-dimensional DEM simulations on cubical assemblies of spheres were performed using the commercial code PFC3D [6]. Periodic boundary conditions were applied in the analyses to create a sample that is effectively of infinite size [27] and is free from the boundary effects associated with rigid boundaries [12], allowing the particle numbers to be kept to a reasonable level.
The input parameters common to all simulations are presented in Table 1. The same parameters were used for sample preparation and in the final state. The Hertz Mindlin contact model was used [7,13]. The particles were initially generated in the central control cell in random, non-overlapping locations. The result was a high porosity system that was not percolating, that is, as there were no contacts in the system, it could not transmit stress. To create a more realistic, dense packing, the radius expansion technique was employed [19], whereby particles were expanded and then allowed to rearrange and come into static equilibrium before another cycle of expansion was undertaken. This process continued until a percolating system at a mean normal stress level of about 200 kPa was obtained. The expansion factor by which particle radii were increased started at a fairly high level and was reduced in a controlled manner with each cycle of expansion. This allowed the target stress level to be reached fairly reliably without simulation times becoming excessive.
For these boundary-free simulations, the average sample stresses (r ij ) were calculated using a weighted average of the particle stresses.
where V is the sample volume, N p is the total number of particles, r ðpÞ ij is the average stress tensor within particle p and V (p) is the volume of particle p. For each particle, the average particle stress was calculated by summing the contact forces acting on it and by averaging the particle volume, as follows: where N c,p is the number of contacts on particle p, x i p is the particle centroidal position and F j c is the force transmitted at contact position x i c . Details on the calculation of r ij and r ðpÞ ij are given by Potyondy and Cundall [19].
Once the mean normal stress reached a value that just exceed 200 kPa, the expansion process was terminated. In order to ensure that the system was stable at the end of the simulation, a further check was carried out to ensure that the coordination number (i.e. the average number of contacts per particle) was stable and remained constant for at least 20,000 cycles. It was difficult to control the simulations to terminate at exactly 200 kPa; however, a comparison of data obtained in a number of trial simulations revealed that there were no measurable changes in the overall material microstructure when samples at 200 and 300 kPa were compared (e.g. for Sample G1-20, the connectivity distributions, that is, the distribution of the number of contacts per particle, were identical at 220, 257 and 290 kPa). Consequently, while the exit stress levels varied between 219 and 255 kPa, as given in Table 3, a valid comparison between the samples can be made. Note that the same coefficient of friction (l = 0.3) was used in all simulations to allow for meaningful comparison. One could argue that it would be appropriate to target the use of a consistent relative density; however, it is difficult to achieve a relative density value in DEM simulation that is directly comparable with the definition of relative density used in laboratory testing. Furthermore, prior experimental studies of internal stability have typically not considered the influence of relative density.
No gravitational force was applied to the particles during the simulations. The process created very ideal systems, made up of spherical particles, which had an isotropic fabric and were subject to an isotropic stress state. The criterion for assessing material stability considers only Particle density q kg/mm 3 2.6 9 10 -6 Coefficient of inter-particle friction l -0.3 the PSD and does not account for particle shape or material anisotropy, and so this simulation approach is appropriate for a fundamental analysis of criteria that are based solely on PSD. The samples were also relatively homogeneous. Note that prior experimental research (e.g. [14]) has also used ideal spherical ballotini to assess internal stability while isolating the effects of particle shape from PSD. Table 3 includes data on the per cent by mass of loose fines (C \ 2) in each sample considered. No clear correlation between these data and the Kezdi criterion is evident, highlighting the need to use particle-scale measures when examining this problem.
For each sample created, the connectivity, C, of the individual particles (i.e. the number of contacts each particle has with its neighbouring particles) was considered. The average connectivity equals the coordination number, Nc. To account for the fact that in widely graded soils there are an overwhelming number of fine particles that make up only a small mass fraction, and the volume-weighted coordination number, Nc vol , proposed by Shaebani et al. [23] was also used. This is defined as is the connectivity of particle p. Note that the form of this definition is consistent with the stress calculation (see Eq. 1).
To analyse the variation in the contact force magnitudes within the sample, normalised contact forces, f, were considered. These are defined as f ¼ N=hNi where N is the contact force and hNi is the average contact force. These data are presented in probability density function (PDF) plots, which have been used in prior research into the fundamental behaviour of granular materials [20,27]. The macroscale parameters void ratio and average stress level were also measured.

PSDs analysed
Normalised PSDs by mass for the 8 samples considered are presented in Fig. 2 Table 2, accompanied by an assessment of the internal stability of each sample according to the Kézdi criterion. Series G1 and G2 are gapgraded samples. The gap ratios (defined by the ratio of the diameter of the smallest coarse particle to the largest fine particle) are the same as those of the G1 and G2 gradings tested by Honjo et al. [5]. The number following the series name signifies the fines content by percentage mass (e.g. G1-20 has 20 % fines by mass). The G1 series has a gap ratio of 2.8. The empirical stability of the G1 samples increases with fines content. G1-10 is classified as unstable, G1-20 as borderline stable and G1-30 and G1-40 are both classified as stable. The grading G2-20 has a small gap ratio of 1.8 and is classified as stable according to the Kézdi [10] criterion. G1-30 and G1-40 have fines contents that are, respectively, close to and above the critical fines content proposed by Kenney and Lau [8].

and in
The grading A is similar to an unstable gap-graded sample that was tested by Skempton and Brogan [25]. However, while the gap ratio (3.6) and fines content (15 %) are similar to Skempton and Brogan's sample, the coarse fraction used in this study is much more uniform, and the largest coarse particles are therefore smaller. This was done to reduce the difference in size between coarse and fine particles, thereby allowing a more representative number of coarse particles to be used while keeping computational costs reasonable. The sample, therefore, has a significantly lower Kézdi (D 0 15 =d 0 85 ) max than Skempton and Brogan's original sample. This highlights that relatively small changes in PSD can lead to large changes in empirical stability. The Dam Filter grading matches that of material collected from an embankment dam currently under construction in southern Europe. This has a high level of empirical stability. The final sample presented has a Uniform grading and a high level of empirical stability.
Due to the large differences in diameter between the fine and coarse soil fractions in an internally unstable soil, a   [25] large number of DEM particles are required to ensure a representative number of coarse particles. For example, for grading A, which has the largest gap ratio, each particle of diameter d 90 has a volume equivalent to 250 particles of diameter d 10 . The number of particles in the samples generated ranged from 10,000 to 62,000 as indicated in Table 2.

Results
The void ratio and stress level of each sample at the end of each simulation are presented in Table 3, and coordination number is presented in Table 4. The minimum void ratio was found at a fines content of approximately 30 %. This is similar to experimental results, for example, [29]. However, as pointed out by Kenney and Lau [9], void ratio alone is not an accurate criterion for assessing internal stability.

Connectivity distributions
A number of samples have coordination numbers that are less than one. This is due to the need to neglect gravity in the periodic cell model, and so unconnected particles can exist when the system is in an equilibrium state. Were gravity to be enabled and a rigid base boundary added in lieu of periodic boundaries, the minimum coordination number would be 1 as the unconnected particles would fall to rest on other particles. However, those particles would still not participate in stress transmission. Neglecting gravity allows easy identification of particles that are not transmitting stress. The coordination number generally increases with increasing empirical stability. The very low coordination numbers show that a large number of particles are loose and potentially mobile. However, this does not necessarily mean that there is a large mass of erodible particles for two reasons: (1) the loose particles are overwhelmingly small, so although large in number, they are a relatively small proportion of the overall mass (as shown in Table 3); (2) the loose particles must also be small enough to pass through the void network, which is not considered in this analysis. The coordination number is also related to the fines content of the grading. All samples had mechanical coordination number values (average number of contacts per particle, considering only particles with two or more contacts) [27] greater than 4 (see Table 4), a requirement for mechanical stability. (This requirement applies only to perfectly spherical, frictional spheres, for example, were rolling resistance to be included in the simulation, the number of contacts required for stability would be reduced.) Sample G1-20 is illustrated in Fig. 3a, b. The coarse and fine fractions are distinguished by colour (the coarse particles are coloured in light grey, while the fine particles are coloured black). Figure 3a illustrates the entire sample, while Fig. 3b illustrates the particles with connectivity of two or more. A connectivity of two or greater was chosen to represent particles that are capable of transferring stress. As is evident from Fig. 3b, the majority of the coarse particles have a connectivity that exceeds 2 and form a continuous primary fabric. The remaining particles are insufficient to fill the void spaces in between and sit suspended and contactless within this void space, that is, the soil has an under-filled grading.
A relative frequency plot of the sample connectivities is illustrated in Fig. 4a. Figure 4b, c, respectively, plot the probabilities of a particle having a connectivity of zero, that is, completely loose particles, and a connectivity greater than or equal to four, the number of contacts required for mechanical stability of frictional spheres [31], against the Kézdi (D 0 15 =d 0 85 ) max criterion. In Fig. 4a, b, it is shown that the three gradings classified as unstable or borderline unstable by the Kézdi criterion (A, G1-10 and G1-20) show the highest probability of a particle in the sample being unconnected and only a very small  probability of a given particle having a connectivity C [ 1.
Here, the primary, stress transmitting fabric tends to be formed by a small number of coarse particles. As the empirical stability of the gradings increases, the probability of a particle being unconnected generally reduces, and there is a corresponding increase in the probability of a particle having a connectivity of C [ 1. The most stable Uniform grading has the highest probability of a particle having a connectivity greater than or equal to 4 (Fig. 4c). This is indicative of a structure that is relatively regular, with many particles participating in stress transmission and therefore not readily susceptible to erosion. A similar pattern of connectivities can be seen for the Dam Filter material, which is also highly stable. Referring to Fig. 4c, the general trend is that the probability of a particle having C C 4 increases systematically with increasing stability. A microscale boundary between the samples that are assessed as empirically stable and unstable can be identified at approximately P = 0.075. G1-40 has a coordination number and a connectivity distribution similar to those obtained for the highly stable samples, despite having a significantly higher (D 0 15 =d 0 85 ) max . The fines content is high enough to ''over fill'' any voids between the coarse particles. In this case, the fine and coarse particles both contribute to form a stable, wellconnected fabric. Figure 3c, d illustrates the particle packing G1-40 and shows that many of the fine particles have a connectivity of at least 2. As they are surrounded by fine particles, the coarse particles in G1-40 have extremely high connectivities (going beyond the scale used in Fig. 4).
Sample G1-30 has an inhomogeneous fabric. This emerges during the initial placement of the particles in the periodic cell. The implementation of the specimen generation algorithm within PFC3D gives a sample with a greater density of fine particles in the voids between the coarse particles close to the boundary. Whereas in some areas of the sample, the fabric is similar to that of an overfilled material, in other areas, an under-filled fabric exists. While this is a limitation of the simulation approach used, similar magnitudes of heterogeneity would be observed in physical tests or in the field. Fig. 4, where it can be seen that the connectivity distribution lies between the most stable samples and the most unstable samples. The loose fines in the under-filled areas would be susceptible to erosion from the individual voids. However, the overall soil mass would be likely to self-filter, as fines eroded from under-filled areas would be retained once they came into contact with the stable, over-filled areas.

G1-30 shows behaviour that is intermediate between the under-and over-filled gradings. This intermediate behaviour is reflected in
A plot of coordination number, Nc, against Kézdi (D 0 15 =d 0 85 ) max is presented in Fig. 5a, and a tentative trend line relating Nc and (D 0 15 =d 0 85 ) max is proposed. There is a general decrease in coordination number with instability.
A reduction in gap size leads to an increase in connectivity, with G2-20 having a higher coordination number than G1-20. There is a significant increase in coordination number at the transition between under-filled to over-filled fabric. The over-filled sample G1-40 has a coordination number similar to those for the highly stable uniform samples. The partially over-filled sample G1-30 has a coordination number between the under-filled and over-filled states.
The volume-weighted coordination number, Nc vol , is plotted against Kézdi (D 0 15 =d 0 85 ) max in Fig. 5b. The stable samples generally have higher values of Nc vol . A tentative dividing line between samples that are stable and unstable at the microscale is drawn on Fig. 5b at Nc vol = 4, which shows agreement with macrostability, defined by (D 0 15 =d 0 85 ) max \ 4. An over-filled fabric, distinct from a uniform fabric, can be easily identified using Nc vol , as shown by the very high value for G1-40 (Nc vol = 16.93) and the intermediate value for the partially over-filled G1-30 (Nc vol = 8.08).

Contact force distributions
From a particle-scale perspective, stress in a granular material is transmitted via a heterogeneous network of particle contacts. ''Chains'' or ''columns'' of highly stressed particles are connected via ''strong'' contacts, that is, contacts where f [ 1. Thornton and Antony [26] demonstrated that these strong contact forces carry almost all the deviatoric load through a sample, whereas ''weak'' contact forces, with f \ 1, generally contribute only to the mean stress. It is therefore reasonable to assume that a load carrying primary matrix would be represented by f [ 1, and poorly confined particles by f \ 1. A PDF plot of the normalised contact force values is presented for each sample in Fig. 6. The diagram only considers active compressive contacts, and so unconnected particles are not included. The most stable samples (Uniform and Dam Filter), in which the stress transfer is distributed well between particles, have PDF plots that show a high probability of contact forces close to the mean and a relatively low probability of very high or very low contact forces.
The most unstable samples, A and G1-10, have PDF plots that are remarkably similar to those obtained for the most stable samples. The fines contents of these samples are sufficiently low for there to be almost no force chains formed entirely of fine particles in the voids between the coarse particles; almost all the fines have connectivities of 0 or 1. As unconnected particles do not contribute to the PDF plots, the data presented represent only the primary matrix. The similarity of the PDF plots for highly stable and highly unstable materials is a drawback to their use for assessing internal stability. However, were gravity to be turned on, the loose fines would come to rest on other particles and would generate a very small contact force. In this case, it is postulated that the PDF plots for highly unstable samples would resemble that of G1-20, but with a more exaggerated peak at f = 0.
Samples G1-20 and G2-20 have an irregular distribution of contact forces. The fines content is high enough that fine particles are able to form a few continuous, low-magnitude force chains. The peak at f = 0 indicates contacts transmitting very low forces of 0 \ f \ 0.5, which form between poorly confined fine particles and leave them susceptible to erosion.
Sample G1-30 is notable as it has the highest probability of a contact transmitting either a very low or very high force. This can be attributed to the heterogeneity of the fabric. Within the sample, there are voids that are overfilled with fine particles, and the stress is distributed across a large number of contacts. In contrast, where there are few fines, a single contact might be transmitting the same overall force. G1-40 has a similar PDF plot to that of the stable Dam Filter material, showing that a stable fabric emerges once the critical fines content is exceeded.
Differences in the force transmission networks of the coarse and fine soil fractions can be analysed quantitatively using PDF plots of normalised contact force, considering only the particles between specific diameter limits. The mean normal contact force for the whole sample, hNi, was used for normalisation. Each contact force is considered twice, once for each contacting particle. A contact force is considered to be part of a diameter range if at least one contacting particle is part of that range. The plot for sample Uniform is presented in Fig. 7a, split into three equal diameter limits by mass passing. It can be seen that the contact forces are well distributed through the diameters, but the particles with larger diameters generally transmitting larger forces. In contrast, in sample G1-20 (Fig. 7b), contacts between coarse particles generally transmit the large forces, whereas fine particles transmit primarily weak forces, f \ 1. The coarse fraction has a PDF curve whose shape is similar to the curve obtained for the stable samples Dam Filter or Uniform (see Fig. 6). This demonstrates that a stable, load-carrying matrix is formed by the coarse particles.
Internally unstable soils have a primary fabric that is formed of a small number of particles transmitting relatively large contact forces. As instability increases, the number of particles forming the primary fabric reduces, leading to a reduction in the probability that a given particle will be acted upon by a strong contact force (f [ 1). The probability of a given particle (including unconnected particles) having at least one strong contact force, P strong , is plotted against Kézdi (D 0 15 =d 0 85 ) max in Fig. 8. It can be seen that P strong generally decreases with increasing (D 0 15 =d 0 85 ) max . The transition between unstable, under-filled fabrics and stable, over-filled fabrics can also be identified. A comparison of G1-20 (under-filled) and G1-40 (overfilled) shows a jump from P strong = 0.02 to P strong = 0.19, with the partially over-filled G1-30 lying in between.
All unstable samples have P strong \ 0.05, whereas all the stable samples have P strong [ 0.05, as shown by the dashed line on Fig. 8. This shows that it is possible to draw a microscale distinction between stable and unstable states that agree well with empirically derived rules.

Discussion and conclusions
This paper presents microscale analyses of samples with PSDs of varying internal stability according to an existing empirical guideline, Kézdi's ratio (D 0 15 =d 0 85 ) max . Idealised samples with isotropic fabric and spherical particles are considered; these idealizations allowed the influence of the PSD on the internal topology of the material to be isolated. Periodic boundaries were used to minimise boundary effects and achieve computational efficiency. Prior experimental research (e.g. [14]) also used ideal spherical ballotini to assess internal stability. As a consequence of the very large differences in particle size, some inhomogeneity was induced during sample preparation in some of the samples; however, in a physical test, a sample with these gradings would also segregate under gravity alone.
Visual observation of the samples and analysis of the contact data confirmed that where samples are unstable according to Kézdi's (D 0 15 =d 0 85 ) max criterion, the internal fabric consists of a load-carrying matrix made up of the coarse particle and loose finer particles occupying the voids between the coarse particles. The micromechanical analysis of the simulation data considered the probability of a particle in each sample attaining a specific connectivity, the coordination number, volume-weighted coordination number and the probability of a normalised compressive contact being transmitted at each contact. Excluding consideration of the ''over-filled'' samples, as the filter stability increases, the overall trend is for the probability of a given particle to have no contacts with other particles to decrease and the probability of a particle having four or more contacts to increase. These data confirm that internally unstable samples, according to Kézdi's criterion, have a large number of loose fines. The average connectivity or coordination number, Nc, reduces with reducing internal stability. The large numbers of loose, fine particles in highly internally unstable soils dominate the coordination number values. The volume-weighted coordination number Nc vol , which accounts for the loose fines, also tends to increase in value with increasing empirical stability. While the probability that a given contact in the material transmits, a compressive force of a given magnitude does not seem to be a good indicator of relative stability, there seems to be a clear correlation between the probability that a particle will participate in the strong force chains (i.e. the stress transmitting contact network) and the stability of the material.
The materials considered here were mainly gap-graded, and the results indicate that a reduction in gap ratio leads to a more stable fabric. The concept of a critical fines content, that is, a fines content that is large enough to facilitate a more internally stable fabric where the fines participate in stress transmission had been proposed by Kenney and Lau [8]. Clear evidence of this phenomenon was observed in the simulations. This supports Kenney and Lau's suggestion that materials with fines contents above the critical level are inherently internally stable [8].
No clear transition was found between stable and unstable structures. Rather, the material fabric became progressively more unstable with increasing empirical instability according to the Kézdi criterion. A similar progressive increase in instability has been shown experimentally [3]. Micromechanical-based criteria that agree with the empirical rule can be identified within this progressive change, for example, all empirically unstable samples had probability of a given particle having a connectivity of C C 4 of less than 0.075, whereas all stable samples had a probability greater than this.
The ''coarse matrix, loose fines'' fabric shown in the unstable gradings is a prerequisite for internal instability. For a soil to be considered internally unstable, it must also have constrictions that are large enough to allow fine material to pass between voids. This may be another reason why clear transitions between stable and unstable gradings are not found. Future micromechanical-based research must address this issue. Note that in principle, DEM simulations could be used to examine clogging as well as erosion. However, such studies would require either a coupled DEM-fluid model (using, for example, the approach proposed by Tsuji et al. [28]) or a detailed analysis of the pore network (using, for example, the algorithm of Reboul et al. [21]).