Variations in Soil Erosion Resistance of Gully Head Along a 25-Year Revegetation Age on the Loess Plateau

: The e ﬀ ects of vegetation restoration on soil erosion resistance of gully head, along a revegetation age gradient, remain poorly understood. Hence, we collected undisturbed soil samples from a slope farmland and four grasslands with di ﬀ erent revegetation ages (3, 10, 18, 25 years) along gully heads. Then, these samples were used to obtain soil detachment rate of gully heads by the hydraulic ﬂume experiment under ﬁve unit width ﬂow discharges (2–6 m 3 h). The results revealed that soil properties were signiﬁcantly ameliorated and root density obviously increased in response to restoration age. Compared with farmland, soil detachment rate of revegetated gully heads decreased 35.5% to 66.5%, and the sensitivity of soil erosion of the gully heads to concentrated ﬂow decreased with revegetation age. The soil detachment rate of gully heads was signiﬁcantly related to the soil bulk density, soil disintegration rate, capillary porosity, saturated soil hydraulic conductivity, organic matter content and water stable aggregate. The roots of 0–0.5 and 0.5–1.0 mm had the highest beneﬁt in reducing soil loss of gully head. After revegetation, soil erodibility of gully heads decreased 31.0% to 78.6%, and critical shear stress was improved by 1.2 to 4.0 times. The soil erodibility and critical shear stress would reach a stable state after an 18-years revegetation age. These results allow us to better evaluate soil vulnerability of gully heads to concentrated ﬂow erosion and the e ﬃ ciency of revegetation. on the Loess Plateau and the location of sampling sites in Nanxiaohegou watershed. The r The revegetation the of of the gully to The D r was and disintegration and conductivity, the in the of the first choice for revegetation to restrain gully headcut erosion. Revegetation erodibility of gully heads by 31% to 78.6% and improve the critical by 1.2 4.0 This study allows us to better evaluate vulnerability of gully head to concentrated erosion during revegetation. Further studies the effect of the different combinations of different root architecture on


Introduction
Soil erosion is recognized as a global environmental problem, which severely damages infrastructure, causes land degradation and water pollution, and threatens the safety of human production and life [1][2][3]. In the past few decades, many scholars have made many efforts to study the process and mechanism of soil erosion, establish many soil erosion prediction models and try to control soil erosion [4][5][6][7]. At present, a set of soil erosion control measures system integrating engineering measures, agricultural measures, and biological measures has been formed [8][9][10], especially vegetation measures play an extremely important role in soil erosion control [11,12].
Previous studies have shown that revegetation can effectively reduce soil erosion. For example, Wang et al. [13] found that soil detachment capacity of abandoned farmland was 1.02 to 2.29 times greater than four restored lands. Li et al. [14] reported that the ratios of the soil detachment capacity of cropland to those of orchard, shrubland, woodland, grassland, and wasteland were 7.14, 12.29, 25.78, 28.45, and 46.43, respectively. The improvement of soil erosion resistance by revegetation is mainly controlled

Sampling Sites Selection
During our investigation of gully heads, some cracks developed near gully heads. Kompani-Zare et al. [39] and Guo et al. [21] stated that soil samples from 0 to 30 cm depth near the gully heads (the distance was less than 5 m) can represent soil properties of the gully heads. Therefore, in consideration of collapsibility and vertical joints development of loess, the sampling plots were established about 1.0 m meters from gully heads to ensure safety. As a result, four natural restoration grasslands with different ages (3, 10, 18, 25 years) were selected ( Figure 1). The natural restoration age was confirmed by consulting the village elders and scientists at the scientific experimental station. The slope aspect and gradient, elevation, soil type, and previous farming practices of the selected sites were similar to minimize the effects of these factors on the experimental results. For comparison, one corn-planted farmland site, with a topography similar to that of the grasslands, was selected as a control. The basic information of the five selected sites is listed in Table 1.

