Review on Numerical Simulation of the Internal Soil Erosion Mechanisms Using the Discrete Element Method

: Internal erosion can trigger severe engineering disasters, such as the failure of embankment dams and uneven settlement of buildings and sinkholes. This paper comprehensively reviewed the mechanisms of soil internal erosion studied by numerical simulation, which can facilitate uncovering the internal erosion mechanism by tracing the movement of particles. The initiation and development of internal erosion are jointly inﬂuenced by the geometric, mechanical, and hydraulic conditions, which determine the pore channels and force chains in soil. The geometric conditions are fundamental to erosion resistance, whereas the mechanical conditions can signiﬁcantly change the soil erosion resistance, and the hydraulic conditions determine whether erosion occurs. The erosion process can be divided into particle detachment, transport, and clogging. The ﬁrst is primarily affected by force chains, whereas the latter two are mostly affected by the pore channels. The stability of the soil is mainly determined by force chains and pore channels, whereas the hydraulic conditions act as external disturbances. The erosion process is accompanied by contact failure, force chain bending, kinetic energy burst of particles, and other processes due to multi-factor coupling.


Introduction
Internal erosion may be induced by the protective layer being stripped by sediment deposits [1] which are piled up by sediment transport [2], causing significant global economic losses. According to Foster [3] and Richards [4], approximately 0.5% of the embankment dams worldwide exhibit failure due to internal erosion, which is the primary failure mode, accounting for about 46% of all failure types. Dam failure induced by internal erosion and even overtopping [5] is subject to early detection, which is important to ensure their safety and security [6]. Besides, internal erosion can cause uneven settlement of buildings [7], extensive ground collapse or sinkholes [8][9][10][11][12][13] (as shown in Figure 1), damages in shield tunnels [14] and in dike foundations [15], and clogging of oil wells [16]. Internal erosion can be classified as suffusion, backward erosion, concentrated erosion, or contact erosion [17][18][19][20][21], of which induced failure was fully surveyed and reviewed by [22,23]. Suffusion is a long-term phenomenon in which fine soil particles are washed away through pore channels between coarse particles [24]. Soils susceptible to suffusion are described as internally unstable [25], and those non-susceptible to suffusion are called internally stable soils. Backward erosion refers to a process by which pipes are formed due to the removal of soil at the exit. As the erosion further develops, pipes extend in the opposite direction of the seepage. Concentrated leak erosion mainly occurs at cracks, whereas contact erosion exists at the interface of two layers of soil with significant differences in particle size.
The pore channels can be simplified as relatively large "chambers" and relatively small "constrictions" connecting the chambers [26]. The mechanism of internal erosion (as shown in Figure 2) is complex due to three coupled processes [27][28][29][30][31]: Particles are detached from the soil skeleton, the detached particles are transported through pore channels, and clogging occurs near constrictions.
Detachment is a process in which the equilibrium state is broken under the action of fluid flow and is the result of the destruction of the force chain [32][33][34]. Therefore, the most easily detached particles are fine particles that are deposited solely due to gravity or have a relatively small contact force. On this basis, Hosn [35] distinguished detached particles according to their unbalanced force. Although gravity-deposited particles cannot transmit effective stress, they provide lateral support for force-chain particles. However, the forcechain particle detachment can lead to stress redistribution and large volume deformation [36,37], significantly changing the soil structure and mechanical properties.
The transport process of the detached particles is typically determined by the constriction size distribution (CSD) [38,39]. According to Hosn [35], a detached particle can be transported through the pore channels and is considered an eroded particle if its diameter is smaller than Dc35 (i.e., 35% of the constrictions are finer than this size). Meanwhile, particle clogging may occur because the particle diameter exceeds the constriction size or arch formation occurs [7,40]. However, clogging can reduce permeability, leading to an increase in excess pore pressure and strength degradation. Finally, the arch can be destroyed by fluid flow, causing permeability reduction and abrupt erosion.
The pore channels can be simplified as relatively large "chambers" and relatively small "constrictions" connecting the chambers [26]. The mechanism of internal erosion (as shown in Figure 2) is complex due to three coupled processes [27][28][29][30][31]: Particles are detached from the soil skeleton, the detached particles are transported through pore channels, and clogging occurs near constrictions. Numerous studies have been carried out to clarify the internal soil erosion mechanism, including theoretical analyses [27,41], experiments [42][43][44], and numerical simulations [34,45,46]. Numerical simulations can accurately model the internal erosion mechanism under fluid flow by tracing the movement of detached particles, visualizing the evolution of the pore channels, and quantifying the force chain variation, unlike experiments. In early research, Hu [47,48], Luo [49], Cividini [50,51], and Sterpi [7] investigated internal erosion using the continuous medium finite element method (FEM), which cannot capture the discrete behavior of geomaterials. The introduction of the discrete element method (DEM) coupled with hydro-mechanical models allowed for modeling the soil deformation characteristic and particle migration, as well as the evolution of the pore channels and force chains at the particle scale.
Therefore, this article first briefly introduces commonly used numerical simulation methods for analyzing internal erosion and then systematically summarizes the internal erosion mechanism from geometric, mechanical, and hydraulic perspectives. Several issues that require further study are highlighted to improve the accuracy and applicability of numerical simulations and improve our understanding of soil erosion mechanisms.

