Next Article in Journal
Economic Impact of Overtopping and Adaptation Measures in Catalan Ports Due to Sea Level Rise
Next Article in Special Issue
Study on Multi-Scale Coupled Ecological Dispatching Model Based on the Decomposition-Coordination Principle
Previous Article in Journal
Batch and Column Scale Removal of Cadmium from Water Using Raw and Acid Activated Wheat Straw Biochar
Previous Article in Special Issue
Local Scour of Armor Layer Processes around the Circular Pier in Non-Uniform Gravel Bed
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Internal Stability Evaluation of Soils

1
State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China
2
Ocean College, Zhejiang University, Hangzhou 310027, China
3
Department of Civil, Environmental, and Geo-Engineering, University of Minnesota, Minneapolis, MN 55455-0116, USA
*
Author to whom correspondence should be addressed.
Water 2019, 11(7), 1439; https://doi.org/10.3390/w11071439
Submission received: 18 May 2019 / Revised: 9 July 2019 / Accepted: 9 July 2019 / Published: 12 July 2019

Abstract

:
Suffusion constitutes a major threat to the foundation of a dam, and the likelihood of suffusion is always determined by the internal stability of soils. It has been verified that internal stability is closely related to the grain size distribution (GSD) of soils. In this study, a numerical model is developed to simulate the suffusion process. The model takes the combined effects of GSD and porosity (n) into account, as well as Wilcock and Crowe’s theory, which is also adopted to quantify the inception and transport of soils. This proposed model is validated with the experimental data and shows satisfactory performance in simulating the process of suffusion. By analyzing the simulation results of the model, the mechanism is disclosed on how soils with specific GSD behaving internally unstable. Moreover, the internal stability of soils can be evaluated through the model. Results show that it is able to distinguish the internal stability of 30 runs out of 36, indicating a 83.33% of accuracy, which is higher than the traditional GSD-based approaches.

1. Introduction

Suffusion is the mass movement of fine particles through the pore space of a coarser matrix driven by seepage forces [1]. It is one of the main mechanisms initiating internal erosion within levees, earth dams and foundations [2,3,4] as well as watershed hillslopes [5,6].
A lot of experimental [7,8,9,10,11,12,13,14] and numerical (CFD) studies [15,16] were launched and expected to reveal the underlying mechanisms corresponding to suffusion. Most of these studies were macroscopic and under the important assumption that the movement of soil is a continuous process. This assumption has never been verified through or connected directly to a particle-scale study. However, some of following studies in granular system start to show their potential to support this assumption. For instance, as DEM models [17,18,19] and CFD+DEM models [20] are developed, the suspension, collision or spin of particles have been described in detail, but still seem to be irrelevant to the internal structure fracture (based on regime transition research [21]) or massive movement (based on rheology research [22,23,24]). Since then, macroscopic studies on suffusion are still mainstream, especially when focusing on engineering practices.
Suffusion can take on many forms, but it is particularly insidious when it occurs under a relatively low hydraulic gradient. Hydraulic engineers determine the internal stability by examining the critical gradient of a soil which directly affects the likelihood of diffusion. Generally, the soil is internally stable if ic < 0.7 [7].
The internal stability of soils is closely related to its grain size distribution (GSD). As pointed out by Wan and Fell [1], internally unstable soils are usually broadly graded soils with particles from silt (or clay) to gravel size, whose GSD curves are concave upward (or known as gap-graded soils). GSD-based assessment approaches, based on experimental evidence, have been used extensively in practice [8,9,10,11,12,13,14]. Most of them distinguish coarse and fine particles within the soil and incorporate characteristic sizes into an assessment index or indices.
Early studies accounted for two particles sizes with an additional division size. For example, Kézdi [25] divided soils into coarse and fine components. The maximum value of D15/d85, which depends on the division size, is used as the index of internal stability. D15 (d85) is the size for which 85% (15%) by weight is finer in the coarse (fine) components, respectively. The soil is evaluated as unstable if D15/d85 > 4. Similarly, Kenney and Lau [26] introduced a “shape curve” to identify an internally unstable soil. They defined F as the percentage of particles that are finer than an adjustable size d and H + F as the percentage of particles finer than the size of 4d. The material is assessed as internally unstable when the minimum value of H/F is less than 1.3 originally, [26] and a modified lower value of 1.0 later [27]. Different from the usage of an adjustable particle size, Burenkova [28] considered three characteristic sizes and evaluated the soils unstable when d90/d60 is greater than 1.86log(d90/d15) + 1 or less than 0.76log(d90/d15), where dx (x = 15, 60 or 90) hereafter is the diameter for which x% by weight of the soil particles are finer.
In recent decades, the idea of a more precise representation of the GSD curve has been further developed. Wan and Fell [1] extracted four sizes from the GSD curve. Two are for coarse particles and the other two for the fine counterpart of the soils. They defined two variables, s1 = 15/log(d20/d05) and s2 = 30/log(d90/d60), for each soil sample. The well-known assessment criterion is that a soil sample is internally unstable when s1 < 22 and s2 > 80. That is, an internally unstable soil is featured as d20/d05 > 4.8 and d90/d60 < 2.4. Also, its accuracy is higher than any other GSD-based approaches mentioned above.
The accuracy of the GSD-based approaches has been improving since more and more GSD information has been added in. In terms of seepage flow dynamics, however, these approaches are still empirical and with ambiguous meaning. Early research suggests that the gap-graded soils are apt to being internally unstable since fine particles easily pass through, but do not fill in the voids between the coarse particles if subjected to seepage [9]. The filtration of detached fine particles has therefore been well accepted as the mechanism responsible for the internal stability. Overall, suffusion appears as the combination of detachment, transport, and filtration of fine particles. They may depend on the excessive shear stress related to particle inception [15,29], sediment transport capacity if definable, and the size of the voids, respectively. In addition to the filtration effect, the other mechanisms may dominate in some cases, but have received less attention so far. Such incomplete understanding hampers the explanation of other important factors such as porosity [30,31] and loading history [32].
This paper presents a numerical simulation of fine particle dynamics subjected to seepage flow, within an extended numerical framework of Fujisawa’s study [15]. We take into account the combined effects of GSD and porosity. Wilcock and Crowe’s theory [33] is adopted to quantify the inception and transport of each size group of soil particles. We reproduce the phenomena observed in seepage tests and explain why the soils with specific GSD behave more internally unstable, and we predict the internal stability of soil through the model.

