An Assessment Model for the Erosion Occurrence of Gap-Graded Sand-Gravel Soils under Variable Seepage Direction

Abstract

The movement of particles caused by erosion is one of the main reasons for the destruction of projects, such as dams, tunnels, and foundation pits. This study highlights a theoretical model to assess the occurrence of erosion in gap-graded, sand-gravel soils under variable seepage direction based on the critical hydraulic conditions of particle initiation. The model introduced the effects of relative exposure degree, relative hidden degree, and seepage direction by considering the difference in particle initiation conditions. On the basis of the variable-section capillary tube model formed by the skeletal pores, the mechanical analysis of the movable particles in the pores was performed, and the formulas for the critical hydraulic conditions were obtained according to the moment balance equation. Subsequently, the coupled CFD-DEM method and the available experimental data were used for validation. The comparison revealed a deviation of 0.0268 for the mean of the ratio between the calculated and simulated values compared to 1, with a covariance (COV) of 0.0344. Further, the mean value of the ratio between the calculated and test values compared to 1 had a maximum deviation of 0.095 and a covariance (COV) of 0.0143. The high degree of agreement between the data proved that the theoretical model can assess the occurrence of erosion more accurately. Finally, based on the theoretical model, the study further explored the effects of seepage direction and relative particle position on the variability in particle initiation conditions, thus finding that, unlike in other studies, the effect of seepage direction was not linear.

1. Introduction

The internal erosion of soils under the effect of groundwater seepage can seriously affect the performance of geotechnical buildings or foundations, which is one of the biggest threats responsible for the destruction of dams, tunnels, and foundation pits [1,2,3]. Most of the related studies have focused on both geometric and hydraulic conditions affecting the occurrence of erosion. When seepage hydraulic conditions exceed a critical threshold, soils that are considered potentially unstable due to their geometric conditions will exhibit particle migration, and the seepage erosion of soils occurs at that time [4].

Geometric conditions include particle shape [5,6], particle gradation [7,8,9], compactness [10], etc. Many scholars have proposed various methods for assessing the internal stability of soils based on different theories and experimental data, such as Kézdi [7], Kenney [8], Aberg [9], and Liu [10], and the method [9,10] of determining internal stability based on soil particle gradation and relative density has been validated by a large number of studies as well as engineering practices. Hydraulic conditions include the magnitude and direction of the seepage velocity or hydraulic gradient, which is the external condition that causes the soil particles to initiate and erosion to occur [11,12,13].

Obviously, a large number of theoretical and experimental studies have been conducted by numerous scholars in order to accurately assess the occurrence of erosion in unstable soils. For example, Kondratchev combined the characteristics of particle movement along the pore, based on the theory of force equilibrium of the particles, to derive the formula for calculating the critical hydraulic gradient [14]:

$$i_{\mathrm{cr}}=\left(G_{\mathrm{s}}-1\right)/\left(1+0.43d_1^2/d^2\right)$$

(1)

In this equation, i_cr is the hydraulic gradient that is critical for erosion to occur. $G_{\mathrm{s}}$ is the relative density of the soil, $G_{\mathrm{s}}={\gamma }_{\mathrm{s}}/{\gamma }_{\mathrm{w}}$, where ${\gamma }_{\mathrm{s}}$ is the particle weight, ${\gamma }_{\mathrm{w}}$ is the water weight, d₁ is the equivalent pore diameter of non-cohesive soils, and d is the particle diameter of mobile particles. In addition, a large number of model test studies of erosion have shown that particle size [15], particle shape [15,16], particle gradation [15,17,18], porosity [19,20], and fines content [21,22] all have an effect on the critical hydraulic conditions of particle initiation.

However, the above studies can only be used to determine the hydraulic conditions for particle initiation in the simple cases of upward or horizontal seepage, and they were not applicable to two-dimensional seepage erosion, such as dams, foundation pits, and tunnels. Therefore, the study of hydraulic conditions started to consider the difference in particle initiation conditions caused by the seepage direction. For example, Indraratna [23] proposed a theoretical model for the migration of fine particles in soil along pore channels based on the assumption that the seepage direction is parallel to the pore channels. Xu [24] performed a force and moment equilibrium analysis of the particles on the basis of this model and discussed the effect of the particle initiation method. Zhang [25] recalculated the drag force of moving particles considering the influence of the surrounding particles and proposed a new calculation method for the critical hydraulic gradient. Huang [26] obtained the calculation formula for the critical hydraulic gradient considering the seepage direction and different exposure positions relative to the fixed particles based on the limiting moment balance condition for rolling around fixed particles. Although the effect of seepage direction was considered in the study, the verification of particle initiation in each seepage direction was not seen, and the accuracy of the critical hydraulic conditions obtained from the calculations needs to be confirmed. In terms of test studies, although Richards [27] and Liang [28] assessed the effect of seepage direction on the critical seepage velocity of erosion occurrence based on self-designed test equipment, the effect of each seepage direction on the critical hydraulic conditions under two-dimensional seepage still needs further study due to the limitations of the test conditions.