Fluid-Particle Coupled Method in Numerical Simulations of Internal Erosion
Internal erosion is a coupling process of fluid flow and particle migration; thus, a Detachment is a process in which the equilibrium state is broken under the action of fluid flow and is the result of the destruction of the force chain [32][33][34]. Therefore, the most easily detached particles are fine particles that are deposited solely due to gravity or have a relatively small contact force. On this basis, Hosn [35] distinguished detached particles according to their unbalanced force. Although gravity-deposited particles cannot transmit effective stress, they provide lateral support for force-chain particles. However, the force-chain particle detachment can lead to stress redistribution and large volume deformation [36,37], significantly changing the soil structure and mechanical properties.
The transport process of the detached particles is typically determined by the constriction size distribution (CSD) [38,39]. According to Hosn [35], a detached particle can be transported through the pore channels and is considered an eroded particle if its diameter is smaller than D c35 (i.e., 35% of the constrictions are finer than this size). Meanwhile, particle clogging may occur because the particle diameter exceeds the constriction size or arch formation occurs [7,40]. However, clogging can reduce permeability, leading to Water 2021, 13, 169 3 of 17 an increase in excess pore pressure and strength degradation. Finally, the arch can be destroyed by fluid flow, causing permeability reduction and abrupt erosion.
Numerous studies have been carried out to clarify the internal soil erosion mechanism, including theoretical analyses [27,41], experiments [42][43][44], and numerical simulations [34,45,46]. Numerical simulations can accurately model the internal erosion mechanism under fluid flow by tracing the movement of detached particles, visualizing the evolution of the pore channels, and quantifying the force chain variation, unlike experiments. In early research, Hu [47,48], Luo [49], Cividini [50,51], and Sterpi [7] investigated internal erosion using the continuous medium finite element method (FEM), which cannot capture the discrete behavior of geomaterials. The introduction of the discrete element method (DEM) coupled with hydro-mechanical models allowed for modeling the soil deformation characteristic and particle migration, as well as the evolution of the pore channels and force chains at the particle scale.
Therefore, this article first briefly introduces commonly used numerical simulation methods for analyzing internal erosion and then systematically summarizes the internal erosion mechanism from geometric, mechanical, and hydraulic perspectives. Several issues that require further study are highlighted to improve the accuracy and applicability of numerical simulations and improve our understanding of soil erosion mechanisms.

Fluid-Particle Coupled Method in Numerical Simulations of Internal Erosion
Internal erosion is a coupling process of fluid flow and particle migration; thus, a suitable coupled model is required for the analysis. Generally, the DEM is used to compute the motion of particles, but there are various solutions for fluid motion.

DEM Model of Particle Migration
The motion of particles in the DEM is computed by Newton's second law that considers the gravity, contact force, and interaction force between the particles and fluid. Commonly used contact models include linear models and the nonlinear Hertz-Mindlin model. However, the contact properties considered by different scholars are substantially different, as listed in Table 1. Hu [52] and Kawano [24] adopted a large modulus which is similar to the natural sandy soil. Due to the negative correlation between the soil modulus and the DEM time step, Cai [53], Guo [54], Zou [55], Abdoulaye [56], and Tao [18,34] used a relatively small modulus to increase the time step. The influence of a relatively small modulus should be further analyzed and corrected to improve the applicability and accuracy of the simulation results.

