Open access peer-reviewed chapter

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

Written By

Francesco Federico

Submitted: 31 October 2016 Reviewed: 08 February 2017 Published: 06 September 2017

DOI: 10.5772/67785

Chapter metrics overview

1,412 Chapter Downloads

View Full Metrics


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.


  • internal erosion
  • contact between base and transition materials
  • suffusion
  • numerical procedure

1. 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):

Figure 1.

Internal erosion mechanisms.

  1. 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].

  2. 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.

  3. 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].

  4. 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.


2. 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):

Figure 2.

Filter formation mechanism.

Figure 3.

Physical scheme of R.E.V.

  • fluid mass balance equation

    [ ( 1 c s ) n v f ] = [ ( 1 c s ) n ] t E1

    being v f the velocity vector of the fluid phase;

  • solid mass balance equation

    ( c s n v s p ) = [ ( 1 c s ) n ] t E2

    being v s p the velocity vector of the transported particles in suspension ( v s p = χ v f , with χ ϵ (0; 1]); cs, 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) [510].

Equation Parameters
Litwiniszyn (1966) n t = c 1 [ n 0 n ] c s c 2 n
  • n0 = initial porosity;

  • c1 and c2 = experimental coefficients (erosion of particles if c1 = 0; c2 > 0; deposition of particles if c1 > 0; c2 = 0).

Sakthivadivel (1966) n t = P 0 [ n 0 n ] n 0 Q c s  
  • P0 = initial probability of deposition of small particles within the voids of the granular material;

  • Q = fluid flow rate.

Vardoulakis (2004) n t = Λ ( 1 n ) c s n v s p - = erosion law constitutive (experimental) parameter with dimension of inverse length.
Saada (2005) n t = a c s - a = positive scalar
Zhang et al. (2013) ( 1 c s ) n t = n c s t k P f c s + k ρ f g c s
  • k = permeability coefficient, depending on n and cs;

  • Pf = fluid pressure;

  • ρf = fluid density.

Indraratna and Vafai (1997) γ m = γ w V w + γ s V s V w + V s = ( 1 c s ) γ w + c s γ s 1  
  • γw = the unit weight of fluid

  • γs = the unit weight of solid phases; Vw = the fluid volume

  • Vs = the volume of particles affected by smaller diameter than the average one (d0) of the voids, defined as

    d 0 = 2.67 n 1 n d h

    with d h = 1 / i Δ P i d ¯ i ; ΔPi = Pi+1Pi, granulometric passing percentage associated with d ¯ i = (di+1 + di)/2.

Table 1.

Some “kinetic” equations available in the literature.

The unknown variables of the problem are n, cs, vf, depending on t and x.

In the proposed numerical procedure, the relationship proposed by Indraratna and Vafai (1997) [10] is applied (Table 1).

If the following hypotheses are assumed: (i) vsp = vf (χ =1) and (ii) unidirectional flow (1-D case), since v D = n v f ( v D , Darcy’s velocity), by combining Eqs. (1) and (2), it is obtained:

v D , x c s x = [ ( 1 c s ) n ] t E3

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.


3. 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.

3.1. Characterization of the granular material

3.1.1. Permeability coefficient

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]:

k = χ γ w μ w n 3 ( 1 n ) 2 d h 2 E4

being dh the equivalent diameter of grains, previously defined (Table 1); χ, a numerical coefficient; and μw, the water viscosity.

3.1.2. 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

F ( V * ) = V min V * f ( V ) d V E5

from which it is possible to write the following “compatibility” equations:

V min V max f ( V ) d V = 1 E6
  V min V max V f ( V ) d V = V ¯ E7

being Vmin and Vmax the minimum and maximum pore volumes, respectively, corresponding to the smallest (dmin) and largest (dmax) particle diameters (Vmin = 0.16855∙dmin3; Vmax = 0.16855∙dmax3, if the pore of maximum volume is composed by three particles; Vmax = 0.476∙dmax3, if the pore of maximum volume is composed by four particles); V ¯ , the expected value of the pores volume:

V ¯ = n 1 n i Δ P i i 6 Δ P i π d ¯ i 3 E8

where dj represents the diameter of a generic particle and ΔPj its corresponding percentage by weight. By maximizing the configurational entropy associated with the distribution of the pore volumes [12, 15]:

E = κ V min V max f ( V ) ln [ f ( V ) ] d V E9

with κ = κ ∙Nvtot (κ is a constant, Nvtot is the total number of pores equal to the total number of particles Nptot), through the Lagrange’s multipliers method, the following probability density function of pore volumes is obtained:

f ( V ) = e β V V min V max e β V E10

The Lagrange’s multiplier coefficient β is numerically deteminated through the following equation, obtained by introducing Eq. (10) in the “compatibilityEq. (7):

[ β ( V min V ¯ ) + 1 ] e β V min [ β ( V max V ¯ ) + 1 ] e β V max = 0 E11

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 β V min e β V max e β V min E12

For an assigned volume V of a pore, the volume of the largest particle (Vcs), able to move through the porous material, satisfies the relation Vcs < 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 (Dcs,min) and largest (Dcs,max) particles passing through the smallest (Vcs,min) and largest (Vcs,max) pores (Figure 4); in other words, the minimum and maximum constriction sizes:

Figure 4.

Constriction size for pores formed by (a) three and (b) four particles.

D c s , min = 2 ( 3 0.5 / 3 1 / 2 ) D V c s , min = 1.94 10 3 D 3 E13

(pores formed by three spherical particles)

D c s , max = ( 2 0.5 1 ) D V c s , max = 3.72 10 2 D 3 E14

(pores formed by four spherical particles)

By defining the coefficient η = Vcs/V, it is obtained:

η min = V c s , min V min = 1.15 10 2 E15
  η max = V c s , max V max = 7.81 10 2 E16

Thus, the volume of constriction sizes and the corresponding diameter (CSD) can be generally evaluated as:

V c s = η ( V ) V E17
D c s = η ( V ) V 6 π 3 E18

If a linear change of η with volume V is simply assumed, it is obtained:

η ( V ) = 1 V max V min [ ( η max η min ) V + ( η min V max η max V min ) ]   E19

Several methods to determine the CSD are available in technical literature [1622]. Some of them take into account the loose (Np = 4, being Np the number of particles forming voids) and dense (Np = 3) soil states or relative densities [1621]; 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).

Figure 5.

Grain size curve [19] and CSD curves obtained by different methods (d and P are the diameter and weight percentage related to grain size distribution, respectively; Dcs and Fcs are the diameter and weight percentage related to CSDs).

Figure 6.

Grain size curve [22] and CSD curves obtained by different methods (d and P are the diameter and weight percentage related to grain size distribution, respectively; Dcs and Fcs are the diameter and weight percentage related to CSDs).

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.

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; Dcs and Fcs are the diameter and weight percentage related to CSDs).

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.

3.2. Problem’s setting and governing equations