Sampling Sites Selection
During our investigation of gully heads, some cracks developed near gully heads. Kompani-Zare et al. [39] and Guo et al. [21] stated that soil samples from 0 to 30 cm depth near the gully heads (the distance was less than 5 m) can represent soil properties of the gully heads. Therefore, in consideration of collapsibility and vertical joints development of loess, the sampling plots were established about 1.0 m meters from gully heads to ensure safety. As a result, four natural restoration grasslands with different ages (3, 10, 18, 25 years) were selected ( Figure 1). The natural restoration age was confirmed by consulting the village elders and scientists at the scientific experimental station. The slope aspect and gradient, elevation, soil type, and previous farming practices of the selected sites were similar to minimize the effects of these factors on the experimental results. For comparison, one corn-planted farmland site, with a topography similar to that of the grasslands, was selected as a control. The basic information of the five selected sites is listed in Table 1.

Sampling and Measurement of Soil and Root
In this study, seven soil and root property parameters including soil bulk density, capillary porosity, soil disintegration rate, soil water-stable aggregate, soil saturated hydraulic conductivity, organic matter content, and root mass density were measured. Firstly, three repeated sampling plots (5 m × 5 m) were established in each of gully head sampling sites with the litter layer removed, and topsoil samples (0-30 cm) were collected. Then, three cutting rings (200 cm 3 ) were used to randomly collect soil samples in each plot, and a total of nine samples were oven-dried at 105 • C for 24 h to determine the soil bulk density of each gully head site. Similarly, the other 9 soil samples were also collected by cutting rings of 200 cm 3 to determine the soil saturated hydraulic conductivity by applying the constant water head test method. Three cutting rings (100 cm 3 ) were used to collect soil samples for the measurement of soil capillary porosity [33]. Three man-made steel cubical boxes (5 cm in length) were used to collect soil samples for measuring soil disintegration rate by using a disintegration box [14,40]. Lastly, the other three samples were randomly collected in each plot to form a mixed sample. A total of 45 mixed samples were obtained and used for laboratory analyses of organic matter content and water-stable aggregate and its stability. These mixed soil samples were air-dried at room temperature, with large roots and organic residues manually removed. Sieves with apertures (0.25, 0.5, 1.0, 2.5, and 5.0 mm) were used to test the water-stable aggregate. The potassium dichromate external heating method was used to measure the soil organic matter content.

Hydraulic Flume Experiments
A hydraulic flume experiment was conducted to determine the soil resistance to concentrated flow upstream gully heads ( Figure 2). The size of the flume was 2.0 m long and 0.15 m wide similar to the one used by De Baets et al. [28,31], which was enough to make water flow along the slope soil. An opening (0.5 m length and 0.1 m wide) was set at the bottom of the flume, and a metal sample box with the same size was used to collect undisturbed soil samples so that the surface of the soil sample was at the same level of the flume surface. The space between sampling box and flume edge sealed with painter's mastic to prevent boundary effects. According to the study of Guo et al. [40], the flume experiment was carried out under five different unit width flow discharges of 2, 3, 4, 5, 6 m 3 s −1 , and thus, a total of 100 samples (5 sites × 5 flow discharges × 4 replications) were collected to measure soil resistance of gully heads. To simulate real flow generation conditions, the soil should be saturated by using a watering pot before experiment. During the experiment, a portable flow meter instrument (LS300-A) with 1.5% accuracy was used to measure flow velocity which was regarded as the flow velocity scoring soil area. Runoff and sediment samples were collected with sampling tanks, and the sampling time was recorded. The measured flow velocity was modified according to flow regime [41]. Sampled sediment was oven-dried at 105 • C for 24 h to determine the soil loss amount (SLA, kg). Thus, the soil detachment rate (D r , kg m −2 s −1 ) could be calculated as follows: please check font size, please check all reference citation where SLA is the oven-dry mass of every sediment sample (kg), T is the experimental period (s) and A is the soil sample area (m 2 ). In addition, the relative soil detachment rate (RD r ) was calculated as the ratio between D r for the root-permeated soil samples and that for the farmland topsoil samples, tested at the same condition [28].
increase the dispersion of soil and then were placed on a 0.25 mm sieve and washed with tap water using low-pressure head. The living roots, plant debris and some pebbles were left on the sieve. Only the living roots were picked out carefully using tweezers one by one [15]. The washed roots were classified into 4 levels (0-0.5, 0.5-1.0, 1.0-2.0, and >2.0 mm) by vernier caliper and then were ovendried for 24 h at 65 °C and weighed to calculate root mass density (RBD, kg m −3 ).