In summary, previous research on the critical hydraulic conditions of particle initiation has ignored the difference between the particle initiation conditions and the verification of the calculation accuracy under the variable seepage conditions. There is no research result for the critical hydraulic conditions of particle initiation that can integrate the influence of the seepage direction and the relative position of the particles. It is important to note that the difference in particle initiation conditions is not only due to the seepage direction but also due to the random distribution of particles in the soil. Therefore, it is of great significance to deeply understand and verify the hydraulic conditions of particle initiation that are affected by the direction of seepage and the relative position of particles in assessing soil erosion under different conditions.

In this paper, based on the critical hydraulic conditions of particle initiation, a theoretical model was proposed to assess the occurrence of erosion in gap-graded sand-gravel soil under variable seepage direction. A capillary tube model formed by skeleton pores was used in the model, and mechanical analysis of the movable particles in the pores was carried out. The influencing factors of seepage direction, relative exposure degree, and relative hidden degree were introduced. According to the moment balance equation, the critical hydraulic condition of particle initiation was obtained to assess the occurrence of erosion. The validity and superiority of the theoretical model for assessing erosion occurrence under variable seepage direction were verified using discrete element numerical simulations and available experimental data. In addition, the effects of seepage direction, relative exposure degree, and relative hidden degree on the critical hydraulic conditions for erosion occurrence were discussed and analyzed. It should be pointed out that the derivation process of the theoretical analysis was based on the assumptions of d₁ pore diameter, spherical particles, and laminar flow in the pore channel.

2. Prediction of Critical Seepage Velocity for Erosion Occurrence under Variable Seepage Direction

2.1. Force Analysis of Movable Particles

According to the difficulty of the movement of non-cohesive soil particles under seepage, the particles can be divided into skeletal particles with larger particle sizes and loose particles with smaller particle sizes [10]. According to the characteristics of particle movement under erosion, a regular variable-section capillary tube model is used to describe the irregular pore system formed by the skeletal particles. The migration of loose particles in the pore can be approximated as the movement of particles in a circular tube of fluid [29]. The pore channel in the capillary tube model is parallel to the flow direction, and the minimum diameter, d₁, and maximum diameter, d₂, of the skeleton pore channel are calculated using Formulas (2) and (3) [30].

$$d_1=2.67\frac{n_{\mathrm{a}}}{1-n_{\mathrm{a}}}\frac{D_{\mathrm{h}}}{\alpha }$$

(2)

$$d_2=1.86d_1$$

(3)

where α is the shape coefficient of soil particles, D_h is the effective particle diameter of skeletal particles, and n_a is the porosity of the hypothetical soil composed of skeletal particles.

When the particle diameter is not larger than the minimum diameter of the pore, d₁, these particles can be taken out of the soil, as long as the hydraulic conditions are met, and the part of the loose particles will be called movable particles. When the particle diameter is between d₁ and d₂, this part of the loose particles only moves within a certain range of the pore channel, and may even block the pore. Therefore, this study assumed that the pore diameter is d₁.

The initiation of non-cohesive soil particles under seepage can be divided into three different modes: rolling, sliding, and lifting [31]. The previous related studies have shown that the rolling mode presents the lowest critical seepage velocity compared to the other modes [24,32]. Moreover, in seepage fields, the flow velocities are usually small and rolling conditions are relatively easy to achieve. Therefore, in this study, the initiation mode of the particles was considered to be rolling.

In addition, the particle diameter determines that the particles do not exist as a single particle in the pore channel, so the particle initiation is influenced by its surrounding particles. Specifically, two key factors—the relative hidden degree and the relative exposure degree—were introduced to consider the difference in particle initiation caused by the different relative positions of particles when analyzing the force on particles.

