Internal Stability Evaluation of Soils

: Su ﬀ usion constitutes a major threat to the foundation of a dam, and the likelihood of su ﬀ usion is always determined by the internal stability of soils. It has been veriﬁed 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 su ﬀ usion process. The model takes the combined e ﬀ ects 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 su ﬀ usion. By analyzing the simulation results of the model, the mechanism is disclosed on how soils with speciﬁc 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. Q.F.; Q.F., and Q.F. and


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 i c < 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 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 E k A k dVdt, where E k denotes the erosion rate of kth size group section defined as the volume of eroded particles per erodible area per time, A k denotes the erodible area per unit volume of kth size group, dV = dx 1 dx 2 dL, and dx 1 , dx 2 , and dL are the dimensions of the unit volume, respectively.
Water 2019, 11, x FOR PEER REVIEW 3 of 19 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 twodimensional 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: 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.

Closures
Closures of items erosion rate Ek, erodible area Ak. and hydraulic conductivity ks are listed below: The governing equations of mass and momentum conservation are as follows: where n denotes the porosity of soil mass, C denotes the solid concentration of seepage flow, v i denotes the seepage velocity in the two directions (i = 1, 2), S r denotes the saturation of soils, h denotes the local water head, k s denotes the hydraulic conductivity, x i 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 E k A k differs from EA in the study of Fujisawa et al. [15] since the dependence of E k on the size of kth group is taken into account. If the dependence is neglected, i.e., E k = E holds for each size group, the formulation by Fujisawa et al. [15] will be recovered.

Closures
Closures of items erosion rate E k , erodible area A k . and hydraulic conductivity k s are listed below: 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 E k is evaluated taking into account GSD, i.e., where τ is shear stress exerted on the soil particles, E k 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). 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., 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: where I = h/L denotes the local hydraulic gradient; g denotes the acceleration of gravity; ρ = Cρ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 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: where I = h/L denotes the local hydraulic gradient; g denotes the acceleration of gravity; ρ = Cρ 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 d k and d k+1 , the critical shear stress τ ck is evaluated based on Equations (7) and (8): where d sm 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; F s 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. "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): 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.

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, πdk 2 , and the number of such particles per unit volume, 6(1 − n)(Pk+1 − Pk)/πdk 3 . Therefore, the Ak reads, where Pk+1 and Pk denote the percentage passing by weight of soil particles for the size finer than dk+1 and dk respectively.

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

The Erodible Area A k
Within a volume of dV, the A k for the kth size group can be quantified with the total projective area of the given size d k , that is, the product of the projective area of a particle, πd k 2 , and the number of such particles per unit volume, 6(1 − n)(P k+1 − P k )/πd k 3 . Therefore, the A k reads, where P k+1 and P k denote the percentage passing by weight of soil particles for the size finer than d k+1 and d k respectively. The hydraulic conductivity k s 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 1.565 10 n 2.3475 (1 − n) 1.565 (10) where d 10 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. (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 k s 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 P k corresponding to the particle size d k , and then to specify d 10 in Equation (11). For this purpose, one can quantify the P k (t + dt) after a time period dt as

Update of d 10 in Equation
A rearrangement of Equation (11) without considering the second-order term results in Accordingly, the size of d 10 corresponding to 10% for which soil particles by weight are fine after the time period dt can be quantified approximately as where d x and d x+1 can be obtained through exhaustive method among all soil size groups, under the condition that item (P x − 10%) + (P x+1 − 10%) reaches the smallest absolute value.

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 v i 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.

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.

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. Variables in closures are updated in every cell of simulation domain, and calculations above are repeated at following time steps.

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.

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.

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

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 gapgraded soils and continuous-graded soils were used in the experiments. (see Figure 6).

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

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. The definition of gap-graded soils varies in different studies. However, among them, most gapgraded 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.

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.

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

Parameter Calibration
During the simulation, the pure water density ρw is set as 1000 kg/m 3 , 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.  (14) Based on Equation (14) 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,

Parameter Calibration
During the simulation, the pure water density ρ w is set as 1000 kg/m 3 , and the viscosity coefficient µ w is set as 10 −6 Pa·s. Besides GSD, other initial conditions of the soils (including initial porosity n 0 , initial permeability coefficient k s0 , and density of soil particles ρ s ) are also reproduced according to the experimental records.
There are two ways to obtain the initial porosity n 0 of the soil: (1) directly from the experimental records, and (2) by calculating n 0 = 1 − V s0 /V, where V s0 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).  (14) Based on Equation (14) Based on these runs, A, B, b and α are calibrated by trial and error. A 0 = 0.021, B 0 = 0.015, and b 0 = 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): 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):  Figure 8 presents the calibration results for runs in Section 3.3: 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].

q-i Relationship
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.

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

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

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

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.  [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].
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 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].
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 k s 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.

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

Determination of Critical Hydraulic Gradient
When i reaches critical hydraulic gradient i c , 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: 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].
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.

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

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 = 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 = i c ), one may notice an abrupt change of q-i curve, and this abrupt change may indicate the failure of soil fabrics.
When i < i c , 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 = i c or i > i c , more particles are lost until the failure of whole structure occurs.
Following the justification, i c 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]. 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]. Figure 11. Comparison between computed critical hydraulic gradient and experimental critical hydraulic gradient.

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 Figures 12 and 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.

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 Figures 12 and 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 < i c = 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. 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]. Figure 11. Comparison between computed critical hydraulic gradient and experimental critical hydraulic gradient.

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 Figures 12 and 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.

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.

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

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 d r ) remaining in soil fabrics. According to Equation (7), τ cr = min(τ ck ) as it has been already known that d r = min(d k ), and it will increase when suffusion proceeds. Meanwhile, as the particles are continuously eroded out of the seepage field, the porosity ρ and hydraulic conductivity k s of the soils become higher, and based on Equation (10), τ will also increase. When i = i c , τ 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 i k < i c , when i = i k , if the loading time is long enough, and all the particles with sizes d < d k will be eroded out of the soil fabrics. Under this assumed situation, d r = d k , and τ cr = τ ck . At the moment, the volume of the soils V s will be decreased by k% (ρ s is a constant), and the porosity of the soil will be n = n 0 + (1 − n 0 )k%.
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.

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%. According to the results shown in Figure 14, the critical hydraulic gradient is i c = 0.28, which is in accordance to the experimental data of Skempton and Brogan [9], and since i c < 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%.

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 i c 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.   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 h l,t in Equations (3) and (4) can be solved numerically based on the following Equation: Figure A1. Geometry of the cells.
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 h l,t in Equations (3) and (4)