The heterogeneous porous material is decomposed into several elements (Figure 8a), each characterized by initial grain size curve (Pi,j,0), porosity (ni,0) and permeability (ki,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 (Vor,i,t), deposited/accumulated particles (Vacc,i,t; Vacc,i,0 = 0), due to migration phenomena, and water saturating the ith element (Vw,i,t). According to the Kozeny-Carman Eq. (4), the permeability ki,t is expressed 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 (d0); and (c) unplugged particle (d < d0); d, particle diameter; d0, average size of a pore channel (adapted from Ref. [2]).

k i , t = χ γ w μ w n i , t 3 ( 1 n i , t ) 2 d h i , t 2 E20

The variables Pi,j,t and ni,t (and then ki,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:

[ k 1 , t l 1 k 2 , t l 2 0 0 0 k 2 , t l 2 k 3 , t l 3 0 0 0 k i , t l i k i + 1 , t l i + 1 0 0 0 k N 1 , t l N 1 k N , t l N 1 1 1 ] ( Δ h 1 , t Δ h 1 , t Δ h 2 , t Δ h N , t ) = ( 0 0 0 Δ H ) E21

Therefore, the suspension rate Qt through the elements of the section Ω, and the volume of the suspension Vout,i, composed by the scoured particles dragged by the seeping fluid, Vs,out,i,t, entered and washed out Vw,out,i,t from each element, is the same during each Δt:

Q t = Ω k i , t Δ h i , t l i ;   V out , t = Q t Δ t E22

Furthemore, Qf,in = Qt and Qs,in = 0 are assumed: in the first element of the system, during each Δt, the incoming fluid volume Vw,in,t = Vs,out,1,t. Eq. (3) is thus discretized as

v D , t c s , i + 1 , t c s , i , t l i = φ f , i , t + 1 φ f , i , t Δ t E23

with ϕ f = (1 - cs)∙n. Eq. (23) can be rewritten as

v D , t Ω Δ t ( c s , i + 1 , t c s , i , t ) = Ω l i ( φ f , i , t + 1 φ f , i , t ) E24

Then, it is possible to evaluate the temporal variation of the fluid volume within the elements which constitute the porous material:

Q t Δ t ( c s , i + 1 , t c s , i , t ) = ( V f , i , t + 1 V f , i , t ) E25

The particles can be scoured if subjected to a flow velocity greater than the local, critical flow rate (vcr) [14, 25]; the analysis of the actions on a movable particle and the dynamic equilibrium along the flow direction [2] allow the estimate of vcr.

In some cases, the drag force FD (Stokes law) may overcome the maximum local shear force related to the effective weight of the particle and the acting confining stresses (FS) (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, vcr is expressed as follows:

v c r = n 3 μ w [ ( γ s γ w ) d 2 6 + α d 2 ( σ z ' + σ y ' ) ] tan ϕ E26

(laminar flow)

[ 24 μ w ρ w d ( v c r / n ) +   5.6 μ w ρ w d ( v c r / n ) + 0.25 ] 1 4 ρ w ( v c r n ) 2 ( γ s γ w ) d 3 ( sin   α + cos   α tan ϕ ) λ   ( σ 1 ' + σ 2 ' )   tan ϕ E26.bis

(turbulent flow)

Eq. (33.bis) (turbulent flow regime) must be numerically solved. The original material is subjected to strong confining actions (frictional forces, geometric hindrances); high flow velocities are needed to mobilize the plugged particles. Conversely, the accumulated, unplugged particles, may be easily scoured during simulation; the corresponding critical flow velocity assumes small values (Eq. (26)).

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 > vcr). To determine the scoured particles (diameter dj) composing the (ith) element, dj is compared with the constriction sizes of the (i + 1)th element. In particular, for each dj, the percentage Fcs,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 (dj) may pass (100 – Fcs,i,j,t) is null, the considered (dj) 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 > vcr).

Figure 9.

Comparison between particles (diameter dj) and CSD.

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 (Sacc,i,t) and original eroded (Sor,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]:

γ m i , t = ( S a c c i , t 100 V acc i , t + S or i , t 100 V or i , t ) γ s i + V w i , t γ w S acc i , t 100 V acc i , t + S or i , t 100 V or i , t + V w i , t E27

The material scoured from the ith element, at time t (Vs,out), depends on the specific weight of the filtering suspension, γm [10]:

V s , out t , i = γ m t , i γ w γ s γ w V out , t E28

The scoured volume of material is composed both by Vor,out and Vacc,out [27]:

V s , out t , i = V or , out t , i + V acc , out t , i ;   V acc , out t , i = V acc t , i S acc i , t V acc t , i S acc i , t + V or t , i S or t , i E29

Vor,out and Vacc,out are decomposed into their granular fractions [10]:

Δ V or , out i , t , j = V or , out t , i P or i , t , j P or i , t , j 1 S or i , t ;   Δ V acc , ou t i , t , j = V acc , out t , i P acc i , t , j P acc i , t , j 1 S acc i , t E30

The element within which each scoured fraction is deposited may be determined by taking into account the corresponding length Lmig,j of the migration path; Lmig,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 ) , mj 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]:

Figure 10.

Proposed numerical procedure flowchart.

Δ V acc i , t + 1 , j = Δ V acc i , t , j + Δ V accin i , t , j Δ V accout i , t , j ;   Δ V or i , t + 1 , j = Δ V or i , t , j Δ V or , out i , t , j E31

3.3. 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].

3.3.1. 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.

Test number 5 6 7 8
Grain size curves of coupled sand and gravel A,3 A,2 B,1 A,1
Sand content of gravel before test (% by weight) 8% 4% 0% 0%
Maximum sand migration length (cm) 75 110 125 215
Amount of sand deposited within gravels (kg) at the end of test 3.17 7.88 9.25 14.5
Negligible migration after (hours) 3.25 3.50 3.00 6.50

Table 2.

Synopsis of Atmazidis experimental tests (adapted from Ref. [28]).

Figure 11.

Grains size curves of involved materials (transition, T; base, B), adapted from Ref. [2].

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.

Figure 12.

Sand content (s.c.) retained in the gravel transition at the end of the filtering process; theoretical vs experimental results: (a) s.c. after 6.5 hours (Sand A-Gravel 1); (b) s.c. after 3 hours (Sand B-Gravel 1); (c) s.c. after 3.5 hours (Sand A-Gravel 2); and (d) s.c. after 3.25 h (Sand A-Gravel 3), adapted from Ref. [2].

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].

