Rainfall Erosion Damage of Residual Soil Slope in Intermittently Frozen Area Based on Discrete Element Method

This study developed a discrete element simulation model based on the 2D particle flow code (PFC2D), of which the mesoscopic parameters were calibrated by the indoor experiments, to investigate the rainfall erosion damage of residual soil slope in the intermittently frozen area. It is to be noted that the runoff scouring action was simulated according to the equivalent rainfall method, the soil particles on the slope were given initial velocity, and the water absorption was considered by increasing the unit weight. The results indicated that the scouring action only caused superficial erosion with the main damage region at the foot, regardless of the FT effect. A splitting phenomenon was observed in the lower part of the steeper slope under the FT effect. Moreover, regardless of the FT effect, the gentler slope tended to incur spalling rather than a splitting phenomenon, where the soil particles slid along the structural plane with strong anti-scouring ability. Besides, the gentler slope maintained higher stability and shorter scouring time. Finally, the scouring velocity increased the erosion damage to a large extent.


INTRODUCTION
The intermittently frozen area was defined as an area where the near-surface soil is frozen annually from one to fifteen days per year (Willmott and Matsuura, 2001;Zhang et al., 2003;Chen, et al., 2020). Intermittently frozen areas are widely distributed in the world, accounting for 6.6% of the exposed lands in the coldest month of the year in the Northern Hemisphere (Baranov, 1964;Zhang et al., 2003;Lu et al., 2017). Taking China, for example, it is mostly distributed in the south of the Yangtze River and the north of the Pearl River Basin, covering an area of about 1.9 million square kilometers and accounting for about 20% of the total national area (Que et al., 2017;. Compared with the permafrost and seasonal frozen soils, the intermittently frozen soil experiences high-frequency freeze-thaw cycles, resulting from the frequently fluctuated ambient temperature around the freezing point (Woo and Winter, 1993;Chen and Li, 2008;Thomas et al., 2009). It should be noted that Que et al. (2017) observed 16 freeze-thaw times a day in December 2013 to mid-February 2014 in Wuyi Mountain, Fujian Province, China. Usually, it will experience superficial spalling occurring on the soil slope, especially under the influence of continuous rainfall or rainstorm. That might develop the instability of soil further, causing potential safety problems and economic loss to the vehicles and pedestrians along with the projects (Que et al., 2017;Li, 2013). Therefore, there is a great need to study further the strength reduction of soil slopes under the shortterm freeze-thaw cycle and the mechanism of superficial spalling induced by rainfall erosion.
Recently, the discrete element technique, like particle flow code (PFC) (Zhang et al., 2013;Liang et al., 2017;Wang et al., 2020;, has been found to be suitable for not only dealing with large deformation discontinuity problems but also illustrating the force and motion pattern from the mesoscopic scale Saman et al., 2014;Shi and Xu, 2015;Xu, 2017;Zhang, 2017;. Hence, it has been currently used in geotechnical engineering, including the superficial freeze-thaw spalling of slope (Evans et al., 2009;Park and Song, 2009;Zhu et al., 2021).
Several studies have been carried out on the erosion of the slope surface caused by rainfall through PFC. Tsuji et al. (1993) introduced computational fluid dynamics to the discrete element method and verified the reliability of the fluid-solid coupling analysis method in simulating two-phase flows using discrete elements. Zeghal and Shamy (2005) implemented the fluid-solid coupling analysis method to analyze saturated sand, which simplified the N-S equation and by which the seepage problem of the soil slope was investigated. Based on this, Bai et al. (2012) determined the permeability coefficient for the PFC method. Xiong et al. (2013) performed a threedimensional simulation of the rainfall slope collapse by fluid-solid coupling analysis, where the rainfall was simulated by applying pressure on the boundary. The results manifested that it was feasible to study the collapse formation mechanism via PFC using the fluid-solid coupling methods. Song (2013) and Wu et al. (2014) analyzed the evolution process of slope erosion and the distribution of slope erodibility by using the aforementioned method. Hu (2014) performed the rainfall equivalent via increasing the gravity and reducing shear strength, but the strength-reduced parameters were not calibrated. Ma (2012) and Wu et al. (2013) carried out an erosion failure simulation of steep slopes via a two-dimensional PFC (PFC 2D ) model, where the effect of runoff scouring was simulated by the particles' movement on the slope surface.
Overall, it can be found that rainfall erosion can be simulated through the fluid-solid coupling analysis method via PFC. However, most studies were limited to the saturated soil. There is no research related to the rainfall equivalent method. Besides, the parameter effects are still unclear either, such as slope gradient, freeze-thaw, and scouring velocity.
To fulfill this gap, this paper carries out an in-depth study about the effects of rainfall erosion damage of residual soil slope in the intermittently frozen area. Specifically, this paper addresses three questions: How to build a slope scouring model? What are the basic effects of scouring on the slope? How the slope gradient, freeze-thaw action, and scouring velocity affect the erosion damage? The objective of the present study is threefold: 1) to build a PFC 2D scouring model, where the mesoscopic parameters are calibrated by indoor experiments, 2) to investigate the basic characteristics of slope under the runoff erosion, and 3) to perform parameter effects on the rainfall erosion, including slope gradient, freeze-thaw action, and scouring velocity.