Parameter Calculation after the Experiments
Soil particle is detached when flow shear stress exceeds the critical shear stress [6]. Soil erodibility parameter (Kr) and critical shear stress (c) were estimated for every natural restoration stage as the slope coefficient and intercept on the abscissa axis of the regression line between soil detachment rate and shear stress as described in the WEPP model as follow: Generally, soil detachment rate can be considered as zero when root reached the infinity. To quantify the relationship between detachment rate and root mass density, the Hill curve was selected to simulate the relationship between them [21,31,42]. The Hill curve is expressed as follows [43]: where RDr is relative soil detachment rate; Xr is root mass density; K, a and b are constants. The parameter a determines the shape of the curve, b determines the steepness of the curve and K is the asymptote of Dr for infinitesimal Xr values. Additionally, the Hill curve can be used to evaluate the ability of roots to increase soil resistance against concentrated flow erosion. According to Li et al. [44], b^(1/a) is plant specific and can be used as an index to compare the effectiveness of different plant roots in reducing soil erosion rates: the lower b^(1/a), the more effective the plant root. When the value of Xr is b^(1/a), the soil detachment rates is reduced by 50%.

Statistical Analysis and Plotting
The analysis of one-way ANOVA followed by multiple comparisons with LSD was applied to assess the differences of soil properties (Soil bulk density, soil capillary porosity, soil disintegration rate, saturated hydraulic conductivity, organic matter content, and water-stable aggregate) and root mass density among the five revegetation ages. All soil and root variables of each revegetation ages were tested whether the data exhibited a normal distribution and variance homogeneity by Shapiro- In addition, the flow depth (h, m) and shear stress (τ, Pa) were also calculated as follows: where h is flow depth (m), q is the flow discharge (m 3 s −1 ), w is the width of the flume (m), v is the mean flow velocity (m s −1 ). ρ is the water mass density (kg m −3 ), g is the gravity constant (m s −2 ), and S is the slope steepness (m m −1 ). After each scouring test, a steel cubical box (10 cm in length) was used to take soil sample in the center of soil area of scouring flume, and the sample was soaked in tap water for about one hour to increase the dispersion of soil and then were placed on a 0.25 mm sieve and washed with tap water using low-pressure head. The living roots, plant debris and some pebbles were left on the sieve. Only the living roots were picked out carefully using tweezers one by one [15]. The washed roots were classified into 4 levels (0-0.5, 0.5-1.0, 1.0-2.0, and >2.0 mm) by vernier caliper and then were oven-dried for 24 h at 65 • C and weighed to calculate root mass density (RBD, kg m −3 ).

Parameter Calculation after the Experiments
Soil particle is detached when flow shear stress exceeds the critical shear stress [6]. Soil erodibility parameter (K r ) and critical shear stress (τ c ) were estimated for every natural restoration stage as the slope coefficient and intercept on the abscissa axis of the regression line between soil detachment rate and shear stress as described in the WEPP model as follow: Generally, soil detachment rate can be considered as zero when root reached the infinity. To quantify the relationship between detachment rate and root mass density, the Hill curve was selected to simulate the relationship between them [21,31,42]. The Hill curve is expressed as follows [43]: where RD r is relative soil detachment rate; X r is root mass density; K, a and b are constants. The parameter a determines the shape of the curve, b determines the steepness of the curve and K is the asymptote of D r for infinitesimal X r values. Additionally, the Hill curve can be used to evaluate the ability of roots to increase soil resistance against concentrated flow erosion. According to Li et al. [44], bˆ(1/a) is plant specific and can be used as an index to compare the effectiveness of different plant roots in reducing soil erosion rates: the lower bˆ(1/a), the more effective the plant root. When the value of X r is bˆ(1/a), the soil detachment rates is reduced by 50%.