3.3.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].

Sample k (m/s) N Gravel (%) Sand (%)
A 0.0045 0.34 85 (G) 15 (S1)
B 0.0084 0.37 85 (G) 15 (S1)

Table 3.

Main properties of tested materials.

Figure 13.

Grain size curves: (a) materials used to prepare the samples and (b) tested samples.

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 ϵ (0, 0.28) and ϵ (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.

Figure 14.

Numerical results vs experimental values: seepage flow (v) vs i, sample A.

Figure 15.

Numerical results vs experimental values: seepage flow (v) vs i, sample B.


4. 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 each B, different transitions T, matching or not the Terzaghi’s criteria [30, 31], have been considered (Figure 16a and b). Results obtained for some B-T combinations are reported in Figures 1721 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.

Figure 16.

Grain size distributions of bases and transitions analyzed: base B1 (a), base B2 (b). Terzaghi’s criteria have been highlighted (piping ratio D15f/D85b < 4, permeability ratio D15f/D15b > 4), adapted from Ref. [32].

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].

Figure 18.

Simulation results for B1-T5 system. Downstream backward erosion is negligible; stabilization of the particle migration phenomenon is reached. Simulations repeated at ∆H = 5 m and ∆H = 20 m showed no remarkable influence of the hydraulic head, adapted from Ref. [32].

Figure 19.

Simulation results for B2-T1 system. The downstream backward erosion is appreciable; stabilization of the particle migration phenomenon is not guaranteed. Simulations repeated at ∆H = 5 m and ∆H = 20 m showed a remarkable influence of the hydraulic head; if ∆H = 20 m, an uncontrolled increase of flow rate can be observed, adapted from Ref. [32].

Figure 20.

Simulation results for B2-T3 system. Downstream backward erosion is negligible; stabilization of the particle migration phenomenon is not attained. Simulations repeated at ∆H = 5 m and ∆H = 20 m showed no remarkable influence of the hydraulic head. Length of the filter is about 70 cm, adapted from Ref. [32].

Figure 21.

Non-dimensional discharge flow Q/Q0 along time t for systems B1-T1, B1-T5, B2-T1 e B2-T3. B1-T1 and B2-T1 systems show the continuous increase of the discharge flow along time; for systems B1-T5 and B2-T3, the flow stabilizes after a few hours at ≈20–30% of their initial value, adapted from Ref. [32].

Figure 22.

Evolution of grain size distribuition in the elements of B1-T3. The CSD of the filter is reported, adapted from Ref. [32].

Figure 23.

Transitions (T)-base (B) systems, with graded (a) or almost uniform (b) grain size curves, adapted from Ref. [32].

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].


5. Analysis of cases

5.1. Suorva dam

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 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 d15 of 1.0 mm, compared to a Sherard recommended size of 0.7 mm. Such “some-erosion” filters would soon seal if erosion initiated.

Figure 24.

Suorva dam: possible leakage paths during 1983 incident (adapted from Ref. [33]).

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.

5.2. 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].

Figure 25.

San Valentino dam: (a) main cross section (A-A), settlements and piezometric measuring points; (b) longitudinal profile: (cc) concrete cutoff (grey colored), (fb) fissured bedrock, (al) alluvium, (fo) fan-outwash; and (c) grain size of materials (adapted from Ref. [34]).

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 dmax = 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 wopt = 6.6% (γd = 20.8 kN/m3).

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].

Figure 26.

San Valentino dam: Ratio and difference between the yearly maximum pie-zometric heads Hpz,max (piezometers paAM2 and paAC, Sec. A-A) and the yearly maximum reservoir elevation Hmax (the piezometric heads are referred to the elevation of the lower point of the valley 1470 m a.s.l.), adapted from Ref. [34].

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.

Figure 27.

San Valentino dam: permeabilities k (m/s) of dam materials and foundation soils, (AC) after construction and (A) current (adapted from Ref. [34]).

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].

Figure 28.

San Valentino dam: measured and theoretical p.h., in 2011 (adapted from Ref. [34]).

5.3. 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 Qt/Q0 versus time t is shown in Figure 29 for the three analyzed case; Q0 is the initial value of the suspension flow rate: Qt/Q0 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].

Materials k0 (m/s) n0 ϕ (°) c’ (kPa)
San Valentino (Base) 1 × 10−6 0.3 25 0
San Valentino (Transition) 1 × 10−6 0.3 25 0
Sourva (Base) 1 × 10−6 0.3 25 0
Sourva (Finer Transition) 2.5×10−6 0.3 25 0
Sourva (Coarser Transition) 5 × 10−6 0.3 25 0