DEM and Fluid Flow Coupled Method
Several scholars have proposed coupled methods to consider the influence of fluid flow in the DEM. Hakuno and Tarumi [57] developed a fluid network model coupled with the DEM. The flow channel between the pores was assumed to follow the Hagen-Poiseuille pipe flow, as shown in Figure 3a. Cui [16], Tang [9], and Goodarzi [58] coupled Darcy's law with the DEM. Tsuji [59] and Xu and Yu [60] proposed the volume-averaged coarse-grid approach to mesh the fluid domain into cells, in which the physical domains overlapped with the soil particles. The locally averaged Navier-Stokes (N-S) equation [61] was solved in each cell, and the impact of the particles was incorporated by considering the porosity and interaction force. This model is called the computational fluid dynamics (CFD)-DEM method (Figure 3b). The Lattice-Boltzmann method (LBM) does not solve the N-S equation directly, but solves the kinetic gas theory equation, which discretizes the fluid domains and particles into fluid nodes and solid nodes, and simulates the fluid flow by tracking the evolution of the density distribution of the fluid nodes [62], as shown in Figure 3c. The LBM-DEM has been widely used to simulate the particle-fluid interaction since it does not require solving the complex N-S equation [63][64][65]. The smoothed particle hydrodynamics (SPH) model was developed using a mesh-free method [66], which discretizes the fluid domains into fluid particles, as shown in Figure 3d Several scholars have proposed coupled methods to consider the influence of fluid flow in the DEM. Hakuno and Tarumi [57] developed a fluid network model coupled with the DEM. The flow channel between the pores was assumed to follow the Hagen-Poiseuille pipe flow, as shown in Figure 3a. Cui [16], Tang [9], and Goodarzi [58] coupled Darcy's law with the DEM. Tsuji [59] and Xu and Yu [60] proposed the volume-averaged coarse-grid approach to mesh the fluid domain into cells, in which the physical domains overlapped with the soil particles. The locally averaged Navier-Stokes (N-S) equation [61] was solved in each cell, and the impact of the particles was incorporated by considering the porosity and interaction force. This model is called the computational fluid dynamics (CFD)-DEM method ( Figure 3b). The Lattice-Boltzmann method (LBM) does not solve the N-S equation directly, but solves the kinetic gas theory equation, which discretizes the fluid domains and particles into fluid nodes and solid nodes, and simulates the fluid flow by tracking the evolution of the density distribution of the fluid nodes [62], as shown in Figure 3c. The LBM-DEM has been widely used to simulate the particle-fluid interaction since it does not require solving the complex N-S equation [63][64][65]. The smoothed particle hydrodynamics (SPH) model was developed using a mesh-free method [66], which discretizes the fluid domains into fluid particles, as shown in Figure 3d. By introducing the equation of state, the fluid pressure can be determined explicitly, and the N-S equation can be solved. The value of the fluid variable at any position can be obtained by interpolating the value of the fluid particles in the kernel.  Figure 4. The fluid network method does not require discretization and meshing, but its parameters are difficult to obtain. The Darcy model, which only needs the permeability, has a straightforward calculation, but the low Reynolds number limits its application. The CFD-DEM solves the fluid flow in cells larger than the particle size, which reduces the accuracy, but the model has relatively high computational efficiency; thus, it is widely used. The LBM-DEM can solve the fluid flow at a pore-scale far smaller than the particle size, capturing the fluid flow characteristics between the pores and achieving high calculation accuracy. However, the computing cost is high, and only less than ten thousand particles can be simulated [29,45,62,68,69]. The SPH does not suffer from accuracy reduction caused by mesh distortion in large deformation and can handle complex boundaries easily; however, its dynamic fluid pressure is not accurate. has relatively high computational efficiency; thus, it is widely used. The LBM-DEM can solve the fluid flow at a pore-scale far smaller than the particle size, capturing the fluid flow characteristics between the pores and achieving high calculation accuracy. However, the computing cost is high, and only less than ten thousand particles can be simulated [29,45,62,68,69]. The SPH does not suffer from accuracy reduction caused by mesh distortion in large deformation and can handle complex boundaries easily; however, its dynamic fluid pressure is not accurate.

Research on Factors Affecting Internal Erosion and Its Mechanism
Many scholars [69][70][71][72] unanimously agree that the initiation and development of internal erosion are jointly influenced by the geometric, mechanical, and hydraulic conditions, as shown in Figure 5. The coupling of these conditions affects the pore channels and force chains, influencing the particle detachment, transport, and clogging, which, in turn, determine the occurrence, initiation time, and the amount of erosion.

