A Theory of Slope Shear Scouring and the Failure Mechanism of PFC3D on a Gangue Slope

: In this paper, scouring shear failure theory is optimized based on the gangue slope near the thermal power ﬁeld in Baiguo Town, Panzhou City, Guizhou Province. Based on the particle ﬂow PFC (particle ﬂow code) 3D ﬂuid–solid coupling method, the scouring failure mechanism of ditch no. 5 of the gangue slope is comprehensively analyzed from the perspectives of the failure mode, displacement, motion track, and stress–strain. We consider the scouring shear theory in respect of (c, ϕ ); this theory is dominated by two types of scouring intensity factors and can effectively explain the internal mechanisms of gully formation. The rainfall scouring failure of gangue slopes can be divided into four stages: (1) integral splash erosion and local pitting at the bottom of the slope; (2) erosion diversion and pitting in the slope; (3) the tributary–slope crest extension schist erosion stage; and (4) integral gully erosion and the landslide stage. The failure process is not only characterized by discontinuous failure but also occurs in the order of bottom–middle–branch–top. A three-section stepped effect is observed during the process in which the gangue is scoured and destroyed, which fully veriﬁes the intermittent characteristics of the scouring and destruction of gangue slopes. During the whole process, the maximum displacement is concentrated at the top of the slope, and its proportions are as follows: top of the slope > tributary > middle of the slope > foot of the slope. The peak displacement of the slope crest in the horizontal Y-direction accounts for 41.76%, and that in the Z-direction accounts for 45.84%. Scouring deposits can be divided into the arc erosion deposit mode and the sector erosion deposit mode. Mainstream gullies primarily control whether deposits are characterized as arc or straight erosion deposits. The later stage of the second phase of scouring is the incubation period of the tributary gully. The large accumulation makes the stress at the bottom of the slope increase sharply, and the ﬂuctuation value is between 2 and 6.8 MPa. The generation of the branch notch is mainly determined by X-direction stress, and 8.6 MPa is the critical stress. In efforts to prevent and control rainfall and landslide, the slope foot area should be preferentially protected, and the soil mass in the slope should be reduced to prevent the maximum energy ﬂuctuation caused by scouring, so as to prevent signiﬁcant displacement damage of the slope top.


Introduction
Gangue is the solid waste discharged during coal mining and washing. It refers to a black-gray rock with a low carbon content and a higher level of textural hardness than the coal associated with the coal seam during the coal-forming process [1]. There are now about 7 billion tons of gangue in China, and this value continues to increase year by year [2,3]. Gangue dumps discharged from mines are essentially formed by piling on the ground. There are more than 1600 gangue dumps formed by the long-term accumulation of gangue, and Guizhou accounts for about 65% of the country's total output. The accumulation of gangue not only takes up a large amount of land resources but it also forms many gangue Since scouring erosion mainly occurs during heavy rainstorms, the splash erosion caused by raindrops at the beginning of rainstorms is not included in this paper, only the scouring and penetration of the water flow in gully no. 5 on the slope surface are considered. The scouring and penetration model of gangue slopes is shown in Figure 1. In deriving the scouring theory, some scholars [37] consider the influence of cohesion c, while others [38] consider the influence of the internal friction angle, but the two factors have not been discussed further. In this paper, the two methods are combined to optimize the shear stress of particles attached to the flow, taking into account the change in the instantaneous velocity of particles in PFC particle flow. It can be seen from the slope scouring phenomenon of viscous gravel soil that during a rainstorm, the flow pattern is turbulent. The flow is affected by soil particles on the slope and forms many vortices of different sizes. Their position, shape, and flow rate are constantly changing as they advance, and the shear stress of flow on the slope surface [39] is expressed as follows: It can be seen from the slope scouring phenomenon of viscous gravel soil that during a rainstorm, the flow pattern is turbulent. The flow is affected by soil particles on the slope and forms many vortices of different sizes. Their position, shape, and flow rate are constantly changing as they advance, and the shear stress of flow on the slope surface [39] is expressed as follows: τ w = τ wa +τ wb (1) τ wa = −ρv x v y (2) τ wa is the viscous shear stress in the laminar flow. τ wb is the turbulent shear stress. ρ is the density of the water flow, and v x and v y are fluctuating velocities in the X-and Y-directions, respectively, as shown in Figure 1. Defined by the pulsation velocity, it can be seen that: where v x0 and v y0 are instantaneous velocities of soil particles and v x and VY are the average velocities of water flow. Due to the fast flow rate of scouring, according to the relationship between the flow rate and seepage rate, the pulsating seepage speed in the Y-direction can be expressed as: where ξ is the effective porosity between particles. From Planter's hypothesis and the Nikuradse hydraulic roughness test, it can be concluded that: where R is the hydraulic radius, and ∆x is the spacing between Section 1 and Section 2 on the slope. Formulas (2)-(5) can be substituted for Formula (1) to obtain: For rock and soil slopes, the scouring capacity of the slope surface should account for the scouring strength of the water flow. At the same time, the influence of the internal friction angle of the soil particles should be given priority in the analysis. Considering the internal friction angle, the bed washout rate g p can be expressed as: q s = h c n R 2 3 J 1 2 (8) τ c = γh c J (9) where ρ c is the density of soil particles on the slope surface and q s is the flow rate per unit width; the Manning coefficient n = 0.025 is used. The hydraulic radius R = h c ; J = sin 45 • = 0.706. τ c is the critical pulling force. γ is the weight of water and is usually taken as 1 × 104 N/m 3 . h c is the slope runoff depth. ϕ is the internal friction angle of the soil and θ is the gradient.
According to the continuous equation of soil particles on the slope, it can be seen that: ∂g p ∂x + ρ c g ∂y ∂t = 0 (10) where ρ c is the dry density of soil particles: We differentiate Equation (10) to obtain the difference equation: It is assumed that the scouring depth obtained by the above method is ∆y 1 , considering the ϕ factor. Furthermore, the total scouring depth ∆y at ∆x spacing is: ∆y = ∆y 1 + ∆y 2 (13) where ∆y 2 is the scouring depth considering cohesion factor c. According to the progressive scouring strength of the entry-exit section, solution (12) can be obtained as follows: We substitute Formulas (7)- (9) into Formula (14) to obtain: (15) ω = 0.015q s J 0.5 τ c 0.6tan ϕ cos θ (16) In this equation, ω is the scouring coefficient, which is related to the water depth of the slope's surface, the inclination angle of the slope body, and the friction angle of the soil body. τ w1 and τ w2 are the flow shear stresses of Section 1 and Section 2, respectively. Equation (15) is simplified to: where k 1 is defined as a type I scouring intensity factor, which represents the thickness of the soil layer eroded by the water flow per unit of time, and is related to the internal friction angle ϕ and the slope position ∆x. The c-value cannot be neglected when considering certain highly cohesive soils. In this situation, the scouring stress of runoff on soil particles on the slope surface is the sum of the impact-induced stress and flow shear stress, which is expressed as follows: k 2 is defined as a type II scouring intensity factor. It is defined as the thickness of the soil layer eroded by the water flow per unit of time, taking into account the cohesion c between soil particles.
Appl. Sci. 2023, 13, 5066 6 of 27 In this equation, c is the cohesion of the soil, σ ω is the impact stress of water flow on soil particles, and ω is the erosion coefficient, which is the rate at which the cohesion of soil particles decreases to the critical stress of scouring under the action of water flow penetration per unit of time; it is a function related to soil characteristics that reflect the sensitivity of soil scouring erosion. By combining the scouring test with the fitting of the first pot global optimization algorithm, the following conclusions can be drawn: ω = 0.01cKB 2.3 sin θ µ + 0.759 (21) where B is the compactness of the slope. The permeability coefficient K is: Then, we consider: By substituting Equation (23) into Equation (20), we obtain: Finally, according to the above formula, an expression for the scouring depth at any position of soil mass with a time variation can be obtained as follows:

Case Data
In this paper, we take the gangue slope near the Panzhou Electric Field in Guizhou Province, China, as an example. According to on-site reconnaissance, the upper lithological layer in the study area is black gangue with a loosely structured Quaternary Holocene artificial accumulation layer (Q4ml). The lower layer is grey-white limestone of the Lower Triassic Yongning Town Formation (T1yn) with an occurrence of 70-75 degrees and 40-45 degrees. Coal gangue piles in the tailings dam reservoir area are mainly located in the left tail slope area of the dam area, with a total area of about 670,000 m 3 (as shown in Figure 2). The top elevation of the coal shale slope is about 1710 m, the height is about 40-60 m, and its slope is 25-40 degrees. It belongs to an unstable slope. It constitutes a major safety hazard to the ash yard and residential buildings under the hill. It can be seen that the lithology of the slope mainly comprises the manual accumulation of gangue with a very loose structure, especially at the corner of the slope (the red area in Figure 2). Under the scouring of rainwater, rainwater collects along the slope on both sides at the natural gully of the corner, forming the largest gully (no. 5) with a width of 3-10 m, a length of about 90 m, and a depth of up to 7 m (see Figure 3). The main focus of this study is the area near gully no. 5 at the corner of the accumulation body. safety hazard to the ash yard and residential buildings under the hill. It can be seen that the lithology of the slope mainly comprises the manual accumulation of gangue with a very loose structure, especially at the corner of the slope (the red area in Figure 2). Under the scouring of rainwater, rainwater collects along the slope on both sides at the natural gully of the corner, forming the largest gully (no. 5) with a width of 3-10 m, a length of about 90 m, and a depth of up to 7 m (see Figure 3). The main focus of this study is the area near gully no. 5 at the corner of the accumulation body.

PFC Principle
The particle flow code (PFC) method is based on the mechanics of discontinuous media and is a mesoscopic discrete element method. This method defines the material as a particle assembly composed of finite particles, which are rigid bodies with a three-dimensional mass and spherical particles. Fine particles and different contact points on each particle have the characteristics of discontinuity. The particle flow method is mainly applied to scientific research on basic problems such as the basic characteristics of geotechnical materials and the material dynamic response. In the PFC3D numerical simulation, the constitutive relationship is the contact force-displacement relationship, the motion equation follows Newton's second law, and the contact between particles is determined by the given contact model. In this paper, the contact model is a parallel bond model. safety hazard to the ash yard and residential buildings under the hill. It can be seen that the lithology of the slope mainly comprises the manual accumulation of gangue with a very loose structure, especially at the corner of the slope (the red area in Figure 2). Under the scouring of rainwater, rainwater collects along the slope on both sides at the natural gully of the corner, forming the largest gully (no. 5) with a width of 3-10 m, a length of about 90 m, and a depth of up to 7 m (see Figure 3). The main focus of this study is the area near gully no. 5 at the corner of the accumulation body.

PFC Principle
The particle flow code (PFC) method is based on the mechanics of discontinuous media and is a mesoscopic discrete element method. This method defines the material as a particle assembly composed of finite particles, which are rigid bodies with a three-dimensional mass and spherical particles. Fine particles and different contact points on each particle have the characteristics of discontinuity. The particle flow method is mainly applied to scientific research on basic problems such as the basic characteristics of geotechnical materials and the material dynamic response. In the PFC3D numerical simulation, the constitutive relationship is the contact force-displacement relationship, the motion equation follows Newton's second law, and the contact between particles is determined by the given contact model. In this paper, the contact model is a parallel bond model.

PFC Principle
The particle flow code (PFC) method is based on the mechanics of discontinuous media and is a mesoscopic discrete element method. This method defines the material as a particle assembly composed of finite particles, which are rigid bodies with a threedimensional mass and spherical particles. Fine particles and different contact points on each particle have the characteristics of discontinuity. The particle flow method is mainly applied to scientific research on basic problems such as the basic characteristics of geotechnical materials and the material dynamic response. In the PFC3D numerical simulation, the constitutive relationship is the contact force-displacement relationship, the motion equation follows Newton's second law, and the contact between particles is determined by the given contact model. In this paper, the contact model is a parallel bond model.

Simplified Model of Gangue Clump
The literature [40] shows that the appearance of gangue blocks is extremely irregular, with certain elongations and flat areas, which leads to an obvious biting action and stress concentration between gangue blocks. In order to simulate this characteristic of gangue as far as possible, and to improve the calculation efficiency, we use the basic unit "Ball" in PFC3D software to simulate particles in bulk materials such as crushed gangue, coarsegrained soil, and non-sticky sand. However, a ball is a rigid body with a rounded surface, which cannot simulate the biting action, stress concentration, and fracture between particles. Therefore, based on the clump command in PFC3D and two real shapes of a gangue block (block and sheet), a three-dimensional numerical model is simulated. First, it is imported into PFC3D using the geometry import command, which provides an enclosed restriction space for the generation of a gangue block in PFC3D. Then, the Pebble unit is used to fill the target restricted-space area. During the filling process, the Pebble continuously adjusts its radius and position in the filling area, according to the BubblePack algorithm, until the error between the filled volume of Pebbles and the target restricted-space volume is less than the set threshold, stopping filling and completing the generation of the clump aggregate. The final generated clump templates of the two shape categories of waste rock blocks and their shapes' characteristic parameters are shown in Figures 4 and 5. "Ratio" and "distance" are the two key words used in the BubblePack algorithm to control the effect of clump generation. "Ratio" controls the ratio of the Pebble's minimum radius to the maximum radius used to generate the clump template, ranging from 0 to 1.0; "distance" controls the smoothness of the clump template generation, ranging from 0 to 180. The larger the value, the smoother the surface of the clump generated. When the value is 0, only point contact, with no overlap, is allowed between the PBBs used to generate the clump. When the value is 160, the PBBs used to generate the clump are allowed to overlap completely. From Figures 4 and 5, it can be seen that when ratio = 0.5 and distance = 160, the clump template generated by the BubblePack algorithm not only restores the original appearance of the gangue block, but it also takes into account the basic external shape and calculation efficiency of the gangue block; therefore, this method was adopted.
which cannot simulate the biting action, stress concentration, and fracture between particles. Therefore, based on the clump command in PFC3D and two real shapes of a gangue block (block and sheet), a three-dimensional numerical model is simulated. First, it is imported into PFC3D using the geometry import command, which provides an enclosed restriction space for the generation of a gangue block in PFC3D. Then, the Pebble unit is used to fill the target restricted-space area. During the filling process, the Pebble continuously adjusts its radius and position in the filling area, according to the BubblePack algorithm, until the error between the filled volume of Pebbles and the target restricted-space volume is less than the set threshold, stopping filling and completing the generation of the clump aggregate. The final generated clump templates of the two shape categories of waste rock blocks and their shapes' characteristic parameters are shown in Figures 4 and 5. "Ratio" and "distance" are the two key words used in the BubblePack algorithm to control the effect of clump generation. "Ratio" controls the ratio of the Pebble's minimum radius to the maximum radius used to generate the clump template, ranging from 0 to 1.0; "distance" controls the smoothness of the clump template generation, ranging from 0 to 180. The larger the value, the smoother the surface of the clump generated. When the value is 0, only point contact, with no overlap, is allowed between the PBBs used to generate the clump. When the value is 160, the PBBs used to generate the clump are allowed to overlap completely. From Figures 4 and 5, it can be seen that when ratio = 0.5 and distance = 160, the clump template generated by the BubblePack algorithm not only restores the original appearance of the gangue block, but it also takes into account the basic external shape and calculation efficiency of the gangue block; therefore, this method was adopted.

Comparative Analysis of Triaxial Tests
In this study, a triaxial test was used to reasonably adjust the mesoscopic parameters of the PFC model. Rock samples were taken from the coal gangue slope near the Panzhou Electric Field in Guizhou Province, China, and an analysis was carried out using an SJ-IA mechanical testing system. The coal gangue sample was cut and polished to ensure that the test piece met the required accuracy. The test piece was a cylinder with a diameter of 50 mm and a height of 100 mm, and the error range was controlled within ±2 mm, as shown on the left of Figure 6.
For the simulation study of particle flow, the key condition to determine the rationality of the model is the selection of the model's meso-parameters. The micro-parameters of a material cannot be directly assigned by real macro-mechanical parameters, so they need to be adjusted using an equivalent comparative test to ensure that the material's behavior is close to the real effect. In this study, a large number of repeated triaxial tests were carried out by continuously changing the micro-parameter settings. At the same time, the numerical model of gangue specimens was established based on the servo principle of the PFC3D triaxial test; additionally, the stress-strain comparison analysis and micro-parameter calibration of the gangue were carried out. Finally, the appropriate micro-parameters were determined, as can be seen on the right of Figure 6. The model generates particles in the cylinder wall with a diameter of 5 m and a height of 10 m according to real triaxial testing. Considering the influence of the particle size effect on the test results [41], for a ratio of model size to average particle size >30, the appropriate change in the particle size does not affect the simulation results but it will improve the calculation efficiency. The particle size distribution of the PFC is shown in Table 1. The model established in this paper was repeatedly checked and adjusted. The average particle size is 1.9 m, the ratio of the minimum slope model size (66 m) to 1.9 m is 37, and the ratio of the maximum model size (160 m) to 1.8 m is 89. These ratios are much larger than the standard value of 30, so the particle size values in this model are reasonable. The three-axle confining pressure was loaded according to the values of 100 KPa, 200 KPa, and 300 KPa from the laboratory test.

Comparative Analysis of Triaxial Tests
In this study, a triaxial test was used to reasonably adjust the mesoscopic parameters of the PFC model. Rock samples were taken from the coal gangue slope near the Panzhou Electric Field in Guizhou Province, China, and an analysis was carried out using an SJ-IA mechanical testing system. The coal gangue sample was cut and polished to ensure that the test piece met the required accuracy. The test piece was a cylinder with a diameter of 50 mm and a height of 100 mm, and the error range was controlled within ±2 mm, as shown on the left of Figure 6.

Comparative Analysis of Triaxial Tests
In this study, a triaxial test was used to reasonably adjust the mesoscopic parameters of the PFC model. Rock samples were taken from the coal gangue slope near the Panzhou Electric Field in Guizhou Province, China, and an analysis was carried out using an SJ-IA mechanical testing system. The coal gangue sample was cut and polished to ensure that the test piece met the required accuracy. The test piece was a cylinder with a diameter of 50 mm and a height of 100 mm, and the error range was controlled within ±2 mm, as shown on the left of Figure 6. For the simulation study of particle flow, the key condition to determine the rationality of the model is the selection of the model's meso-parameters. The micro-parameters of a material cannot be directly assigned by real macro-mechanical parameters, so they need to be adjusted using an equivalent comparative test to ensure that the material's behavior is close to the real effect. In this study, a large number of repeated triaxial tests

Calibration of the Meso-Parameters of Gangue
The trend in the stress-strain curves is the key to determining the relationship between the meso-parameters and macro-parameters; this is also the micro-parameter calibration method adopted by most scholars. On this basis, the stress-strain curves were obtained, as shown in Figure 7. Compared with the indoor triaxial test, the trends of the two stressstrain curves obtained from the numerical triaxial test under the same ambient pressure are essentially the same as those of the indoor test, and the peak deviation stress error of the gangue is largely maintained within the range 0~1 KPa. maximum model size (160 m) to 1.8 m is 89. These ratios are much larger than the stand value of 30, so the particle size values in this model are reasonable. The three-axle con ing pressure was loaded according to the values of 100 KPa, 200 KPa, and 300 KPa f the laboratory test. The trend in the stress-strain curves is the key to determining the relationship tween the meso-parameters and macro-parameters; this is also the micro-parameter bration method adopted by most scholars. On this basis, the stress-strain curves w obtained, as shown in Figure 7. Compared with the indoor triaxial test, the trends of two stress-strain curves obtained from the numerical triaxial test under the same amb pressure are essentially the same as those of the indoor test, and the peak deviation st error of the gangue is largely maintained within the range 0~1 KPa. In terms of the mechanical properties, we added the attenuation value and atten tion coefficient of two kinds of coal gangue as reference factors to enrich the mesosc parameter calibration index so that the macro effect reflected by the final paramete closer to that of the real material. This further reflects the brittleness of the coal gan material. In this paper, we quote a mechanical index and a rock strength attenuation efficient, which were proposed by Peng Jun et al. [42], to characterize the post-p strength attenuation behaviors of rock samples. The specific algorithm is: In terms of the mechanical properties, we added the attenuation value and attenuation coefficient of two kinds of coal gangue as reference factors to enrich the mesoscopic parameter calibration index so that the macro effect reflected by the final parameters is closer to that of the real material. This further reflects the brittleness of the coal gangue material. In this paper, we quote a mechanical index and a rock strength attenuation coefficient, which were proposed by Peng Jun et al. [42], to characterize the post-peak strength attenuation behaviors of rock samples. The specific algorithm is: where D s is the strength attenuation coefficient with a value range of [0, 1], ∆σ is the strength attenuation value, σ P 1 is the peak strength, and σ r is the residual strength. According to Formulas (27) and (28), the attenuation value and attenuation coefficient of the gangue were calculated as shown in Figure 8. It can be seen that with the increase in the confining pressure, the attenuation value of gangue strength will gradually increase, but the attenuation coefficient will decrease with the increase in the confining pressure, so the brittleness will gradually decrease. It is determined that the confining pressure of 100 KPa~300 KPa characterizes a process of gradually weakening ductility development for the gangue. Although brittleness gradually decreases, from the trend of the increasing intensity attenuation value, the brittleness of the gangue increases. This is because the internal particle structure distribution of the gangue is looser than that of general rocks, and there are differently sized gangue particles. Small gangue particles are tightened first, which weakens the brittleness. The brittleness is enhanced by contact within a large-sized gangue. Between 100 KPa and 300 KPa, the inner part of the gangue is characterized by the contact shrinkage process of a small-sized gangue, so its brittleness is gradually smaller. The error of the attenuation value is controlled within 10 KPa, and the error of the attenuation coefficient is controlled within 0.025; both parameters are within a reasonable error range.
According to Formulas (27) and (28), the attenuation value and attenuation cient of the gangue were calculated as shown in Figure 8. It can be seen that w increase in the confining pressure, the attenuation value of gangue strength will gr increase, but the attenuation coefficient will decrease with the increase in the co pressure, so the brittleness will gradually decrease. It is determined that the co pressure of 100 KPa~300 KPa characterizes a process of gradually weakening duct velopment for the gangue. Although brittleness gradually decreases, from the tren increasing intensity attenuation value, the brittleness of the gangue increases. Th cause the internal particle structure distribution of the gangue is looser than that eral rocks, and there are differently sized gangue particles. Small gangue parti tightened first, which weakens the brittleness. The brittleness is enhanced by within a large-sized gangue. Between 100 KPa and 300 KPa, the inner part of the is characterized by the contact shrinkage process of a small-sized gangue, so its bri is gradually smaller. The error of the attenuation value is controlled within 10 K the error of the attenuation coefficient is controlled within 0.025; both parame within a reasonable error range. In conjunction with Table 2, which shows the macro-mechanical parameters o from numerical triaxial tests, it can be concluded that the macro-mechanical char tics of the micro-mechanical parameters in this numerical triaxial test are essentia sistent with those of the laboratory test. In conclusion, the micro-parameters tested numerical triaxial test show that the macro-mechanical properties have a certain of reliability and that the PFC materials are available. The calibration values of t micro-parameters are shown in Table 3. In conjunction with Table 2, which shows the macro-mechanical parameters obtained from numerical triaxial tests, it can be concluded that the macro-mechanical characteristics of the micro-mechanical parameters in this numerical triaxial test are essentially consistent with those of the laboratory test. In conclusion, the micro-parameters tested via the numerical triaxial test show that the macro-mechanical properties have a certain degree of reliability and that the PFC materials are available. The calibration values of the PFC micro-parameters are shown in Table 3.  Table 3. Mesoscopic parameters of the PFC model of gangue.(D is particle density; E * is contact elastic modulus; K * is stiffness ratio; µ is friction coefficient; λ is parallel bonding radius product factor; K c * is parallel bonding stiffness ratio; E c * is elastic modulus of parallel bonding; σ c is parallel bonding tensile strength; c is parallel bonding force; ϕ is parallel bonding friction angle; d w is fluid density; v s is fluid viscosity coefficient.) (1) Monitoring points and measuring balls. Based on the above research, we designed a PFC gangue slope model with a size of 160 × 66 × 67 m (length × width × height); the total number of clump models was 71,971, as shown in Figure 9. The monitoring of the calculated data shown in this paper was mainly realized using a measuring ball and a monitoring point. A measuring sphere with a radius of 5.5 m was set and evenly distributed over the entire slope surface to monitor the real-time changes in the stress and strain values of the soil mass across the entire slope surface. Regarding the monitoring points, based on the damage at the no. 5 gully site, monitoring points were set at the top, middle, and lower surfaces (monitoring points 1, 2, 3, and 4) and side surfaces (monitoring point 5). In order to explore the scouring depth, one vertical monitoring point was set at 8 m below each surface-level monitoring point. There were 10 monitoring points in total, so the layout can accurately monitor the displacement of the gully position and the change in the movement track.
(2) Fluid-solid coupling model. When setting up the fluid-solid coupling model, the flow field was simulated as a fluid module, as shown in Figure 10. Many stability tests were carried out on the model to ensure that the imbalance force between particles was less than 10-5 N before the application of fluid. For the convenience of a later analysis, the slope area was divided into four areas: A, B, C, and D (catchment areas A and C; scouring area B; and permeation area D). In order to achieve a more realistic flushing simulation effect, based on previous scholars' models of flow field action [43], the flow field was divided into two layers: the flushing layer and the permeable layer. The thickness of both layers is set at 1 m and moves downwards with the destruction depth of soil particles until the desired depth and a new equilibrium state are reached. In addition, the penetration layer weakens the bond strength between particles and accelerates the formation of gullies, assuming the flow rate is 0. Many scholars pay attention to the simplicity of the flow field and cannot perfectly restore the complex turbulent process in the flow field when they study the fluid-solid coupling of the PFC. If the flow field is not specifically laid out, the particles on the slope will disperse everywhere, affecting the final failure mode and disturbing the clarity of the data. In order to invert the shape of gully no. 5 as far as possible, and to achieve the aims of our research, the scouring layer was divided into the scouring area and the catchment area. The scouring area was set according to the area of gully no. 5; see the dark area in Figure 10. According to this logic, it is also necessary to reduce the intensity of the fluid meso-parameters in the catchment area. At the same time, the flow rate in the PFC will deviate from the actual slope runoff speed. If the flow rate is too high, the unevenness of the slope's surface will cause particles to collide with each other and bounce off of the slope. If the velocity is too low, the expected effect will not be achieved due to the cohesion between the particles. After repeated trial calculations, the scouring speed was determined to be 2.3 m/s. When the slope particles are scoured, the original equilibrium state is broken, and the particles begin to gradually slide off the slope along with the fluid movement and accumulate at the foot of the slope until a new equilibrium state is reached. At this time, the particles are essentially stationary, and this is considered the end of the destruction process. After undertaking many trial calculations, it was found that when the calculation reaches 80 min, the particles are basically stationary, and the maximum velocity is less than 10 −4 m/s. (2) Fluid-solid coupling model. When setting up the fluid-solid coupling model, the flow field was simulated as a fluid module, as shown in Figure 10. Many stability tests were carried out on the model to ensure that the imbalance force between particles was less than 10-5 N before the application of fluid. For the convenience of a later analysis, the slope area was divided into four areas: A, B, C, and D (catchment areas A and C; scouring area B; and permeation area D). In order to achieve a more realistic flushing simulation effect, based on previous scholars' models of flow field action [43], the flow field was divided into two layers: the flushing layer and the permeable layer. The thickness of both layers is set at 1 m and moves downwards with the destruction depth of soil particles until the desired depth and a new equilibrium state are reached. In addition, the penetration layer weakens the bond strength between particles and accelerates the formation of gullies, assuming the flow rate is 0. Many scholars pay attention to the simplicity of the flow field and cannot perfectly restore the complex turbulent process in the flow field when they study the fluid-solid coupling of the PFC. If the flow field is not specifically laid out, the particles on the slope will disperse everywhere, affecting the final failure mode and disturbing the clarity of the data. In order to invert the shape of gully no. 5 as far as possible, and to achieve the aims of our research, the scouring layer was divided into the scouring area and the catchment area. The scouring area was set according to the area of gully no. 5; see the dark area in Figure 10. According to this logic, it is also necessary to reduce the intensity of the fluid meso-parameters in the catchment area. At the same time, the flow rate in the PFC will deviate from the actual slope runoff speed. If the flow rate is too high, the unevenness of the slope's surface will cause particles to collide with each other and bounce off of the slope. If the velocity is too low, the expected effect will not be achieved due to the cohesion between the particles. After repeated trial calculations, the scouring speed was determined to be 2.3 m/s. When the slope particles are scoured, the original equilibrium state is broken, and the particles begin to gradually slide off the slope along with the fluid movement and accumulate at the foot of the slope until a new equilibrium state is reached. At this time, the particles are essentially stationary, and this is considered the end of the destruction process. After undertaking many trial calculations, it was found that when the calculation reaches 80 min, the particles are basically stationary, and the maximum velocity is less than 10 −4 m/s.

Failure Mode
Under the action of fluid scouring and penetration, weak bonds between the surface particles of the gangue slope are lost, and the slope's mechanical strength characteristics are weakened, forming a gully, which affects the stability of the slope. Therefore, it is important to study the failure mode of the PFC gangue slope model of scouring under rainfall conditions to determine the gully formation mechanisms of such slopes and methods for slope protection and treatment.
We recorded the failure mode of the PFC gangue slope within 80 min of rainfall

Failure Mode
Under the action of fluid scouring and penetration, weak bonds between the surface particles of the gangue slope are lost, and the slope's mechanical strength characteristics are weakened, forming a gully, which affects the stability of the slope. Therefore, it is important to study the failure mode of the PFC gangue slope model of scouring under rainfall conditions to determine the gully formation mechanisms of such slopes and methods for slope protection and treatment.
We recorded the failure mode of the PFC gangue slope within 80 min of rainfall scouring; see Figure 11. The failure types of the whole process can be divided into four stages: (1) Integral splashing and local etching of the slope bottom. This phase corresponds to the 10-20 min period shown in Figure 11. The early stag of rainfall occurs from 0 to 10 min and no obvious slope erosion occurs. Clump particl do not disengage from the slope. However, a displacement shape similar to that of gul (1) Integral splashing and local etching of the slope bottom. This phase corresponds to the 10-20 min period shown in Figure 11. The early stage of rainfall occurs from 0 to 10 min and no obvious slope erosion occurs. Clump particles do not disengage from the slope. However, a displacement shape similar to that of gully 5 is formed in area B of the slope's surface. At this time, the scouring condition of the slope's surface is still in the early stage of catchment. The effect of hydraulic force on the clump particles is mainly concentrated in the valley. An arc-shaped radial flow channel is formed in the whole area (B), with the displacement of clump particles at 0-0.4 m and the maximum displacement concentrated at the bottom of the slope. As the hydraulic downward action destroys the particle contact model, the interaction between the particles at the bottom of the slope increases and the displacement increases. At 20 min, the clump particles on the surface of the slope begin to be destroyed and they disperse and disengage from the slope. The particles are taken away from the suspension module, and the maximum displacement is 9.79 m. This process is called sheet erosion, where rills appear on the bottom surface of the slope. At the same time, the catchment on the slope's surface gradually transfers from Zone C to Zone A, and the particle displacement can reach 1.5 m due to catchment erosion.
(2) Erosion diversion and sheet erosion stages in the middle section of the slope. The rainfall continues to progress, and the scouring time reaches 30 min. At this time, the denuded particles at the bottom of the slope gradually increase in number, forming a muddy shape and accumulating downward. At the same time, denudation develops upward and reaches the middle of the slope, leading to the start of the sheet erosion. It is not difficult to determine that the erosion resistance of the soil particles on the slope is less than the runoff scouring force, and the particles are in a local unstable state. At this time, the shear force between the runoff particles washes away the adjacent particles and forms a small erosion hole in the middle of the slope. Since the inter-particle hydraulic action gradually shifts to Zone A during the catchment stage, the main gully begins to flow in the direction of Zone A in the middle of the slope. It is interesting to note that the particles are eroded transversely in the direction of Zone A and form a Z-shaped scouring ditch with the rill at the bottom of the slope. After the rainfall scouring reaches 40 min, the downward scouring capacity of the main stream is greater than the lateral damage capacity of the catchment tributary. This results in the upward denudation speed of the main gully, which exceeds that of the lateral tributaries. With the progression of the rainfall, the runoff energy continues to increase. At the same time, due to the infiltration of the rainfall, the strength of the soil on the slope decreases, and the small erosion holes formed in the slope gradually deepen and widen towards the slope's top, forming headward erosion. The accumulation at the bottom of the slope also gradually diffuses; the shape of the gully is not limited to the strip gully, and the width of the gully gradually expands outward. At this time, the direction and structure of the mainstream and tributaries are formed, and the runoff channel gradually becomes visible.
With the continuous development of rill erosion, the main gully continues to widen under the action of hydraulic force. Up to 50 min, the displacement of the accumulated particles can reach 98.32 m. At this time, the main stream's erosion capacity is less than that of the tributaries. The upper part of the main stream gully is suspended from erosion, and the lower part is further accumulated downwards. The gully section is exposed (a fault phenomenon), forming a large erosion pit with a depth of about 8 m. At the same time, the denudation of the tributaries is ongoing, and the trend of the slope gully changes from two paths to one, forming a large closed Z-type gully. The flow path is diverted, so the scouring capacity is mainly concentrated in the direction of the branch in area A. Based on the tension between the hydraulic and particle forces, the tributary valleys extend and increase continuously. When the time point of 60 min is reached, the peak displacement of the accumulated particles reaches 122.23 m, the side valleys are penetrated, and the transverse particle tension is reduced. At this time, due to continuous rainfall, the force chain between particles near the section of the erosion pit in area B reaches the limit state. The direction of denudation changes from the branch direction in area A to the main direction in area B and continues to produce schistosity upwards, forming headward erosion at the top of the slope.
(4) Overall Gully Erosion and the Collapse Stage According to the phenomenon whereby two scouring diversions occurred in the previous stage, the effect of hydraulic pressure on particles reaches its peak at 70 min, and the scouring intensity is the most intense at this point. When the top erosion pit section is gradually destroyed and collapses, the energy at the top of the slope is gradually released, the particles begin to accelerate the erosion, and a large number of particles are washed to the bottom of the slope and accumulate outward. At this time, the maximum displacement of the accumulated particles is 122.24 m, which is similar to the displacement observed in the previous stage, indicating that the whole model is in a state of gradual convergence at this time. Finally, the model tends to converge at 80 min, and, at this time, the gully of the lateral tributary in Zone A also reaches the maximum depth. The particles on the slope top of Zone B are completely connected through the pipe, forming a complete and continuous multidirectional gully with a peak displacement of 149.72 m.
In contrast to the failure mode, the process of the rain-erosion-induced failure of the coal gangue slope has the characteristics of discontinuous failure. Moreover, the failure sequence is bottom of the slope-middle of the slope-lateral tributary-top of the slope. The tributary gully in area A on the slope is determined by the water flow trend in the early catchment stage, and the whole gully trend is diverted in the middle and late stages of rainfall. However, the main stream's hydraulic action soon causes the secondary diversion in area B and completes the penetration damage at the top. This phenomenon provides a favorable reference for proving the shear erosion theory and a further detailed analysis in the later stage.

Comparison between Theory and Simulation
After the PFC simulation calculations, the final gully failure pattern was obtained. The model angle was adjusted and the failure form of the coal gangue slope near Panzhou Electric Field in Guizhou Province was compared and observed, as shown in Figure 12. It can be seen from the figure that the shape of gully no. 5 simulated by the PFC is similar to the on-site failure shape of the slope, and the failure conditions of the main gully and the tributary gully are basically consistent. We achieved the preliminary aim of a PFC simulation. In order to further verify the rationality of the model and the theory, based on the rainfall scouring shear failure theory previously developed by our team, two types of scouring intensity factors were considered in depth, and the scouring depth optimization theory was further developed in combination with the characteristics of the PFC time effect calculations. Through a field survey, the data of one point every 10 m were input into Formula (22) to obtain the failure depth of each point after rainfall scouring and we compared it with the site failure depth and PFC simulation depth, as shown in Figure 13. Figure 14 shows the final failure critical line under three scouring failure conditions. Combining the two figures, we can see that the depth curve trends of the PFC calculation, the optimization theoretical calculations, and the actual data are consistent. The maximum damage point of gully no. 5 is located within the 30~50 m area in the middle of the slope, and the actual damage depth can reach 12 m. However, the peak depth calculated via the PFC and the optimization theory is relatively large and the curves are relatively close to each other. It can be seen from the profile that the scouring has formed a three-stage ladder effect, which is concentrated at the top of the slope, the middle of the slope, and the bottom of the slope (0~25 m, 25~65 m, and 65~110 m). This step effect fully verifies the intermittent characteristics of the erosion failure of the coal gangue slope. Since the diversion and fault phenomena of rainfall scouring are concentrated at the position of 30~60 m in the middle of the slope, the erosion pit here is the deepest. Finally, it can be seen that the depth calculated via the depth optimization scouring theory, the gully depth under the PFC simulation, and the actual failure depth have a good consistency, which proves the rationality of the simulation and the theoretical deduction. At the same time, the combination of these methods provides a valuable reference for early warnings and the hydrodynamic theory of coal gangue slope engineering disasters.

Displacement Analysis
The change in the displacement with the time of each monitoring point of the model is shown in Figure 15. The displacement curve can further elucidate the scouring and damage laws of the slope's bottom, middle, top, and tributary areas. From the change in the surface displacement (Figure 15a,c), the trends for the displacement curves of each point in the horizontal Y-direction and the vertical Z-direction are essentially consistent, indicating that the surface soil particles in the scouring area are subject to more hydraulic scouring when they follow the flow and the constraints between soil particles are lower, and when the movement direction is consistent with the slope direction. From the perspective of deep particles, the particles at the bottom of the slope and at the top of the slope show similar curve trends in the process of comparing the horizontal displacement curve with the vertical displacement curve; however, the particle displacement curve trends in the middle of the slope and in the tributary region are quite different. This shows that when the rainfall infiltration reaches the deep soil particles, the internal disturbance of the particles in the middle of the slope is greatest, which causes the water flow to change in the middle of the slope, thus forming a gully tributary.

Displacement Analysis
The change in the displacement with the time of each monitoring point of the model is shown in Figure 15. The displacement curve can further elucidate the scouring and damage laws of the slope's bottom, middle, top, and tributary areas. From the change in the surface displacement (Figure 15a,c), the trends for the displacement curves of each point in the horizontal Y-direction and the vertical Z-direction are essentially consistent, indicating that the surface soil particles in the scouring area are subject to more hydraulic scouring when they follow the flow and the constraints between soil particles are lower, and when the movement direction is consistent with the slope direction. From the perspective of deep particles, the particles at the bottom of the slope and at the top of the slope show similar curve trends in the process of comparing the horizontal displacement curve with the vertical displacement curve; however, the particle displacement curve trends in the middle of the slope and in the tributary region are quite different. This shows that when the rainfall infiltration reaches the deep soil particles, the internal disturbance of the particles in the middle of the slope is greatest, which causes the water flow to change in the middle of the slope, thus forming a gully tributary.

Displacement Analysis
The change in the displacement with the time of each monitoring point of the model is shown in Figure 15. The displacement curve can further elucidate the scouring and damage laws of the slope's bottom, middle, top, and tributary areas. From the change in the surface displacement (Figure 15a,c), the trends for the displacement curves of each point in the horizontal Y-direction and the vertical Z-direction are essentially consistent, indicating that the surface soil particles in the scouring area are subject to more hydraulic scouring when they follow the flow and the constraints between soil particles are lower, and when the movement direction is consistent with the slope direction. From the perspective of deep particles, the particles at the bottom of the slope and at the top of the slope show similar curve trends in the process of comparing the horizontal displacement curve with the vertical displacement curve; however, the particle displacement curve trends in the middle of the slope and in the tributary region are quite different. This shows that when the rainfall infiltration reaches the deep soil particles, the internal disturbance of the particles in the middle of the slope is greatest, which causes the water flow to change in the middle of the slope, thus forming a gully tributary. According to the analysis of the failure mode, we deduced that the coal gangue slope model is divided into four stages (0-20 min, 20-40 min, 40-60 min, and 60-80 min) in the process of scouring failure. Within the period 0-20 min, the surface particle displacement (MP1, MP2) at the bottom and middle of the slope increases linearly, reaches the maximum value, and tends to be flat at 10 min. The maximum horizontal displacement of the surface particle MP1 at the bottom of the slope during this period is 12 m, and the maximum vertical displacement is 2 m; the horizontal peak displacement of surface particle MP2 in the middle of the slope is 40 m, and the vertical peak displacement is 13 m. During this period, there is no change in displacement at the top of the slope and at the tributaries, which indicates that the first stage of scouring damage impacts the slope bottom first, but the damage trend is obviously more violent in the middle of the slope. When reaching the second stage (20-40 min), there is still no displacement change at the slope bottom, the slope top, or the tributaries. The particle displacement of the surface soil in the middle of the slope increases slowly, but the displacement of its deep particles (Figure 14d) fluctuates, indicating that the erosion damage is concentrated in the middle of the slope at this time. The surface particles are washed with the water flow, and the deep particles are damaged by turbulence. When the scouring reaches the third stage (40-60 min), the surface particle displacement curve between the slope top and the tributary area begins to increase sharply. The maximum horizontal displacement of the slope top can reach 102 m, accounting for 41.76% of the overall horizontal displacement, and the maximum vertical displacement of the slope top can reach 42 m, accounting for 45.85% (see Figure 16). It can be seen from the displacement curve that the displacement at the top of the slope changes prior to the displacement of the tributaries. This shows that after the scouring and accumulation in the middle of the slope, the tributary gap is formed, which is also the reason for the formation of multiple gullies caused by scouring. However, at this time, the stress is not concentrated in the direction of the tributary, but at the top of the slope, so the direction of the gully is not completely changed at this time. It can be seen from Figure 15b-d that the particle displacement curve in the deep part of the tributary changes slowly and is unstable, indicating that the scouring force of the tributary is not strong before the main gully is connected; this again verifies that the stress in the third stage is more highly concentrated at the top of the slope. In the fourth stage (60-80 min), except for the deep particles in the tributary, the particle displacement in other areas tends to be stable, and the particles in the monitoring area have largely been washed to the bottom of the slope with the rain and completed their accumulation. The deep particles in the tributary are still involved in the process of the continuous growth of displacement, which According to the analysis of the failure mode, we deduced that the coal gangue slope model is divided into four stages (0-20 min, 20-40 min, 40-60 min, and 60-80 min) in the process of scouring failure. Within the period 0-20 min, the surface particle displacement (MP1, MP2) at the bottom and middle of the slope increases linearly, reaches the maximum value, and tends to be flat at 10 min. The maximum horizontal displacement of the surface particle MP1 at the bottom of the slope during this period is 12 m, and the maximum vertical displacement is 2 m; the horizontal peak displacement of surface particle MP2 in the middle of the slope is 40 m, and the vertical peak displacement is 13 m. During this period, there is no change in displacement at the top of the slope and at the tributaries, which indicates that the first stage of scouring damage impacts the slope bottom first, but the damage trend is obviously more violent in the middle of the slope. When reaching the second stage (20-40 min), there is still no displacement change at the slope bottom, the slope top, or the tributaries. The particle displacement of the surface soil in the middle of the slope increases slowly, but the displacement of its deep particles (Figure 14d) fluctuates, indicating that the erosion damage is concentrated in the middle of the slope at this time. The surface particles are washed with the water flow, and the deep particles are damaged by turbulence. When the scouring reaches the third stage (40-60 min), the surface particle displacement curve between the slope top and the tributary area begins to increase sharply. The maximum horizontal displacement of the slope top can reach 102 m, accounting for 41.76% of the overall horizontal displacement, and the maximum vertical displacement of the slope top can reach 42 m, accounting for 45.85% (see Figure 16). It can be seen from the displacement curve that the displacement at the top of the slope changes prior to the displacement of the tributaries. This shows that after the scouring and accumulation in the middle of the slope, the tributary gap is formed, which is also the reason for the formation of multiple gullies caused by scouring. However, at this time, the stress is not concentrated in the direction of the tributary, but at the top of the slope, so the direction of the gully is not completely changed at this time. It can be seen from Figure 15b-d that the particle displacement curve in the deep part of the tributary changes slowly and is unstable, indicating that the scouring force of the tributary is not strong before the main gully is connected; this again verifies that the stress in the third stage is more highly concentrated at the top of the slope. In the fourth stage (60-80 min), except for the deep particles in the tributary, the particle displacement in other areas tends to be stable, and the particles in the monitoring area have largely been washed to the bottom of the slope with the rain and completed their accumulation. The deep particles in the tributary are still involved in the process of the continuous growth of displacement, which means that after the completion of the main gully's penetration failure, the energy is completely transferred to the tributary and continues to complete the erosion failure of the tributary. The surface soil is damaged by laminar flow, the displacement of particles is large, and the curve changes in a stepwise manner. The deep soil is greatly affected by turbulence, and the displacement is relatively small, but the displacement curve fluctuates greatly; particles collide with each other and even rebound. During the whole process, the maximum displacement is concentrated at the top of the slope, and its proportions are as follows: top of the slope > tributaries > middle of the slope > bottom of the slope. The peak displacement of the slope top in the horizontal Y-direction accounts for 41.76%, and that in the Z-direction accounts for 45.84% ( Figure 16). It is surmised that in the process of preventing and controlling rainfall-induced landslides, priority should be given to protecting the sliding soil mass at the bottom of the slope, and the soil mass in the middle of the slope should be reduced to prevent the maximum energy fluctuation caused by scouring, so as to prevent extensive displacement damage at the top of the slope.
Appl. Sci. 2023, 13, x FOR PEER REVIEW 21 of 28 means that after the completion of the main gully's penetration failure, the energy is completely transferred to the tributary and continues to complete the erosion failure of the tributary. The surface soil is damaged by laminar flow, the displacement of particles is large, and the curve changes in a stepwise manner. The deep soil is greatly affected by turbulence, and the displacement is relatively small, but the displacement curve fluctuates greatly; particles collide with each other and even rebound. During the whole process, the maximum displacement is concentrated at the top of the slope, and its proportions are as follows: top of the slope > tributaries > middle of the slope > bottom of the slope. The peak displacement of the slope top in the horizontal Y-direction accounts for 41.76%, and that in the Z-direction accounts for 45.84% (Figure 16). It is surmised that in the process of preventing and controlling rainfall-induced landslides, priority should be given to protecting the sliding soil mass at the bottom of the slope, and the soil mass in the middle of the slope should be reduced to prevent the maximum energy fluctuation caused by scouring, so as to prevent extensive displacement damage at the top of the slope.

Scour Tracking and Accumulation
During the numerical calculations, the positions of eight typical monitoring particles were selected, as shown in Figure 17, and the tracked particles were analyzed in combination with the results of the failure mode. At the bottom of the slope, MP1 is preferentially denuded by the thin-layer water flow, and MP6 is greatly affected by turbulence, but this effect is limited to the large particle resistance at the bottom of the slope; the movement tracks of the two are not obvious, and they are characterized by a slow creeping phenomenon, which is defined as a slow accumulation at the bottom of the slope. With the development of the main gully, particles MP2 and MP7 in the middle of the slope move along the main gully to the bottom of the slope under the scouring action, and then accumulate in the Y-direction along with the particles at the bottom of the slope. The mixing effect of the huge laminar flow and turbulent flow in the middle of the slope causes the development of tributaries. It can be seen that the MP5 and MP10 particles move from the side to the middle of the slope, then enter the mainstream gully and accumulate in the positive Y-direction at the bottom of the slope. We can clearly see that the particles at different depths in the tributary area have different stacking heights at the bottom of the slope, and the deep particles accumulate on the surface particles due to the later denudation sequence. Therefore, in the deposits at the bottom of the slope, the soil particles washed out from the deep layer cover the soil particles that originally accumulated on the surface layer. Finally, the two particles MP3 and MP4 are scoured at the top of the slope. They all move along the direction of the main gully from the top of the slope to the middle

Scour Tracking and Accumulation
During the numerical calculations, the positions of eight typical monitoring particles were selected, as shown in Figure 17, and the tracked particles were analyzed in combination with the results of the failure mode. At the bottom of the slope, MP1 is preferentially denuded by the thin-layer water flow, and MP6 is greatly affected by turbulence, but this effect is limited to the large particle resistance at the bottom of the slope; the movement tracks of the two are not obvious, and they are characterized by a slow creeping phenomenon, which is defined as a slow accumulation at the bottom of the slope. With the development of the main gully, particles MP2 and MP7 in the middle of the slope move along the main gully to the bottom of the slope under the scouring action, and then accumulate in the Y-direction along with the particles at the bottom of the slope. The mixing effect of the huge laminar flow and turbulent flow in the middle of the slope causes the development of tributaries. It can be seen that the MP5 and MP10 particles move from the side to the middle of the slope, then enter the mainstream gully and accumulate in the positive Y-direction at the bottom of the slope. We can clearly see that the particles at different depths in the tributary area have different stacking heights at the bottom of the slope, and the deep particles accumulate on the surface particles due to the later denudation sequence. Therefore, in the deposits at the bottom of the slope, the soil particles washed out from the deep layer cover the soil particles that originally accumulated on the surface layer. Finally, the two particles MP3 and MP4 are scoured at the top of the slope. They all move along the direction of the main gully from the top of the slope to the middle of the slope. However, due to the influence of height difference accumulation, when moving to the bottom of the slope, the two particles change their motion trajectory: the MP3 particles move diagonally in the XY-direction, while the MP4 particles move in the positive Y-direction. Therefore, we know that MP6 and MP1 tracks belong to slow accumulation movement, MP7 and MP2 and MP3 and MP4 tracks belong to fast scouring movement, and MP5 and MP10 belong to DC steering movement. Depending on the scouring location, the resulting particle trajectories vary considerably.
Appl. Sci. 2023, 13, x FOR PEER REVIEW 22 of 28 of the slope. However, due to the influence of height difference accumulation, when moving to the bottom of the slope, the two particles change their motion trajectory: the MP3 particles move diagonally in the XY-direction, while the MP4 particles move in the positive Y-direction. Therefore, we know that MP6 and MP1 tracks belong to slow accumulation movement, MP7 and MP2 and MP3 and MP4 tracks belong to fast scouring movement, and MP5 and MP10 belong to DC steering movement. Depending on the scouring location, the resulting particle trajectories vary considerably.  Interestingly, we combined several sets of motion trajectory curves and changed the perspective to obtain Figure 18. It can be seen from Figure 18a that there are two gullies in the whole scouring process, which are defined as the first and second gullies (i.e., tributary gullies and mainstream gullies) in order of completion. A further analysis shows that there are two modes of scouring pile integration, namely, the arc erosion accumulation mode (Figure 18b) and the fan erosion accumulation mode (Figure 18c). Arc erosion accumulation is mainly controlled by deep particles, namely, MP3, MP6, MP7, and MP10. Fan erosion accumulation is mainly affected by shallow particles, namely, MP1, MP2, MP4, and MP5. In the actual rainfall scouring process, the mainstream gullies rarely exist alone, and most of the slope erosion failure modes are characterized by the combination of multiple tributaries and the main stream, which is the effect of multiple gullies (only two gullies are considered in this model). The main stream gully mainly controls arc-shaped or linear erosion and accumulation. Due to their complex role in the slope, the interaction between laminar flow and turbulent flow is pronounced, which provides conditions for the development of tributaries and further forms the fan erosion and accumulation model. This provides favorable conditions for us to further study the internal mechanisms of coal gangue scouring and accumulation. Interestingly, we combined several sets of motion trajectory curves and changed the perspective to obtain Figure 18. It can be seen from Figure 18a that there are two gullies in the whole scouring process, which are defined as the first and second gullies (i.e., tributary gullies and mainstream gullies) in order of completion. A further analysis shows that there are two modes of scouring pile integration, namely, the arc erosion accumulation mode (Figure 18b) and the fan erosion accumulation mode (Figure 18c). Arc erosion accumulation is mainly controlled by deep particles, namely, MP3, MP6, MP7, and MP10. Fan erosion accumulation is mainly affected by shallow particles, namely, MP1, MP2, MP4, and MP5. In the actual rainfall scouring process, the mainstream gullies rarely exist alone, and most of the slope erosion failure modes are characterized by the combination of multiple tributaries and the main stream, which is the effect of multiple gullies (only two gullies are considered in this model). The main stream gully mainly controls arcshaped or linear erosion and accumulation. Due to their complex role in the slope, the interaction between laminar flow and turbulent flow is pronounced, which provides conditions for the development of tributaries and further forms the fan erosion and accumulation model. This provides favorable conditions for us to further study the internal mechanisms of coal gangue scouring and accumulation.

Stress-Strain Analysis
Finally, in order to elucidate the internal failure mechanism of the model, the four representative measuring balls in the model were selected, their stress-strain data were extracted, and the real-time stress-strain curves of the PFC scouring model were obtained, as shown in Figure 19. At the beginning of the scouring process, due to the initial contact between rainwater and the coal gangue particles, the overall stress on the slope surface increased, with a maximum of 4.9 MPa located at the bottom of the slope. In the later part of the first stage of scouring, the stress change in each area is not obvious, but the strain in the slope changes sharply, fluctuating between 5 and 8. When reaching the second stage (20-40 min), there are large stress fluctuation values at the bottom and middle of the slope, and the stress increases sharply at 20 min. At this time, a large number of particles in the middle of the slope accumulate towards the bottom of the slope. The large accumulation causes the stress at the bottom of the slope to increase sharply, with a fluctuation value between 2 and 6.8 MPa. It is worth noting that the stress generated by accumulation is

Stress-Strain Analysis
Finally, in order to elucidate the internal failure mechanism of the model, the four representative measuring balls in the model were selected, their stress-strain data were extracted, and the real-time stress-strain curves of the PFC scouring model were obtained, as shown in Figure 19. At the beginning of the scouring process, due to the initial contact between rainwater and the coal gangue particles, the overall stress on the slope surface increased, with a maximum of 4.9 MPa located at the bottom of the slope. In the later part of the first stage of scouring, the stress change in each area is not obvious, but the strain in the slope changes sharply, fluctuating between 5 and 8. When reaching the second stage (20-40 min), there are large stress fluctuation values at the bottom and middle of the slope, and the stress increases sharply at 20 min. At this time, a large number of particles in the middle of the slope accumulate towards the bottom of the slope. The large accumulation causes the stress at the bottom of the slope to increase sharply, with a fluctuation value between 2 and 6.8 MPa. It is worth noting that the stress generated by accumulation is greater than the stress caused by erosion and denudation. At this time, the stress and strain in the tributary area also begin to increase slightly, with the stress reaching 1.7 MPa and the strain reaching 5.1, indicating that the later part of the second stage is the incubation period of the tributary gully. After that, in the third stage (40-60 min), the bottom of the slope continues to accumulate, and the stress continues to fluctuate. However, due to the dense particles at the bottom of the slope, the particles are denuded less, and there is no obvious strain. At this time, the stress of the tributary increases sharply to 3.8 MPa, the tributary notch begins to form, and the strain reaches six. Then, the stress concentration is transferred from the tributary to the top of the slope, the stress of the tributary decreases to 0.4 MPa, and the stress at the top of the slope increases to 4.4 MPa. The trend of the strain curve is similar: the maximum strain at the top of the slope reaches 7.5, and then decreases sharply. At the fourth stage (60-80 min), the stress in the tributary area increases again, reaching 3.9 MPa, and the final scouring damage is completed. An alternation can be seen in the change in stress and strain at the top of the slope and in the tributary after the second stage of scouring. During the whole scouring process, the stress generated by the accumulation at the bottom of the slope is always the largest, but the damage it causes is the smallest. The strain fluctuation caused by the scouring in the middle of the slope is the largest, and the damage it causes is the strongest. greater than the stress caused by erosion and denudation. At this time, the stress and strain in the tributary area also begin to increase slightly, with the stress reaching 1.7 MPa and the strain reaching 5.1, indicating that the later part of the second stage is the incubation period of the tributary gully. After that, in the third stage (40-60 min), the bottom of the slope continues to accumulate, and the stress continues to fluctuate. However, due to the dense particles at the bottom of the slope, the particles are denuded less, and there is no obvious strain. At this time, the stress of the tributary increases sharply to 3.8 MPa, the tributary notch begins to form, and the strain reaches six. Then, the stress concentration is transferred from the tributary to the top of the slope, the stress of the tributary decreases to 0.4 MPa, and the stress at the top of the slope increases to 4.4 MPa. The trend of the strain curve is similar: the maximum strain at the top of the slope reaches 7.5, and then decreases sharply. At the fourth stage (60-80 min), the stress in the tributary area increases again, reaching 3.9 MPa, and the final scouring damage is completed. An alternation can be seen in the change in stress and strain at the top of the slope and in the tributary after the second stage of scouring. During the whole scouring process, the stress generated by the accumulation at the bottom of the slope is always the largest, but the damage it causes is the smallest. The strain fluctuation caused by the scouring in the middle of the slope is the largest, and the damage it causes is the strongest. The peak values of the measurement ball data of the whole model were extracted, and the stress and strain were divided into the three directions of the XYZ data; the peak stress and strain curve is shown in Figure 20. It can be seen from the figure that the higher peak stress and strain values of the coal gangue slope scouring failure are distributed in the horizontal Y-direction of the middle of the slope, the horizontal X-direction of the tributary, and the vertical Z-direction of the slope bottom. From the peak data in the X-direction, it can be seen that the transverse peak stress increases from the top of the slope to the bottom of the slope, and the maximum value can reach 8.6 MPa in the tributary area. At this time, the peak strain also reaches the maximum value of 10.3 in the tributary region, and then decreases at the bottom of the slope (because the bottom of the slope is dense). This shows that the generation of the branch notch is mainly determined by the X-direction stress, and 8.6 MPa is the critical stress. According to the Y-direction curve, the scouring damage in the middle of the slope is the most severe, with a peak stress of 8.9 MPa and a peak strain of 14.3. The Y-direction peak stress at the top of the slope is close to that of the tributary, and both remain between 5.5 and 6.2 MPa. It can be observed that in the distance from the middle of the slope to the tributary to the bottom of the slope, the peak stress-strain curve in the Y-direction shows a "U"-type distribution feature. However, The peak values of the measurement ball data of the whole model were extracted, and the stress and strain were divided into the three directions of the XYZ data; the peak stress and strain curve is shown in Figure 20. It can be seen from the figure that the higher peak stress and strain values of the coal gangue slope scouring failure are distributed in the horizontal Y-direction of the middle of the slope, the horizontal X-direction of the tributary, and the vertical Z-direction of the slope bottom. From the peak data in the X-direction, it can be seen that the transverse peak stress increases from the top of the slope to the bottom of the slope, and the maximum value can reach 8.6 MPa in the tributary area. At this time, the peak strain also reaches the maximum value of 10.3 in the tributary region, and then decreases at the bottom of the slope (because the bottom of the slope is dense). This shows that the generation of the branch notch is mainly determined by the X-direction stress, and 8.6 MPa is the critical stress. According to the Y-direction curve, the scouring damage in the middle of the slope is the most severe, with a peak stress of 8.9 MPa and a peak strain of 14.3. The Y-direction peak stress at the top of the slope is close to that of the tributary, and both remain between 5.5 and 6.2 MPa. It can be observed that in the distance from the middle of the slope to the tributary to the bottom of the slope, the peak stress-strain curve in the Y-direction shows a "U"-type distribution feature. However, combined with the overall analysis, the stress response in the middle of the slope is stronger; this area can be described as the hub area for the development of scouring and destruction. The contact force chain of its internal particles is more likely to break and is more destructive. It can be seen from the Z-direction data that due to gravity and scouring, the energy reaching the bottom of the slope gradually increases; moreover, the peak stress gradually increases from the top of the slope to the bottom of the slope. The maximum Z-direction peak stress at the bottom of the slope can reach 10.3 MPa, but the peak strain here is only 4.6. According to the law of the stress-strain curve, the internal mechanisms of erosion failure in coal gangue slopes that were emphasized in the previous analysis are more clearly defined. combined with the overall analysis, the stress response in the middle of the slope is stronger; this area can be described as the hub area for the development of scouring and destruction. The contact force chain of its internal particles is more likely to break and is more destructive. It can be seen from the Z-direction data that due to gravity and scouring, the energy reaching the bottom of the slope gradually increases; moreover, the peak stress gradually increases from the top of the slope to the bottom of the slope. The maximum Zdirection peak stress at the bottom of the slope can reach 10.3 MPa, but the peak strain here is only 4.6. According to the law of the stress-strain curve, the internal mechanisms of erosion failure in coal gangue slopes that were emphasized in the previous analysis are more clearly defined.

Discussion
This paper optimizes the wash shear failure theory and combines it with the discrete element PFC3D method to verify its rationality, which provides a new idea for the PFC flu-structure interaction method. In addition, the article enriches the modeling method of the PFC coal gangue slope model and provides a valuable reference for the early warning and simulation of coal gangue slope engineering disasters. However, when building the PFC3D model, this paper assumes that the model is an ideal boundary model, and the types of clusters are limited to three, and the relationship between the ratio of various clusters and the model parameters is not considered in depth. At the same time, the setting of the fluid module is too small, and the fluid-structure coupling design requires further improvement. Although only one kind of soil material (coal gangue) is considered in this study, the clump modeling method, the test process of particle flow meso-parameters and the fluid-structure coupling model of erosion damage can be further used for the simulation of clay and sand. Interestingly, we happened to have conducted a project on the improvement of red clay in Sichuan, China, recently, so the simulation method will be further improved. This paper is based on the real geometric shape of the scene to model, which is representative but not comparative. In the next research work, the rainfall erosion mechanism of different slopes with different angles and shapes will be the focus. The influence of pore water pressure on the failure mechanism of a coal gangue slope is not discussed in this paper, which requires further study.

Conclusions
(1) Based on the previously derived rainfall scouring shear failure theory, a scouring depth optimization theory was further developed by taking into account the two types of scouring intensity factors (c and φ) and combining the characteristics of the

Discussion
This paper optimizes the wash shear failure theory and combines it with the discrete element PFC3D method to verify its rationality, which provides a new idea for the PFC flu-structure interaction method. In addition, the article enriches the modeling method of the PFC coal gangue slope model and provides a valuable reference for the early warning and simulation of coal gangue slope engineering disasters. However, when building the PFC3D model, this paper assumes that the model is an ideal boundary model, and the types of clusters are limited to three, and the relationship between the ratio of various clusters and the model parameters is not considered in depth. At the same time, the setting of the fluid module is too small, and the fluid-structure coupling design requires further improvement. Although only one kind of soil material (coal gangue) is considered in this study, the clump modeling method, the test process of particle flow meso-parameters and the fluid-structure coupling model of erosion damage can be further used for the simulation of clay and sand. Interestingly, we happened to have conducted a project on the improvement of red clay in Sichuan, China, recently, so the simulation method will be further improved. This paper is based on the real geometric shape of the scene to model, which is representative but not comparative. In the next research work, the rainfall erosion mechanism of different slopes with different angles and shapes will be the focus. The influence of pore water pressure on the failure mechanism of a coal gangue slope is not discussed in this paper, which requires further study.

Conclusions
(1) Based on the previously derived rainfall scouring shear failure theory, a scouring depth optimization theory was further developed by taking into account the two types of scouring intensity factors (c and ϕ) and combining the characteristics of the PFC time effect calculation. Verification shows that the curve trends of the PFC calculation, the optimization theory calculation, and the actual slope depth data are consistent.
The combination of these methods provides a valuable reference for early warnings of coal gangue slope engineering disasters and for the hydrodynamic theory. (2) The rain erosion failure of the coal gangue slope can be divided into four stages: (1) overall splash erosion and local schist erosion at the bottom of the slope; (2) erosion diversion and sheeting in the middle of the slope; (3) the tributary-slope crest extension and etching stage; and (4) the integral gully erosion and landslide stage. The scouring failure process has the characteristics of a discontinuous failure. Moreover, the failure sequence is as follows: slope bottom-slope middle-lateral tributary-slope top. The tributary gully in area A on the slope is determined by the water flow trend in the early catchment stage. The process of coal gangue scouring and destruction forms a three-stage ladder effect, which is concentrated at the top of the slope, the middle of the slope, and the bottom of the slope (0~25 m, 25~65 m, and 65~110 m). This step effect fully verifies the intermittent characteristics of the erosion failure of the coal gangue slope. (3) During the whole process, the maximum displacement is concentrated at the top of the slope, according to the following proportions: top of the slope > tributaries > middle of the slope > bottom of the slope. The peak displacement of the slope top in the horizontal Y-direction accounts for 41.76%, and that in the Z-direction accounts for 45.84%. In efforts to prevent and control rainfall-induced landslides, priority should be given to the protection of the sliding soil mass at the bottom of the slope, and the focus should be on reducing the soil mass in the middle of the slope to prevent the maximum energy fluctuation caused by scouring, so as to prevent extensive displacement damage at the top of the slope. The latter part of the second stage of scouring is the incubation period for the tributary gully. During the whole scouring process, the strain fluctuation caused by scouring in the middle of the slope is the largest; this is the hub area for the development of scouring damage. The peak stress is up to 8.9 MPa, the peak strain is up to 14.3, and the damage caused by scouring is the strongest. The larger peak stress and strain values of coal gangue slope scouring failure are distributed in the horizontal Y-direction of the middle of the slope, the horizontal X-direction for the tributary, and the vertical Z-direction for the slope bottom. The generation of the branch notch is mainly determined by the X-direction stress, and 8.6 MPa is the critical stress. (4) There are two modes of erosion accumulation integration, namely, the arc erosion accumulation mode and the fan erosion accumulation mode. Arc-shaped erosion and accumulation are mainly controlled by deep particles, and fan-shaped erosion and accumulation are mainly affected by shallow particles. In the actual process of rainfall erosion, most of the slope erosion and failure modes show the effect of multi-ditch collection. The main stream gully mainly controls arc-shaped or linear erosion and accumulation. Due to their complex roles in the slope, the interaction between laminar flow and turbulent flow is intense, providing conditions for the development of tributaries and further forming the fan erosion and accumulation model. This provides favorable conditions for a more effective study of the internal mechanisms of coal gangue scouring and accumulation.

Data Availability Statement:
The datasets generated for this study are available from the corresponding author upon reasonable request.