2. Numerical Model

Numerical models have been successfully applied to analyze seepage flow through embankments, dikes, and soil foundations. With a reasonable model for particle movement within soil matrix, numerical models may be able to evaluate internal stability of soils with a certain GSD. Our new model is an alternative approach for exploring internal stability.

2.1. Governing Equations

In terms of representing the dynamics of fine particles subjected to seepage forces, numerical simulation at the microscopic scale provides a potential approach. In 2010, Fujisawa et al. described the process of suffusion (or the process of porosity n change with time t) within a unit volume of seepage field by using ∂n/∂t = EA [15]; The involved parameters E and A are the erosion rate defined as the volume of eroded particles per erodible area per time, and the erodible area per unit volume.
Taking the GSD into consideration, as shown in Figure 1, we assume that the particles from the soil skeleton would be eroded into the pore flow and increase its concentration, considering the two-dimensional seepage flow in the directions 1 and 2. During the time period of dt, the volume of eroded soil particles is ∑EkAkdVdt, where Ek denotes the erosion rate of kth size group section defined as the volume of eroded particles per erodible area per time, Ak denotes the erodible area per unit volume of kth size group, dV = dx1dx2dL, and dx1, dx2, and dL are the dimensions of the unit volume, respectively.
The governing equations of mass and momentum conservation are as follows:
n t = E k A k
n C t + v i C x i = E k A k
n S r t + v i x i = 0
v i = k s h x i
where n denotes the porosity of soil mass, C denotes the solid concentration of seepage flow, vi denotes the seepage velocity in the two directions (i = 1, 2), Sr denotes the saturation of soils, h denotes the local water head, ks denotes the hydraulic conductivity, xi denotes the distance along the ith direction (i = 1, 2), and t denotes the time passing by. Details about the closure equations will be introduced in Section 2.2.
Item ∑EkAk differs from EA in the study of Fujisawa et al. [15] since the dependence of Ek on the size of kth group is taken into account. If the dependence is neglected, i.e., Ek = E holds for each size group, the formulation by Fujisawa et al. [15] will be recovered.

2.2. Closures

Closures of items erosion rate Ek, erodible area Ak. and hydraulic conductivity ks are listed below:

2.2.1. The Erosion Rate Ek