Research on Factors Affecting Internal Erosion and Its Mechanism
Many scholars [69][70][71][72] unanimously agree that the initiation and development of internal erosion are jointly influenced by the geometric, mechanical, and hydraulic conditions, as shown in Figure 5. The coupling of these conditions affects the pore channels and force chains, influencing the particle detachment, transport, and clogging, which, in turn, determine the occurrence, initiation time, and the amount of erosion.

Geometric Conditions
Geometric conditions are inherent to the soil, representing a key factor determining erosion resistance. Compared with widely graded soils, gap-graded soils are more prone to internal erosion. In gap-graded soils, the fines content (F), gap ratio (G, defined by the ratio of the diameter of the smallest coarse particle to the largest fine particle), and density

Geometric Conditions
Geometric conditions are inherent to the soil, representing a key factor determining erosion resistance. Compared with widely graded soils, gap-graded soils are more prone to internal erosion. In gap-graded soils, the fines content (F), gap ratio (G, defined by the ratio of the diameter of the smallest coarse particle to the largest fine particle), and density (D r ) have significant impacts on erosion. Besides, the particle shape is crucial for erosion resistance due to its impact on pore channels and interlocking. However, few studies have investigated the particle shape due to the difficulty of simulating irregular particles.

Effect of Geometric Conditions on Stress Distribution
Gap-graded soils can be divided into three fabric states based on two critical fines contents (S* and S max ): (1) "Underfilled" fabric when F < S* and the fines do not fill the pores; (2) "filled" fabric when S* < F < S max and the fines fill the pores; (3) "overfilled" fabric when F > S max and the coarse particles are separated from the fine particles [74], as shown in Figure 6. Skempton [75] proposed that S max = 35%, and S* was between 24% and 29% due to the density. The fabric determines the stress distribution in the soil and influences the particle detachment process. Meanwhile, it also affects the pore channels, which are the transport paths of the detached particles. The stress distribution can be quantified by the stress reduction factor α, which is the ratio of the average effective stress in the fines to the stress in the sample [75]. As the fines' content and density increase and the gap ratio decreases, the contact force between the particles strengthens, and the connectivity increases, as reported by Li [76], Ahmadi [28], and Shire [71,74]. Shire [71,74] investigated the relationship between α and the fines content, the density, and the gap ratio, as shown in Figure 7. It was found that α decreased rapidly with the gap ratio in the underfilled state since more particles were deposited due to gravity, and there were weak contacts. α remained at a relatively high level, and the contribution of the coarse and fine particles to the stress transmission was equivalent in the overfilled state. α decreased linearly with the gap ratio in the filled state ( Figure 7a). Moreover, the density had the largest impact in the sample with the largest gap ratio sample since the pore channels were wider in this sample, as shown in Figure 7b. This can confirm that the soil fabric is closely related to fines content, gap ratio, and density. The stress distribution can be quantified by the stress reduction factor α, which is the ratio of the average effective stress in the fines to the stress in the sample [75]. As the fines' content and density increase and the gap ratio decreases, the contact force between the particles strengthens, and the connectivity increases, as reported by Li [76], Ahmadi [28], and Shire [71,74]. Shire [71,74] investigated the relationship between α and the fines content, the density, and the gap ratio, as shown in Figure 7. It was found that α decreased rapidly with the gap ratio in the underfilled state since more particles were deposited due to gravity, and there were weak contacts. α remained at a relatively high level, and the contribution of the coarse and fine particles to the stress transmission was equivalent in the overfilled state. α decreased linearly with the gap ratio in the filled state (Figure 7a). Moreover, the density had the largest impact in the sample with the largest gap ratio sample since the pore channels were wider in this sample, as shown in Figure 7b. This can confirm that the soil fabric is closely related to fines content, gap ratio, and density. rapidly with the gap ratio in the underfilled state since more particles were deposited due to gravity, and there were weak contacts. α remained at a relatively high level, and the contribution of the coarse and fine particles to the stress transmission was equivalent in the overfilled state. α decreased linearly with the gap ratio in the filled state (Figure 7a). Moreover, the density had the largest impact in the sample with the largest gap ratio sample since the pore channels were wider in this sample, as shown in Figure 7b. This can confirm that the soil fabric is closely related to fines content, gap ratio, and density.
(a) (b) Figure 7. The relationship between the stress reduction factor α and the fines content, gap ratio, and density (reprinted from [71,74]). (a) The relationship between the stress reduction factor α and the fines content and gap ratio in dense samples. (b) The relationship between the stress reduction factor α and the density and gap ratio at a fines content of 30%.