As shown in Figure 1, the forces on the movable particles A under erosion include the drag force, F_D, of the flow; the underwater gravity, G′; the hydrostatic pressure, F_P; the support force, N, of the surrounding particles; and the frictional force, F_f, between the particles. The drag force can be calculated by the resistance that the particles overcome when moving in the viscous laminar flow [33], its magnitude is related to the flow velocity of the pore, and the direction is along the theoretical bed formed by the surrounding particles, which is determined by Formula (4). The hydrostatic pressure is the water pressure difference acting on the projected area of particles [29]; its magnitude is related to the hydraulic gradient. Combined with Darcy’s law, its direction is the same as the flow velocity of the pore, which is determined by Formula (5). The underwater gravity is determined by Formula (6). When particle A rolls around the contact point with the surrounding B particles, the moment generated by the rolling friction force, F_f, between the particles and the supporting force, N, of the surrounding particles is zero, which is ignored in the calculation.

$$F_{\mathrm{D}}=3\mathit{\pi }{\mu }_{\mathrm{w}}dv$$

(4)

$$F_{\mathrm{P}}=\frac{\pi d^2}{4}{\delta }_{\mathrm{s}}i{\gamma }_{\mathrm{w}}$$

(5)

$$G^{\prime }=\frac{1}{6}\mathit{\pi }({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}})d^3$$

(6)

In the above equations, d is the diameter of movable A particles, ${\mu }_{\mathrm{w}}$ is the viscosity of water, i is the hydraulic gradient at both ends of the pipe, and ${\delta }_{\mathrm{s}}$ is the length of the pore channel, which can be taken as ${\delta }_{\mathrm{s}}=d$.

According to the rolling critical state of the movable particle A, the moment balance equation is established with the contact point, O, as the center.

$$F_{\mathrm{P}}\sin \beta \cdot \frac{d}{2}+F_{\mathrm{D}}\sin \psi \cdot \frac{d}{2}-G\sin \left(\theta -\beta \right)\cdot \frac{d}{2}=0$$

(7)

where $\theta$ represents the angle between the pore channel direction and the vertical direction, and $\beta$ and $\psi$ are the angles associated with the position of the movable particle.

2.2. Analysis of the Relative Exposure Degree and Relative Hidden Degree

The distribution position of movable particles directly affects the critical conditions of particle initiation. The relative position of the movable particle A is represented by the exposure degree and the hidden degree [34]. ${\mathsf{\Delta }}_1$ is the vertical distance along the pore channel between the lowest point of the particle and the point of contact with the following particle, and ${\mathsf{\Delta }}_2$ is the vertical distance between the lowest point of the particle and the contact point with the surrounding particle along the bed surface. According to its correlation with the radius, the relative exposure degree of movable particle A is expressed as ${\mathsf{\Delta }}_1^{\prime }$, and the relative hidden degree is expressed as ${\mathsf{\Delta }}_2^{\prime }$. From Figure 1, it can be seen that their geometric relationship with the angle associated with the particle position can be expressed as:

$$\sin \beta =\frac{R-{\mathsf{\Delta }}_1}{R}=1-{\mathsf{\Delta }}_1^{\prime },$$

(8)

$$\sin \psi =\frac{R-{\mathsf{\Delta }}_2}{R}=1-{\mathsf{\Delta }}_2^{\prime }.$$

(9)

Assuming that the distribution law of relative exposure obeys a skew-normal distribution, the mean value of relative exposure on the slope is calculated according to Formula (10) [35].

$$\overline{{\mathsf{\Delta }}_1^{\prime }}=0.215-\frac{0.04}{{\left[0.04{\cos }^2\theta +0.25{\left(\cos \theta -1\right)}^2\right]}^{1.5}}\left[\frac{1}{24}-{\left(\frac{1-\sin \theta }{2}\right)}^2\left(\frac{1}{6}+\frac{\sin \theta }{3}\right)\right]$$

(10)

Considering that both the relative hidden degree and the relative exposure degree represent the positional relationship between the movable particles and the surrounding particles, the relative hidden degree is the same as the relative exposure degree in the same seepage direction based on the same assumption.

2.3. Critical Seepage Velocity for Particle Initiation

Assuming that the pore diameter is d₁, according to Poiseuille theory [23], and assuming the number of pores per unit area [30], the flow velocity, v, of the pore system per unit area is:

$$v=n\left(\frac{{\gamma }_{\mathrm{w}}}{{\mu }_{\mathrm{w}}}\right)\left(\frac{d_1^2}{32}\right)i,$$

(11)

where n is the porosity of the soil.

Further, considering the relationship between the flow velocity, $v$, of the pore system and the seepage velocity, $v_s$, the seepage velocity at which erosion occurs is obtained as:

$$v_s=n^2\left(\frac{{\gamma }_{\mathrm{w}}}{{\mu }_{\mathrm{w}}}\right)\left(\frac{d_1^2}{32}\right)i.$$

(12)

Combining Darcy’s law, v = ki, the permeability coefficient, k, of the pore system is:

$$k=n^2\left(\frac{{\gamma }_{\mathrm{w}}}{{\mu }_{\mathrm{w}}}\right)\left(\frac{d_1^2}{32}\right).$$

(13)

Therefore, according to the moment balance equation, the critical seepage velocity of particle initiation can be calculated using the following formula.

$$v_{\mathrm{s}}=n^2\frac{\left({\gamma }_{\mathrm{s}}-{\gamma }_{\mathrm{w}}\right)d^2{d}_1^2\sin \left(\theta -\beta \right)}{\left(48d^2\sin \beta +18d_1^2\sin \psi \right){\mu }_{\mathrm{w}}}$$

(14)

It can be seen from the calculation formula that the critical seepage velocity is affected by the density, the diameter of the movable particles, soil gradation, porosity, the seepage direction, and the relative position of the particles. When the soil parameters are constant, the critical seepage velocity only varies with ${\mathsf{\Delta }}_1$, ${\mathsf{\Delta }}_2$, and $\theta$.

3. Simulation of Erosion Occurrence under Variable Seepage Direction

3.1. Numerical Model

The discrete element method (DEM), based on Newton’s second law of motion and the force-displacement law, is a widely used analytical method to simulate the motion behavior of large-mass dispersed particles. According to this feature, the coupled DEM–CFD method is applied to simulate the erosion of gap-graded sand-gravel soils in different seepage directions, and the proposed formula for calculating the critical velocity of fine particle initiation is verified.

Similar to the forces on the particles in the theoretical analysis, in the discrete element simulation, the forces exerted by the fluid on a particular particle, i, include the drag force ($F_{\mathrm{D}}$), the fluid pressure gradient force ($F_{\mathrm{P}}$), and the buoyancy force ($F_{\mathrm{B}}$),

$$F_i^f=F_{\mathrm{D}}+F_{\mathrm{P}}+F_{\mathrm{B}},$$

(15)

where the drag force is calculated using the expression proposed by Di Felice [36], which is valid for dense particle states and a wide range of Reynolds numbers.

