Particle Migration Phenomena Related to Hydromechanical Effects at Contact between Different Materials in Embankment Dams

The compatibility between a fine grained base (B) material and a downstream coarser particulate transition (T), under the seepage forces related to hydraulic gradients, plays a key role in safety of earthfill dams. This aspect, that is the possible migration of fine grains through the voids larger than their size, is analyzed according to a numerical procedure simulating the 1D-coupled particle migration and seepage unsteady states. The procedure accounts for grain size curve, constriction sizes and porosity of the materials as well as the rate of the suspension and drag forces associated with the seepage flow; friction triggered by normal contact forces induced by confining pressure is considered too. The procedure has been systematically applied to: (i) simulate the newly formed filter (F) at the contact of different B-T systems, (ii) review the criteria proposed by Terzaghi, and (iii) analyze the particle migration phenomena that affected some embankment dams. and time) variability of granulometric properties, voids volume, porosity ( n ), permeability ( k ), flow velocity, local piezometric flow


Introduction
Internal erosion phenomena are an important safety issue for embankment dams, dikes, and levees as shown by historical failures and incidents [1]. These phenomena [1] occur when soil particles are dragged by seepage in embankments (e.g., at the contact between core and transition or downstream materials), foundations, from embankment to foundation, around and into conduits through embankments and adjacent walls supporting embankments.
Four mechanisms of internal erosion have been recently distinguished [1] (Figure 1): a. Concentrated leaks. The cracks are induced by differential settlement during construction or in operation of the embankment, by hydraulic fracture in case of low stresses, compared to the internal pressure, around erosion channels as well as by the desiccation at high levels in the fill. The sides of cracks may be eroded by the leaking water [1].
b. Backward erosion. The erosion process begins at a free surface on the downstream side of a dam and evolves beneath or within the embankment, through the development of erosion channels or piping.
c. Contact erosion. (i) Seepage flow parallel to the contact between coarse grained materials and fine soils may erode the smaller particles; (ii) particle migration of fine grained materials (base B) through the voids of a coarser materials (transition T), under seepage flow oriented along a direction almost orthogonal to the interface of the contacting materials [2].
d. Suffusion. The small particles of soil are dragged by the seepage flow through the pores of the coarser particles.
Concentrated leaks, suffusion, backward erosion (e.g., at the downstream face of the embankment or of the embankment core) may trigger piping erosion phenomena within the embankment or its foundation. The presence of adequate filters that autonomously generate within transition materials can prevent these triggering mechanisms from evolving at a larger scale and forming pipes [3]. Therefore, the problems of granulometric stability at the contact between materials characterized by different grain size curves, under seepage forces normally oriented with respect to the average contact surface, assume particular importance. The assessment of transition zones and downstream zones to act as filters is well understood and in many dams it can be demonstrated that these will provide adequate filter protection even if they do not fully meet modern "filter" design criteria [2]. The control of the granulometric stability of base (B) requires a correctly designed protective transition (T) [2,4]; its voids, related to the grain size distribution and porosity, must be small enough to stop the migrating B particles within short distances, thus avoiding limit states related to the backward erosion possibly appearing under the shape of flow pipes ( Figure 2); T must also allow a safe drainage of B, avoid clogging and blinding, inducing in turn uncontrolled increase of interstitial pressure.