Statistical Analysis and Plotting
The analysis of one-way ANOVA followed by multiple comparisons with LSD was applied to assess the differences of soil properties (Soil bulk density, soil capillary porosity, soil disintegration rate, saturated hydraulic conductivity, organic matter content, and water-stable aggregate) and root mass density among the five revegetation ages. All soil and root variables of each revegetation ages were tested whether the data exhibited a normal distribution and variance homogeneity by Shapiro-Wilk test and Levene test, respectively. If the data failed to meet the two conditions, the Kruskal-Wallis test was performed for the above analysis. The interaction effect of flow discharge and revegetation age was detected using a two-way ANOVA. Pearson's correlation analysis was used to determine linkages among soil properties, root mass density, and soil detachment rate. Relationships among soil detachment rate, soil properties, flow shear stress and restoration age were analyzed by the regression method. The data analyses were conducted in SPSS v. 16.0 statistical software (IBM Corp., Armonk, NY, USA). The figure plotting was conducted by Origin v. 2020 (OriginLab Corp., Northampton, MA, USA). Figure 3 illustrates that the six soil properties of gully heads exhibited a significant increase or decrease with revegetation age. Compared with slope farmland, the soil bulk density (SBD) and soil disintegration rate (SDR) of revegetated gully heads significantly decreased by 5.7-18.6% and 28.8-80.5%, respectively (p < 0.05, Figure 3a,c), while the soil capillary porosity (SCP), saturated soil hydraulic conductivity (SHC), organic matter content (OMC) and water-stable aggregate (WSA) significantly increased by 3.9-13.8%, 17.4-236.2%, 34.2-221.8%, and 27.7-64.4%, respectively (p < 0.05, Figure 3b,d-f). Figure 4 illustrates that the roots of 1-2 mm in slope farmland had the relatively higher root mass density (RMD, 0.20 kg m −3 ) and accounted for 39% of total RMD. Notably, after revegetation, the RMD of >2.0 mm was significantly greater than those of the other three root diameters, and it can account for 40-61% of total RMD. When revegetation age was greater than 3 years, there was a significant difference in RMD among four root diameter levels (p < 0.05). In addition, we found that the RMD of four root diameters (except for >2.0 mm) showed a non-significant increase in the first three-years and then significantly increased.

Effect of Revegetation on Soil and Root Properties of Gully Heads
These results were similar to previous findings regarding the effects of revegetation on soil properties [45][46][47]. In fact, the improvement of soil properties of gully heads with revegetation age can be attributed to the accumulation of fresh plant residues in surface soil as well as roots and decomposed root residues in subsurface soil [48]. These materials can be directly transformed into soil organic matter and thus provide energy/carbon sources and nutrients for soil microorganisms [49,50], further promoting the development of soil aggregation and enhancing the cohesion of soil particles [51]. Hence, vegetation restoration would induce the formation of macroaggregates and increase the water stability of aggregates [19,52].
Water 2020, 12, 3301 7 of 16 root mass density (RMD, 0.20 kg m −3 ) and accounted for 39% of total RMD. Notably, after revegetation, the RMD of >2.0 mm was significantly greater than those of the other three root diameters, and it can account for 40-61% of total RMD. When revegetation age was greater than 3 years, there was a significant difference in RMD among four root diameter levels (p < 0.05). In addition, we found that the RMD of four root diameters (except for >2.0 mm) showed a nonsignificant increase in the first three-years and then significantly increased.   These results were similar to previous findings regarding the effects of revegetation on soil properties [45][46][47]. In fact, the improvement of soil properties of gully heads with revegetation age can be attributed to the accumulation of fresh plant residues in surface soil as well as roots and decomposed root residues in subsurface soil [48]. These materials can be directly transformed into soil organic matter and thus provide energy/carbon sources and nutrients for soil microorganisms