Effect of Soil Fabric on Pore Channels and Force Chains
In addition to affecting the stress distribution, the soil fabric also has a substantial impact on the pore channels and soil skeleton. The pore channels are wide and have high connectivity, making the transport of detached particles less tortuous in the underfilled soil fabric. Therefore, the soil is susceptible to erosion, and the higher the fines content, Figure 7. The relationship between the stress reduction factor α and the fines content, gap ratio, and density (reprinted from [71,74]). (a) The relationship between the stress reduction factor α and the fines content and gap ratio in dense samples. (b) The relationship between the stress reduction factor α and the density and gap ratio at a fines content of 30%.

Effect of Soil Fabric on Pore Channels and Force Chains
In addition to affecting the stress distribution, the soil fabric also has a substantial impact on the pore channels and soil skeleton. The pore channels are wide and have high connectivity, making the transport of detached particles less tortuous in the underfilled soil fabric. Therefore, the soil is susceptible to erosion, and the higher the fines content, the more intense the erosion is. This phenomenon was confirmed by the numerical simulation results of Zou [55], Hu [46], and Cai [53], who found that the erosion ratio gradually increased with an increase in the fines content when the fines content was less than 30%. However, the filled and overfilled soil fabrics show more complex behavior. On the one hand, the pore channels become narrower, and the pore connectivity worsens with an increase in the fines content, resulting in a more tortuous transport path of the detached particles and increased clogging probability. On the other hand, the coarse particles, the main components of the soil skeleton, are separated by the fine particles which weaken the stability of the skeleton, as the fines content increases. In addition, the erosion of fine particles, bearing effective stress, leads to disturbance of the soil skeleton. The soil skeleton failure will happen when the soil skeleton can't resist the drag force.
Therefore, in the underfilled fabric, with an increase in the fines content, more fine particles with low stress lead to a more intense erosion. However, it is not clear in the filled and overfilled fabric due to the complex effects, which deserve careful investigation.

Effect of Gap Ratio on Pore Channels and Flow Velocity
In addition to the fines content, the gap ratio also determines the pore channels. Tao [77], Cai [53], and Zhang [68] analyzed the effect of the gap ratio on the erosion ratio. Tao [77] found that the erosion ratio decreased gradually as the gap ratio increased from 1.25 to 2.75. However, Cai [53] found that the erosion ratio increased as the gap ratio increased from 2 to 6; these results were similar to those of Zhang [68].
Tao [77] obtained results opposite to those of Cai [53] and Zhang [68] under constant flow rate boundary conditions since the size of pore channels enlarges with the gap ratio, resulting in two effects. First, the expansion of pore channels reduces the local flow velocity, which lowers the driving force and decreases the number of detached particles under constant flow rate conditions. Second, the expansion of the pore channels allows for easier migration of the detached particles through the constriction, thus enhancing particle transport and weakening clogging. The first effect dominates the erosion process at low gap ratios since the fine particles readily form arches, causing clogging and resulting in fewer detached particles. In contrast, the second effect dominates the erosion process at high gap ratios since few arches form, and clogging is less common. Therefore, there should be a critical gap ratio (G cr ) where the two effects are equivalent. It is meaningful to fully investigate the relationship between G cr and soil fabric and mechanical conditions.

Mechanical Conditions
Mechanical conditions, including the confining pressure and stress state, have substantial influences on the mechanical behavior of soil. To date, several studies were conducted on the influence of the mechanical conditions on contact anisotropy [78,79] and soil shear behavior [80,81], whereas few studies focused on soil erosion.

Effect of Confining Pressure
As the confining pressure increases, the contact force increases, which strengthens the friction between particles. Therefore, a strong driving force is needed to detach the forcechain particles from the soil skeleton, which significantly reduces the quantity of erodible soil particles. However, the increase in the confining pressure has no inhibitory effect on the detachment of gravity-deposited particles since these particles hardly transfer any effective stress [68]. Moreover, the pore channels are compacted at a high confining pressure, slowing the transport of detached particles. In addition, both numerical simulations [69] and experiments [82] have proved the existence of compression arches (here, the compression arch is part of the force chain, unlike the arch that forms clogging), which bear most of the driving force and form a protective zone behind them [69], as shown in Figure 8. However, negative effects also occur with an increase in confining pressure. The erosion of the force-chain particles can destroy the equilibrium state of the soil skeleton, leading to a large unbalanced force in the remainder of the force chain due to the loss of the large contact force, as shown in Figure 9. As a result, force-chain bending and local soil kinetic energy bursts may occur, reducing the soil stability. A suitable confining pressure can strengthen the friction, enhance the compression arch, and narrow the pore channels, improving the soil stability. However, excessive confining pressure may destroy the compression arches and soil skeleton [82] due to a relatively large unbalanced force, resulting in an unstable state.

