## 2 Erosion Theory
### Shear Erosion Principle  
We subdivided the coal gangetic slope near the No. 3 gully and established erosion shear and rainwater infiltration models (Figures 3 and 4). The slope model had a height of about 66 m, a width of 66 m, and the incline of the slope was 25-38\({}^{\circ}\). As erosion mainly occurs during heavy rainstorms, the erosion parts are mainly concentrated in the natural gully (near the No. 3 gully). We ignored the raindrop splash erosion in the early stage of rainfall, and only considered the water flow erosion and infiltration of the No. 3 gully on the slope surface. We used the study by [PERSON] et al. [28], on the erosion shear theory of gravel soil, in order to calculate the erosion depth of the soil over time for derivation of the erosion theory. However, the considered speed factor includes a certain error range, which cannot match the particle mesoscopic motion speed in discrete element PFC and, therefore, needed to be refined. During heavy rain, the erosion phenomenon of the viscous gravel soil slope presents a turbulent water flow pattern. The water flow is affected by the soil particles on the slope, forming many vortices of different sizes, with positions, shapes, and flow rates constantly changing as it progresses. At first, we express the shear stress of the slope surface flow [29] using the following equations:  
\[\tau_{uv}=\tau_{vu}+\tau_{vvb}, \tag{1}\]  
\[\tau_{vu}=-\rho v_{x}^{\prime}v_{y}^{\prime}, \tag{2}\]  
where \(\tau_{vu}\) is the viscous shear stress in laminar motion, \(\tau_{vib}\) is the turbulent shear stress, \(\rho\) is the flow density, and \(v_{x}^{\prime}\) and \(v_{y}^{\prime}\) as the pulsating flow velocities in the \(x\) and \(y\) coordinate directions, respectively (Figure 4). Defined by the pulsation speed, the velocities are calculated from the following equations:  
\[ \begin{cases}v_{x}^{\prime}=v_{x0}-\overline{v}_{x},\\ v_{y}^{\prime}=v_{y0}-\overline{v}_{y},\end{cases}  \tag{3}\]  
where \(v_{x0}\) and \(v_{y0}\) are the instantaneous velocities of soil particles, and \(\overline{v}_{x}\) and \(\overline{v}_{y}\) are the average water flow velocities. According to the definitions of water flow velocity and infiltration velocity, the pulsating infiltration velocity in the y-direction is expressed by the equation:  
\[v_{y}^{\prime}=v_{x}^{\prime}\cdot\xi=(v_{x0}-\overline{v}_{x})\cdot\xi, \tag{4}\]  
where \(\xi\) is the effective porosity between particles. According to the assumption of Prandtl and the hydraulic roughness test of Nikuradse, the following equation can be derived [30]:  
\[\tau_{vub}=\left(\frac{\overline{v}_{x}}{2.5!n\frac{R}{\overline{\mu}}+4.75} \right)^{2}\rho, \tag{5}\]  
where \(R\) is the hydraulic radius and \(\mu\) is the slope roughness. According to previous experimental studies [31; 32], the surface roughness can be described by the average particle size of soil particles (i.e., \(\mu=\overline{d}\), where \(\overline{d}\) is the average particle size of PFC particles). Substituting Equations (2)-(5) into Equation (1) yields:  
\[\tau_{vv}=\left(\frac{\overline{v}_{x}}{2.5!n\frac{R}{\overline{d}}+4.75} \right)^{2}\rho-(v_{x0}-\overline{v}_{x})^{2}\cdot\xi\rho=\rho\left[\left( \frac{\overline{v}_{x}}{2.5!n\frac{R}{\overline{d}}+4.75}\right)^{2}-\xi(v_{ x0}-\overline{v}_{x})^{2}\right]. \tag{6}\]  
Figure 4: Shear erosion theory model. The instantaneous velocity of soil particles can be monitored in the discrete element PFC calculation, while the average velocity of the fluid and the effective porosity of the particles are known quantities. Therefore, the optimized shear stress of water flow (Equation (6)) can be applied to the simulation calculation in PFC. For rock and soil slopes, the erosion ability of the slope surface should consider the erosion intensity of water flow. Coal gangue is a relatively coarse scattered rock and soil particle, which generally appears as agglomerates, lumps, or flakes during the erosion process.