Effect of Revegetation Age on Soil Detachment of Gully Heads
As illustrated in Figure 5a, the D r of gully heads showed a significant decrease during the 25-year revegetation. This result was not agreed with the conclusion of Wang et al. [16] who stated that soil detachment capacity of sloped lands fluctuated with abandonment time, and the soil detachment capacity of the slope farmland was significantly greater than those of the abandoned farmlands. The difference was mainly attributed to the great difference in erosion environment (e.g., plant type, geomorphological feature, climate) significantly affecting the succession process [36,47]. The mean D r of slope farmland was 1.6 to 3.0 times greater than those of revegetated gully heads, which indicated that the revegetation played a role in enhancing the soil resistance of gully heads to concentrated flow erosion. of linear functions (y = p × x + q), and the p-value decreased with revegetation age, indicating that the sensitivity of Dr of the gully heads to concentrated flow erosion gradually decreased with increasing restoration age. Besides, the interacted effect of revegetation age and unit width discharge significantly affected Dr (p < 0.001) ( Table 2).   Figure 6 showed that Dr was positively correlated with soil bulk density and soil disintegration rate (p < 0.01), but negatively correlated with capillary porosity, saturated hydraulic conductivity, organic matter and water-stable aggregate of >0.25 mm (p < 0.01). Regression analysis showed that Dr increased with soil bulk density as a power function (Figure 7a), which showed an opposite trend with the Wang et al. [13] and Yu et al. [15]. Lower soil bulk density was caused by greater root physical and soil organisms' activities, and thus a soil with lower bulk density was harder to be detached. Additionally, Dr decreased with capillary porosity as a logarithmic function (Figure 7b), which was caused by physically binding and chemically bonding effect of root improving soil structure and porosity and hence increasing soil resistance to erosion [28,53]. Soil disintegration rate referred to the dispersion speed of soil contacting with water, which is an important factor determining soil resistance to erosion [14]. In this study, the soil disintegration rate decreased with the revegetation time ( Figure 3c) and Dr increased linearly with an increase in soil disintegration rate (Figure 7c). This is attributed to root wedging mechanism preventing soil from detaching that roots can bind soil and tie surface soil layer into strong and stable subsurface soil layer [14,54]. Saturated  Figure 5b shows the D r of gully heads of slope farmland and four restored grasslands varied with flow discharge. The optimal relationships between D r and flow discharge were fitted, which can reflect the response of D r of gully heads to concentrated flow induced by rainstorms of different recurrence intervals. It was found that the response of D r of slope farmland to flow discharge could be expressed as a power function (y = m × x n ), and the n-value was greater than 1, indicating the soil loss of gully heads increases at an increased speed with increasing flow. However, for the restored gully heads, the optimal relationships between D r and flow discharge could be described by a series of linear functions (y = p × x + q), and the p-value decreased with revegetation age, indicating that the sensitivity of D r of the gully heads to concentrated flow erosion gradually decreased with increasing restoration age. Besides, the interacted effect of revegetation age and unit width discharge significantly affected D r (p < 0.001) ( Table 2).  Figure 6 showed that D r was positively correlated with soil bulk density and soil disintegration rate (p < 0.01), but negatively correlated with capillary porosity, saturated hydraulic conductivity, organic matter and water-stable aggregate of >0.25 mm (p < 0.01). Regression analysis showed that D r increased with soil bulk density as a power function (Figure 7a), which showed an opposite trend with the Wang et al. [13] and Yu et al. [15]. Lower soil bulk density was caused by greater root physical and soil organisms' activities, and thus a soil with lower bulk density was harder to be detached. Additionally, D r decreased with capillary porosity as a logarithmic function (Figure 7b), which was caused by physically binding and chemically bonding effect of root improving soil structure and porosity and hence increasing soil resistance to erosion [28,53]. Soil disintegration rate referred to the dispersion speed of soil contacting with water, which is an important factor determining soil resistance to erosion [14]. In this study, the soil disintegration rate decreased with the revegetation time ( Figure 3c) and D r increased linearly with an increase in soil disintegration rate (Figure 7c). This is attributed to root wedging mechanism preventing soil from detaching that roots can bind soil and tie surface soil layer into strong and stable subsurface soil layer [14,54]. Saturated soil hydraulic conductivity is an integrating parameter for several physical characteristics such as bulk density, porosity, and mechanical composition. The conclusion that the D r decreased with increasing soil hydraulic conductivity by a power function is reasonable (Figure 7d) because this study and previous research findings have also indicated that changes in soil bulk density and porosity influenced soil detachment and also were affected by revegetation (e.g., Neves et al. [55]; Zhang et al. [56]). A negative power function was found between D r and soil organic matter content (Figure 7e). The accumulation of soil organic matter in soil could promote the formation of aggregate and enhance the cohesion of soil particles [57]. Hence, water-stable aggregate also was an indicator determining soil resistance to flow erosion [19]. The D r decreased as a linear function of water-stable aggregate of >0.25 mm (Figure 7f). The results were in agreement with the findings of Li et al. [14]. However, in Wang et al. [13,16] studies, no significant relationships were found between D r and organic matter and water-stable aggregate of >0.25 mm, probably caused by small variations of the two factors in their studies and difference in land use between their studies and this study (Podwojewski et al. [58]).

