The Seepage-Destruction Mechanism of Water Inrush Channel of Sandstone Fault Filling Using the EDEM-Fluent Method

By the EDEM-Fluent coupling calculation method, the formation mechanism of the seepage failure water inrush channel and the migration of particles in the sandstone fault filling body are studied. Under the condition of variable hydraulic gradient, the whole seepage process can be divided into three stages: slow seepage stage, sudden seepage stage, and stable seepage stage. In the stage of slow seepage, the mass of lost particles is small. In the stage of sudden seepage, particles are lost on a large scale. In the stable seepage stage, the model is basically in a stable state. During the seepage process, the particles in the outlet zone will move before the particles in the inlet zone under the action of the seepage force. The overall movement trend of particles can be predicted, while the movement trajectory of a single particle is irregular. The change trend of the contact quantity between particles is basically consistent with the change in the quality of the lost particles. Moreover, the change in the contact quantity of particles is caused by the loss of the filler particles.


Introduction
During the construction of mountain tunnels, geological disasters of di erent degrees will be su ered and water inrush is common [1][2][3]. Its frequent occurrence and harmfulness bring great resistance to tunnel survey, design, and construction [4][5][6][7][8]. When crossing karst caves and fault fracture zones under water-rich conditions, karst-type water inrush usually occurs [9][10][11]. Besides, the ssure type water inrush occurs under the in uence of water storage structure and rainfall, and the water inrush pressure even exceeds 2 MPa [12][13][14]. Fault-type water inrush occurred during the excavation of two regional rock fault planes under severe weathering conditions, and the water ow was very huge, even reaching 1800 m 3 /h [15].
With the increasing demand for safety in the process of tunnel construction, a series of empirical standards and analysis models are proposed for the occurrence mechanism of water inrush, the prediction of water inrush volume, water inrush location, and water inrush characteristics. Guo et al. [16] used UDEC software based on the discrete element method, which was combined with laboratory tests, and analyzed the occurrence mechanism of water inrush in karst tunnels from macroperspectives and microperspectives. It was clari ed that the secondary stress redistribution of surrounding rock and karst water pressure with a high head easily caused water inrush in karst tunnels. Using the nite di erence program, Jiang et al. [17] calculated the instantaneous and stable water in ow during deep and long tunnel excavation in bedrock ssure aquifers and established a reliable and practical hydrogeological model. Based on the nite element method, Li [18] explored the seepage characteristics of groundwater and identi ed the location of inrush water as the vault and the arch waist on both sides and proposed the cement-water glass double slurry waterblocking measures. Zhao and Yang [19] analyzed the characteristics of water in ow when the tunnel passes through the fault fracture zone. When the tunnel collapses, the water in ow pressure will drop sharply and the existence of surface water will increase the water in ow and grouting thickness of the tunnel.
However, the traditional analysis method cannot be suitable for the complex fault fracture zone conditions, and the simulation of fluid movement is also relatively effective. e SPH method can be applied to solve complex fluid dynamics [20]. Bai et al. [21] proposed an SPH-FDM boundary method to analyze heat transfer in homogeneous media with discontinuous interfaces [22]. It should be noted that hydrogeological conditions are not the only determinant of water inrush type and seepage path [5,[23][24][25]. Bai et al. [26,27] established a nonlinear coupled transport equation of water movement driven by temperature in unsaturated soil. e water distribution under the action of the temperature gradient is uneven, which verified the hysteresis effect of water movement caused by the thermal cycle. A thermo-mechanical composition model is established for the particle rearrangement problem caused by the thermodynamics of porous granular materials, which can be used to describe the path problem in the process of water migration [28]. Bai et al. [29] explored the dynamic characteristics of suspended particles in porous media under variable temperature and proposed that water viscosity, adsorption effect, electric double layer effect, and particle kinetic energy will all affect particle motion.
e Zhaitang Mountain Tunnel of the G109 expressway has sandstone geological conditions. e fault fracture zone will be exposed during the construction process, the fault runs through the Zhaitang Tunnel and Zhaitang reservoir, and the shortest horizontal distance is only 320 meters. Based on the geotechnical dominant surface theory, it belongs to the dangerous water-conducting section [30,31]. With the tunnel excavation, the permeability coefficient of the fault fracture zone increases and the water source in the Zhaitang reservoir may carry the filling material in the fault to gush water into the tunnel along the fault fracture zone, thus forming a relatively through seepage channel, which is prone to seepage damage type water inrush of the filling body [9,[32][33][34][35]. Based on the DEM-CFD coupling calculation method, Zhou et al. [36] used a 100 mm × 50 mm × 100 mm rock mass model to simulate the engineering scale water inrush process with a real time of 1.8 s, showing the formation and expansion process of the dominant channel of water inrush. rough the selfdesigned test system, Wang et al. [37] made samples with a diameter of 100 mm and a height of 500 mm to qualitatively study the mechanism of water inrush caused by the migration of fault-filling particles. Xu et al. [38] used PFC3D software to build a fault fractured rock mass model with a diameter of 60 mm and a height of 120 mm and simulated the disaster process of water inrush in a microscopic view.
rough the laboratory test method, Chen and Lyu [39] made a 30 mm × 30 mm × 110 mm sample of the broken rock mass to discuss the change law and influencing factors of the particle migration on the hydraulic properties of the sandstone in the fault fracture zone during the occurrence of water inrush.
During the formation of seepage failure type water inrush channel in sandstone fault filling body, there is the interaction between particle phase and fluid phase. Not only the independent migration and stress of solid particles and fluid but also the exchange of energy, momentum, and mass between solid particles and fluid should be considered. Based on the fluid-solid coupling theory, this paper adopts the coupling method of discrete element software EDEM 2018 and computational fluid dynamics software Fluent 18.0, which is more accurate than the traditional method (i.e., PFC method), to establish a small-scale broken sandstone model and carry out three-dimensional numerical simulation. e formation mechanism of the water inrush channel and the migration law of particles in the seepage process are explored.

EDEM Geometric Model.
e fault fracture zone of tunnel excavation is a dangerous section for water diversion. Using EDEM software and the discrete element method, the model of broken sandstone inside the fault is established, which is expressed by several particle combinations. However, it is difficult to establish a stratum model which is equivalent to the actual engineering size in the process of numerical simulation due to the limitation of computer capabilities. erefore, combined with actual working conditions, a 40 mm × 40 mm × 100 mm rectangular broken sandstone millimeter scale model [40,41] with an initial porosity of 18.27% was studied to discuss the formation mechanism of seepage failure type water inrush channel of sandstone fault backfill under the action of reservoir seepage force.
As shown in Figure 1, the dark gray area is the geometric area where particles are generated, and the red border is the calculation domain. at is, if the particle movement exceeds the calculation domain, it will be deleted automatically. e highest water level of the Zhaitang reservoir is lower than the design elevation of the tunnel, so the flow in this model is seepage from bottom to top. e six faces of the broken sandstone model are all walls. e bottom and top faces are, respectively, set as fluid inlets and outlets. e physical properties of the surrounding walls are consistent with those of the sandstone, and none of them has water permeability. e top face of the model was deleted so that the particles inside the model can gush outward under the action of single   Advances in Civil Engineering bottom-up water seepage. Among them, the initial hydraulic gradient at the inlet is 0.1, increasing by 0.1 every 0.1 s in a stepped manner, and the fluid outlet pressure is constant at 0. Broken rock mass is considered to be composed of rock mass skeleton and filling medium. In combination with the actual conditions of the tunnel, two particles with different particle size ranges are used to represent the skeleton particles (2-15 mm) and filling particles (0.25-2 mm), respectively, to simulate the broken sandstone. e ratio of filling particles to skeleton particles is 1 : 3, all particle shapes are simplified to spheres, and they are rigid elements that cannot be deformed. EDEM software based on the discrete element principle is used to generate broken sandstone particles. e particle sizes of filling particles and skeleton particles are randomly generated within their desirable range, and a total of 18985 particles were generated. Figure 2 shows the broken sandstone model generated by EDEM with a gravity of 9.8 KN/m 3 and a negative direction along the Z axis.
It should be noted that fault fractured rock mass is formed by bonding a series of discontinuous particles, and different particle contact models will produce different contact forces. Based on this, the contact between skeleton particles is set as Hertz-Mindlin with the bonding model to simulate the adhesion between particles. Skeleton particles with certain adhesion constitute a skeleton network. Once the adhesion between particles is destroyed, the particles will change to Hertz-Mindlin (No Slip) model when they contact again. According to field measurement and theoretical analysis, the physical parameters of particles are shown in Table 1.
Due to the difference between the mesh generation in EDEM software and the mesh size of fluid computing, when using EDEM-Fluent for coupling calculation, it is necessary to divide the mesh of EDEM and fluid computing, respectively. e mesh size of EDEM is directly related to the minimum particle size, which is usually about twice the minimum particle radius. e selected EDEM grid size is 2.5 times of the minimum particle radius, and the minimum radius of generated particles is 0.125 mm. Hence, the grid size is 0.3125 mm, including a total of 5611647 grids. To ensure the grid quality and the accuracy of calculation results, a single fluid calculation grid must contain at least three solid particles, and the size of the fluent fluid grid is 10 mm × 10 mm × 10 mm, including a total of 160 fluid computing grid cells.

EDEM Calculation Principle.
In the discrete element method, the particle motion follows Newton's second law of motion [42], and the interaction between particles is a transient process. Within a sufficiently small time step, it can be assumed that the velocity and acceleration of a single particle are fixed values and only the particles or boundaries in direct contact with it play a role without being affected by the indirect action of other particles. e basic equation of motion for a single particle is as follows: where m is the particle mass, u is the linear displacement, θ is the angular displacement, f d , f c , an d f are the damping force contact force and volume force on the particles, respectively, I is the moment of inertia, and r c is the force arm of the contact force on the particles. During EDEM calculation, the particle and boundary position will be obtained. Taking the relative displacement between particles as the basic physical variable, the normal force and tangential force of particles will be obtained according to the relationship between force and displacement in the contact model, so as to obtain the resultant force and resultant moment. en, the physical quantities such as displacement, velocity, and acceleration of particles will be calculated according to Newton's second law of motion, and the continuous cycle calculation will be carried out. Fluent calculation ensures that the fluid motion follows the continuity equation and momentum conservation equation.

Particle Movement.
e total duration of seepage failure of the simulated sandstone fault filling body is 1 s, the time step of EDEM is set as 1 × 10 − 7 s, and the data saving interval   Figure 3 shows the particle state in the model when t � 0, t � 0.5 and t � 1.0 s in the seepage process [22]. e dark green particles are filling particles and the yellow-green particles are skeleton particles. Different time points correspond to different hydraulic gradients [27]. By observing the particle state at different time points, based on the loss speed and amount of particles, the movement law of particles under the action of fluid seepage force with gradient increase can be obtained.
When the seepage force is small, the filling particles in the outlet area along the direction of seepage force and some skeleton particles outside the skeleton network will gush out before the particles in the inlet area [43], and the particle position in the inlet area has no obvious change and is in a stable state. With the increase of hydraulic gradient, the resultant force of fluid drag force, gravity, and indirect contact force of particles in the middle section occurs gradually upward and the contact relationship between particles also changes. e flow will carry a large number of filling particles and some small skeleton particles from the high water pressure to the low water pressure area. At this stage, the filling particles gush out on a large scale, some small skeleton particles will also be lost, and the large skeleton particles will move unsteadily with the change of contact relationship. As most of the filling particles have gushed out, the remaining filling particles will gradually stabilize, but there are still a small part of filling particles and small skeleton particles slowly gushing out. Due to the change of the contact relationship between particles, the particles moving from the middle and lower part to the outlet area under the action of seepage force will also gush out. e filling particles in the inlet area cannot gush out because of the large skeleton particles and form a stable structure with the remaining skeleton particles.
When t � 1 s, the seepage simulation process ends, and most of the filling particles and a small part of the skeleton particles have gushed out, which further increases the porosity, permeability coefficient, and water velocity of the fault fracture zone, forming a relatively smooth water inrush channel. In addition, the loss of particles is obviously nonuniform. Although the particles migrate first in the outlet area, the particles do not gush out layer by layer, indicating that the water velocity in the particle pores is different in the same horizontal plane, and the seepage path is irregular [44]. However, under the action of the drag force of the water flow, it drives the filling particles and some skeleton particles to gush out, thus forming a relatively smooth water inrush channel.
According to the motion law of a single particle, the motion trajectories of three particles in the outlet area, the middle section, and the inlet area are randomly selected (Figure 4). In the process of gushing, a single particle moves under the combined action of gravity, water drag force, normal contact force, and tangential contact force between different particles. e upward migration and gushing of filling particles and some small skeleton particles along the direction of seepage force need to pass through the pores between skeleton networks. e particle sizes of single particles randomly selected may be different, and the gravity of particles with different sizes is different [45,46]. Because the pore size and the flow velocity and direction are different between the skeleton particles, the size and direction of the flow drag force on different particles are different. e contact relationship between a single particle and other particles is not exactly the same [47][48][49], and the sizes and directions of normal, tangential, and angular velocities of particles affected by water flow are also different. Similarly, there are differences in the moment of inertia of particles with different particle sizes [50]. erefore, the sizes and directions of normal and tangential contact forces on a single particle at different times are different, which leads to different resultant forces of a single particle under the action of water flow.
Actually, the movement of particles is irregular and will not strictly follow the action direction of seepage force. For example, the single particle movement tracks in the middle section selected in Figure 4 move downward first and then upward. Figure 5 more intuitively shows the trajectory of the particle, in which the purple particle is the characteristic particle, the dark green particle is the filling particle, and the yellow-green particle is the skeleton particle. e motion track of the characteristic particle starts because the particle is a skeleton particle free from the skeleton network. Under the action of initial water pressure, the force of seepage force on the skeleton particle is not enough to make it move upward and gush out, but the drag force of water flow can carry the filling particles around the particle upward, thus changing the contact relationship between the particle and other filling particles. With the increase of seepage force, the surrounding filling particles continuously gush out and the contact force between particles decreases. e gravity of the characteristic particles gradually plays a leading role in the resultant force, and the resultant force direction of the particles is gradually downward, so the particles move downward first. As the seepage force continues to increase, the drag force of the water flow on the characteristic particles continues to increase and the resultant force on the particles gradually moves upward. e drag force begins to dominate To clarify the influence of the particle indirect contact force on the motion trajectory of the characteristic particle, the change of the particle indirect contact force of particle 1 in the whole seepage process is output in the EDEM postprocessing ( Figure 6). Within 0-0.2 s, the filling particles begin to move under the action of the drag force of water flow, and the contact relationship between the characteristic particles and the surrounding filling particles changes, so the contact force on the particles begins to decrease gradually. At this time, particle 1 does not move, and at the same time, the pores around the particle increase. After 0.2 s, gravity plays a leading role, and the characteristic particles begin to move downward, but the drag force of the water flow on the particles increases with the increase of the hydraulic gradient, and the particles begin to move upward until they gush out of the channel at 0.6 s. During the movement of the particle, it is subjected to a short and slight contact force at some time because the particle collides with other particles during the movement. e trajectory of a single particle is also different in the middle section shown (Figure 4). In the process of particle emission, some particles do not move upward along the direction of seepage force but remain stationary after a period of upward movement until the end of seepage. e reason for the particle trajectory is that at the beginning, the filling particles move upward under the action of the water drag force, get stuck between the skeleton network formed by the large skeleton particles after a certain distance, and then reach a stable state. Because the characteristic particle is located inside the model and it is difficult to observe directly, other particles are hidden and only the characteristic particle is displayed (Figure 7). e characteristic particle moves upward slowly under the initial water pressure, and it moves violently at 0.8 s and finally stops the upward trend.

Mass Change of Lost
Particles. To explore the particle emission during the seepage process, using the postprocessing function of EDEM software, the change of residual particle mass with time is shown in Figure 8, which can directly reflect the particle loss speed and mutation point at different time points and the particle emission in a certain period of time. Figure 9 shows the variation of lost particle mass with time during the whole seepage process. It can be seen that within 0-0.3 s, the slope of the curve is relatively flat, and the particle loss speed is relatively slow. At 0.3 s, the slope of the curve increases sharply, the particles begin to lose on a large scale and the particle loss speed and amount increase. After 0.8 s, the curve slope tends to be flat again, the particle loss speed gradually decreases, and the particle loss quality tends to be stable. e particle loss process can be divided into three stages. In the first stage (slow seepage stage, 0-0.3 s), the particle loss mass in this stage is 5.0 g, accounting for 5.41% of the total mass loss in the whole seepage process. e initial hydraulic gradient is 0.1, the drag capacity of the flow is low, the filling particles in the outlet area will move upward slowly, and the skeleton particles will not move because of their own gravity and the large contact force between the particles. e loss mass of particles in 0-0.1 s is 2.8 g, accounting for 56% of the total mass of particles lost in the slow seepage stage. is is because the upward movement of particles in the outlet area is not hindered by other particles. After this part of the particles gushes out, the rest of the filling particles and some of the small skeleton particles free from the skeleton network begin to move upward. However, due to the drag force of water flow, the outflow velocity and outflow volume of particles are relatively small. In practical engineering, the seepage volume is small, and the water flow is relatively clear.
is stage is the best time to take corresponding measures. Grouting and water plugging at appropriate positions can avoid the formation of water inrush channels. In the second stage (sudden seepage stage, 0.3 s-0.8 s), the critical hydraulic gradient corresponding to the starting point is 0.4 and the total mass loss of particles is 85.3 g, accounting for 92.32% of the total mass loss in the whole seepage process. With the emission of particles gathered in the outlet area in the slow seepage stage, the pores between the skeleton particles gradually increase and the particle loss speed and amount suddenly change, forming a relatively through seepage channel. is stage is the main stage of particle loss in the whole seepage process. In fact, the seepage volume increases and the water quality is turbid. In the third stage (stable seepage stage, 0.8-1 s), the total mass loss of particles in this stage is 2.1 g, accounting for 2.27% of the total mass loss in the whole seepage stage. As most of the filling particles and some of the skeleton particles that are free from the skeleton network have gushed out, only a few particles will slowly gush out after 0.8 s. e remaining filling particles interact with the large skeleton particles to form a new skeleton. e seepage reaches a stable state and forms a smooth seepage channel.

Change of Particle Contact Quantity.
In the process of seepage, particles move upward and gush out, because the contact relationship between particles is changing and the Advances in Civil Engineering normal contact force and tangential contact force are changing. e contact number between particles is described by the change of coordination number in EDEM with time. e coordination number is greater than the contact number between particles, but the larger the coordination number, the more the contact number between particles. As shown in Figure 10, the particle coordination number versus the time curve is shown. e red curve represents the change of total coordination number between particles, and the blue curve represents the change in coordination number of skeleton particles.
Within 0-0.3 s, the coordination number of skeleton particles did not change much, while the total coordination number of particles decreased slightly to a certain extent. After 0.3 s, the coordination number of skeleton particles began to decrease because the continuous loss of filling particles led to the change of contact force on skeleton particles under the action of seepage force, which made skeleton particles free from the skeleton network move, However, the contact relationship between the skeleton particles in the skeleton network has not changed. e change of contact quantity between particles is consistent with the seepage stage, and the decrease rate of contact quantity between particles is relatively slow in the process of slow seepage. In the stage of sudden seepage, the contact quantity between filling particles and between filling particles and skeleton particles changes abruptly. During the steady seepage stage, the contact quantity between particles reaches a steady state. After the seepage process, the particles in the model are basically in a stable state, and only a small number of particles will continue to lose. erefore, the contact number between particles will not change significantly.

Change of the Porosity.
e increasing porosity in the process of seepage is the fundamental reason for the formation of water inrush channels. Because the particle loss in the inlet area is smaller than the total loss, the porosity changes in the inlet area and the total are studied, respectively. By adding "Geometry Bin" in EDEM postprocessing, the change of porosity in the inlet area is analyzed, where the size is set to 40 mm × 40 mm × 50 mm. e evolution curve of porosity with time cannot be directly generated by EDEM. erefore, the change of porosity is calculated through the change of particle volume or mass in the region. e calculation is shown as follows: where P is the porosity, V is the total volume of the model, m s is the particle mass in the model,ρ 0 is the bulk density of the material, and ρ s is the material density. Figure 11 (red curve) shows that the porosity in the entrance region decreases first and then increases rapidly and finally becomes stable with time. In the stage of slow seepage, the hydraulic gradient is relatively small. For larger particles in the outlet area, gravity plays a leading role. Under the downward condition of the resultant force, it moves slowly downward and continuously flows into the "Geometry Bin." erefore, in this stage of seepage, the porosity in the inlet area decreases. As the hydraulic gradient reaches the critical hydraulic gradient and enters the sudden seepage stage, the filling particles in the inlet area will move upward rapidly under the drag force of water flow, so the porosity will increase rapidly. After entering the stable seepage stage, there is still a small amount of filling particles moving to the outlet area, but most of the particles have been in a relatively stable state, so the porosity will slowly increase and tend to be stable.
For the overall porosity, the porosity versus time curve (Figure 11, blue curve) can be obtained according to Figure 8. e change of overall porosity is a direct manifestation of particle loss. e time change curve of overall porosity is consistent with the time change curve of lost particle mass in Figure 9. In the slow seepage stage, because a small amount of particles gush out, the overall porosity increases slowly. In the stage of sudden seepage, the hydraulic gradient reaches the critical hydraulic gradient and the overall porosity    Advances in Civil Engineering increases rapidly with the large amount of particles gushing out. After reaching the stable seepage stage, as the particles have basically remained stable, the overall porosity of the model also tends to be stable.

Conclusions
In the process of seepage, the particles will migrate under the action of the flow seepage force and the particles in the outlet area will move before the particles in the inlet area. As the particles in the outlet area gush out, the particles in the inlet area will move and gush out, forming a relatively smooth water inrush channel. Under the condition of variable hydraulic gradient, based on the loss speed and quality of particles, the whole seepage process can be divided into three stages: slow seepage stage, sudden seepage stage, and stable seepage stage. ere are critical hydraulic gradients in the slow seepage stage and sudden seepage stage. e particle loss speed is slow under the condition of the low hydraulic gradient. After reaching the critical hydraulic gradient, the particle loss speed increases significantly. After most of the filling particles gush out, the particle loss speed slows down and tends to be stable. In the process of water inrush, the overall movement trend of particles can be roughly predicted, but the movement trajectory of a single particle is irregular. e change of contact quantity is mainly caused by the loss of filling particles. e change trend of contact quantity between particles is basically consistent with the change in particle loss quality. In the stage of slow seepage, the contact quantity between particles decreases slowly. e contact quantity between particles decreases rapidly after entering the stage of sudden seepage. In the steady seepage stage, the contact quantity between particles is basically in a stable state.
Data Availability e original contributions presented in the study are included in the article, and further inquiries can be directed to the corresponding author.

Authors' Contributions
e main contribution of SS in this paper is methodology, and the other authors (WC, YW, FH, and PZ) contributed to investigation and analysis.