Particle migration phenomena modeling
By referring to the representative elementary volume (R.E.V.) of a granular material composed by the volume of particles in suspension, the fluid phase volume and the solid phase (stationary solids) volume (Figure 3), the particle migration process can be described, along space x and time t, by the following system of governing partial differential equations (PDEs):  • fluid mass balance equation being v ⇀ f the velocity vector of the fluid phase; • solid mass balance equation being v ⇀ sp the velocity vector of the transported particles in suspension (v with χ e (0; 1]); c s , the particles in suspension concentration; n, the volumetric porosity; the contribution of diffusion term has been neglected; • "kinetic equation," describing the deposition and erosion processes of the deposited/accumulated particles within voids [5]; to this purpose, different formulations have been proposed in the past ( Table 1) [5][6][7][8][9][10].
The unknown variables of the problem are n, c s , v f , depending on t and x.
If the following hypotheses are assumed: Darcy's velocity), by combining Eqs. (1) and (2), it is obtained: For a complete simulation of the particles 1D migration and its evolution towards limit (clogging, blinding, complete erosion) granulometric or stability conditions, the (space and time) variability of granulometric properties, voids volume, porosity (n), permeability (k), flow velocity, local piezometric gradients, flow direction, as well as the particles erodibility should be taken into account. However, the variability of the previously defined physical-mechanical variables, especially of the voids volume distribution (VVD) [11,12], cannot be appropriately described through a simple equation.
Thus, the above formulations may be safely applied to real cases only if the limits deriving from extreme schematization of the analyzed process are removed.
To solve this extremely complex problem, a numerical procedure has been developed. The procedure takes into account grain size curve, constriction sizes distribution (CSD), porosity of the particulate materials, the rate of the suspension, piezometric gradients, drag forces associated with the seepage flow; the friction triggered by normal contact forces induced by confining pressure related to the effective stress state, is considered too.
Through the simulation of the 1D coupled particle migration and seepage unsteady states, this procedure aims: (i) to preliminarily evaluate if design protective transition zones, between core and downstream materials ( Figure 2), are adequate to the "natural" filters formation (time of formation, capability of new-formed filter to limit the transport/erosion of smaller particles); (ii) to understand if, in existing earth dams not provided with transitions, the formation of a "natural" transition zone, due to deposition of smaller particles within the voids of the downstream material (typical of widely graded morainic materials [2]), occurs and allows to interrupt/stop particle erosion or migration.

Proposed numerical procedure
To model the particle migration processes and evaluate the safety of earth structures against serviceability or ultimate limit states, an advanced characterization of the granular material as well as the simulation of the (space and time) variability of its properties must worked out.

Equation Parameters
Litwiniszyn (1966) ∂n ∂t ¼ c 1 ½n 0 À n Á c s À c 2 Á n • n 0 = initial porosity; • c 1 and c 2 = experimental coefficients (erosion of particles if c 1 = 0; c 2 > 0; deposition of particles if c 1 > 0; c 2 = 0).  The permeability coefficient (k) represents the fundamental parameter figuring in the equation that describes the seepage of a fluid through a porous medium; it mainly depends on grain size properties and porosity. k can be evaluated through the Kozeny-Carman relationship [13,14]: being d h the equivalent diameter of grains, previously defined ( Table 1); χ, a numerical coefficient; and μ w , the water viscosity.