21, 13, x FOR PEER REVIEW 9 of 18
leading to a large unbalanced force in the remainder of the force chain due to the loss of the large contact force, as shown in Figure 9. As a result, force-chain bending and local soil kinetic energy bursts may occur, reducing the soil stability. A suitable confining pressure can strengthen the friction, enhance the compression arch, and narrow the pore channels, improving the soil stability. However, excessive confining pressure may destroy the compression arches and soil skeleton [82] due to a relatively large unbalanced force, resulting in an unstable state.

Effect of Anisotropic Stress State
Numerous experiments [42,43,83,84] have proved that the erosion process varied significantly under different stress states, such as isotropy, triaxial compression, and the triaxial extension stress path. Since the soil is uniformly compressed in the isotropic stress state, the pore channels become narrower; thus, some constrictions are smaller than the diameter of the fine particles, blocking the transport path of the detached particles and weakening the erosion process. The anisotropic stress state is relatively complex. Since the force chains evolve towards the direction of the major principal stress in anisotropic stress

Effect of Anisotropic Stress State
Numerous experiments [42,43,83,84] have proved that the erosion process varied significantly under different stress states, such as isotropy, triaxial compression, and the triaxial extension stress path. Since the soil is uniformly compressed in the isotropic stress state, the pore channels become narrower; thus, some constrictions are smaller than the diameter of the fine particles, blocking the transport path of the detached particles and weakening the erosion process. The anisotropic stress state is relatively complex. Since the force chains evolve towards the direction of the major principal stress in anisotropic stress states [29,45,85], the pore distribution shows anisotropy, which also follows the direction of the major principal stress. Therefore, the transport path of the detached particles is relatively unobstructed, and the erosion develops rapidly when the seepage flows in the direction of major principal stress. In contrast, the seepage direction is perpendicular to the pore channel when the anisotropy is perpendicular to the major principal stress; thus, the detached particles are prone to be redeposited on the soil skeleton. Moreover, the transport path is relatively tortuous, and clogging occurs readily. It deserves attention that the anisotropic stress state may have a relatively larger impact on the pore channels in underfilled soils and on the force chains in filled and overfilled soils [76], emphasizing the importance of the soil fabric. Besides, in anisotropic stress states, the soil may dilate or contract, leading to the widening or narrowing of the pore channels, which affects the transport rate of the detached particles [29,45]. Foroutan [85] assessed the influence of the medium stress ratio b = (σ 1 − σ 2 )/(σ 1 − σ 3 ) on the macro-and micro-features of undrained granular samples. It was found that the effective mobilized friction angle reached the maximum when b was 0.5, indicating that the stress state has a significant influence on soil friction.
The influence of the mechanical conditions is relatively complex because it affects the friction and compression arches and causes contraction or dilatation of the soil. Thus, the mechanical conditions affect particle detachment, transport, and clogging by changing the force chains and pore channels. Additional research is required to understand this process comprehensively. For example, research should be conducted to (1) quantify the impact of friction on erosion, (2) explore the impact of the unbalanced force on compression arch failure due to the detachment of the force-chain particles, (3) investigate the geometry and mechanical boundaries of soil contraction and dilation, and (4) clarify the relationship between the soil fabric and the mechanical conditions.

The Process of Soil Structure Reorganization
Soil may lose stability, causing changes in the soil structure and mechanical properties when suffering external disturbance. The hydraulic conditions are external disturbance factors. Wautier [86] used the second-order work criterion to prove that the soil stability varied with the direction. When the stress increment exceeds the threshold value, which depends on the loading rate, the direction of loading, the soil properties (the geometric conditions and mechanical conditions), and the microstructural reorganization will be induced, i.e., soil instability is triggered. Therefore, the direction of hydraulic loading, the hydraulic gradient, and the hydraulic loading rate are crucial factors affecting soil stability. As the fluid flows through the soil, the additional driving force destroys several contacts, reaching the Mohr-Coulomb limit sliding condition, which indicates the occurrence of local particle detachment. Subsequently, the contact quantity and contact force decrease [34,62], leading to the continuous loss of lateral support and resulting in unjamming and bending of a few local force chains associated with a local burst of kinetic energy. This burst of kinetic energy propagates to the entire soil sample and provokes a generalized unjamming of the force chains [86], leading to microstructural reorganization until a new stable state is reached, as shown in Figure 10.