Table 4.

Physical and mechanical parameters characterizing the analyzed materials.

Figure 29.

Ratio Qt/Q0 vs time t (hours) ; Q0, initial value of the suspension flow rate Qt (adapted from Ref. [2]).

After 12 hours, Qt/Q0 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 Qt/Q0 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; Qt slightly increases compared with the Suorva considered cases. Stabilization occurs after 6 hours and Qt becomes only 1.2 times greater than Q0.

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].

Figure 30.

Particles accumulation in the protective transitions (d = 0.002 mm) (adapted from Ref. [2]).


6. 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.


  1. 1. ICOLD. Internal Erosion of Existing Dams, Levees and Dykes, and Their Foundations. In: Bridle, R. and Fell, R., Eds., Bulletin 164, 2013, Volume 1: Internal Erosion Processes and Engineering Assessment, International Commission on Large Dams, Paris.
  2. 2. Federico F, Montanaro A. Granulometric stability of moraine embankment dam materials: theoretical procedure and back-analysis of cases. In: Proceedings of 6th International Conference on Dam Engineering, C. Pina, E. Portela, J. Gomes (ed.), 15-17 February 2011, Lisbon, Portugal, pp. 423–439
  3. 3. Tran DK, Prime N, Froiio F, Callari C, Vincens E. Numerical modeling of backward front propagation in piping erosion by DEM-LBM coupling. European Journal of Environmental and Civil Engineering 2016, Taylor & Francis ed., DOI: 10.1080/1964189.2016.1248794
  4. 4. Moraci N, Mandaglio MC, Cazzuffi D. I geotessili con funzione di filtro a contatto con terreni granulari: criteri e parametri di progetto. Italian Geotechnical Journal. 2010;2:47–71
  5. 5. Litwniszyn. Colmatage-scouring kinetics in the light of stochastic birth-death process. Bull. Acad. Pol. Sci., Ser. Sci. Techn. 1966;14 (9):561
  6. 6. Sakthivadivel R. Theory and mechanism of filtration of non-colloidal fines through a porous medium. Tech. Rep. HEL 15-5, 1966, University of California, Berkeley
  7. 7. Vardoulakis I. Fluidization in artesian flow conditions: Hydromechanically unstable granula media. Geotechnique. 2004, 54(3):165–177
  8. 8. Saada Z, Canou J, Dormieux L, Dupla JC, Maghous S. Modelling of cement suspension flow in granular porous media. International Journal for Numerical and Analytical Methods in Geomechanics. 2005;29:691–711
  9. 9. Zhang XS, Wong H, Leo CJ, Bui TA, Wang JX, Sun WH, Huang ZQ. A thermodynamics-based model on the internal erosion of earth structures. Geotechnical and Geological Engineering. 2013;31:479-492. DOI: 10.1007/s10706-012-9600-8
  10. 10. Indraratna B, Vafai F. Analytical model for particle migration within base soil - filter system. Journal of Geotechnical Engineering. 1997;123(2):100–109
  11. 11. Musso A, Federico F. Geometrical probabilistic approach to the design of filters. Italian Geotechnical Journal. 1983;17(4):173–193
  12. 12. Musso A, Federico F. Pore size distribution in filtration analyses. In: Proceedings of XI ICSMFE, S. Francisco. 1985;1:1207–1212
  13. 13. Harr ME. Groundwater and Seepage. New York: Dover Publications Inc., 1962
  14. 14. Kovacs G. Seepage Hydraulics. Elsevier Scientific Publishing Company, Amsterdam; Oxford, New York, 1981
  15. 15. Federico F, Musso A. Some advances in the geometric-probabilistic method for filter design. In: Proceedings of International Conference on “Filters and Filtration Phenomena in Geotechnical Engineering”, Karlsruhe, October 1992, pp. 75–82
  16. 16. Silveira A. An analysis of the problem of washing through in protective filters. In: Proceedings of the 6th International Conference on “Soil Mechanics and Foundation Engineering”, ICSMFE, Montréal. 1965;2:551–555
  17. 17. Silveira A, de Lorena Peixoto T, Nogueira J. On void size distribution of granular materials. In: Proceedings of the 5th Pan-Am Conference on “Soil Mechanics and Foundation Engineering”, Buenos Aires, Argentina, 1975, pp. 161–176
  18. 18. Indraratna B, Raut AK, Hhabbaz H. Constriction-based retention criterion for granular filter design. Journal of Geotechnical and Geoenviromental Engineering (ASCE). 2007;133(3):266–276
  19. 19. Dallo YAH, Wang Y. Discussion of ‘A new theoretical method to evaluate the internal stability of granular soils’. Canadian Geotechnical Journal. 2012;49: 866–868. DOI:10.1139/T2012-036
  20. 20. Moraci N, Mandaglio MC, Ielo D. Reply to the discussion by Dallo and Wang on ‘A new theoretical method to evaluate the internal stability of granular soils’. Canadian Geotechnical Journal. 2012;49:869–874. DOI: 10.11.39/T2012-047
  21. 21. Wang Y, Dallo YAH. On estimation of the constriction size distribution curve for cohesionless soils. European Journal of Environmental and Civil Engineering. 2014;18(6):683–698, DOI: 10.1080/19648189.2014.909335
  22. 22. To HD, Scheuermann A, Williams DJ. A new simple model for the determination of the pore constriction size distribution. In: Proceedings of the 6th International Conference on Scour and Erosion, 31 August 2012, Paris, France, pp. 295–303
  23. 23. Witt K. Filtrationsverhalten und Bemessung von Erdstoff-Filtern. Technical Report 104, D. Inst. F. Bodenmechanik und Felsmechanik, Karlsruhe, 1986
  24. 24. Federico F, Iannucci M. Effects of granulometric stability of cohesionless materials on safety of geotechnical structures. In: Proceedings of International Symposium on Geomechanics from “Micro to Macro” 2014, Cambridge, UK. 2014; 2, pp. 913–918
  25. 25. Perzlmaier S. Hydraulic criteria for internal erosion in cohesionless soil. Internal Erosion of Dams and their Foundations - Fell & Fry (eds.). London Taylor & Francis Group, ISBN 978-0-415, 2007, pp. 179–190
  26. 26. Biswas S. Study of cohesive soil granular filter interaction incorporating critical hydraulic gradient and clogging. Engineering-Research Master, University Of Wollongong, NSW, Australia, 2005
  27. 27. Federico F, Montanaro A. Numerical analyses of granulometric stability of moraine dam cores. In: Proceedings of 7th European Conference on Numerical Methods in Geotechnical Engineering, Trondheim, Norway, 2010
  28. 28. Atmatzidis DK. A study of sand migration in gravel. In: Proceedings of 12th International Conference on Soil Mechanism and Foundation Engineering., Session 8/3, Rio De Janeiro, Brazil, 1989, pp. 683–686
  29. 29. Skempton AW, Brogan JM. Experiments on piping in sandy gravels. Geotechnique. 1994;44(3):449–460
  30. 30. Terzaghi K. Der Grundbruch and Stauwerken and Scine Verhutung. Reprinted in “From Theory to Practice in Soil Mechanics”. New York: John Wiley and Sons, Inc.; 1922. pp. 450–476
  31. 31. Terzaghi K. Soil mechanics: a new chapter in engineering science. Journal of the Institution of Civil Engineers. 1939;12:106–141
  32. 32. Federico F, Jappelli R, Iannucci M. The Progress of Internal Erosion from Typical Hydromechanical Response of Contacts in Embankment Dams. In: Proceedings of XVI International Technical Conference on Dam Monitoring 2015 – TKZ 2015 - Session 1 “Internal Erosion in Earth Structures – Poland, 29 September – 2 October
  33. 33. Nilsson A. The susceptibility of internal erosion in the Suorva dam. Internal Erosion of Dams and their Foundations, Editors R Fell and J-J Fry, Taylor and Francis; 2007. pp. 167–172.
  34. 34. Pinamonti P, Scienza M, Catalano A, Del Gizzi F, Jappelli R, Federico F, Montanaro A. Safety of the San Valentino Earth Dam after 60 Years of Operation. In: Proceedings of 8th ICOLD European Club Symposium 2010 - DAM SAFETY, Sustainability in a Changing Environment, 22–23 September 2010, Innsbruck, Austria
  35. 35. Federico F, Jappelli R, Montanaro A, Scienza M. Self-protection of a quasi-homogeneous embankment dam revealed by advanced analyses and monitoring. In: Proceedings of 8th European Conference on Numerical Methods in Geotechnical Engineering, June 18-20, 2014, Delft, The Netherlands

Written By

Francesco Federico

Submitted: 31 October 2016 Reviewed: 08 February 2017 Published: 06 September 2017