Response of Soil Detachment to Soil Properties
Water 2020, 12, x FOR PEER REVIEW 9 of 15 study and previous research findings have also indicated that changes in soil bulk density and porosity influenced soil detachment and also were affected by revegetation (e.g., Neves et al. [55]; Zhang et al. [56]). A negative power function was found between Dr and soil organic matter content (Figure 7e). The accumulation of soil organic matter in soil could promote the formation of aggregate and enhance the cohesion of soil particles [57]. Hence, water-stable aggregate also was an indicator determining soil resistance to flow erosion [19]. The Dr decreased as a linear function of water-stable aggregate of >0.25 mm (Figure 7f). The results were in agreement with the findings of Li et al. [14]. However, in Wang et al. [13,16] studies, no significant relationships were found between Dr and organic matter and water-stable aggregate of >0.25 mm, probably caused by small variations of the two factors in their studies and difference in land use between their studies and this study (Podwojewski et al. [58]).   . Correlation matrix among soil detachment rate, soil properties, and root mass density. Note: Dr, SBD, SCP, SDR, SHC, OMC, WAS, RMD < 0.5, RMD 0.5-1.0, RMD 1.0-2.0, and RMD > 2.0 refers to the soil detachment rate, soil bulk density, soil capillary porosity, soil disintegration rate, saturated soil hydraulic conductivity, organic matter content, water-stable aggregate, root mass density of <0.5 mm, root mass density of 0.5-1.0 mm, root mass density of 1.0-2.0 mm, and root mass density of >2.0 mm, respectively.