Three Effects of Hydraulic Loading Direction
The importance of the hydraulic loading direction, i.e., the seepage direction, is reflected in three aspects: (a) The soil stability is directional; (b) the seepage direction determines the effect of gravity; (c) the pore channel's direction affects the seepage action. Gravity has a pronounced influence on the number of detached particles since it is the driving force when the seepage is in the direction of gravity. In contrast, gravity is the resistant force when the seepage is opposite to it. Unlike stress disturbance which is transmitted through the force chain, the hydraulic disturbance is affected by the direction of the pore channels. When the pore channel direction is close to the direction of the force chain [42] or the same as the seepage direction, the seepage flow is unobstructed, and the migration of the detached particles is significantly enhanced [45]. On the contrary, when the pore channel direction is inconsistent with the seepage direction, the contact hindering the flow will cause continuous failure and restructuring under the impact of fluid flow. Therefore, as erosion continues, the contact direction will evolve toward the seepage direction [87,88]. Xiong [36] analyzed the effect of the seepage direction and found that the erosion ratio decreased rapidly as the angle of the seepage direction and gravity direction increased from 0° to 90°. It deserves attention that the directivity of the soil stability is likely related to the directivity of the force chains and pore channels, although the quantitative relationship is unclear. Further, the stability index of the soil in different directions may change during the erosion process since the erosion changes the pore channels and force chain.

Three Effects of Hydraulic Loading Direction
The importance of the hydraulic loading direction, i.e., the seepage direction, is reflected in three aspects: (a) The soil stability is directional; (b) the seepage direction determines the effect of gravity; (c) the pore channel's direction affects the seepage action. Gravity has a pronounced influence on the number of detached particles since it is the driving force when the seepage is in the direction of gravity. In contrast, gravity is the resistant force when the seepage is opposite to it. Unlike stress disturbance which is transmitted through the force chain, the hydraulic disturbance is affected by the direction of the pore channels. When the pore channel direction is close to the direction of the force chain [42] or the same as the seepage direction, the seepage flow is unobstructed, and the migration of the detached particles is significantly enhanced [45]. On the contrary, when the pore channel direction is inconsistent with the seepage direction, the contact hindering the flow will cause continuous failure and restructuring under the impact of fluid flow. Therefore, as erosion continues, the contact direction will evolve toward the seepage direction [87,88]. Xiong [36] analyzed the effect of the seepage direction and found that the erosion ratio decreased rapidly as the angle of the seepage direction and gravity direction increased from 0 • to 90 • . It deserves attention that the directivity of the soil stability is likely related to the directivity of the force chains and pore channels, although the quantitative relationship is unclear. Further, the stability index of the soil in different directions may change during the erosion process since the erosion changes the pore channels and force chain.