Distribution of the volume of voids and corresponding distribution of constriction sizes
According to a "geometric-probabilistic" model [11,12,15], the most reliable distribution of the volume of voids within porous media corresponds to a situation of maximum "disorder" of the granular material (particles and voids). The probability that a generic void volume V assumes smaller or equal value than V * (cumulative function of probability) is from which it is possible to write the following "compatibility" equations: being V min and V max the minimum and maximum pore volumes, respectively, corresponding to the smallest (d min ) and largest (d max ) particle diameters (V min = 0.16855•d min , if the pore of maximum volume is composed by three particles; V max = 0.476•d max 3 , if the pore of maximum volume is composed by four particles); V, the expected value of the pores volume: where d j represents the diameter of a generic particle and ΔP j its corresponding percentage by weight. By maximizing the configurational entropy associated with the distribution of the pore volumes [12,15]: with κ 0 = κ •N vtot (κ is a constant, N vtot is the total number of pores equal to the total number of particles N ptot ), through the Lagrange's multipliers method, the following probability density function of pore volumes is obtained: The Lagrange's multiplier coefficient β is numerically deteminated through the following equation, obtained by introducing Eq. (10) in the "compatibility" Eq. (7): Thus, β depends on the porosity n through the expected value V of the pores volume.
Through Eq. (10), it is possible to obtain the cumulative probability function of the volume of voids: FðVÞ ¼ e ÀβV À e ÀβVmin e ÀβVmax À e ÀβVmin ð12Þ For an assigned volume V of a pore, the volume of the largest particle (V cs ), able to move through the porous material, satisfies the relation V cs < V. By assuming spherical particles (D, diameter), on the basis of geometric observations [11,15], it is possible to determine the diameters of the smallest (D cs,min ) and largest (D cs,max ) particles passing through the smallest (V cs,min ) and largest (V cs,max ) pores ( Figure 4); in other words, the minimum and maximum constriction sizes: (pores formed by three spherical particles) (pores formed by four spherical particles) By defining the coefficient η = V cs /V, it is obtained: Thus, the volume of constriction sizes and the corresponding diameter (CSD) can be generally evaluated as: If a linear change of η with volume V is simply assumed, it is obtained: Several methods to determine the CSD are available in technical literature [16][17][18][19][20][21][22]. Some of them take into account the loose (N p = 4, being N p the number of particles forming voids) and dense (N p = 3) soil states or relative densities [16][17][18][19][20][21]; other approaches enable the consideration of the porosity (n) resulting in a density dependent CSD [22,23].
The comparison between some of methods available in literature and the method proposed and adopted by the author, previously described, is carried out ( Figures 5 and 6).
It is observed that the author's method provides CSD curves ranged between the ones obtained through the models proposed by the authors of Refs. [16,17,19], corresponding to the dense soil state ( Figure 5), and by To et al. [22], for assigned values of porosity n ( Figure 6).
The porosity n, neglected by the authors of Refs. [16,17,19], may greatly influence the CSD; as a result, there might be theoretically the possibility of particle movements under seepage forces even for fairly poorly graded soils if the porosity would be high enough [22].
The comparison between the author's method and well-known experimental results [23] is reported in Figure 7. In the considered lab experiment [23], the constrictions sizes of a sample of gravel with grain sizes between 8 and 63 mm and porosity n = 0.34 were manually measured.    It is worth observing that the best interpretation of experimental CSD, through the model proposed by the author, is obtained for the porosity value n = 0.34 characterizing the tested material.