$$\left\{\begin{array}{l}{F}_{\mathrm{D}}=\frac{1}{8}{C}_{\mathrm{d}}{\rho }_f\pi d^2(u-v)|u-v|n^{-\chi }\\ C_{\mathrm{d}}={\left(0.63+\frac{4.8}{\sqrt{{\mathrm{Re}}_{\mathrm{p}}}}\right)}^2\\ R{e}_{\mathrm{p}}=\frac{n{\rho }_{f}d|u-v|}{{\mu }_{\mathrm{w}}}\\ \chi =3.7-0.65\exp \left(-\frac{{\left(1.5-{\log }_{10}{\mathrm{Re}}_{\mathrm{p}}\right)}^2}{2}\right)\end{array}\right.,$$

(16)

where $d$ is the particle diameter, $C_{\mathrm{d}}$ is the drag force coefficient of a single spherical particle, ${\rho }_f$ is the fluid density, u is the particle velocity, and $\chi$ is the correlation function that modifies the drag force coefficient by considering the presence of other particles in the system.

The fluid pressure gradient force ($F_{\mathrm{P}}$) is represented by the following formula.

$$F_{\mathrm{P}}=\frac{1}{6}\pi d_i^3{\gamma }_{\mathrm{w}}i$$

(17)

The buoyancy force (F_B) is represented by the following formula.

$$F_{\mathrm{B}}=\frac{1}{6}\pi d_i^3{\rho }_{f}g$$

(18)

In addition, the contact force for particle, i, is calculated as [34]:

$$F_{ij}^c=-F_{\mathrm{n}}\cdot n_{\mathrm{c}}+{\mathit{F}}_{\mathrm{s}},$$

(19)

where $F_{\mathrm{n}}$ is the normal force acting on the contact surface, ${\mathit{F}}_{\mathrm{s}}$ is the shear force acting on the contact surface, and ${\mathit{n}}_{\mathrm{c}}$ is the unit normal vector on the contact surface.

The linear model uses a linear spring with constant normal stiffness, ${\mathit{k}}_{\mathrm{n}}$, and shear stiffness, ${\mathit{k}}_{\mathrm{s}}$, to describe the intrinsic mechanical behavior of the interparticle contact. $F_{\mathrm{n}}$ and ${\mathit{F}}_{\mathrm{s}}$ are calculated as follows [37].

$${\mathit{F}}_{\mathrm{n}}=\left\{\begin{array}{l}{k}_{\mathrm{n}}{g}_{\mathrm{s}},g_{\mathrm{s}}<0\\ 0,\mathrm{otherwise}\end{array}\right.$$

(20)

$${\mathit{F}}_{\mathrm{s}}=\left\{\begin{array}{l}{\mathit{F}}_{\mathrm{s}}^{\mathrm{c}},‖{\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}‖<{\mathit{F}}_{\mathrm{s}}^{\mathsf{\mu }}\\ {\mathit{F}}_{\mathrm{s}}^{\mathsf{\mu }}\left(\frac{{\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}}{‖{\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}‖}\right),\mathrm{otherwise}\end{array}\right.$$

(21)

$${\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}={\left({\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}\right)}_0-k_{\mathrm{s}}\mathsf{\Delta }{\mathsf{\delta }}_{\mathrm{s}}$$

(22)

$${\mathit{F}}_{\mathrm{s}}^{\mathsf{\mu }}=-\mathsf{\mu }{F}_{\mathrm{n}}$$

(23)

where $g_{\mathrm{s}}$ denotes the spacing between the contact surfaces; ${\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}$ is the calculated shear force; ${\left({\mathit{F}}_{\mathrm{s}}^{\mathrm{c}}\right)}_0$ is the initial shear force; $\mathsf{\Delta }{\mathsf{\delta }}_{\mathrm{s}}$ is the relative shear displacement increment; ${\mathit{F}}_{\mathrm{s}}^{\mathsf{\mu }}$ is the friction force at the contact, i.e., the maximum shear force; and $\mathsf{\mu }$ denotes the friction coefficient between the friction surfaces.

So, at any moment, t, the equation of motion governing particle, i, is

$$\left\{\begin{array}{l}{m}_i\frac{\partial {\mathit{u}}_i}{\partial t}={\displaystyle \sum _{i=1}^{n_i^c}{F}_{ij}^c}+F_i^f+m_i\mathit{g}\\ I_i\frac{\partial {\omega }_i}{\partial t}={\displaystyle \sum _{i=1}^{n_i^c}{M}_{ij}}\end{array}\right.,$$

(24)

where m_i and I_i are the mass and moment of inertia of particle i, respectively; $u_i$ and ${\omega }_i$ are the linear and angular velocities of particle, respectively; $\mathit{g}$ is the acceleration of gravity. The forces involved include the contact force, $F_{ij}^c$, and torque, $M_{ij}$, acting on particle i by particle j or the wall and the particle–fluid interaction forces, $F_i^f$, acting on particle i; $n_i^c$ is the number of contacts of particle i.

In order to make the analysis feasible and effective, the gradation of the studied sand-gravel soil sample was simplified to some extent, and the key characteristics of the missing fraction were retained. At the same time, considering the influence of calculation speed, the maximum and minimum grain diameter ratios were controlled to be within 10. The coarse-grain group used gravels with grain diameters from 2.0 to 5.0 mm, and the fine-grain group consisted of only fine particles with a grain diameter of 0.6 mm. The shape of the particles used in the simulation was spherical. Three kinds of gap-graded sand-gravel soils were selected, and the gradation is shown in Figure 2. The main difference in the sample gradation was the content of fine particles, fc, the mass fractions of the fine particles were 15%, 20%, and 25%, respectively. In the process of sample preparation and testing, the acceleration of gravity was set to 9.81 m/s², and the direction was vertical downward. The sample was prepared using the layered under-compaction method [38], and the sample size was 20 × 20 × 20 mm. The linear model was selected for the contact model, and

Table 1. Table of model parameters.

Model	Parameter Type	Value
Particle model	Fine particle diameter/mm	0.6
	Coarse particle diameter/mm	2.0~5.0
	Particle density/(kg/m3)	2650
	Friction coefficient	0.5
	Effective modulus/Pa	1 × 109
	Normal-to-shear stiffness ratio	1.25
	Initial porosity	0.3
Fluid model	Fluid density/(kg/m3)	1000
	Dynamic viscosity/($\mathrm{Pa}\cdot \mathrm{s}$)	0.001

shows the parameters related to DEM particles and contacts.

Taking the upward seepage of Sample 1 as an example, after the soil sample was prepared, a series of equal-sized hexahedral grids were created to describe the fluid region with the number of 2 × 2 × 2. At the same time, the upper wall was replaced by a grid wall with a uniform square, thus forming a particle flow numerical model with only the upper grid wall as a free opening boundary, as shown in Figure 3.

In addition, to facilitate the simulation analysis, 8 equal-volume measurement regions were established to monitor and record the porosity of the sample. The measuring region was a sphere with a radius of 5 mm, and the center of the measuring region was at the center of each fluid element.

3.2. Erosion Simulation Test

Firstly, the parameter of the fine particle initiation rate was defined to describe the erosion effects of different seepage velocities, and its value is the percentage of the number of slightly moving particles along the seepage direction relative to the number of all particles. In the erosion simulation of the sample, the maximum velocity of particle initiation was sought by observing the migration state of fine particles at different seepage velocities and monitoring the initiation rate of fine particles. At a certain seepage velocity, most of the fine particles in the sample had obvious displacement. With the application of the graded seepage velocity, the initiation rate of fine particles gradually increased and stabilized at a certain value, and the initiation rate of particles was not affected by the continuous increase in the seepage velocity; that is to say, the seepage velocity was the critical velocity of particle initiation.

During the simulation process, the increasing vertical seepage velocity was applied, in turn, and the velocity change rate was 0.05 cm/s. The initiation states of fine particles under different seepage velocities are shown in Figure 4.

It can be seen from Figure 4 that when the seepage velocity was $v=0.00 \mathrm{cm}/\mathrm{s}$, the fine particles in the soil sample basically remained motionless; when the seepage velocity was $v=0.40 \mathrm{cm}/\mathrm{s}$, the particles in the fluid element at the lower right corner of the selected soil sample section had obvious vertical displacement. However, with the increase in seepage velocity, the number of moving particles in the cross-section increased until the seepage velocity was $v=0.7 \mathrm{cm}/\mathrm{s}$, and all the fine particles in the cross-section were almost moving and had obvious vertical displacement.

According to the initiation state of fine particles at each level of the seepage velocity, it was found that the particle initiation was regionalized, and the fine particles in each region at the seepage velocity were initiated, in turn. It is believed that this is caused by the difference in the local porosity of each fluid element because the sample is not uniform, thus resulting in differences in the local porosity of each region. The smaller the local porosity of the region is, the greater the pore flow velocity will be under constant seepage velocity and the stronger the erosion will be, and the fine particles will be the first to move. For the whole soil sample, only when most of the fine particles in the whole cross-section area are eroded can it be considered that the velocity reaches the critical velocity of fine-particle initiation under the overall porosity.

At the same time, the fine-particle displacement and the rate of fine-particle initiation were monitored during the simulation. Based on the principle of discrete element simulation, when the fluid–particle interaction force is applied to all particles, the fine particles satisfying the equation of motion will move slightly. Therefore, for all samples, the rate of fine-particle initiation will increase rapidly and tend to be stable in a very short time. This phenomenon was also confirmed by the monitoring results in the simulation. The variation curves of the initiation rate of fine particles at each level of the seepage velocity are shown in Figure 5.

It can be seen from the figure that the maximum initiation rate of fine particles was maintained at about 80%, and the minimum seepage velocity to reach the maximum initiation rate was 0.65 cm/s. That is, at this seepage velocity, the fine particles with a 0.6 mm particle diameter were eroded. Compared with the theoretical calculation results, the validity of the critical velocity calculation formula was initially verified.

The prediction of the critical velocity for particle initiation under variable seepage direction was achieved by changing the flow direction of the fluid element. The simulations were divided into 5 groups according to the magnitude of the seepage direction angle—namely 90°, 120°, 135°, 150°, and 180°—and 15 groups of simulations were completed for 3 kinds of samples. The critical velocity of particle initiation obtained from the formula calculation and numerical simulation under the different conditions is shown in

Table 2. The critical velocity, calculated by formula and numerical simulation under different conditions.

Sample	Direction of Seepage/°	Calculation Results/(cm/s)	Simulation Results/(cm/s)
Sample 1	90	0.557	0.55
	120	0.622	0.60
	135	0.678	0.65
	150	0.815	0.75
	180	0.706	0.65
Sample 2	90	0.664	0.65
	120	0.742	0.75
	135	0.808	0.80
	150	0.971	0.95
	180	0.841	0.85
Sample 3	90	0.864	0.85
	120	0.965	0.95
	135	1.050	1.05
	150	1.263	1.20
	180	1.094	1.05

, and the comparative analysis showed that the formula calculation results were basically consistent with the numerical simulation results.

The comparison of the critical velocity obtained from the calculation formula and numerical simulation is shown in Figure 6. If the position of the scatter distribution is closer to the diagonal line in the figure, it means that the data are more consistent and the formula predicts the critical velocity of particle initiation better. As shown in Figure 6, the data points are basically distributed along the contour line, and the mean of the ratio between test and predicted values compared to 1 had a deviation of 0.0268, a covariance (COV) of 0.0344, and a maximum deviation of 6.145%, which is within ±7%, overall.

4. Verification of Existing Test Data

Soil samples from existing erosion tests [18,20] were used to verify the prediction formula for the critical velocity. One of the data sources was the indoor erosion test conducted by Skempton based on four kinds of sand-gravel materials prepared by themselves, in which the test soil samples C and D were internally stable. The other is the vertical and horizontal seepage tests using a specially developed test device on five different kinds of non-cohesive soils prepared by Ahlinhan. In the tests, the initial relative density of the soil samples was changed and soil samples A1 and A2 were internally stabilized soils. Therefore, soil samples A, B, E1, E2, and E3 were selected for verification, respectively. It should be noted that different tests have different critical conditions for identifying the occurrence of piping, and the formula cannot select a uniform particle size as the identified particle diameter for the occurrence of piping. Therefore, based on the analysis of existing studies and for safety reasons, for the Skempton piping test, it is assumed that when the soil particle loss amount reaches 5% of the soil mass, the corresponding seepage velocity is taken as the critical velocity of soil particle movement. That is, the particle diameter was taken as d = d₅. Similarly, d = d₁₀ was chosen for the Ahlinhan piping test. The parameters of test soil samples are shown in

Table 3. Parameters of test soil samples.

Sample	Initial Relative Density	Natural Void Ratio e	Limiting Particle Diameter xa [9]/mm	Effective Particle Diameter of Skeletal Particles Dh/mm	Porosity of the Hypothetical Soil na	Minimum Pore Diameter d1/mm	Particle Diameter d
A	-	0.52	1.000	3.235	0.44	0.568	$d_5$
B	-	0.59	0.431	2.488	0.43	0.417	$d_5$
E1	0.5	0.62	0.537	1.279	0.47	0.500	$d_{10}$
	0.6	0.60			0.46	0.486	$d_{10}$
	0.85	0.55			0.44	0.452	$d_{10}$
E2	0.41	0.55	0.439	1.191	0.46	0.443	$d_{10}$
	0.5	0.52			0.45	0.427	$d_{10}$
	0.85	0.41			0.41	0.361	$d_{10}$
E3	0.75	0.52	0.96	2.166	0.44	0.386	$d_{10}$
	0.95	0.46			0.42	0.355	$d_{10}$

.

The critical velocity of the test soil sample was predicted using the formula, and the scatter points of the predicted data and the test data were plotted in Figure 7. As can be seen from Figure 7a, under the horizontal seepage, the scatter points are located roughly in the area on both sides of the diagonal. The mean of the ratio between test and predicted values compared to 1 had a deviation of 0.095, and the covariance (COV) was 0.0143, so the two sets of data had a high degree of agreement. The same good agreement was found between the two sets of data under vertical seepage, as shown in Figure 7b, with a deviation of 0.0152 from 1 for the mean of the ratio between the test and predicted values and a covariance (COV) of 0.0199. After comparison with the results calculated by the literature method [39], it was clear that the prediction of the critical velocity for particle initiation using the formula was highly accurate for both horizontal and vertical seepage, and the validity of the formula was verified.

In addition, under vertical seepage, because the Kondratchev method considered the same forces on the particles as this study, the data points predicted by it were plotted in Figure 7b, and it can be seen that some of the data points of the critical velocity of particle movement predicted by the Kondratchev method are far from the diagonal region; the data error was larger and the mean of the ratio between test and predicted values compared to 1 had a deviation of 0.370. It can be seen that the accuracy of the critical hydraulic condition formula proposed in this paper was significantly better than the prediction results provided by the Kondratchev method.

5. Analysis of the Differences in Particle Initiation Conditions

5.1. Seepage Direction

The simulated sample was used to explore the relationship between the influence of seepage direction on the critical velocity of particle initiation, and the variation curves of critical velocity under different seepage directions were plotted according to the comparative data of critical velocity in Table 2, as shown in Figure 8. It can be seen from Figure 8 that, for seepage erosion in different directions—although the critical velocity of the vertical seepage was greater than that of horizontal seepage, and the prediction result was the same as the Ahlinhan [20] seepage test result—the trend was not linear. With the increase in seepage direction angle, the critical velocity showed a trend of increasing and then decreasing. When the seepage direction angle is smaller than a certain angle, the critical flow velocity increases with the increase in the seepage direction angle, which is the same as the current universal research conclusion [26,27,28]. When the seepage direction angle is larger than the specific angle, the critical velocity tends to decrease with the increase in the seepage direction angle. This phenomenon is a new finding in this article due to the lack of relevant angular conditions in the test study.

It is well-explained in the theoretical study from the perspective of particle force analysis: under the influence of the relative exposure degree, the gravity of the particles in different seepage directions contributes differently to the critical driving force for initiation. Both Richards [27] and Zhou [40] emphasized the effect of particle gravity in their studies. When the movable particles in the pore roll around the contact point with the downstream particles, the effect of the relative exposure position makes the gravity of the particles act as the driving force for particle initiation when the seepage direction angle is smaller than the exposure angle of the particles, and the particles can be moved in the pore by their own gravity. With the increase in the seepage direction angle, the driving role of the weight of the particles becomes smaller and smaller. When the seepage direction angle is equal to the exposure angle of the particles, the driving effect of the particles’ gravity disappears. When the seepage direction angle is greater than the exposure angle, the particles’ gravity becomes the force that hinders the movement of the particles, and as the seepage direction continues to increase, the hindering effect of the particles’ gravity becomes larger and larger until the particles’ gravity becomes the hindering force completely, and the critical hydraulic condition required for initiation reaches the maximum. The hindering effect of the particles’ gravity decreases again, gradually, as the seepage direction angle continues to increase, and the hydraulic conditions required for initiation are gradually reduced. As a result, the critical velocity is influenced by the relative exposure degree, with the increase in the seepage direction angle, and shows a trend of increasing and then decreasing.

5.2. The Relative Exposure Degree and the Relative Hidden Degree

Skempton test soil sample A was used to explore the influence of the relative position of the particles on the critical velocity of particle initiation. The variation curves of the critical velocity with the relative exposure degree and the relative hidden degree, under different seepage directions, are shown in Figure 9 and Figure 10.

It can be seen from Figure 9 that there is a curvilinear correlation between the relative exposure and the critical velocity due to the influence of the seepage direction. Under an increasing seepage direction angle, the critical velocity increases first and then decreases with the increase in relative exposure. When $\theta =90°$ or $\theta =180°$, the critical velocity changes obviously with the relative exposure degree. At a certain seepage direction angle in the middle, the influence curve slope of relative exposure trends toward 0. With the increase in relative exposure degree, the increment of critical velocity is very small and, at this time, the relative exposure degree has little influence on the critical velocity. As for the effect of the relative hidden degree, it can be seen from Figure 10 that the critical hydraulic condition increases with the increase in relative hidden degree.

Combined with the force analysis of particle initiation, as the relative hidden degree of particles increases, the drag force driving the particles to roll around the contact point becomes smaller and smaller. Correspondingly, a larger force is needed to overcome the gravity moment, and the critical velocity is also larger at this time. When the relative exposure degree of the particles increases on the one hand, the hydrostatic pressure driving the particles to initiate gets smaller and smaller on the other hand. With the increase in the seepage direction angle, the driving effect of the gravity on the particles rolling around the contact point will gradually weaken or the hindering effect will gradually strengthen, which will lead to the increase in the critical hydraulic condition of particle initiation. When the seepage direction angle gradually increases to the extent that the downstream particles no longer bear the gravity of the moving particles, the particle initiation only needs to overcome the gravity, not the gravity moment, and the influence of the relative exposure degree will be very small. When continuing to increase the seepage direction angle, the hindering effect of the gravity on the particle rolling around the contact point will gradually weaken with the increase in the relative exposure degree, and the critical velocity required for particle initiation will gradually become smaller. As a result, the critical velocity of particle initiation increases and then decreases with the increase in the relative exposure degree and the increase in the seepage direction angle.

6. Conclusions

Aiming at examining the differences in particle initiation conditions, an erosion occurrence assessment model, based on the critical hydraulic conditions of particle initiation under variable seepage, was proposed. The following conclusions were obtained:

• The model comprehensively considered the effects of the seepage direction and the relative position of particles. Based on the rolling theory of soil particles, the formula for calculating the critical velocity of particles under the variable seepage condition was derived. Further, the hydraulic conditions of seepage erosion in different directions were verified using discrete element numerical simulation and existing test data, which created conditions for subsequent application to two-dimensional seepage.
• In the analysis of the influencing factors, the relative hidden degree was positively correlated with the critical hydraulic conditions of particle initiation, while the relative exposure degree and the seepage direction interact and jointly influence each other. Different from the existing conclusions, the effect of seepage direction on the critical hydraulic conditions was not linear; however, the critical velocity showed a trend of increasing and then decreasing with the increase in seepage direction, due to the influence of relative exposure. Meanwhile, the critical velocity also showed a trend of increasing and then decreasing with the increase in relative exposure in different seepage directions.