The Critical Hydraulic Gradient
When the hydraulic disturbance is insufficient to cause contact failure and forcechain bending, a small amount of erosion due to detached particles does not result in the reorganization of the microstructure. Therefore, similar to the stress disturbance, suffusion occurs when the hydraulic gradient exceeds the threshold value [62,[89][90][91], i.e., the critical hydraulic gradient. As the hydraulic gradient increases, the number of detached particles increases because more and more particles overcome the resistance force. In addition, with the increase in the hydraulic gradient, an increased number of force chains are destroyed and recombined, leading to periodic clogging and blowout of the clogging clusters [46]. The clogging-unclogging events lead to abrupt permeability and porosity changes, increasing soil instability [46,92] and releasing some of the fine particles blocking the pores, thereby increasing the transport distance of the detached particles. Therefore, the final erosion amount and erosion rate increase with an increase in the hydraulic gradient [89]. Moreover, Zhang [68] found that the erosion ratio increased roughly linearly with the hydraulic gradient. Since the erosion process changes the microstructure of soil, the mechanical properties of the soil will change after erosion [36,87]; this topic requires additional in-depth research.
The critical hydraulic gradient decreases with an increase in the hydraulic loading rate [34]. One interpretation is that the soil does not have sufficient time to restructure in response to the failure of the force chains since the rapidly increasing hydraulic loading rate will lead to a rapid increase in soil kinetic energy, which will be quickly transferred to the entire sample. In addition, experimental research showed that the hydraulic loading path and salt concentration had a significant and complex impact on erosion [93][94][95]. It is necessary to carry out in-depth research on the relevant hydraulic conditions to understand the impact of hydraulic conditions at the microscopic scale to deal with complex engineering applications.

Comprehensive Mechanism of Internal Erosion at the Microscopic Level
The previous analysis revealed the micro-mechanism of internal soil erosion, as depicted in the conceptual model in Figure 11. It was found that the geometric conditions and mechanical conditions determine the number and morphology of the pore channels and force chains in the soil, whereas the hydraulic conditions act as an external disturbance. The soil stability has directivity, which may result from the direction of the pore channels and force chains. Particle detachment can occur when the driving force overcomes the resistant force, indicating that the process is primarily affected by the hydraulic gradient and the force chains. The erosion is initiated when numerous particles are detached from the skeleton. Further, force-chain bending occurs due to lateral support loss, which can trigger a local kinetic energy burst in the soil. This effect is transferred to the entire sample, resulting in large-scale force-chain bending. Some detached particles are redeposited in the pores or become part of the force chain due to microstructural reorganization. Meanwhile, some of the remaining detached particles are transported through the constrictions and are removed from the soil sample; this process is influenced by the pore channels and seepage direction. The other detached particles are clogged near the constrictions due to their size or the arch formation. These particles will be transported after microstructural reorganization. The processes of redeposition, transport, and clogging determine the erosion amount. The contributory factors to internal erosion at the macroscopic level (see Figure 5) are considered at the microscopic level in the conceptual model in Figure 11.

Conclusions
Soils exhibit complex erosion characteristics due to the coupled effects of the geometric, mechanical, and hydraulic conditions. This article provides an in-depth review of the erosion mechanism at the microscopic level. The main findings and recommendations for future work are as follows: (1) Soil fabric is the primary factor affecting erosion. In the underfilled fabric, with an increase in the fines content, although the pore channels become narrower, more particles with low stress lead to a more intense erosion. In the filled and overfilled fabric, fine particles bear a higher stress than the underfilled fabric. However, a large unbalanced force induced by the detachment of force-chain particles results in a local kinetic energy burst in the soil. The combination of these effects causes the complex behavior of soils with filled and overfilled fabrics, which deserves careful investigation. (2) A critical gap ratio Gcr exists since the gap ratio impacts the local flow velocity and pore channel size, below which the erosion weakens when the gap ratio increases; otherwise, the erosion significantly intensifies. It is meaningful to fully research the relationship between Gcr and soil fabric and mechanical conditions.

Conclusions
Soils exhibit complex erosion characteristics due to the coupled effects of the geometric, mechanical, and hydraulic conditions. This article provides an in-depth review of the erosion mechanism at the microscopic level. The main findings and recommendations for future work are as follows: (1) Soil fabric is the primary factor affecting erosion. In the underfilled fabric, with an increase in the fines content, although the pore channels become narrower, more particles with low stress lead to a more intense erosion. In the filled and overfilled fabric, fine particles bear a higher stress than the underfilled fabric. However, a large unbalanced force induced by the detachment of force-chain particles results in a local kinetic energy burst in the soil. The combination of these effects causes the complex behavior of soils with filled and overfilled fabrics, which deserves careful investigation. (2) A critical gap ratio G cr exists since the gap ratio impacts the local flow velocity and pore channel size, below which the erosion weakens when the gap ratio increases; otherwise, the erosion significantly intensifies. It is meaningful to fully research the relationship between G cr and soil fabric and mechanical conditions.
(3) Like the response to an increase in the fines content, the pore channels become narrower, and the contact force increases with an increase in the confining pressure. However, a larger confining pressure, resulting in a relatively large unbalanced force, can destroy the compression arches in the soil. In addition, the anisotropic stress state causes directionality of the pore channels and force chains. Therefore, it is necessary to quantify the effects of soil anisotropy coupled with seepage direction on erosion characteristics. (4) The coupling of geometric, mechanical, and hydraulic conditions affects the pore channels and force chains, influencing particle detachment, transport, and clogging, which, in turn, determine the occurrence, initiation time, and the amount of erosion, which are hard to quantify from available studies. The following mechanical effects deserve further in-depth investigations: The impact of unbalanced force on compression arch failure, the boundary of soil contraction and dilation, the relationship between soil fabric and mechanical conditions. Besides, the critical hydraulic gradient depends on the hydraulic loading rate and direction, and the soil properties, which is not quantitatively clear.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. Intermediate principal stress σ 3 Minor principal stress