METHODOLOGY Experimental Design
As shown in Figure 1, in order to investigate the effects of freeze-thaw on erosion damage of the residual soil slope induced by rainfall, the indoor test and PFC 2D simulation are performed in depth. First, the freeze-thaw experiments are conducted to verify the feasibility of the PFC simulating models by comparing the stress-strain curves. At the same time, the details of the two aforementioned experiments will be fully presented, especially the numeric scouring model. Then, based on the PFC 2D models of which the mesoscopic parameters were calibrated, both the basic behaviors and the parameter effects of scouring action on the slope are going to be thoroughly analyzed, including the slope morphology, contact, and denudation amount. Finally, a discussion will be made on the mechanism of freeze-thaw on erosion damage of the residual soil slope induced by rainfall.

Experimental Scheme
The water content of samples is set at 21.5%, and the dry density is 1.49 g/cm 3 . The testing temperature range is −7-15°C. The FT cycles were set under 0, 2, 4, 6, and 8 times. The water supplement was conducted at the sealed condition. It is to be noted that the triaxial shear test was performed under the consolidated drainage condition. The temperature range was determined according to the field detection and the historical record.

Specimen Preparation and Freeze-Thaw Triaxial Shear Test Specimen Preparation
The specimens were prepared according to JTG E40-2007(Research Institute of Highway Science, 2010. It is noted that the soil was screened by a 2 mm sieve in order to maintain the homogeneity and the similarity of the initial soil structure. The mass percentage of each soil grain size (gradation) 1-2 mm, 0.5-1 mm, 0.25-0.5 mm, 0.075-0.25 mm, 0.03-0.075 mm, 0.075-0.03 mm, and <0.075 mm is 4. 45, 14.43, 12.66, 15.17, and 53.29, respectively.

Freeze-Thaw Triaxial Shear Test
The equipment utilized in this research is TCK-1, which was produced by Nanjing Soil Instrument Co., Ltd. To avoid misoperation, five specimens were prepared for each group. Three of them were used for shear tests under different confining pressures (100, 150, 200 kPa), and the other two were set as a standby.

Model Generation and Initial Setting of Calibration Modes Porosity Conversion
In PFC 2D simulation, the porosity can be converted from the practical soil model by an empirical equation (He et al., 2014). Wang believed that the parabolic modes seemed to be more appropriate after performing the comparison between several methods . By using the same approach in this research, shown in the following equation, the converted porosity is 18.65%: where p2D and p3D are the porosity in 2D and 3D, respectively.

Geometry
The dimension of the specimen is 32 mm × 15.64 mm, the soil particle radius is set between 0.65 and 1.3 mm, and the density is 2,655 kg/m 3 . The model consists of two pairs of rigid walls: the vertical panel that is applying the load in the vertical direction and the servo panel that is controlling pressure on the lateral side of the specimen.

Contact Models and Boundary Conditions
Contact models. Two models were utilized to consider the bonding and debonding conditions of the contacts between particles. For contacts that have not yet broken, it is widely accepted that the contact model is much more appropriate in simulating the clay material in this case. However, the opinions about which model is more appropriate are diverse, although the accuracy of the stress and strain curves from the default model is extensively acknowledged in the debonding case.
In most studies, the linear contact model is selected as the default model. However, it was found that increasing the interparticle friction cannot change the peak strength of the specimen, which was similar to that in the studies performed by Skinner (1969) and Oda et al. (1982). However, Suiker and Fleck (2004)  and Thornton (2000) believed that it could be explained by particle rolling, which is also supported in this paper. It has been acknowledged that the rotational velocity will be enormous while the material reaches its yield point. This shows a rolling effect on the surface, which undoubtedly reduces the friction coefficient. Thus, the established rotation resistance linear model (rrlinear) was selected as the default model (Iwashita and Oda, 1998;Iwashita and Oda, 2000;Jiang et al., 2005;Jiang et al., 2009;Liu et al., 2013). Boundary. The flexible boundary condition was achieved by setting the lateral wall stiffness 1/10 of the particle stiffness, while the vertical wall stiffness was the same as that of the particles.

Calibration of the Stress and Strain Curves
The stress-strain curves of the simulated and experimental results are shown in Figures 2A,B. It can be found that the correlation between these two curves is very high before the peak shear strength. Although both curves tend to decrease first and then remain unchanged after the peak value, the simulated curves seem to decrease rapidly at the initial stage after the peak strength (generally around 4% axial strain). And then, they remain basically unchanged, and the final residual strength is also lower than that in the testing curves. This is mainly due to the complicated nature of the size, shape, distribution, and stress characteristics of the practical soil particles, and the residual strength is provided by friction and structural occlusion force. Additionally, the particles are two-dimensional discs in the 2D simulation, where the point-type inter-particle contacts tended to slip under high force, which resulted in the attenuation of the residual strength.

Sample Generation
As shown in Figure 3, two slope models with a gradient of 45°and 60°were established. According to the study performed by Ni et al. (2000), changing the particle size appropriately has little effect on the simulation results. Besides, by integrating the operation ability of the computer and the running timing of the model, the radius of the soil particles was determined in a range from 0.1625 to 0.325 mm. Finally, 9,124 and 8,057 particle elements were generated in the 60°and 45°slopes, respectively.

Scouring Simulation
Generally, the simulation of equivalent rainfall is realized by changing the macroscopic and mesoscopic parameters of soil. Therefore, in order to simulate the equivalent rainfall scouring process more precisely, this study improves the method proposed by Ma (2012) in two aspects. On the one hand, the soil particles on the slope are given an initial velocity; on the other hand, the water absorption was considered by increasing the unit weight of soil particles. Specifically, the scouring process can be generally divided into the following four steps. Firstly, the generation and gravity balance of the soil particles are performed. Then, the allocation of the contact models is carried out based on the freeze-thaw zoning. Here, it was assumed that the width of the freeze-thaw zone (Zone B) accounts for 1/10 of the total width of the top of the slope. After that, the identification of the surface  particles is made by the programmed codes. Finally, the walls are going to be removed, and the velocity and the extra gravity will be imposed on the surface particles.

Scouring Velocity and Terminal Condition
Scouring velocity. When setting the scouring velocity, the particles will impact each other and bump away from the slope under high velocity; on the contrary, it does not work due to the cohesion between the particles. After repeated trials, a scouring velocity of 0.1 m/s was determined. In this study, the velocity is also variable to discuss the influence of flow velocity on the slope stability. Due to the change in direction of particle motion after collision, it is impossible to apply a constant flow velocity to the slope particles at present. Hence, only the influence of single scouring action on the slope stability will be discussed here.
Terminal condition. The original balance of superficial soil particles is broken under scouring, and the particles begin to move down the slope with water flow until a new balance is reached, where the movement of the particle is nearly stationary. After several trials, it is found that when the operation reaches 2.5 million cycles, the particles are stationary, and the maximum velocity is less than 1e-3 m/s. Therefore, the terminal condition was set when the total cycles reached 2.5 million times.

Simulation Scenarios and Parameters
The simulation scheme and main parameters are set, as shown in Table 1.

Basic Behaviors of Scouring Action
Morphology As shown in Figure 4, in order to distinguish easily, the slope is evenly divided into six layers by colors from top to bottom according to height. It can be found that only the shallow layer of the slope was destroyed after being scoured, while the other parts of the slope were not deformed. It is to be noted that the gray background represents the initial slope area, which does not participate in the calculation.
From Figures 4A-D, it can be found that the runoff only caused damage to the slope surface, leaving a little slide of soil particles, candle holes, and cavity prototypes at the top and foot of the slope. Besides, the particles at the slope foot are first washed and denuded and then accumulated on the slope feet, thereby making the soil in the upper layer in an unbalanced state to move downward. And then, under the continuous action of runoff shearing force, the shallow soil particles on the upper slope are gradually eroded, resulting in the particles in the deeper layer being disturbed and slid down along the structural surface with higher erosion resistance and accumulation at the foot of the slope, which finally can be called the spalling phenomenon.

Contacts
The contacts between particles can be illustrated in Figures 5A-D. As noted, the contacts in blue and red represent the initial and the default anti-rotation contact model, respectively. It can be found that most of the particles slide and pile up at the foot of the slope, which is in good agreement with the aforementioned slope morphology. Besides, the rotation resistance model is the major contact model of spalling fragments, which indicates that the spalling objects are non-sticky and loose soil granules. This is in line with the practical field investigation. However, part of them still obey the contact-bonding model, which is mainly due to the scouring effect that destroyed the soil structure at the weakest point of anti-scouring ability, avoiding the contact model there, which can be classified as fragmentary spalling.

Denudation Amount
The rotation angle of the particles is near zero when the soil structure is stable, while it increases when being spalled. Here, it is assumed that the contacts of soil particles have been denuded when the finial absolute rotation angle is greater than 0.5°. Therefore, by monitoring the particle rotation angle, the damage degree of the superficial slope can be quantified. It can be seen from Table 2 that, under runoff scouring, the amount of each layer is sorted as the foot (fifth layer) > top layer (first layer) > fourth layer > third layer > second layer, indicating that the most of the spalling occurred at the lower part of the slope.

Slope Gradient Morphology
The erosion failure process of the 45°slope before the FT effect is shown in Figures 6A-D. As noted, the scouring time of the 45°slope is shorter than that of the steeper one (60°slope). For example, the shape of the 45°slope changes slightly after 0.8 M steps, while the shape change of the 60°slope can still be observed at 2.5 M steps. This is mainly because the scouring force was counteracted by adhesion and friction between the particles, leaving the particles of the 45°slope to maintain a stable state more easily. This indicates that the slope gradient has a significant influence on stability under the action of runoff scouring, and the smaller slope has a stronger antiscouring ability.

Denudation Amount
The denudation amount of 45°and 60°slopes is shown in Table 2. As noted, it is evident that the amount of spalling particles was decreased significantly under the lower gradient slope, while the amount of each layer was the same as that of the steeper slope.

Freeze-Thaw Action Morphology
The erosion failure process of 45°and 60°slopes after the FT effect is shown in Figures 7A-H from a to d and e to h, respectively. Compared with Figures 4, 6, it can be found that, at the initial stage of erosion (within 0.8 M steps), the runoff only caused damage to the slope surface both with and without the FT effect. Moreover, in the middle and late stages (after 0.80 M steps), considering the FT effects, the particles cut into the lower layer of soil rather than rolling down under scouring action, resulting in a splitting effect on the shallow particles in the lower slope. The split blocks are dispersed into smaller blocks and particles under   the scouring effect and then piling up at the foot of the slope. The underlying reason for splitting was the poor cohesion among soil particles in the freeze-thaw zone. When the runoff was wrapping up the soil particles and moving downward, the soil in that zone was likely to be punched, and the phenomenon tended to occur.
Besides, in the final stage (above 2.5 M steps), the slope foot was damaged more significantly with the action of FT, which can be ascribed to the stress concentration at the foot of the slope. So, the particles deformed more largely under frequent FT disturbances. Finally, compared to the 60°slope, the impacts of runoff scouring on the 45°slope are quite limited on the top, regardless of the FT effect. During this process, the splitting from the upper layer soil to the lower does not occur.

Denudation Amount
The denudation of 45°and 60°slopes is shown in Table 2. As noted, it is evident that the number of spalling particles was increased significantly after considering the FT effect, which can be seen from the increased total number of 45°and 60°slopes by 147 and 156, rising 83.75 and 43.1%, respectively. Besides, the FT effect has a larger impact on the lower parts of the slope. For example, the denudation amount in the fourth and fifth layers of the slope considering the FT effect is much larger than that without.

Particle Movement
In order to detect the FT effect on the particle movement at different locations during the scouring process, three measurement spheres were selected on the 60°slope to record the position changes under different conditions, as shown in Figure 8, and the corresponding trajectories are presented in Figures 9A-C. It can be found that the soil particles on the lower layer of the slope surface reach the steady state first than those on the upper layer. This is mainly due to the shorter movement distance of the former particles, resulting in smaller potential energy, which is faster dissipated by the bonding and friction between the particles. Besides, when considering the FT effect, the spheres (5,247 and 649) located on the upper and middle layers reach a steady state faster than those without. The reason is that the adhesion between the particles of the shallow soil after the FT effect is reduced, leaving less kinetic energy needed to be consumed to destroy the contacts between the particles. Finally, the time of the sphere (5,529) located on a lower layer to reach the steady state is the same in both scenarios. However, when considering the FT effect, the particles move first downward and then upward, which indicates particles at the foot move more actively.

Denudation Amount
The denudation amount of each layer under the FT effect is shown in Table 3. It can be found that the denudation amount    increases significantly as the velocity is increased. With the increase in velocity, the denudation amount of particles in the top layer (the first layer) of the slope is unchanged, while that in the lower layer increases gradually. The reason for that is the flow was only loaded once rather than continuously, leaving short disturbance time to the top layers and producing shallow superficial spalling in the lower layer. Increasing the velocity resulted in a larger erosion load on the lower surface. So, the erosion amount increases gradually.

CONCLUSIONS
A PFC 2D scouring model, of which the mesoscopic parameters were calibrated by the indoor experiments, was developed to investigate the rainfall erosion damage of residual soil slope in the intermittently frozen area. The main conclusions are addressed in the following three aspects: • At the beginning of erosion, only the slope surface was damaged under both with and without the FT effect. And then, the splitting phenomenon was observed rather than rolling down in the middle and late stages. The most damaged part was at the foot, as was more evident under the action of FT in the final stage. Moreover, the FT effect caused less damage to the gentler slope.
• The denudation amount of spalling particles was increased significantly after considering the FT effect. Finally, regardless of the FT effect, the particles in the lower layer reached the stable state first, and the time was the same under these two conditions. • The pit size at the foot and the denudation amount were both increased under higher scouring velocity.
The future study will focus on the development of more realistic contact models and computing effective methods, which facilitate revealing the freeze-thaw effect on the slope surface. Besides, the coupling effects of water flow and soil particle movement on slope stability will also be explored.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, and further inquiries can be directed to the corresponding author.