Response of Soil Loss of Gully-Head to Root Traits
Significant correlation was found between D r and RMD of 0-0.5 mm, 0.5-1.0 mm, 1.0-2.0 mm and >2.0 mm (p < 0.01, Figure 6), of which the RMD of 0-0.5 mm had the highest correlation with D r , indicating that roots of each diameter level had the significant impact on soil erosion of gully heads, especially the fibrous root of 0-0.5 mm. Furthermore, the Hill curve could well simulate the relationships between D r and RMD of different root diameters with R 2 varying from 0.42 to 0.57 ( Figure 8). As illustrated in Figure 8, the D r showed a rapid decrease when RMD of 0-0.5, 0.5-1.0, 1.0-2.0 and >2.0 mm ranged from 0 to 0.25 kg m −3 , 0 to 0.3 kg m −3 , 0 to 0.5 kg m −3 , and 0 to 1.0 kg m −3 , respectively, implying that soil erosion of gully heads could be controlled once vegetation restoration or root growth in soil. Although the roots were limited in density and flexible in early revegetation stage, whereas, roots can contribute to soil cohesion and additional strength, and be crucial in reduction of soil erosion [28,59]. Additionally, root system can bind soil and tie surface soil layer into strong and stable subsurface soil layer [54]. Well-developed root system had great physical binding and chemical bonding effect that could well bind soil particles and soil aggregates together and enhance soil resistance to erosion [16,28,42].
In addition, judged by fitted efficiency (R 2 ), the optimal results were found in RMD of 0-0.5 mm (Figure 8a), indicating that fibrous root of 0-0.5 mm is the optimal root system reducing soil loss of gully heads. However, Li et al. [44] reported that the ability of plant roots to decrease soil erosion mainly depended on the number of fibrous roots <1.0 mm. Shangguan et al. [53] also found a similar result but recommended root surface area density as the root variable. The reason may be that plant species with contrasting root architectures have a different erosion reduction effect [25]. Additionally, Amezketa [60] and Gyssels et al. [61] reported that monocotyledonous plants are superior to dicotyledonous plants and grasses are better than cereals in stabilizing soil aggregates.
According to Li et al. [44], bˆ(1/a) can be used as an indicator to compare the effectiveness of different diameter roots in reducing soil erosion. The lower bˆ(1/a), the more effective the diameter root. especially the fibrous root of 0-0.5 mm. Furthermore, the Hill curve could well simulate the relationships between Dr and RMD of different root diameters with R 2 varying from 0.42 to 0.57 ( Figure 8). As illustrated in Figure 8, the Dr showed a rapid decrease when RMD of 0-0.5, 0.5-1.0, 1.0-2.0 and >2.0 mm ranged from 0 to 0.25 kg m −3 , 0 to 0.3 kg m −3 , 0 to 0.5 kg m −3 , and 0 to 1.0 kg m −3 , respectively, implying that soil erosion of gully heads could be controlled once vegetation restoration or root growth in soil. Although the roots were limited in density and flexible in early revegetation stage, whereas, roots can contribute to soil cohesion and additional strength, and be crucial in reduction of soil erosion [28,59]. Additionally, root system can bind soil and tie surface soil layer into strong and stable subsurface soil layer [54]. Well-developed root system had great physical binding and chemical bonding effect that could well bind soil particles and soil aggregates together and enhance soil resistance to erosion [16,28,42]. In addition, judged by fitted efficiency (R 2 ), the optimal results were found in RMD of 0-0.5 mm (Figure 8a), indicating that fibrous root of 0-0.5 mm is the optimal root system reducing soil loss of gully heads. However, Li et al. [44] reported that the ability of plant roots to decrease soil erosion mainly depended on the number of fibrous roots <1.0 mm. Shangguan et al. [53] also found a similar result but recommended root surface area density as the root variable. The reason may be that plant species with contrasting root architectures have a different erosion reduction effect [25]. Additionally, Amezketa [60] and Gyssels et al. [61] reported that monocotyledonous plants are superior to dicotyledonous plants and grasses are better than cereals in stabilizing soil aggregates.
According to Li et al. [44], b^(1/a) can be used as an indicator to compare the effectiveness of different diameter roots in reducing soil erosion. The lower b^(1/a), the more effective the diameter

Effect of Revegetation Age on Soil Erosion Resistance of Gully Head
Rill soil erodibility parameter (K r ) and critical shear stress (τ c ) were employed to characterize the soil erosion resistance of gully heads [13,21], and were determined by the WEPP model (Equation (4)). The fitted linear function between D r and shear stress was illustrated in Figure 9. The slope of the fitted line is equal to the K r , and the K r of the restored grasslands were 31% to 78.6% less than that of slope farmland. In addition, we found that the soil erodibility of 3-year restored grassland rapidly declined by 31% compared with the slope farmland, indicating the short-term revegetation can rapidly reduce soil erodibility of gully heads. The K r of grasslands in this study ranged from 0.0009 to 0.0029 s m −1 , which were less than those reported by Li et al. [14]. Wang et al. [13] found that averaged K r of restored lands of abandoned farmland was 0.0024 s m −1 that was close to those of this study. The difference was mainly caused by differences in land use, plant species and restoration time. The soil samples were taken from different vegetation restoration models (korshinsk peashrub, black locust, Chinese pine and mixed forest of amorpha and Chinese pine) in the study of Wang et al. [13], and the restoration age (37 years) was greater than that of this study (3 to 25 years). Regression analysis found that the K r decreased with restoration time in an exponential function and showed a slight change when restored age was greater than 18 years (Figure 10). study. The difference was mainly caused by differences in land use, plant species and restoration time. The soil samples were taken from different vegetation restoration models (korshinsk peashrub, black locust, Chinese pine and mixed forest of amorpha and Chinese pine) in the study of Wang et al. [13], and the restoration age (37 years) was greater than that of this study (3 to 25 years). Regression analysis found that the Kr decreased with restoration time in an exponential function and showed a slight change when restored age was greater than 18 years ( Figure 10).  In addition, the c increased with restoration age by a power function ( Figure 10). However, the result was inconsistent with the finding of Wang et al. [16] that critical shear stress varied with restoration age in a nonlinear pattern, reaching the minimum at the restored age of 18. The difference in the temporal change of critical shear stress between Wang et al. [13] and this study was caused probably by differences in soil properties and vegetation characteristics of the sampling sites.
Compared with slope farmland, the c of the grasslands was improved by 1.2 to 4.0 times, while c of restored land had a little decrease when restored time was more than 18 years ( Figure 10). The result further indicated that revegetation can effectively improve the soil erosion resistance of gully head to concentrated flow, and the critical shear stress would reach a stable state after a 18-year revegetation.

Conclusions
This study was carried out to explore the effect of revegetation age on soil erosion resistance of gully heads in the gully region of the Loess Plateau. The results showed that revegetation significantly improved soil properties and promoted root accumulation of gully heads. The mean Dr In addition, the τ c increased with restoration age by a power function ( Figure 10). However, the result was inconsistent with the finding of Wang et al. [16] that critical shear stress varied with restoration age in a nonlinear pattern, reaching the minimum at the restored age of 18. The difference in the temporal change of critical shear stress between Wang et al. [13] and this study was caused probably by differences in soil properties and vegetation characteristics of the sampling sites. Compared with slope farmland, the τ c of the grasslands was improved by 1.2 to 4.0 times, while τ c of restored land had a little decrease when restored time was more than 18 years ( Figure 10). The result further indicated that revegetation can effectively improve the soil erosion resistance of gully head to concentrated flow, and the critical shear stress would reach a stable state after a 18-year revegetation.

Conclusions
This study was carried out to explore the effect of revegetation age on soil erosion resistance of gully heads in the gully region of the Loess Plateau. The results showed that revegetation significantly improved soil properties and promoted root accumulation of gully heads. The mean D r of slope farmland was 1.6 to 3.0 times greater than those of revegetated gully heads. The revegetation can effectively weaken the sensitivity of soil erosion of the gully heads to concentrated flow. The D r of gully heads was positively related to bulk density and disintegration rate and negatively related to soil capillary porosity, saturated soil hydraulic conductivity, organic matter content, and water-stable aggregate. Roots of 0-0.5 mm and 0.5-1.0 mm were the most effective roots in reducing soil erosion of gully head, and the native plant species with rich root of 0.5-1.0 mm and 0-0.5 mm were recommended as the first choice for revegetation to restrain gully headcut erosion. Revegetation can reduce soil erodibility of gully heads by 31% to 78.6% and improve the critical shear stress by 1.2 to 4.0 times. This study allows us to better evaluate soil vulnerability of gully head to concentrated flow erosion during revegetation. Further studies are needed to quantify the effect of the different combinations of vegetation types with different root architecture types on soil erosion resistance of gully heads.