Problem's setting and governing equations
The heterogeneous porous material is decomposed into several elements (Figure 8a), each characterized by initial grain size curve (P i,j,0 ), porosity (n i,0 ) and permeability (k i,0 ); i, j, and t define the system element, materials granular fractions and the elapsed time, respectively [24]. Each element is schematically composed by original material (V or,i,t ), deposited/accumulated particles (V acc,i,t ; V acc,i,0 = 0), due to migration phenomena, and water saturating the i th element (V w,i,t ). According to the Kozeny-Carman Eq. (4), the permeability k i,t is expressed as The variables P i,j,t and n i,t (and then k i,t ) evolve because of erosion-deposition processes, associated with particle migration, causing an unsteady seepage flow. The unsteady state is simply analyzed by considering a sequence of steady states (time interval, Δt; "successive steady states" method); for each Δt, the continuity equation holds: Figure 7. Grain size distribution (GSD) [23], experimental CSD curve [23] and CSD curves obtained by the proposed method (d and P are the diameter and weight percentage related to grain size distribution, respectively; D cs and F cs are the diameter and weight percentage related to CSDs). Á Therefore, the suspension rate Q t through the elements of the section Ω, and the volume of the suspension V out,i , composed by the scoured particles dragged by the seeping fluid, V s,out,i,t , entered and washed out V w,out,i,t from each element, is the same during each Δt: Furthemore, Q f,in = Q t and Q s,in = 0 are assumed: in the first element of the system, during each Δt, the incoming fluid volume V w,in,t = V s,out,1,t . Eq. (3) is thus discretized as Figure 8. (a) Problem setting. One-dimensional unsteady seepage flow through a heterogeneous base (B)-transition (T) system; B and T are divided into elements; a constant head difference ΔH is imposed. Forces acting on a migrating particle: (b) plugged particle (d 0 ); and (c) unplugged particle (d < d 0 ); d, particle diameter; d 0 , average size of a pore channel (adapted from Ref.
Δt ð23Þ Then, it is possible to evaluate the temporal variation of the fluid volume within the elements which constitute the porous material: The particles can be scoured if subjected to a flow velocity greater than the local, critical flow rate (v cr ) [14,25]; the analysis of the actions on a movable particle and the dynamic equilibrium along the flow direction [2] allow the estimate of v cr.
In some cases, the drag force F D (Stokes law) may overcome the maximum local shear force related to the effective weight of the particle and the acting confining stresses (F S ) (Figure 8b and c); therefore, the particle can be eroded : α is a coefficient allowing to consider the density of the granular matrix (0 < α ≤ 4/π); α = 4/π for granular matrix composed by spherical particles arranged in hexagonal configuration, most dense state [14,26].
For a horizontal flow path, v cr is expressed as follows: (laminar flow) 24μ w ρ w dðv cr =nÞ þ 5:6 ffiffiffiffiffiffi μ w p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρ w dðv cr =nÞ p þ 0:25 The hydraulic conditions allowing the migration of movable particles are first considered; the analysis of the geometric conditions follows, as long as the previous ones is verified (v > v cr ). To determine the scoured particles (diameter d j ) composing the (ith) element, d j is compared with the constriction sizes of the (i + 1)th element. In particular, for each d j , the percentage F cs,i,j,t of smaller constriction sizes of the (i + 1)th element is defined (Figure 9); if the percentage of openings through which particles (d j ) may pass (100 -F cs,i,j,t ) is null, the considered (d j ) and higher remaining granulometric fractions are not able to cross through the voids of the (i+1)th element and then deposit; conversely, the smaller granulometric fractions can be scoured (if v > v cr ).
The balance of the eroded and deposited granular fractions allows to redefine the grain size distribution, the porosity and the permeability of each element and, then, to determine the piezometric gradients and flow velocity associated with each element of system. Once the fractions of accumulated (S acc,i,t ) and original eroded (S or,i,t ) material are determinated through hydraulic and geometric methods, it is possible to evaluate the specific weight of filtering suspension, by rearranging the equation proposed by Indraratna and Vafai [10]: The material scoured from the ith element, at time t (V s,out ), depends on the specific weight of the filtering suspension, γ m [10]: The scoured volume of material is composed both by V or,out and V acc,out [27]: V or,out and V acc,out are decomposed into their granular fractions [10]:  ΔV or, outi,t,j ¼ V or, outt,i Á P ori,t,j À P ori,t,jÀ1 S ori,t ; ΔV acc, outi,t,j ¼ V acc, outt,i Á P acci,t,j À P acci,t,jÀ1 S acci,t ð30Þ The element within which each scoured fraction is deposited may be determined by taking into account the corresponding length L mig,j of the migration path; L mig,j depends on the probability of a particle not encountering a smaller constriction size. The length covered by an assigned particle up to its arrest is based on concepts of stereology; it depends on the constriction sizes distribution (CSD) related to the PSD (pores size distribution) as well as on the thickness of the filter [12,15]. The length of the migration path is then compared to the length that the particles can cross during each step Δt: L mig, j ¼ minðs Á m j ; U t Á ΔtÞ, m j being the number of constrictions greater than the particle size encountered by the particle along its path; s is the unit step assigned to each comparison (Figure 10). At time t + Δt, the accumulated and original volume fractions within each element become [28]: ΔV acci, tþ1, j ¼ ΔV acci, t, j þ ΔV accini, t, j À ΔV accouti, t, j ; ΔV ori, tþ1, j ¼ ΔV ori,t,j À ΔV or, outi,t,j ð31Þ

Validation
Theoretical analyses have been carried out to validate the proposed numerical procedure by simulating laboratory tests carried out by different authors on selected materials [2,24].

Results by Atmazidis
Atmazidis [28] analyzed the particle migration phenomena at the contact between sands (base, B) and gravels (transition, T), under horizontal seepage flow and constant hydraulic gradient (i = 0.25) conditions. The amount of sand deposited within the gravel pores was measured after each test.
The main features of the experiments and the properties of the tested materials (e.g., grain size curves) are shown in Table 2 and Figure 11. Figure 12 shows the sand content retained in each gravel T material, migrated from the B material, at the end of each test, in function of the distance from sand-gravel interface.
Through the proposed numerical procedure, the experimental results [28] have been back analyzed. The following input parameters have been considered: T material (gravel) thickness = 2.5 m (experimental value); number of elements in which T material has been divided = 20 (each 12.5 cm length); number of elements in which B material has been divided = 1 (1 m length). The B material erosion resistance has been neglected [2].
The results, in terms of particles migration distance, obtained through the proposed procedure, appreciably approximate the values measured by Atmazidis [28]. Particularly, for "clean" gravel (Gravel 1, 0% initial sand content), particles migration distances are larger than those ones obtained for gravels initially containing small sand content (Gravel 2 and Gravel 3) ( Figure 12).
Furthermore, the numerical procedure allows to simulate the development of a filter at the sand (B)-gravel (T) interface, able to contrast the sand particle migration: high values of the sand content in the upstream layer of the gravel (15 ÷ 25% by weight, Figure 12) mean that a filter is formed [2].

Results by Skempton and Brogan
The effects of an upward seepage flow through a "gap graded" material, under increasing flow rate (Q) conditions, were measured [29]. The progressive increase of Q was obtained by increasing the hydraulic head difference between the lower and the upper surface of the material sample, until the piping process was triggered.    The main properties and the grain size curves of the considered tested material samples are reported in Table 3 and Figure 13, respectively [29].
In all experiments, it is observed an initial linear increase of the average seepage rate (q), following the rise of the average hydraulic gradient (i), until a critical value is achieved; afterwards, the ratio q/i progressively grows if i increases; the erosion of fine particles occurs; while, coarse (gravel) particles are not scoured [24].
To interpret the lab results, each sample (length L = 0.16 m) has been divided in eight elements (length 0.02 m); piezometric head ΔH = i*L (with i e (0, 0.28) and e (0, 0.37) for sample A and B, respectively) has been imposed. Under laminar flow condition, the numerical results seem to well interpret the experimental ones, especially for i < 0.10 (sample A) and i < 0.17 (sample B); for higher values of i, the numerical values follow a linear trend, underestimating the experimental ones. Under turbulent flow condition, for high values of i, the seepage velocity exponentially increase and the numerical results well fit the experimental ones (Figures 14 and 15). The difference between two flow regimes (laminar and turbulent) is due to higher drag forces in the turbulent regime inducing particles migration and, consequentially, the increase of seepage velocity.

Analysis of the empirical Terzaghi's criteria
Simulations of the erosion process of B material at the interface B-T have been carried out [30] for two B materials (extended (B1) and homogeneous (B2) grain size distribution GSD). For   [30,31], have been considered (Figure 16a and b). Results obtained for some B-T combinations are reported in Figures 17-21 in terms of evolution of permeability, porosity, average hydraulic gradient, flow velocity, and flow rate. To better understand the phenomenon, results for the combination B1-T3 are reported in Figures 22 and 23: the base has been divided into 10 elements; the transition has been divided in 30 elements. At the end of the analysis (t = 120 hours) particles of the eighth, ninth, and 10th element of the base are migrated within the first element of the transition forming the filter.  In Figure 22, the CSD (constriction size distribution) of the filter at the interface B1-T3 is shown. Filter's voids, related to the GSD of the 10th element of the base, are small enough to stop the migrating particles, according with Terzaghi's criteria. Analysis of the results, moreover, allows to identify areas for grain size distributions of the transitions that show a similar behaviour (Figure 23a and b). For B1, a first zone which incorporates the permeability ratio (finer transitions) is identified: the filter developes rapidly, in reason of the GSD similar to the B material, which minimizes the opening of the voids. The small percentage of fine particles of the T fosters backward erosion, which developes downstream (no confinement) back to the interface B-T. The stabilization of the phenomenon, therefore, is not guaranteed. The second zone has its upper limit in the piping ratio: the GSD of the transitions is different from the B one. The permeability of T is high while hydraulic gradient is small. The weight increase of the matrix particles close to the drain and, as a consequence, the increase of the resistance force opposite to the particles displacemnts strongly reduce the possibility of downstream backward erosion. The B-F-T system quickly reaches a configuration of granulometric equilibrium. A third zone, with coarser T, is delimited by a limit T, beyond which it is observed the filter (F) is not formated at the interface. Migrating particles, in fact, cross the transition to the drain without being intercepted. For B2, the first zone, with the finer transitions, is located before the permeability ratio: the filter quickly forms at the interface; however, backward erosion phenomena downstream are observed because of the small weight of the finer particles of the transition close to the drain (no confinement). The stabilization of the phenomenon is not ensured; for high values of ΔH, a rise of the flow rate is observed. The second zone is defined by the permeability ratio as lower limit and the piping ratio as upper limit. The influence of the hydraulic load is negligible and no backward erosion downstream is observed. The stabilization of migration is ensured. The third zone has features similar to those identified for the base B1 [32]. Sinkholes and leaks occurred through the moraine core of Suorva dam (Sweden) during periods of floods and high water level in 1983 [33]. In one case, in October 1983, up to 100 l/s of turbid leakage was seen but reduced after 10 days to about half this amount. Figure 24 shows the dam section and the possible leakage routes. The upper part of the dam was constituted by a coarse material, unable to contrast the internal erosion. The intermediate parts were constituted by materials with seal capacity only after an excessive erosion that may induce damages and sinkholes, as occurred at Suorva. Particularly, the erosion occurred in the core, through the entire width of the filter. This avoided the free downward drainage in the Figure 17. Simulation results for B1-T1 system. It is worth observing a downstream backward erosion; stabilization of the particle migration phenomenon is not attained. Simulations repeated at ΔH = 5 m and ΔH = 20 m showed a remarkable influence of the hydraulic head, adapted from Ref. [32]. filter, maintaining the water level in it to only 2 m below reservoir level for a relatively long time. The dam was not provided with any drain along the filter or into the downstream rockfill. In the lower part of the dam, where there was a fine and a coarse filter, the filter was finer, with a maximum d 15 of 1.0 mm, compared to a Sherard recommended size of 0.7 mm. Such "some-erosion" filters would soon seal if erosion initiated.
There were indications of erosion at the core-foundation interface, but no extensive damage.
The grouting records along the dam at the position where the leaks and sinkholes occurred show large grout takes at the possible leakage positions high in the core and lower in the core; and at the base of the core and into the upper parts of the foundation. To protect the dam against internal erosion and provide it to a drainage system to discharge the leakage flows caused by future incidents, a rockfill berm was built [33] on the downstream slope.

San Valentino dam
The 32 m high San Valentino "zoned" embankment dam was built in the 1940s in the basin of Adige river, North Italy, near the small town S.Valentino alla Muta (Bolzano); the dam is 447 m long and 7 m wide; the maximum reservoir elevation is 1498.10 m a.s.l. (Figure 25a). The dam project was modified several times to take into account the expected settlements and flow rates; in 1942, an advanced project was elaborated; the construction works were soon interrupted during the Second World War; construction ended in 1950 [34,35].
Fan-outwash materials (moraine debris with minor amounts of sand and silt) mainly affect the left side of the valley; the subsoil is mainly composed by fine gravel, sand and silt layers, with    peaty lenses, of alluvial and lacustrine deposition, characterized by significant thickness and low permeability [35].
On the right side, a fissured bedrock outcrops or lies at shallow depth. In the central part, the bedrock dips steeply and at the bottom of the valley it disappears under the sediments mantle  ( Figure 25b). The fan-outwash material was selected for the dam construction. The grain size of the quarried soils uniformly extends into the range of sand, gravel and silt with negligible amounts of clay. The core grains had d max = 50-60 mm; to reduce the core per-meability, a small amount of bentonite was added, without appreciable effects (Figure 25c). The embankment material was compacted at w opt = 6.6% (γ d = 20.8 kN/m 3 ).   Piezometers and assestimeters were installed during construction and operation. In the central part of the dam, small and slow increments of the piezometric heads (p.h.) were measured. In the period 1950-1985, the ratio between the measured maximum value of p.h. and the maximum reservoir level increased about 15% (sect. A-A, Figure 26). Afterwards, it became almost constant [35].
Complex FEAs have been deviced to evaluate the properties (geometry, permeability and stiffness heterogeneity) of foundation soils and dam body materials, concrete cutoff, permeability defects of the drain and of the "nominal core", unsaturated hydraulic behavior of materials [34].
Seepage analyses results show that the measured free surface, piezometric heads and leakages cannot been numerically simulated if a homogeneous dam is considered; it is necessary to impose a k-zoning of the dam body. Results also show that the current k distribution differs from the end of construction k distribution (Figure 27), due to suffusion phenomena that affected the core material, mainly the finer fractions. The interpretation of the piezometric heads measured near the cutoff allows to hypothesize a local permeability defect (local increase of k) immediately above the diaphragm wall, at the contact with the embankment materials.
According to this hypothesis, properly modeled, an appreciable agreement (both max values and time evolution) among measured and simulated piezometric heads ( Figure 28) and leakages values, is obtained [34].

Results
Migration phenomena regarding the moraine materials of the S. Valentino and Suorva dams [2,27] have been simulated through the proposed numerical procedure ( Table 4). The following input parameters are assigned: total length of the B-T systems = 3 m (for B, 1 m; for T, 2 m); number of elements dividing the B-T system = 60 (each 5 cm length); piezometric head difference H = 6 m [27]. To favor their migration, B material particles are considered unplugged (i.e., confining stresses are neglected). Numerical results point out that the finer fractions (d ≈ 0.002 mm) of the analyzed core materials are eroded [27]. The ratio Q t /Q 0 versus time t is shown in Figure 29 for the three analyzed case; Q 0 is the initial value of the suspension flow rate: Q t /Q 0 rapidly increases if the Suorva core material is protected by the coarser transition, due to the intense erosion of the finer fractions of B and the corresponding increase of permeability [27].
After 12 hours, Q t /Q 0 still increases: both granulometric and hydraulic stabilizations are not yet occurred. If the Suorva core material is protected by the finer transition, after 4 hours Q t /Q 0 slowly reduces: the erosion of B particles is not still completed; the clogging of the voids of T progressively occurs and the phenomenon seems to achieve a stable state [2]. In San Valentino dam, the washout of the finer particles of B is controlled by the finer protective transition; Q t  slightly increases compared with the Suorva considered cases. Stabilization occurs after 6 hours and Q t becomes only 1.2 times greater than Q 0 .
Very different values of the lengths of the path crossed by the scoured particles, through the analyzed granular transitions (T), are observed ( Figure 30). Within the T material of the S. Valentino dam, the particles mainly deposit and accumulate just a few cm, after the B-T interface. In the finer T material of Suorva dam, the particles mainly accumulate at about 70 cm after the B-T interface; while, in the coarser transition material of Suorva dam, at about 150 cm after the B-T interface [2].

Concluding remarks
An original numerical procedure to simulate the particle migration phenomena in granular media due to seepage flows has been developed. The procedure takes into account voids, constriction sizes, and porosities of the particulate materials (geometric-probabilistic models) as well as the rate of the seeping suspension and piezometric gradients (hydraulic models). The continuous variations (in the space, 1D, and along time, t) of the physical (porosity n) and hydraulic (permeability k) properties of the granular medium, during the coupled seepage flow, deposition and scouring particles processes, are simulated. First, validation of the proposed method has been carried out by simulating selected experimental tests, referred to materials whose diameters range from clay to sand (base materials, B) and sand to gravel (granular transition, T). Then, the proposed numerical procedure has been applied to interpret the empirical Terzaghi's criteria. Finally, the analysis of the granulometric stability of cohesionless moraine materials, constituting the core of earth dams, has been carried out. Through the proposed numerical procedure, the particle migration phenomena occurring at the interface of the contact core-protective transition has been simulated. Results show that the proposed procedure is able to simulate the analyzed deposition-erosion processes; particularly, "not negligible" erosion phenomena, involving the finer fractions of these materials, may cause anomalies and malfunctions in dams whose core is constituted by broadly graded cohesionless materials protected by coarse granular transitions.