For suffusion, a linear relationship has been used between the erosion rate E and the excessive shear stress τ acting on eroded particles. However, this assumption does not distinguish the difference in the inception of particles with different sizes. Here, we adopt Wilcock and Crowe ‘s theory [33] for mixed sediments which found wide applications in geomorphic science and river engineering. That is, the erosion rate Ek is evaluated taking into account GSD, i.e.,
E k = { α ( τ τ c k )   for   τ τ c k 0   for   τ < τ c k
where τ is shear stress exerted on the soil particles, Ek and τck are the erosion rate and critical shear stress for the kth size group, respectively, provided that the soil particles can be divided into several size groups (group number ≥k) and the particles within the same size group have the same critical shear stress.
The selection of α is based on calibration. More details will be introduced in Section 3.3.
The shear stress τ that is exerted on soil particles, is generated by the seepage flow with the soil fabric. In 1980, Hillel [34] assumed that the seepage field within a unit volume (see Figure 2a) can be treated as equivalent tubes in a length of L and radius of R = D/2 (see Figure 2b).
As shown in Figure 2, the equivalence guarantees the same discharge due to seepage flow under the same hydraulic gradient as well as the same porosity in average. The shear stress τ is then defined at the wall of tube. Using the Poiseuille law, it can be obtained by the following:
τ = I 2 ρ g k s μ n
where I = h/L denotes the local hydraulic gradient; g denotes the acceleration of gravity; ρ = s + (1 − C)ρw is the density of seepage flow; ρs denotes the density of soil particles; ρw denotes the density of pure water; μ = μw(1 + 2.5C) is the dynamic viscosity of seepage flow with suspended sediments, which is corrected by the suspension concentration; and μw denotes the viscosity of pure water. It is noted that Equation (6) quantifies an equivalent shear stress within the unit volume, which balances the pressure drop driving the seepage flow.
The equivalent shear stress τ in the model is a macroscopic measure of the microscopic seepage shear stress exerted on the soil particles in the unit volume, and the potential difference between these two items is compensated through the calibration of erodible coefficient α in Equation (5). Meanwhile, the erosion of the fine particles is assumed to be equivalent to the process of sediment transportation, such an idea is also introduced when calculating the critical shear stress τc.
In previous studies [15,35], the critical shear stress τc was assumed to be a constant. However, since the GSD cannot be ignored, coarse particles will protect the fine ones from being eroded (known as the “Hiding Function”). Thus, the critical shear stress varies when particle size differs (see Figure 3). By integrating Wilcock and Crowe’s theory [33] for the inception and transport of non-uniform sediments, for the kth group containing the soil particles of size between dk and dk+1, the critical shear stress τck is evaluated based on Equations (7) and (8):
τ c k = τ c s m ( d k d s m ) b
τ c s m ( ρ s ρ w ) g d s m = A + B exp ( 20 F s )
where dsm denotes the medium size of soil particles; τcsm is the corresponding critical shear stress; b denotes an exponent that may differ from the erosion in open channel flows and is calibrated in the simulation; Fs denotes the content ratio of fine particles (diameter d < 2 mm); and A and B are the two empirical parameters, and can also be calibrated during the simulation.
Parameter A in Equation (8) is usually related to the cohesion of the soil. When there are cohesive particles in the soil, the value of A will be relatively higher [33]; parameter B in Equation (8) is usually related to fluid properties [33]. Due to the nature differences between the Poiseuille flow in our model and the open channel flow in the sediment transport process in Wilcock and Crowe’s study [33], the value of B may also be different. In our simulation, these two parameters are obtained through calibration. More details will be introduced in Section 3.3.

2.2.2. The Erodible Area Ak

Within a volume of dV, the Ak for the kth size group can be quantified with the total projective area of the given size dk, that is, the product of the projective area of a particle, πdk2, and the number of such particles per unit volume, 6(1 − n)(Pk+1Pk)/πdk3. Therefore, the Ak reads,
A k = 6 ( 1 n ) ( P k + 1 P k ) / d k
where Pk+1 and Pk denote the percentage passing by weight of soil particles for the size finer than dk+1 and dk respectively.

2.2.3. The Hydraulic Conductivity ks

The hydraulic conductivity ks is related to the saturation, porosity, and GSD of the soil. A relationship for the hydraulic conductivity of the saturated soil is proposed, and it reads
k s = Γ d 10 1.565 n 2.3475 ( 1 n ) 1.565
where d10 denotes the diameter of 10% for which soil particles by weight are finer, and Γ is a constant (calibrated in Section 3.3). Such a relationship is based on Chapuis’ study [10], as well as the selection of exponents of each item.

2.2.4. Update of d10 in Equation (10)

In the case of suffusion driven by a seepage flow, the moving out of eroded particles will change the local porosity and GSD. As a result, the hydraulic conductivity ks as well as the critical shear stress τck has a feature of dynamic adjustment during the erosion process. In order to capture the dynamic features, it is necessary to formulate the temporal change of Pk corresponding to the particle size dk, and then to specify d10 in Equation (11). For this purpose, one can quantify the Pk(t + dt) after a time period dt as
P k ( t + d t ) = P k ( t ) ( 1 n ) i = 1 k E i A i d t ( 1 n ) E A d t
A rearrangement of Equation (11) without considering the second-order term results in
d P k ( t ) d t = P k ( t ) E A i = 1 k E i A i 1 n
Accordingly, the size of d10 corresponding to 10% for which soil particles by weight are fine after the time period dt can be quantified approximately as
d 10 = ( d x + d x + 1 ) / 2
where dx and dx+1 can be obtained through exhaustive method among all soil size groups, under the condition that item (Px − 10%) + (Px+1 − 10%) reaches the smallest absolute value.

2.3. Numerical Method and Solution Procedure

The numerical model is solved by using a central scheme of the finite volume method (FVM) built on the OpenFOAM platform (version 2.2.2, https://openfoam.org/release/2-2-2/). More details when assigning basic variables (seepage velocity vi and porosity n) and governing equations are shown in Appendix A.
The flow chart indicating the solution procedure is illustrated in Figure 4.
Variables in closures are updated in every cell of simulation domain, and calculations above are repeated at following time steps.

3. Model Verification

In this section, the proposed model is examined by previous experiments (studies by Skempton and Brogan in 1994 [9] and Wan and Fell in 2008 [1]). The experimental works have detailed records of hydraulic gradient (or eroded soil mass) at different stage of seepage flow and are therefore used for model calibration and validation.
Furthermore, in the next section, more experiments (studies by Lau in 1984 [36]; Lafleur et al. in 1989 [8]; and Chapuis in 1996 [10]) are also included when examining the model’s capability in predicting internal stability.

3.1. Overview of Experiments

3.1.1. Devices

The designs in experimental devices between these studies were similar: their experiments were conducted using acrylic cylinders filled up with fully mixed soils (see Figure 5). The upper end of the cylinder was the outlet of the seepage field, which allowed free outflow (also known as being downstream of the seepage field), and the discharge of the seepage flow was measured when running the experiment. The lower end was the inlet of the seepage field (also known as being upstream of seepage field), which connected the filter layer and the water pipe. The filter layer was designed to make the inflow distribution more stable, and the water pipe was used to connect the soil device with the water tank.

3.1.2. Observed Stages of Suffusion Process

With the increase of hydraulic gradient, Skempton and Brogan [9] observed the seepage flow and identified three stages of suffusion: (1) the “dancing-like movement of particles” stage, during which initial signal demonstrating that the suffusion had begun; (2) the “slight but general movement of particles” stage, during which suffusion kept developing and formed the shape of the “pipe” started stretching towards the upstream, and (3) the “strong general piping of particles” stage, where the soil particles were strongly moving out of the soil fabric.
Wan and Fell [1] gave similar characterizations. They identified the flow at the outlet as “clear”, “cloudy”, and “very cloudy” by manually observing and recording, because the degree of turbidity indicated the degree of suffusion, just like the judgment shown in Skempton and Brogan [9]. The “clear” flow indicated no suffusion and other conditions meant erosion at different stages. The soil samples were evaluated as internally unstable when the stage of “slight but general movement of particles” or “cloudy” flow was reached below the hydraulic gradient of a manually fixed value which was set as the critical one (e.g., 0.7 mentioned in Section 1).

3.1.3. Soil Samples

Skempton and Brogan [9] used three different samples of non-cohesive soils and reported three runs of data, namely Runs A, B, C, respectively. In comparison, Wan and Fell [1] presented 18 runs of experimental data with 14 soil samples (upward flow tests only), their soils contained a portion of cohesive particles (diameter d < 0.02 mm). The maximum soil diameter was d = 30 mm. Both gap-graded soils and continuous-graded soils were used in the experiments. (see Figure 6).
The definition of gap-graded soils varies in different studies. However, among them, most gap-graded soils have irregular shapes. For instance, a horizontal part of the curve could be graded from 0.02 mm to 2 mm. Otherwise the soils are defined as continuous-graded.
For the runs that lacked complete curves (e.g., Runs 2 and 7), we assume that the diameter of particles in the missing part is the same size as the finest soil particles shown in curves, so that we can perform simulations for these runs.

3.2. Simulation Domain and Boundary Settings

The simulation domain is proposed considering the circular symmetry of the seepage field, and assuming a two-dimensional flow at the vertical cross-section plane. As shown in Figure 7a, the domain covers a rectangle cross-sectional area. The size of the area is corresponding to the actual size of experimental device. Therefore, for the runs in Skempton and Brogan [9] the height of the domain is set as 155 mm and the width is set as 139 mm; and for the runs in Wan and Fell [1] the height is set as 300 mm and the width is set as 300 mm. The inlet boundary is imposed with a stable, consistent water head at each stage of hydraulic gradient. Taking Run B in Skempton and Brogan [9] as an example, the loading history is shown in Figure 7b.
In 2017, Rochim et al. [32] noted that the loading history might affect the experimental results of critical hydraulic gradient as well as the eroded mass. They found that if they set a longer duration time under a specific hydraulic gradient, they would record a lower critical hydraulic gradient and more eroded mass.
Such a phenomenon can be explained: when the duration time is not long enough under a specific hydraulic gradient, the erosion process may not reach a steady/stable status, which will lead to the unreliable judgment of critical hydraulic gradient.
To eliminate the effect of different loading histories, in our simulation, a long enough duration time and corresponding stable results are used to validate our model.

3.3. Parameter Calibration

During the simulation, the pure water density ρw is set as 1000 kg/m3, and the viscosity coefficient μw is set as 10−6 Pa·s. Besides GSD, other initial conditions of the soils (including initial porosity n0, initial permeability coefficient ks0, and density of soil particles ρs) are also reproduced according to the experimental records.
There are two ways to obtain the initial porosity n0 of the soil: (1) directly from the experimental records, and (2) by calculating n0 = 1 − Vs0/V, where Vs0 denotes the initial soil volume and V denotes the initial volume of seepage field.
Parameters need to be calibrated include α in Equation (5), b in Equation (7), A and B in Equation (8), and Γ in Equation (10). The calibrated parameters are presented in Table 1. These parameters are calibrated using a part of experimental data, i.e., Runs A and B from Skempton and Brogan [9] (for non-cohesive soils) and Runs 2 and 3 from Wan and Fell [1] (for cohesive cases).
Based on these runs, A, B, b and α are calibrated by trial and error. A0 = 0.021, B0 = 0.015, and b0 = 0.067, as suggested by Wilcock and Crowe [33] are used as initial values for A, B, and b, and α0 = 0.0000021, as suggested by Fujisawa et al. [15], is used as the initial value of α.
It should be noted again that Wan and Fell [1] did not record the amount of soil erosion, and α ranges from 0.0001 to 0.01 would not affect the calculation results in Section 3.4.1. However, Skempton and Brogan [9] did record the amount of soil erosion during the experiment, the model can reasonably reproduce the q-i relationship as well as the eroded mass when α = 0.005 is subtracted. Therefore, α is fixed as 0.005 in our model. Value of Γ is calculated based on Equation (14):
Γ = k s 0 ( d 10 1.565 n 0 2.3475 ( 1 n 0 ) 1.565 )

3.4. Calibration Results

3.4.1. q-i Relationship

Figure 8 presents the calibration results for runs in Section 3.3:
It shows that the computed q-i relationships agrees well with measured ones. It should be noted that ±10% change of these parameters will not affect the calibration results.

3.4.2. Eroded Mass

The computed eroded mass in Run A (duration time, 80 min; last hydraulic gradient, i = 0.28) and Run B (duration time, 90 min; last hydraulic gradient, i = 0.34) in Skempton and Brogan [9] are 220.40 g and 77.16 g, which are close to the measured mass 225 g and 75 g, with relative errors less than 3%.

3.5. Model Verification

3.5.1. q-i Relationship in Other Runs

Using the calibrated parameters, we conducted the simulation with the data from other runs, the results of representative runs are shown in Figure 9.
The internally unstable runs (e.g., Runs A and B) and internally stable runs (e.g., Run C) in Skempton and Brogan [9] show obvious difference in shape of q-i curves. For example, when i = 0.28 and 0.34, curves of Runs A and B show a sudden rise, respectively, whereas the curve of Run C shows a stable trend all over the experiment. The explanation is as follows: for Runs A and B, as particles are massively eroded by the seepage flow under a specific hydraulic gradient i, the porosity n and the hydraulic conductivity ks increases, as well as the discharge q; while for Run C, less soil particles are eroded under such a hydraulic gradient which triggers massive suffusion in Runs A and B, or an even higher hydraulic gradient. The simulation result of the data by Wan and Fell [1] show a similar phenomenon.

3.5.2. Eroded Mass for Run C in Skempton and Brogan

The computed eroded mass for Run C in Skempton and Brogan [9] is 57.4 g, close to the measurement of 55 g. The relative error is less than 5%.

3.6. Determination of Critical Hydraulic Gradient

When i reaches critical hydraulic gradient ic, it can be inferred that the soil has undergone a sudden, irreversible, and largely destructive adjustment due to suffusion. Run 10 in Wan and Fell [1] is used as an example to reveal the reason. In Figure 10, the suffusion process can be divided into three different stages according to the different slopes of the q-i curve:
In stage I (i < 0.55), discharge of seepage flow q grows linearly with the hydraulic gradient i, during this stage the out flow is “clear” according to the corresponding experimental record; and in stage II (0.55 < i < 0.76), q grows in a nonlinear, but controlled, trend with i, during this stage the outflow is “cloudy” according to the corresponding experimental record; and in stage III (i > 0.76 = ic), one may notice an abrupt change of q-i curve, and this abrupt change may indicate the failure of soil fabrics.
When i < ic, the suffusion behavior has already taken place in the formation of losses of fine particles in limited scale, and such losses will not affect the stability of the whole soil fabrics. When i = ic or i > ic, more particles are lost until the failure of whole structure occurs.
Following the justification, ic is evaluated by observing the boundary between stage II and stage III for all the runs in the three studies. Figure 11 compares the computed critical gradients with the measured gradients. 18 of the total 21 runs are reproduced satisfactorily (relative error <10%). The three exceptions are Runs 5, 11, 13 in Wan and Fell [1].

3.7. Porosity Distribution in Different Suffusion Stages

In this study, taking Run A in the study by Skempton and Brogan [9] as an example in order to describe the suffusion process more directly and concretely, we present the spatial distribution of porosity corresponding to stage I, II, and III in 3.6 (see Figure 12 and Figure 13), respectively:
Figure 12a denotes the initial porosity distribution (corresponding to stage I), Figure 12b,c denotes the porosity distribution when it reaches a stable one under hydraulic gradient i = 0.1 and i = 0.2 respectively (corresponding to stage II since i < ic = 0.28). When i = 0.1, the high-porosity zone only appears limitedly in the local region near the outlet, which means that there is no significant erosion in the seepage field. When i = 0.2, the high-porosity zone expands for a while, however, it is still limited to the local region near the outlet.
Figure 13a–c denotes porosity distribution when i = ic = 0.28 (corresponding to stage III), during the loading time t = 3000 s, 4000 s and 5000 s. The high-porosity zone grows so rapidly and intensely that it connects the upstream to the downstream. Such a phenomenon indicates an unsteady status under the massive suffusion behavior.

4. Predicting Internal Stability

In Section 3, we describe the suffusion process numerically and evaluate the critical hydraulic gradient ic, and in this section we attempt to give a further investigation in predicting internal stability.

4.1. Shear Stress Comparison Approach

Averagely, suffusion will be triggered when τ > τcr, where τcr denotes the critical shear stress of the finest soil particles (particle size dr) remaining in soil fabrics. According to Equation (7), τcr = min(τck) as it has been already known that dr = min(dk), and it will increase when suffusion proceeds. Meanwhile, as the particles are continuously eroded out of the seepage field, the porosity ρ and hydraulic conductivity ks of the soils become higher, and based on Equation (10), τ will also increase. When i = ic, τ will keep exceeding τcr as the suffusion continues. The suffusion will develop inevitably until the whole fabrics are destroyed.
For each of the soils, there must exist a hydraulic gradient ik < ic, when i = ik, if the loading time is long enough, and all the particles with sizes d < dk will be eroded out of the soil fabrics. Under this assumed situation, dr = dk, and τcr = τck. At the moment, the volume of the soils Vs will be decreased by k% (ρs is a constant), and the porosity of the soil will be n = n0 + (1 − n0)k%.
Therefore, when n = n0 + (1 − n0)k% is given, τcr = τck; where τck is determined by Equations (7) and (8). Since then, a quantitative τcr-n relationship can be computed. A τ-n relationship can also be obtained through Equation (6).
Here, an approach to predict S (stable) and U (unstable) based on the comparison of the relationship between τcr-n and τ-n is illustrated: when i is given, by comparing the relative locations of the τcr-n curve and the τ-n curve in the same coordinate system, we can judge whether suffusion starts, stops or proceeds irreversibly.

4.2. Evaluation Results

Taking Run A in Skempton and Brogan [9] as an example (see Figure 14), under the hydraulic gradient i = 0.28, when n < 0.405, τ > τc, which means suffusion is triggered and developed; when n > 0.405, the τ-n curve almost coincides with the τc-n curve, which means that the suffusion process is irreversible. When i = 0.30, the τ-n curve is always above the τc-n curve and there is no intersection between the two curves. This indicates that under this hydraulic gradient, the suffusion process of the soil is still irreversible.
According to the results shown in Figure 14, the critical hydraulic gradient is ic = 0.28, which is in accordance to the experimental data of Skempton and Brogan [9], and since ic < 0.7, according to the study by Hillel [34], the internal stability of soils in this run is U. It should be noted that α does not involve with the analysis process, which indicates that it does not affect the results of internal stability prediction.
To verify the approach, experimental data from Lau [36], Lafleur et al. [8], and Chapuis et al. [10] are also studied. The GSD curves are presented in Appendix B. This approach is able to distinguish the internal stability of 30 runs out of 36, indicating 83.33% of accuracy, while the most accepted traditional GSD-based approach based on the study by Wan and Fell [1] is 75%.

5. Conclusions

In this study, a new numerical model taking combined effects of GSD and porosity into account is developed to describe the suffusion process. The theory by Wilcock and Crowe [33] in field of sediment transportation is adopted to quantitively describe the inception and transport of finer soil particles. During the simulation, all of the related parameters are calibrated in a solid way. The numerical results showed satisfactory performance when validating with the experimental data of Skempton and Brogan [9] and Wan and Fell [1].
The model can reproduce the q-i relationships of experiments as well as critical hydraulic gradients ic of runs, and describe the suffusion process by showing porosity spatial distribution at different stages.
Inspired by the modeling work, an approach to predicting internal stability is also proposed. Disclosing the relationships among the factors, we found that the soils with specific GSD and porosity are prone to induce the suffusion. The model demonstrates better performance in the judgments of the internal stability (83.33% of correctness) of 36 experimental runs. The accuracy is higher than most of traditional GSD-based approaches. The results provide useful guidelines in macroscopic engineering practice.

Author Contributions

Conceptualization, Q.F. and X.F.; methodology, Q.F. and H.-C.H.; software, Q.F. and T.M.; validation, Q.F., Y.J. and J.W.; numerical runs, Q.F. and X.F.; investigation, Q.F.; data curation, H.-C.H. and J.W.; writing—original draft preparation, Q.F.; writing—review and editing, Q.F., Y.J. and X.F.; visualization, Q.F. and T.M.; supervision, X.F.; project administration, X.F.; funding acquisition, X.F. and Y.J.

Funding

This research is funded by The National Key Research and Development Program of China (Grant No. 2017YFC0804602), and The National Natural Science Foundation of China (Grant No. 51525901). The second author acknowledges the support of open fund of State Key Laboratory of Hydroscience and Engineering, Tsinghua University.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The numerical method for assigning basic variables (seepage velocity vi and porosity n) within each cell are listed as:
v i l , t = Ω l v i t d V l V l
n l , t = Ω l n t d V l V l
where Ωl and Vl denote the region and volume of the lth cell, respectively. vil,t and nl,t denotes the seepage velocity and porosity within lth during tth time step, respectively.
Based on Equations (A1) and (A2), the unknown variable porosity nl,t+1 in governing Equation (1) can be solved numerically as:
n l , t + 1 = n l , t + d t E k t A k t
where Ekt denotes the erosion rate of kth size group section at tth time step, and Akt denotes the erodible area per unit volume of kth size group at tth time step.
Adopting the derivation by Fujisawa et al. [21], the unknown variable concentration Cl,t+1 in governing Equation (2) can be solved numerically as:
V l n l , t + 1 C l , t + 1 n l , t C l , t d t + q i t t l i j Δ s j l = V l n l , t + 1 n l , t d t
As shown in Figure A1, tijl and Δsjl denotes the unit normal vector and length of the jth (j = 1, 2, 3 and 4) boundary segment of the lth cell. Also, item qit tijl can be calculated through Equation (A5) using the FVS scheme:
q i t t l i j = C l , t v l , t t i j l + C r , t v r , t t i j l
where superscript r = r1, r2, r3, and r4 denotes the neighboring cells (shearing the boundary) of lth cell.
Figure A1. Geometry of the cells.
Figure A1. Geometry of the cells.
Water 11 01439 g0a1
In the simulation, item ∂Sr/∂t is always 0 since the soil are fully saturated all the time. Similar to Equations (A4) and (A5), the unknown variable water head hl,t in Equations (3) and (4) can be solved numerically based on the following Equation:
k s l , t grad h 1 r 1 , t 2 grad h 1 l , t + grad h 1 r 2 , t Δ s 1 l + k s l , t grad h 2 r 3 , t 2 grad h 2 l , t + grad h 2 r 4 , t Δ s 1 l = 0
where,
grad h 1 l = h r 1 , t 2 h l , t + h r 2 , t Δ s 1 l
grad h 2 l = h r 3 , t 2 h l , t + h r 4 , t Δ s 2 l

Appendix B

The GSD curve of soils in Lau [36], Lafleur et al. [8], and Chapuis et al. [10] are shown in Figure A2:
Figure A2. GSD curves of soils in (a) internal stable runs; (b) internal unstable runs in studies by Lau [36], Lafleur et al. [8], and Chapuis et al. [10].
Figure A2. GSD curves of soils in (a) internal stable runs; (b) internal unstable runs in studies by Lau [36], Lafleur et al. [8], and Chapuis et al. [10].
Water 11 01439 g0a2

References

  1. Wan, C.F.; Fell, R. Assessing the potential of internal instability and suffusion in embankment dams and their foundations. J. Geotech. Geoenviron. 2008, 134, 401–407. [Google Scholar] [CrossRef]
  2. Wan, C.F.; Fell, R. Investigation of rate of erosion of soils in embankment dams. J. Geotech. Geoenviron. 2004, 130, 373–380. [Google Scholar] [CrossRef]
  3. Fell, R.; Fry, J.J. State of the art on the likelihood of internal erosion of dams and levees by means of testing. In Erosion in Geomechanics Applied to Dams and Levees, 1st ed.; Bonelli, S., Nicot, F., Eds.; John Wiley & Sons Inc: London, UK, 2013; Volume 1. [Google Scholar]
  4. Fleshman, M.S.; Rice, J.D. Laboratory Modeling of the mechanisms of piping erosion initiation. J. Geotech. Geoenviron. 2014, 140, 04014017. [Google Scholar] [CrossRef]
  5. Abrams, D.M.; Lobkovsky, A.E. Growth laws for channel networks incised by groundwater flow. Nat. Geosci. 2009, 2, 193–196. [Google Scholar] [CrossRef]
  6. Howard, A.D.; McLane, C.F. Erosion of cohesionless sediment by groundwater seepage. Water Resour. Res. 1998, 24, 1659–1674. [Google Scholar] [CrossRef]
  7. Adel, H.D.; Bakker, K.J. Internal stability of minestone. In Proceedings of the International Symposium on Modelling Soil–Water-Structure Interactions, Rotterdam, The Netherlands, 29 August–2 September 1988; pp. 225–231. [Google Scholar]
  8. Lafleur, J.; Mlynarek, J. Filtration of broadly graded cohesionless soils. J. Geotech. Eng. 1989, 115, 1747–1768. [Google Scholar] [CrossRef]
  9. Skempton, A.W.; Brogan, J.M. Experiments on piping in sandy gravels. Geotechnique 1994, 44, 449–460. [Google Scholar] [CrossRef]
  10. Chapuis, R.P.; Contant, A. Migration of fines in 0–20 mm crushed base during placement, compaction, and seepage under laboratory conditions. Can. Geotech. J. 1996, 33, 168–176. [Google Scholar] [CrossRef]
  11. Chang, D.S.; Zhang, L.M. A Stress-controlled erosion apparatus for studying internal erosion in soils. Geotech. Test. J. 2011, 34, 579–589. [Google Scholar]
  12. Chang, D.S.; Zhang, L.M. Critical hydraulic gradients of internal erosion under complex stress states. J. Geotech. Geoenviron. 2013, 139, 1454–1467. [Google Scholar] [CrossRef]
  13. Fleshman, M.S.; Rice, J.D. Constant gradient piping test apparatus for evaluation of critical hydraulic conditions for the initiation of piping. Geotech. Test. J. 2013, 36, 834–846. [Google Scholar] [CrossRef]
  14. Rönnqvist, H.; Viklander, P. On the Kenney-Lau approach to internal stability evaluation of soils. Geomaterials 2014, 4, 129. [Google Scholar] [CrossRef]
  15. Fujisawa, K.; Murakami, A. Numerical analysis of the erosion and the transport of fine particles within soils leading to the piping phenomenon. Soils Found. 2010, 50, 471–482. [Google Scholar] [CrossRef]
  16. Li, M.X.; Fannin, R.J. Capillary tube model for internal stability of cohesionless soil. J. Geotech. Geoenviron. 2013, 139, 831–834. [Google Scholar] [CrossRef]
  17. Guo, N.; Zhao, J. A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media. Int. J. Numer. Methods Eng. 2014, 99, 789–818. [Google Scholar] [CrossRef] [Green Version]
  18. Wagner, G.J.; Gregory, J. Coupling of atomistic and continuum simulations using a bridging scale decomposition. J. Comput. Phys. 2003, 190, 249–274. [Google Scholar] [CrossRef]
  19. Li, X.; Ke, W. A bridging scale method for granular materials with discrete particle assembly—Cosserat continuum modeling. Comput. Geotech. 2011, 38, 1052–1068. [Google Scholar] [CrossRef]
  20. Lominé, F.; Scholtes, L. Modeling of fluid–solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: Application to piping erosion. Int. J. Numer. Anal. Methods Geomech. 2013, 37, 577–596. [Google Scholar] [CrossRef]
  21. Bi, D.; Zhang, J. Jamming by shear. Nature 2011, 480, 355–358. [Google Scholar] [CrossRef]
  22. Boyer, F.; Élisabeth, G. Unifying suspension and granular rheology. Phys. Rev. Lett. 2011, 107, 188301. [Google Scholar] [CrossRef]
  23. Cassar, C.; Nicolas, M. Submarine granular flows down inclined planes. Phys. Fluids 2005, 17, 103301. [Google Scholar] [CrossRef]
  24. Yohannes, B.; Hill, K.M. Rheology of dense granular mixtures: Particle-size distributions, boundary conditions, and collisional time scales. Phys. Rev. E 2010, 82, 061301. [Google Scholar] [CrossRef] [PubMed]
  25. Kézdi, A. Soil Physics-Selected Topics-Developments in Geotechnical Engineering-25, 1st ed.; Elsevier Publishing Company, Limited: Essex, UK, 1979; Volume 1. [Google Scholar]
  26. Kenney, T.C.; Lau, D. Internal stability of granular filters. Can. Geotech. J. 1985, 22, 215–225. [Google Scholar] [CrossRef]
  27. Kenney, T.C.; Lau, D. Internal stability of granular filters: Reply. Can. Geotech. J. 1986, 23, 420–423. [Google Scholar] [CrossRef]
  28. Burenkova, V. Assessment of suffusion in non-cohesive and graded soils. In Title of Filters in Geotechnical and Hydraulic Engineering, Proceedings of Geo-Filters, Karlsruhe, Germany, 20–22 October 1992; A A Balkema: Rotterdam, The Netherlands, 1993; pp. 357–360. [Google Scholar]
  29. Reddi, L.N.; Lee, I.M. Comparison of internal and surface erosion using flow pump tests on a sand-kaolinite mixture. Geotech. Test. J. 2000, 23, 116–122. [Google Scholar]
  30. Marot, D.; Rochim, A. Assessing the susceptibility of gap-graded soils to internal erosion: Proposition of a new experimental methodology. Nat. Hazards 2016, 83, 365–388. [Google Scholar] [CrossRef]
  31. Hieu, D.M.; Kawamura, S. Internal erosion of volcanic coarse grained soils and its evaluation. Int. J. Geomate 2017, 13, 165–172. [Google Scholar] [CrossRef]
  32. Rochim, A.; Marot, D. Effects of hydraulic loading history on suffusion susceptibility of cohesionless soils. J. Geotech. Geoenviron. 2017, 143, 04017025. [Google Scholar] [CrossRef]
  33. Wilcock, P.R.; Crowe, J.C. Surface-based transport model for mixed-size sediment. J. Hydraul. Eng. 2003, 129, 120–128. [Google Scholar] [CrossRef]
  34. Hillel, D. Fundamentals of Soil Physics, 1st ed.; Academic Press: San Diego, CA, USA, 1980. [Google Scholar]
  35. Khilar, K.C.; Fogler, H.S. Model for piping-plugging in earthen structures. J. Geotech. Eng. 1985, 111, 833–846. [Google Scholar] [CrossRef]
  36. Lau, D. Stability of Particle Grading of Compacted Granular Filters. Ph.D. Thesis, University of Toronto, Toronto, ON, Canada, 1984. [Google Scholar]
Figure 1. Unit volume of seepage field as the modeling framework.
Figure 1. Unit volume of seepage field as the modeling framework.
Water 11 01439 g001
Figure 2. (a) Seepage field and (b) its equivalent tubes within a unit volume.
Figure 2. (a) Seepage field and (b) its equivalent tubes within a unit volume.
Water 11 01439 g002
Figure 3. Effect of GSD in suffusion process.
Figure 3. Effect of GSD in suffusion process.
Water 11 01439 g003
Figure 4. Flow chart of simulation program (Eq(s): Equation(s)).
Figure 4. Flow chart of simulation program (Eq(s): Equation(s)).
Water 11 01439 g004
Figure 5. Designs of experimental devices.
Figure 5. Designs of experimental devices.
Water 11 01439 g005
Figure 6. GSD curves of soils in (a) internal stable runs; (b) internal unstable runs in the studies of Skempton and Brogan [9] (dashed lines) and Wan and Fell [1] (solid lines).
Figure 6. GSD curves of soils in (a) internal stable runs; (b) internal unstable runs in the studies of Skempton and Brogan [9] (dashed lines) and Wan and Fell [1] (solid lines).
Water 11 01439 g006
Figure 7. (a) Simulation domain and (b) loading history of hydraulic gradient of Run B by Skempton and Brogan [9].
Figure 7. (a) Simulation domain and (b) loading history of hydraulic gradient of Run B by Skempton and Brogan [9].
Water 11 01439 g007
Figure 8. Calibration: dashed lines for computation results of runs in Skempton and Brogan [9]; solid lines for computation results of runs in Wan and Fell [1]; Solid points for test results of runs in Skempton and Brogan [9]; and hollow points for runs in Wan and Fell [1].
Figure 8. Calibration: dashed lines for computation results of runs in Skempton and Brogan [9]; solid lines for computation results of runs in Wan and Fell [1]; Solid points for test results of runs in Skempton and Brogan [9]; and hollow points for runs in Wan and Fell [1].
Water 11 01439 g008
Figure 9. Verification: dashed lines for computation results of runs in Skempton and Brogan [9]; solid lines for computation results of runs in Wan and Fell [1]; solid points for test results of runs in Skempton and Brogan [9]; and hollow points for runs in Wan and Fell [1].
Figure 9. Verification: dashed lines for computation results of runs in Skempton and Brogan [9]; solid lines for computation results of runs in Wan and Fell [1]; solid points for test results of runs in Skempton and Brogan [9]; and hollow points for runs in Wan and Fell [1].
Water 11 01439 g009
Figure 10. Three stages of suffusion process for Run 10 in Wan and Fell [1].
Figure 10. Three stages of suffusion process for Run 10 in Wan and Fell [1].
Water 11 01439 g010
Figure 11. Comparison between computed critical hydraulic gradient and experimental critical hydraulic gradient.
Figure 11. Comparison between computed critical hydraulic gradient and experimental critical hydraulic gradient.
Water 11 01439 g011
Figure 12. Porosity distribution for Run A in Skempton and Brogan [9] when i < ic = 0.28.
Figure 12. Porosity distribution for Run A in Skempton and Brogan [9] when i < ic = 0.28.
Water 11 01439 g012
Figure 13. Porosity distribution for Run A in Skempton and Brogan [9] when i = ic = 0.28.
Figure 13. Porosity distribution for Run A in Skempton and Brogan [9] when i = ic = 0.28.
Water 11 01439 g013
Figure 14. Change of shear stress and critical shear stress with porosity (τ-n and τc-n curves).
Figure 14. Change of shear stress and critical shear stress with porosity (τ-n and τc-n curves).
Water 11 01439 g014
Table 1. Details of parameter calibration.
Table 1. Details of parameter calibration.
ParametersCalibrated ValuesNotes
Runs in Skempton and Brogan [9] (for Non-Cohesive Soils)Runs in Wan and Fell [1] (for Cohesive Soils)
A0.0210.053Following Wilcock and Crowe [25]
B0.120.12
b0.60.6
α0.0050.005Following Fujisawa et al. [20]
ΓBased on Equation (14)Based on Equation (14)

Share and Cite

MDPI and ACS Style

Feng, Q.; Ho, H.-C.; Man, T.; Wen, J.; Jie, Y.; Fu, X. Internal Stability Evaluation of Soils. Water 2019, 11, 1439. https://doi.org/10.3390/w11071439

AMA Style

Feng Q, Ho H-C, Man T, Wen J, Jie Y, Fu X. Internal Stability Evaluation of Soils. Water. 2019; 11(7):1439. https://doi.org/10.3390/w11071439

Chicago/Turabian Style

Feng, Qingfeng, Hao-Che Ho, Teng Man, Jiaming Wen, Yuxin Jie, and Xudong Fu. 2019. "Internal Stability Evaluation of Soils" Water 11, no. 7: 1439. https://doi.org/10.3390/w11071439

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop