Erosion Process and Temporal Variations in the Soil Surface Roughness of Spoil Heaps under Multi-Day Rainfall Simulation

The extensive artificially accelerated erosion of spoil heaps on newly engineered landforms is a key ecological management point requiring better understanding. Soil surface roughness is a crucial factor influencing erosion processes; however, study on spoil heap erosion with a view of surface roughness is lacking. This study investigated the erosion processes and the spatiotemporal variation of surface roughness on spoil heaps, and then, analyzed how the roughness affected the hydrological and sediment yield characteristics. Sequences of four artificial rainstorms with constant rainfall intensity (90 mm/h) were applied to cone-shaped spoil heaps (ground radius 3.5 m, height 2.3 m) of a loess soil containing 30 mass percent rock fragments. The surface elevation was sampled by a laser scanner. For the surface roughness indicators, the root mean square height (rmsh) and the correlation length (cl) increased sharply during the first rainfall event, and in the last three rainfall events, rmsh increased slightly and cl showed a relative decrease. The initial rmsh/cl of the whole slope surface ranged from 0.063 to 0.135, and increased with the rainfall sequence, thus, indicating that the spoil heap surface became rougher. Increasing soil roughness in the rainfall sequence delayed the initial runoff time and increased the runoff yield. The average runoff coefficient of the spoil heaps was 0.658. The average erosion rate of each rainfall event can be simulated by a regression equation of the corresponding average runoff rate and median cl (R-square of 0.816). Soil slumping with an average volume of 0.014 m3 occurred in the first two rainfall events, thus, significantly changing the roughness and peak instant erosion rate. Together, the results revealed the effects of surface roughness on the erosion of spoil heaps and would provide a useful reference for soil loss prediction and control.


Introduction
The accelerated erosion of artificially disturbed land surfaces has received widespread attention in recent years [1][2][3][4]. The soil erosion rates of areas disturbed by construction activity have been shown to be 2 to 40,000 times greater than those in pre-construction conditions [5]. Most studies on erosion have concentrated on industrial, mining, and road construction sites [6][7][8][9][10][11], whereas fewer studies have investigated the erosion of spoil heaps during infrastructure construction. Spoil heaps consist of excavated earth materials that contain rock fragments and have been piled up on disturbed land surfaces [12]. Their steep slopes, loose structures, and low levels of vegetation cover lead to the much more severe soil erosion compared to that of undisturbed surfaces or tilled surfaces [11][12][13].

Experimental Materials and Equipment
The indoor simulated rainfall experiment was carried out in the Simulated Rainfall Hall at the State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau. The rainfall height was 18 m, with an adjustable rainfall intensity ranging from 30 to 350 mm/h [37] and a rainfall uniformity of over 80%.
A homemade platform was used to simulate the conical spoil heap (Figure 1). According to the rainfall characteristics of multiple years on the Loess Plateau [38,39], the rainfall intensity and duration were set to be 90 mm/h and 40 min, respectively. A total of four consecutive rainfalls with intervals of 24 h were performed, and three replicates were conducted.
The area of the steel experimental platform was 3.5 × 3.5 m, and a 2.5 m high right-triangle steel baffle was placed on two sides of the platform. Inclined guiding flumes with slopes of 5 • were provided on the other two sides of the platform, which converged at the water outlet. To facilitate the discharge of runoff and sediment, the bottom of the platform was also on an incline at 3 • towards the water outlet. To reduce the marginal effect of the interface between the steel baffle and the soil after loading, gauze was pasted on the inside of the two upright steel baffles. The soil used in the experiment was collected from the spoil ground of construction projects around Yangling, Shaanxi, and its mechanical composition was as follows: 27.6% clay, 48.2% silt, and 23.3% sand. According to the classification of the IUSS Working Group WRB (2015) [40], the soil belongs to the urbic technosol group. Rock fragments are a common component in spoil heaps. The rock fragments used in this study were collected from the natural landslide massif near a highway in Zhouzhi County, Yangling, Shaanxi. The rock fragments are siliceous limestones with diameters from 1 to 15 cm. Based on the composition of local landslide rock fragments, they were classified into four categories according to the particle diameter, including 1-3.5, 3.5-7, 7-10, and 10-15 cm. The groups were mixed in a ratio of 2:3:3:2, respectively.
Based on the experimental platform, the height of the spoil heap was set to 2.2 m, and the slope was based on the natural angle of repose of the soil. Before stacking, the soil was screened with a 10 mm sieve to remove debris, and the soil moisture was maintained at 12-14%. The resulting soil was then mixed with the pre-prepared rock fragments. To restore the process of spoil heap formation, the stacking materials were transported to 2.8 m above the platform by a conveyor belt and were allowed to fall freely to simulate the dumping and stacking process under the action of gravity. After stacking, the surface was flattened with a board and the final spoil heap was in the shape of a 1/4 conical body ( Figure 1) with a topsoil bulk density of 1.290 ± 0.011 g/cm 3 , a bottom radius of 3.5 m, a height of 2.2 m, a slope length of 4.2 m, and a slope (i.e., the angle of natural repose) of 32° to 35° (62.5-70.0%). The morphological parameters of the simulated spoil heap were consistent with the field data of Zhao et al. [41], in which the slope length and a slope of 80% of the spoil heaps were in the range of 2.5-5.5 m and 33° ± 2°, respectively.

Experimental Procedure
To ensure the consistency of the initial surface in three replicates, 10 min of pre-watering with an intensity of 30 mm/h was performed after stacking, following which the body was allowed to sit for 24 h under a plastic cover. A ring cutter was used to measure the soil bulk density before each rainfall on the upper, middle, and lower parts of the slope. Before the experiment, a rain gauge was placed on the upper, middle, and lower parts of the slope to measure the rainfall per unit time per unit area. The error between the measured and designated rainfall intensity was controlled within 5% to ensure the uniformity of rainfall. The time it took for runoff to be observed at the water outlet was recorded, and runoff sediment samples were collected every minute for the first 4 minutes, and then, every 2 minutes thereafter. The samples were then dried and weighed. After each simulated rainfall, the platform was covered with plastic for 24 h before the next simulation.

The Acquisition of Surface Elevation
The slope was scanned after every simulated rainfall using a Leica MS50 laser scanner to obtain the point-cloud data of the surface elevation. Five point-cloud data files were collected from each The soil used in the experiment was collected from the spoil ground of construction projects around Yangling, Shaanxi, and its mechanical composition was as follows: 27.6% clay, 48.2% silt, and 23.3% sand. According to the classification of the IUSS Working Group WRB (2015) [40], the soil belongs to the urbic technosol group. Rock fragments are a common component in spoil heaps. The rock fragments used in this study were collected from the natural landslide massif near a highway in Zhouzhi County, Yangling, Shaanxi. The rock fragments are siliceous limestones with diameters from 1 to 15 cm. Based on the composition of local landslide rock fragments, they were classified into four categories according to the particle diameter, including 1-3.5, 3.5-7, 7-10, and 10-15 cm. The groups were mixed in a ratio of 2:3:3:2, respectively.
Based on the experimental platform, the height of the spoil heap was set to 2.2 m, and the slope was based on the natural angle of repose of the soil. Before stacking, the soil was screened with a 10 mm sieve to remove debris, and the soil moisture was maintained at 12-14%. The resulting soil was then mixed with the pre-prepared rock fragments. To restore the process of spoil heap formation, the stacking materials were transported to 2.8 m above the platform by a conveyor belt and were allowed to fall freely to simulate the dumping and stacking process under the action of gravity. After stacking, the surface was flattened with a board and the final spoil heap was in the shape of a 1/4 conical body ( Figure 1) with a topsoil bulk density of 1.290 ± 0.011 g/cm 3 , a bottom radius of 3.5 m, a height of 2.2 m, a slope length of 4.2 m, and a slope (i.e., the angle of natural repose) of 32 • to 35 • (62.5-70.0%). The morphological parameters of the simulated spoil heap were consistent with the field data of Zhao et al. [41], in which the slope length and a slope of 80% of the spoil heaps were in the range of 2.5-5.5 m and 33 • ± 2 • , respectively.

Experimental Procedure
To ensure the consistency of the initial surface in three replicates, 10 min of pre-watering with an intensity of 30 mm/h was performed after stacking, following which the body was allowed to sit for 24 h under a plastic cover. A ring cutter was used to measure the soil bulk density before each rainfall on the upper, middle, and lower parts of the slope. Before the experiment, a rain gauge was placed on the upper, middle, and lower parts of the slope to measure the rainfall per unit time per unit area. The error between the measured and designated rainfall intensity was controlled within 5% to ensure the uniformity of rainfall. The time it took for runoff to be observed at the water outlet was recorded, and runoff sediment samples were collected every minute for the first 4 min, and then, every 2 min thereafter. The samples were then dried and weighed. After each simulated rainfall, the platform was covered with plastic for 24 h before the next simulation.

The Acquisition of Surface Elevation
The slope was scanned after every simulated rainfall using a Leica MS50 laser scanner to obtain the point-cloud data of the surface elevation. Five point-cloud data files were collected from each complete experiment. The vertical and horizontal errors of the scanner were less than 1 mm if the scanning distance was less than 5 m. The scanner was set up 3-4 m vertically away from the slope (Figure 1), and three calibration points were marked to determine the relative position of the spoil heap body in the coordinate system before the first scan. Subsequent scans would then initiate only if the error of the three points was less than 3 . This procedure was equivalent to the calibration step in computer data processing [42] and ensured that the five point-cloud data files were in the same coordinate system. The slope area was~11.54 m 2 , the scanning point spacing was set at 5 mm, and~240,000 elevation points were obtained during each scan.

Data Processing of The Point Cloud
In ArcGIS (version 10.4, environmental systems research institute, USA) software, the point cloud was generated into grid files (digital elevation model, DEM) with a resolution of 5 × 5 mm using TIN files as a transition. Five DEMs (A1-A5) were obtained for each complete experiment.
The surface of the spoil heap was a curved surface with a steep slope, and the variation in the surface elevation caused by the curvature and the slope was not caused by micro-topographical changes, which would not be analyzed in this study. To eliminate the effect of curvature and slope on the analysis of surface elevation, A0 was created based on A1. The surface of A0 exhibited the general shape (slope, slope length, and curvature) of the spoil heap, but it did not depict the slight difference in elevation caused by the surface micro-topography. Thus, A0 was A1 with a smooth surface. The difference between A1-A5 and A0, and DEM (A01-A05) were then determined. In this way, the obstruction of slope and curvature on the analysis of surface micro-topography was eliminated by removing the contribution of slope and curvature in elevation. The A0 DEM was created as follows: (1) The A1-dependent DEM with a resolution of 100 × 100 mm was generated using the "re-sampling" tool. (2) Using the "re-sampling" tool again, the DEM obtained from (1) was converted to a DEM with a resolution of 5 × 5 mm. The interpolation method was bilinear interpolation [11]. Finally, the DEMs of A01-A05 were converted to the point files, and the point coordinates were used to evaluate the micro-topography factor. The calculation was carried out for the upper, middle, and lower parts of the slope.

Calculation of Soil Surface Roughness Parameters
The root mean square height (rmsh) and the correlation length (cl) were extracted from the correlation function computed from the measured soil surfaces [43]. The parameters rmsh, cl, and rmsh/cl can characterize the surface features in vertical, horizontal and combined directions [44]. rmsh is one of the random roughness parameters proposed by Allmaras et al. [45]; for the discrete one-dimensional case, the rmsh is given by: where N is the number of height records, z i is the height corresponding to record i, and z is the mean height of all records. The correlation length (cl) represents the horizontal component of roughness, describing the relative locations of heights or the way in which the heights vary along the surface [46]. Once the horizontal distance between two points is greater than cl, their heights may be considered to be statistically independent of one another. cl was calculated from exponential fits of the normalized autocorrelation function [47]: where ρ(h) is the autocorrelation function, which represents the correlation existing between height z of point i (z i ) and that of another point located at a lag distance h (z i+h ), and N(h) is the number of pairs considered in each lag h. cl is defined as the critical length where the heights of two points are regarded as independent, i.e., ρ(h) is equal to 1/e, so that ρ(cl) = 1/e. rmsh and cl explain the roughness in one direction (vertical or horizontal), but a combined roughness parameter is needed to interpret the emissivity of a given rough soil surface. Zribi and Dechambre [47] have proposed parameter z s as a combination of rmsh and cl, thus, accounting for both vertical and horizontal roughness components:

Runoff Characteristics
The runoff is related to the infiltration amount and stronger soil infiltration capacity leads to less runoff yield on the slope. In each rainfall event, the runoff rate increased and then, stabilized with time ( Figure 2). Lower initial runoff rates resulted from a lower water content and a greater infiltration rate of the underlying surface. The infiltration rate gradually decreased with increasing topsoil moisture, and the runoff yield increased correspondingly. When the infiltration rate increased to a stable value, the runoff rate tended to stabilize. The times when the runoff rates stabilized were 30 and 20 min for the first and second rainfall events, respectively, whereas 10 min was required to achieve runoff for the third and fourth events (dotted lines in Figure 2). The shorter increasing period and longer stabilizing period with the rainfall sequence were related to the initial moisture content, which were 18.3%, 22.3%, 25.4%, and 26.5%, respectively, in the first to the fourth rainfall events. Higher initial soil water content resulted in a shorter time to reach a stable flow rate on spoil heaps.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 16 pairs considered in each lag h. cl is defined as the critical length where the heights of two points are regarded as independent, i.e., ρ(h) is equal to 1/e, so that ρ(cl) = 1/e. rmsh and cl explain the roughness in one direction (vertical or horizontal), but a combined roughness parameter is needed to interpret the emissivity of a given rough soil surface. Zribi and Dechambre [47] have proposed parameter zs as a combination of rmsh and cl, thus, accounting for both vertical and horizontal roughness components:

Runoff Characteristics
The runoff is related to the infiltration amount and stronger soil infiltration capacity leads to less runoff yield on the slope. In each rainfall event, the runoff rate increased and then, stabilized with time ( Figure 2). Lower initial runoff rates resulted from a lower water content and a greater infiltration rate of the underlying surface. The infiltration rate gradually decreased with increasing topsoil moisture, and the runoff yield increased correspondingly. When the infiltration rate increased to a stable value, the runoff rate tended to stabilize. The times when the runoff rates stabilized were 30 and 20 min for the first and second rainfall events, respectively, whereas 10 min was required to achieve runoff for the third and fourth events (dotted lines in Figure 2). The shorter increasing period and longer stabilizing period with the rainfall sequence were related to the initial moisture content, which were 18.3%, 22.3%, 25.4%, and 26.5%, respectively, in the first to the fourth rainfall events. Higher initial soil water content resulted in a shorter time to reach a stable flow rate on spoil heaps. As shown in Figure 3, the runoff yield increased with the rainfall sequence, but little differences existed in the third and the fourth rainfall events because of the nearly constant infiltration capacity of the underlying surface. The runoff yield in the fourth rainfall event increased by 29.0% compared with that in the first rainfall event. The runoff coefficient is the ratio of the runoff yield to rainfall during the same period, which increased from 0.526-0.580 in the first rainfall event to 0.696-0.717 in the fourth rainfall event, with a mean value of 0.658. As shown in Figure 3, the runoff yield increased with the rainfall sequence, but little differences existed in the third and the fourth rainfall events because of the nearly constant infiltration capacity of the underlying surface. The runoff yield in the fourth rainfall event increased by 29.0% compared with that in the first rainfall event. The runoff coefficient is the ratio of the runoff yield to rainfall during the same period, which increased from 0.526-0.580 in the first rainfall event to 0.696-0.717 in the fourth rainfall event, with a mean value of 0.658.

Soil Erosion Characteristics
As shown in Figure 4, the soil loss amount of each rainfall event first increased, then decreased with the rainfall sequence. The total soil loss was 17.319, 19.282, and 18.829 kg, respectively, for each repetition, and the maximum soil loss occurred in the third rainfall event for Repetition 1 (4.829 kg) and in the second rainfall event for Repetitions 2 (5.477 kg) and 3 (5.177 kg). The evolution of the erosion rates is shown in Figure 5. During the first two rainfall events, the erosion rate fluctuated dramatically because of the gravitational soil slump ( Figure 6). The topsoil was softened after absorbed water and the shear strength was reduced. Together with the steepness of the slope and the loose structure, the spoil heap slopes had slumping proneness. The peak points of the changing lines of erosion rate ( Figure 5) indicated the time nodes when the soil slump occurred. The erosion rate changed relatively smoothly during the last two rainfall events, showing a trend of an initial increase followed by stabilization after the slumped soil had been eroded.

Soil Erosion Characteristics
As shown in Figure 4, the soil loss amount of each rainfall event first increased, then decreased with the rainfall sequence. The total soil loss was 17.319, 19.282, and 18.829 kg, respectively, for each repetition, and the maximum soil loss occurred in the third rainfall event for Repetition 1 (4.829 kg) and in the second rainfall event for Repetitions 2 (5.477 kg) and 3 (5.177 kg).

Soil Erosion Characteristics
As shown in Figure 4, the soil loss amount of each rainfall event first increased, then decreased with the rainfall sequence. The total soil loss was 17.319, 19.282, and 18.829 kg, respectively, for each repetition, and the maximum soil loss occurred in the third rainfall event for Repetition 1 (4.829 kg) and in the second rainfall event for Repetitions 2 (5.477 kg) and 3 (5.177 kg). The evolution of the erosion rates is shown in Figure 5. During the first two rainfall events, the erosion rate fluctuated dramatically because of the gravitational soil slump ( Figure 6). The topsoil was softened after absorbed water and the shear strength was reduced. Together with the steepness of the slope and the loose structure, the spoil heap slopes had slumping proneness. The peak points of the changing lines of erosion rate ( Figure 5) indicated the time nodes when the soil slump occurred. The erosion rate changed relatively smoothly during the last two rainfall events, showing a trend of an initial increase followed by stabilization after the slumped soil had been eroded. The evolution of the erosion rates is shown in Figure 5. During the first two rainfall events, the erosion rate fluctuated dramatically because of the gravitational soil slump ( Figure 6). The topsoil was softened after absorbed water and the shear strength was reduced. Together with the steepness of the slope and the loose structure, the spoil heap slopes had slumping proneness. The peak points of the changing lines of erosion rate ( Figure 5) indicated the time nodes when the soil slump occurred. The erosion rate changed relatively smoothly during the last two rainfall events, showing a trend of an initial increase followed by stabilization after the slumped soil had been eroded. Soil slump occurred on the upslope position in the three repeated treatments, with a mean slump volume of 0.014 m 3 (Table 1); the corresponding DEM models are displayed in Figure 6. Soil shrinks when it is gradually wetted and the collapse cracks developed on the upper slope [48,49]. The runoff enters the interior of the cracks and exerts a downward thrust on the soil, thus, causing a soil slump to occur on the upper slope.  Figure 7 shows that rmsh first increased, and then, stabilized with cumulative rainfall (CR) on different parts of the slope. The rmsh of initial surfaces before rainfall was 0.209-0.271 cm, with little difference in the three slope parts. However, the final rmsh differed significantly among various slope Soil slump occurred on the upslope position in the three repeated treatments, with a mean slump volume of 0.014 m 3 (Table 1); the corresponding DEM models are displayed in Figure 6. Soil shrinks when it is gradually wetted and the collapse cracks developed on the upper slope [48,49]. The runoff enters the interior of the cracks and exerts a downward thrust on the soil, thus, causing a soil slump to occur on the upper slope.  Figure 7 shows that rmsh first increased, and then, stabilized with cumulative rainfall (CR) on different parts of the slope. The rmsh of initial surfaces before rainfall was 0.209-0.271 cm, with little difference in the three slope parts. However, the final rmsh differed significantly among various slope parts (* P < 0.01; 1.713 cm for the upper slope, 1.384 cm for the middle slope, and 0.669 cm for the lower slope). The greatest rmsh change for the upper slope reflected the dramatic change in elevation as a result of the soil slump. The smallest increase in rmsh for the lower slope indicated the least change in topography.

Spatiotemporal Variation of The Roughness Parameter with Rainfall Accumulation
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 16 parts (*P < 0.01; 1.713 cm for the upper slope, 1.384 cm for the middle slope, and 0.669 cm for the lower slope). The greatest rmsh change for the upper slope reflected the dramatic change in elevation as a result of the soil slump. The smallest increase in rmsh for the lower slope indicated the least change in topography. The first rainfall event had the strongest effect on surface roughness, and rmsh increased the most during the first rainfall event (by 1.318, 0.990, and 0.298 cm for the three slope parts, respectively). The three subsequent rainfall events induced only slight increases in rmsh of 0.186, 0.151, and 0.101 cm for the upper, middle, and lower slopes, respectively. The relationship between rmsh and CR for the different slope parts can be fitted by the exponential function y=a−bc x (a=1.710/1.350/0.636, b=1.501/1.107/0.365, c=0.967/0.966/0.974). The R 2 values were all > 0.9, thus, implying that CR is good to indicate the change in rmsh for soil surfaces of spoil heaps.  The changes in cl show an initial increase followed by a slight decrease for the three slope positions (Figure 8). The initial flat surface corresponded to the shortest cl because of the random variability of heights values of each point record, and no significant difference existed in three slope parts (being 3.142, 4.262, and 4.271 cm for the three slope parts, respectively). The cl value had increments of 7.023, 8.199, and 3.336 cm in the first rainfall event for the three slopes parts, respectively. The first rainfall led to the development of erosional features based on the initial surface, on which the height values of each point are randomly distributed with little correlation between each other. In the first rainfall, the surface self-adjusted to primeval erosion terrain under the effect The first rainfall event had the strongest effect on surface roughness, and rmsh increased the most during the first rainfall event (by 1.318, 0.990, and 0.298 cm for the three slope parts, respectively). The three subsequent rainfall events induced only slight increases in rmsh of 0.186, 0.151, and 0.101 cm for the upper, middle, and lower slopes, respectively. The relationship between rmsh and CR for the different slope parts can be fitted by the exponential function y = a − bc x (a = 1.710/1.350/0.636, b = 1.501/1.107/0.365, c = 0.967/0.966/0.974). The R 2 values were all >0.9, thus, implying that CR is good to indicate the change in rmsh for soil surfaces of spoil heaps.
The changes in cl show an initial increase followed by a slight decrease for the three slope positions (Figure 8). The initial flat surface corresponded to the shortest cl because of the random variability of heights values of each point record, and no significant difference existed in three slope parts (being 3.142, 4.262, and 4.271 cm for the three slope parts, respectively). The cl value had increments of 7.023, 8.199, and 3.336 cm in the first rainfall event for the three slopes parts, respectively. The first rainfall led to the development of erosional features based on the initial surface, on which the height values of each point are randomly distributed with little correlation between each other. In the first rainfall, the surface self-adjusted to primeval erosion terrain under the effect of runoff and erosion and the correlation between the elevation points was enhanced.
Erosion continued to develop based on the primeval erosion terrain in the subsequent rainfall events. The soil surfaces became increasingly discontinuous because of the development of potholes and rills. Consequently, the rate of elevation change increased, and cl decreased. During the subsequent three rainfalls, cl decreased by 10.0%, 17.8%, and 21.3% for the three slope parts, respectively. However, they were still greater than the initial cl values, which indicated that once the surface was shaped by the runoff and erosion processes, the correlation between surface elevation points was greater than the original surface. Figure 7. Change of rmsh with cumulative rainfall (CR) for different parts of the slope. The changes in cl show an initial increase followed by a slight decrease for the three slope positions (Figure 8). The initial flat surface corresponded to the shortest cl because of the random variability of heights values of each point record, and no significant difference existed in three slope parts (being 3.142, 4.262, and 4.271 cm for the three slope parts, respectively). The cl value had increments of 7.023, 8.199, and 3.336 cm in the first rainfall event for the three slopes parts, respectively. The first rainfall led to the development of erosional features based on the initial surface, on which the height values of each point are randomly distributed with little correlation between each other. In the first rainfall, the surface self-adjusted to primeval erosion terrain under the effect of runoff and erosion and the correlation between the elevation points was enhanced. A significant difference existed in cl values after rainfall between the three slope parts (middle slope>upper slope>lower slope). The lower slope was the flattest with less microtopography change in the rainfall sequence, with a smaller cl than those of the other two slope parts. This may be attributed to the widening flow width of the lower slope that shallowed the runoff, resulting in a decreased erosive force and a slighter erosion degree. Although rainfall shaped the upper and middle slopes more obviously, the cl of the upper slope was smaller than that of the middle slope, because the soil slump led to a more broken surfaces and severe changes in elevation.
rmsh/cl indicated the change in surface roughness, combining the horizontal and vertical characteristics. It increased with CR for the three slope parts (Figure 9), and the average rmsh/cl of the whole slope increased from 0.063 to 0.145, verifying the increasing roughness due to the soil erosion process. The initial rmsh/cl was 0.068, 0.057 and 0.0635 for the three slope parts, respectively, with no significant difference. After rainfall, the rmsh/cl differed significantly among the three slopes, occurring in the order upper > middle > lower slope. The decrease in the roughness from the top to the bottom was associated with the erosion degree. The difference in rmsh/cl between the upper and middle slope was 0.047 on average, which was twice that between middle and lower slope. The significant larger rmsh/cl values of upper slope were mainly because of soil slumping ( Figure 6).
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 16 Erosion continued to develop based on the primeval erosion terrain in the subsequent rainfall events. The soil surfaces became increasingly discontinuous because of the development of potholes and rills. Consequently, the rate of elevation change increased, and cl decreased. During the subsequent three rainfalls, cl decreased by 10.0%, 17.8%, and 21.3% for the three slope parts, respectively. However, they were still greater than the initial cl values, which indicated that once the surface was shaped by the runoff and erosion processes, the correlation between surface elevation points was greater than the original surface.
A significant difference existed in cl values after rainfall between the three slope parts (middle slope＞upper slope＞lower slope). The lower slope was the flattest with less microtopography change in the rainfall sequence, with a smaller cl than those of the other two slope parts. This may be attributed to the widening flow width of the lower slope that shallowed the runoff, resulting in a decreased erosive force and a slighter erosion degree. Although rainfall shaped the upper and middle slopes more obviously, the cl of the upper slope was smaller than that of the middle slope, because the soil slump led to a more broken surfaces and severe changes in elevation.
rmsh/cl indicated the change in surface roughness, combining the horizontal and vertical characteristics. It increased with CR for the three slope parts (Figure 9), and the average rmsh/cl of the whole slope increased from 0.063 to 0.145, verifying the increasing roughness due to the soil erosion process. The initial rmsh/cl was 0.068, 0.057 and 0.0635 for the three slope parts, respectively, with no significant difference. After rainfall, the rmsh/cl differed significantly among the three slopes, occurring in the order upper > middle > lower slope. The decrease in the roughness from the top to the bottom was associated with the erosion degree. The difference in rmsh/cl between the upper and middle slope was 0.047 on average, which was twice that between middle and lower slope. The significant larger rmsh/cl values of upper slope were mainly because of soil slumping ( Figure 6).

Relationships between Runoff and Soil Roughness
The runoff initial time (T0) is a comprehensive reflection of the soil texture, slope, and surface roughness [50]. As shown in Figure 10, T0 decreased rapidly, and then, tended to be stable with the rainfall sequence. The exponential equation y=37.238−1497.635*0.046 x (R 2 =0.9) described this

Relationships between Runoff and Soil Roughness
The runoff initial time (T 0 ) is a comprehensive reflection of the soil texture, slope, and surface roughness [50]. As shown in Figure 10, T 0 decreased rapidly, and then, tended to be stable with the rainfall sequence. The exponential equation y = 37.238 − 1497.635*0.046 x (R 2 = 0.9) described this relationship.

Relationships between Sediment Yield Characteristics and Soil Roughness
Runoff is the exogenic force of soil erosion. In this study, the average runoff rate and erosion rate were significantly positively correlated, with a partial correlation coefficient of 0.772 (Table 2). Excluding the dominant role of runoff, the effect of roughness on soil loss seemed to be relatively small. Only the cl in three indicators correlated with the erosion rate positively, among which, the correlation coefficient of the median cl was larger, being 0.701. The median cl represented the average of the initial and final cl of a rainfall event, taking into account the change in cl during rainfall. Regression analysis of the average erosion rate, runoff rate, and initial and median cl (cli, clm), is listed in Table 3. According to the test parameters of the two fitting equations, the combination of runoff rate and clm more reliably predicted the erosion rate, with a satisfactory R 2 of 0.816. The standardized regression coefficients of R and clm were 0.585 and 0.473, respectively, showing the runoff influenced the soil loss rates more than the roughness. Table 3. Regression analysis results of average runoff rate (R, mm/min), median cl (clm , cm), initial cl (cli , cm) and average erosion rate (E, g/min m −2 ). Correlation analysis of the initial soil moisture, soil roughness parameters (average of the entire slope), T 0 , and average runoff rate is shown in Table 2. The initial soil moisture was significantly correlated with both T 0 and runoff rate, with correlation coefficients >0.9. Greater topsoil moisture was associated with a shorter T 0 and greater runoff rate. Therefore, surface roughness delayed runoff initiation and the correlation coefficient of the indicator rmsh/cl was the largest, being −0.873.
Surface roughness still influenced the runoff rate with continuing rainfall, showing that greater roughness led to more runoff outflow. Three initial roughness indicators all had significant correlation coefficients from 0.770 to 0.883, and the median rmsh/cl also provided a verification of this relationship with a coefficient of 0.788 (* P < 0.01).

Relationships between Sediment Yield Characteristics and Soil Roughness
Runoff is the exogenic force of soil erosion. In this study, the average runoff rate and erosion rate were significantly positively correlated, with a partial correlation coefficient of 0.772 (Table 2). Excluding the dominant role of runoff, the effect of roughness on soil loss seemed to be relatively small. Only the cl in three indicators correlated with the erosion rate positively, among which, the correlation coefficient of the median cl was larger, being 0.701. The median cl represented the average of the initial and final cl of a rainfall event, taking into account the change in cl during rainfall. Regression analysis of the average erosion rate, runoff rate, and initial and median cl (cl i , cl m ), is listed in Table 3. According to the test parameters of the two fitting equations, the combination of runoff rate and cl m more reliably predicted the erosion rate, with a satisfactory R 2 of 0.816. The standardized regression coefficients of R and cl m were 0.585 and 0.473, respectively, showing the runoff influenced the soil loss rates more than the roughness. Table 3. Regression analysis results of average runoff rate (R, mm/min), median cl (cl m , cm), initial cl (cl i , cm) and average erosion rate (E, g/min m −2 ).

Discussion
Our results indicate that the runoff coefficients of spoil heaps are about 0.658, however, the runoff coefficients of the undisturbed or tillage surfaces containing rock fragments are generally >0.7-0.8 [51,52], thus, indicating a stronger infiltration capacity of spoil heaps. For spoil heaps, the loose soil structure plus rock fragment fraction led to the development of continuous preferential flow paths or macropore flows, which improved the solution infiltration rates [53].
In addition to the development of cracks on the upper slope, the greater infiltration capacity of spoil heaps is also one of the reasons why soil slump is prone to happen. During the first two rainfalls, the rain sank into the soil that made the topsoil soften thus, reducing the shear strength. Together with the steepness of the slope and the loose structure, the spoil heap slopes had slumping proneness. On the other hand, the safety factors for rocky soil slopes are lower than those without rock fragments, because rock fragments obstruct groundwater flow along the slope; therefore, the groundwater levels must be raised to provide a sufficient gradient to facilitate flow, thereby increasing the pore water pressure and decreasing the slope stability [54,55]. Soil slump in our results echoed the opinion of Peng et al. [14] that rock fragments affected the erosion types and styles of spoil heaps.
Our results show an increasing rmsh and a decreasing cl with CR, which is in contrast to those for tillage surfaces [47,56,57]. Generally, tilled surfaces have larger initial rmsh values due to the artificial soil ridges or excavation. At the same time, the artificial land preparation measures made the initial cl of farmland surface smaller than that of the primeval erosion terrain (after the first rainfall) of spoil heap [21,56,58]. As the cumulative rainfall increased, the soil at bulge was eroded and transported to low-lying places for sedimentation due to the relatively gentle slope. Hence, the undulating soil surface gradually becomes gentler and less rough under the flattening effect of runoff, which is associated with a decreasing rmsh and increasing cl values [56,58]. However, the initial surface of the spoil heap was a naturally smooth surface after dumping, and resulted in a smaller rmsh and larger cl values. Furthermore, the steep slope promoted the eroded soil to move out of the slope rather than sedimentate, therefore, the spoil heap surface became increasingly uneven, and rills developed [59]. As a result, rmsh gradually increased and cl decreased with CR [60].
Increasing soil roughness in the rainfall sequence delayed the initial runoff time. This finding is in agreement with the results reported by Luo et al. [36], Vermang et al. [22] and Zhao et al. [58]. Rough surface retained more water in depression storage; hence, water had a longer residence time on the surface and was more likely to infiltrate [61]. However, more runoff yield was gained on rougher surfaces, which was observed in our study. It is different from the conclusions in previous studies reporting a decrease in total runoff with increasing roughness [28], although the effect on final runoff rate is more variable with both a decreasing [22] and an increasing effect of increasing roughness [18,61]. A higher degree of surface roughness can help flow concentration on the slope [18,62]. On the spoil heaps surface, due to the upper soil slump, rills developed maturely along the lower edge of the slump.
Hence the increasing surface roughness indicated a higher developed drainage network. Furthermore, the steeper slope of spoil heaps made the overland flow tend not to be stored in potholes after the formation of the drainage network [17]. Therefore, an increasing runoff yield with rougher surface was observed.
The close correlation between erosion rate and cl m in our results was also different from some previous studies establishing a relationship between sediment yield and some vertical parameters (rmsh, RR) on tillage land [25,29,63]. For spoil heaps, the changes in horizontal parameter may be closer to the changes in the erosion process, such as slope slumps, which reflect the peculiarity of the soil erosion process of the spoil heaps. In the present study, increasing cl m corresponded with a larger erosion rate. cl m can implicitly indicate the accessibility of the erodible soil to some extent. In the first rainfall event, cl m increased. Correspondingly, the initial smooth surface became uneven because of slump, and the slump body was easily erodible; or, the surface self-adjusted to primeval erosion terrain based on the erosion of surface floating soil, which is also a preferential eroded source. In response, average erosion rate increased. As cl m gradually decreased in the subsequent rainfalls, a more fragmentized surface developed due to flow path incision and the exposure of rock fragments [63,64]. The supply of detachable sediment was reduced, hence, the erosion rate decreased.
Limitations also existed in our study. Soil slump was important for the changes in erosion process and micro-topography for spoil heaps. The cracks on the surface that induced the soil slump were observed, but how deep the cracks were and how rainwater entered the cracks and seeped along the slope to form a sliding zone cannot be revealed. Future work should involve the evolution of the surface topography, the infiltration route and the change of mechanics characteristics in soil matrix, under different rainfall conditions. By means of combination monitoring of internal and external aspects, the critical conditions for soil slump can be explored. The findings will make a significant contribution to the risk assessment of spoil heaps and the accuracy of soil loss prediction will also be improved.

Conclusions
A laboratory study of the effects of multi-day simulated rainfall on spoil heaps revealed the evolution of soil surface roughness and the erosion processes; the influence of varying soil surface roughness on runoff and soil loss was then analyzed.
The average runoff rate of sub-rainfall increased until the third rainfall event and then, became relatively stable. The average erosion rates first increased and then, decreased slightly with the rainfall sequence. The peak instant erosion rate was associated with soil slumping on the upslope in the first two rainfall events, leading to significant changes in the roughness of the upslope.
The first rainfall induced the greatest change in soil roughness. The rmsh and cl increased sharply during the first rainfall event, whereas rmsh increased only slightly in the last three rainfalls, and a relative decrease was observed for cl. The rmsh/cl, describing the roughness characteristics combining the vertical and horizontal features, increased with the rainfall sequence, thus, indicating that the spoil heap surface was becoming rougher.
A larger soil roughness with the rainfall sequence delayed the initial runoff time and increased the runoff yield. The average erosion rate of each rainfall event can be simulated by a regression equation of the corresponding average runoff rate and median cl (R-square of 0.816).
Our results firstly revealed the interaction between surface roughness and soil erosion on spoil heap. The obtain relationship between soil loss and runoff, roughness helps to predict the soil loss more conveniently from the topographic characteristics, which provides a reference to policy formulation and planning during infrastructure construction.