Soil erodibility and its influencing factors on the Loess Plateau of China: a case study in the Ansai watershed

The objectives of this work were to identify the best possible method to estimate soil erodibility (K) and understand the influencing factors of soil erodibility. In this study, 151 soil samples were collected during soil surveys in the Ansai watershed of the Loess Plateau of China. The K values were estimated by five methods: erosion-productivity impact model (EPIC), nomograph equation (NOMO), modified nomograph equation (M-NOMO), Torri model and Shirazi model. The main conclusions of this paper are (1) K values in the Ansai watershed ranged between 0.009 and 0.092 t · hm2 · h/(MJ ·mm · hm2), and the maximum values were 1.9–7.3 times larger than the corresponding minimum values, and the Shirazi and Torri models were considered the optimal models for the Ansai watershed. (2) Different land use types had different levels of importance; the principal components (PCs) accounted for 100 % (native grassland), 48.88 % (sea buckthorn), 62.05 % (Caragana korshinskii), and 53.61 % (pasture grassland) of the variance in soil erodibility. (3) The correlations between soil erodibility and the selected environmental variables differed among different vegetation types. For native grasslands, soil erodibility had significant correlations with terrain factors. For the most artificially managed vegetation types (e.g., apple orchards) and artificially restored vegetation types (e.g., sea buckthorn), soil erodibility had significant correlations with the growing conditions of vegetation. Soil erodibility had indirect relationships with both environmental factors (e.g., elevation and slope) and human activities, which potentially altered soil erodibility.


Introduction
Soil erodibility (K), one of the key factors of soil erosion (Igwe, 2003;Fu et al., 2005;Ferreira et al., 2015), is defined as the susceptibility of soil to erosional processes (Bagarello et al., 2012;Bryan et al., 1989).It has been extensively used in both theoretical and practical approaches to measure soil erosion.However, it is a complex concept affected by many factors, including soil properties (Chen et al., 2013;Wang et al., 2015;Manmohan et al., 2012), terrain (Wang et al., 2012;Mwaniki et al., 2015;Parajuli et al., 2015), climate (Hussein et al., 2013;Sanchis et al., 2012), vegetation (Sepúlveda-Lozada et al., 2009), and land use (Cerdà et al., 1998;Tang et al., 2016).To calculate soil erodibility, many strategies have been researched to understand soil erodibility, including measurements of physical and chemical soil properties, instrumental measurements, mathematical models, and graphical methods (Wei et al., 2017a).Although the direct measurement of soil erosion in large plots under natural rainfall over long periods can provide accurate estimates of soil erodibility, this is a time-consuming and costly method (Bonilla et al., 2012;Vaezi et al., 2016a, b).Therefore, mathematical models are more commonly used to estimate soil erodibility.Some of the most common estimation models are the nomogram model (NOMO) and the modified nomogram model (M-NOMO), which were established by Wischmeier et al. (1971Wischmeier et al. ( , 1978)); the erosion-productivity impact model (EPIC), which was developed by Williams et al. (1990); the best nonlinear fitting formula using the physical and chemical properties of the soil, which was developed by Torri et al. (1997); and the estimation model that uses the average size of the soil geometry developed by Shirazi et al. (1988).Each estimation method differs in terms of applicability, even within the same area, because the different estimation methods include different physical and chemical soil properties (Lin et al., 2017;Wang et al., 2013b;Kiani et al., 2016).Consequently, the estimated results can significantly differ among methods because soil conditions vary by region (Lin et al., 2017;Wang et al., 2013b).Selecting the optimal estimation method of soil erodibility is therefore critical to estimate the amount of soil erosion.
Soil erosion on the Loess Plateau of China is among the highest in the world (Fu et al., 2009;Huang et al., 2016).The area affected by soil and water loss is as large as 4.5 × 10 5 km 2 (∼ 71 % of the local land area), and the longterm average sediment loss is up to 1.6 × 10 9 t (Fu et al., 2017).To maintain water quality and control soil erosion (Fu et al., 2011), the Chinese government has implemented a large-scale policy to convert farmlands to forests and grasslands since the 20th century (Lü et al., 2012;X. M. Feng et al., 2013;Wu et al., 2016).Although the large-scale introduction of vegetation is expected to have reduced soil erosion, the extent of the reduction remains unclear.Therefore, different estimation methods should be used to calculate erosion factors, including the soil erodibility factor.In this study, the Ansai watershed of the Loess Plateau of China was chosen as a case study, and the five abovementioned estimation methods of estimating K value were applied.The objectives of this study were (1) to estimate the soil erodibility factor with different methods, (2) to select the optional method to estimate K, and (3) to understand the influencing factors of soil erodibility for the local area.

Study area
The Ansai watershed (between 36 • 30 45 -37 • 19 3 N and 108 • 5 44 -109 • 26 18 E) is located around the upper reaches of the Yanhe River, in the inland hinterland of the northwestern Loess Plateau.This watershed lies in the northern part of Shanxi Province and borders the Ordos basin.It belongs to the typical loess hilly and gully region and covers an area of approximately 1334 km 2 .The soil type in the study area is loess soil, with low fertility and high vulnerability to erosion (Zhao et al., 2012;Yu et al., 2015).The topography is complex and varied, and the land surface is fragmented into different land uses, dominated by rain-fed farmland, apple orchard, native grassland, pasture grassland, shrubland, and forest (Q.Feng et al., 2013).The elevations within the watershed are high in the northwest and low in the southeast, ranging between 997 and 1731 m above sea level.The watershed belongs to the mid-temperate continental semiarid monsoon climate region.The average annual precipitation is 505.3 mm, and 74 % of the rainfall occurs from June to September.

Sample point setting
The soil data used in this study came from 151 typical sample data sets that were obtained during soil surveys conducted from July to September 2014.The soil type of all 151 sample points is loess soil.Representative vegetation types were selected: (1) natural vegetation: native grasslands (NG); (2) artificially managed vegetation types: apple orchards (AO) and farmland (FL); and (3) artificially restored vegetation types: pasture grasslands (PG), sea buckthorn (SB), Caragana korshinskii (CK), David's peach (DP), and black locust (BL).The distance between each vegetation site sampled was at least 2 km, and the size of each vegetation type was greater than 30 m by 30 m.The selected sample plots were evenly distributed within the study area.The sample plots within the farmland and grassland had a size of 2 m × 2 m, whereas the corresponding dimensions for the sample plots within the shrubland and forest areas were 5 m × 5 m and 10 m × 10 m, respectively.Each sample plot was replicated three times.The locations of the sampling points were determined using a GPS unit (Garmin eTrex 309X, Garmin Ltd. subsidiary in Shanghai, China).The collected soil samples were taken to the laboratory, dried naturally, ground, and sieved with a 2 mm sieve.The soil particle size distributions of the soil samples were evaluated using the hydrometer method.The size classes of soil particles in this study were based on USDA classes and were as follows: sand (0.005-2.0 mm), silt (0.002-0.05 mm), and clay (< 0.002 mm) (Wang et al., 2012).
To fully explore the primary factors influencing soil erodibility in the Ansai watershed, we chose four types of environmental factors: physicochemical soil properties, topographic factors, climate factors, and vegetation factors.Although soil erodibility does not directly depend on environmental factors, soil properties such as soil particle size distribution and soil organic matter can be affected by environmental factors; thus, environmental factors have indirect relationships with soil erodibility.These environmental factors covered 20 independent variables: elevation (Ele), slope position (SP), slope aspect (SA), slope gradient (SG), slope shape (SS), clay content (Cla), silt content (Sil), sand content (San), organic matter (OM) content, soil bulk density (SBD), porosity (Por), average annual rainfall (AAR), vegetation coverage (VC), aboveground biomass (AB), vegetation height (VH), litter biomass (LB), plant density (PD), crown width (Cro), basal diameter (BD), and branch number (BN).All of the data on environmental factors were derived from the field surveys.The main characteristics and sampling numbers for the study area are shown in Table 1, and the sampling points are shown in Fig. 1.Based on the results of the Spearman correlation analysis, we retained some environmental variables that displayed significant correlations (P < 0.05) with soil erodibil-ity to perform a principal component analysis (PCA) and obtain the minimum data set (MDS) (Xu et al., 2008).Only those principal components (PCs) with eigenvalues N > 1.0 and only those variables with highly weighted factor loadings (i.e., those with absolute values within 10 % of the highest value) were retained for the MDS (Mandal et al., 2008).

Research methods
Soil erodibility indicates the degree of difficulty with which soil becomes separated, eroded, and transported by rainfall erosivity (Wang et al., 2013a;Cerdà et al., 2017).The soil erodibility factor, which is commonly known as the K factor in models, is defined as the average rate of soil loss per unit of rainfall erosivity index from a cultivated continuous fallow plot that is 22.1 m long with 9 % slope in the universal soil loss equation (Zhang et al., 2008).To minimize bias from any single estimation method, we estimated the K values using five estimation models (i.e., EPIC, NOMO, M-NOMO, Torri, and Shirazi), which have been widely applied in research on soil erodibility (Wischmeier et al., 1971(Wischmeier et al., , 1978;;Williams et al., 1990;Torri et al., 1997;Shirazi et al., 1988).

K-value estimation using the EPIC model
The erosion-productivity impact model (EPIC) developed by Williams et al. (1990) is as follows: where SAN is percent sand content, SIL is percent silt content, CLA is percent clay content, C is percent organic carbon content, and SN 1 = 1 − SAN/100.The resulting K value is reported in United States customary units of short ton 2.3.2K-value estimation using the NOMO model Wischmeier et al. (1971) proposed this model after analyzing the relationships between soil erosion and five soil characteristic indicators,: percent silt + very fine sand fraction (0.05-0.1 mm), percent sand fraction, soil organic matter content, a code for soil structure, and a code for soil permeability: where M is the product of the percent of silt + very fine sand and the percent of all soil fractions other than clay, OM is soil organic matter content (%), S is soil structure code, and P is soil permeability code.The resulting K value is reported in United States customary units of short ton

K-value estimation using the M-NOMO model
On the basis of the universal soil loss equation (USLE) model, the RUSLE model was modified for calculating soil erodibility; the revised nomograph equation was modified from the previous nomograph equation (Wischmeier et al., 1978).The revised nomograph equation is as follows: where M is the product of the percent of silt + very fine sand and the percent of all soil fractions other than clay, OM is soil organic matter content (%), S is soil structure code, and P is soil permeability code.The resulting K value is reported in United States customary units of short ton

K-value estimation using the Torri model
Torri et al. (1997) established this model in 1997 using data describing soil particle size and soil organic matter content.
The model has few parameters and simple data acquisition.
The formula used for this model is as follows: where OM and c are percent soil organic matter and clay content, respectively.D g can be calculated by using the following formula: where D g is the Napierian logarithm of the geometric mean of the particle size distribution, d i (mm) is the maximum diameter of the ith class, d i−1 (mm) is the minimum diameter and f i is the mass fraction of the corresponding particle size class.We calculated D g based on three particle-size classes: sand, silt, and clay.The resulting K values are reported in the international units of (t 2.3.5 K-value estimation using the Shirazi model   Meanwhile, D g in this model can be calculated by using the following formula: where f i is the weight percentage of the ith particle size fraction (%), m i is the arithmetic mean of the particle size limits for the ith fraction (mm), and n is the number of particle size fractions.The resulting K value is reported in United States customary units of short ton • ac • h/(100 ft • short ton • ac • in).

K-value comparisons
To increase the comparability of the results from the different estimation models, our research adopted the international units for the K values, t • hm 2 • hr/(MJ • mm • hm 2 ).The in-ternational K value is equal to the K value reported in the United States customary units multiplied by 0.1317.To clarify the form of the distribution, we collected the frequency distribution figures of soil erodibility for each model (Wei et al., 2017a, b).The K values obtained using the five methods were normally distributed (P > 0.05).Therefore, the soil erodibility K values measured within the study area were statistically analyzed directly, without the need for data conversion (Fang et al., 2016).To discuss the best possible texturebased method to estimate K, related research on K estimation, especially that involving measured values of K on the Loess Plateau of China, was consulted.A Taylor diagram was also used to compare the models.

Soil erodibility in the Ansai watershed based on five different models
We obtained different values when calculating descriptive statistics of the K value in the Ansai watershed among the different models (Table 2).The range of K values based on the five methods were between 0.032 and 0.060, 0.046 and 0.092, 0.047 and 0.088, 0.009 and 0.066, and 0.018 and 0.044 t • hm 2 • h/(MJ • mm • hm 2 ) for K EPIC , K NOMO , K M-NOMO , K Torri , and K Shirazi , respectively.The maximum values were 1.875, 2.000, 1.872, 7.333, and 2.444 times larger than the corresponding minimum values (Table 2).The differences between the mean and median values were 0.001, −0.001, 0.000, 0.000, and 0.000 t • hm 2 • h/(MJ • mm • hm 2 ) for K EPIC , K NOMO , K M-NOMO , K Torri , and K Shirazi , respectively.The standard deviations (SDs) of the K values were 0.408, −0.447, −1.079, −2.639, and 0.059 for K EPIC , K NOMO , K M-NOMO , K Torri , and K Shirazi , respectively.The skewness values of the K values were 0.946, 0.956, 4.353, 16.872, and 0.009 for K EPIC , K NOMO , K M-NOMO , K Torri , and K Shirazi , respectively.The Cv value of K M-NOMO was 0.067 < 10 %, and the Cv values of K EPIC , K NOMO , K Torri , and K Shirazi were 0.109, 0.110, 0.113, and 0.182, respectively, all of which corresponded to between 10 % and 100 %.
In the Taylor diagrams (Taylor, 2001; Fig. 2), the K values based on the EPIC model were used as the reference objects.The K values based on the Torri, NOMO, and Shirazi models were similar and located close to each other.In contrast, the K values estimated by the M-NOMO and EPIC models were inconsistent with the other K values.

Spearman correlation coefficients of soil erodibility and environmental variables in the Ansai watershed
The correlations between soil erodibility and the environmental variables varied among the different vegetation types (Tables S1-S4 in the Supplement).In general, soil erodibility in artificially managed vegetation types (apple orchards and David's peach) and artificially restored vegetation types (e.g., sea buckthorn and black locust) had significant correlations with vegetation properties.For example, soil erodibility in areas planted with apple orchards had a significant positive correlation with plant density (P < 0.05, Table S1).Soil erodibility in areas with sea buckthorn had significant negative correlations with slope gradient and plant density and significant positive correlations with average annual rainfall and aboveground biomass (P < 0.05, Table S3).Soil erodibility of areas with David's peach had significant positive correlation with aboveground biomass and significant negative correlations with slope gradient, vegetation coverage, vegetation height, crown width, and basal diameter (P < 0.05, Table S4).Soil erodibility in areas with black locust had significant negative correlation with elevation and significant positive correlations with slope position, slope gradient, soil bulk density, vegetation coverage, litter biomass, and branch number (P < 0.05, Table S4).Soil erodibility in areas under other vegetation types, such as grassland or farmland, was more strongly correlated with soil or landscape properties than other impact factors.The results of the analyses of correlations between estimated K values and the selected environmental variables showed that soil erodibility in farmlands had significant positive correlations with slope shape and a significant negative correlation with slope gradient (P < 0.05, Table S1).The soil erodibility of areas with native grasslands had a significant negative correlation with elevation and significant positive correlations with average annual rainfall and slope gradient (P < 0.05, Table S2).The soil erodibility of areas with pasture grasslands did not have significant correlations with environmental variables other than soil organic matter content and soil particle size (P < 0.05, Table S2).
Soil erodibility in areas with Caragana korshinskii had a sig-   nificant positive correlation with elevation and a significant negative correlation with average annual rainfall (P < 0.05, Table S3).

Principal component analysis of soil erodibility under different vegetation types
The PCA identified one PC each for apple orchards, native grasslands, sea buckthorn, Caragana korshinskii, and pasture grasslands, which accounted for 100 %, 48.88 %, 62.05 %, and 53.61 % of the variances, respectively (Table S5).For apple orchards, plant density was the primary contributor to the high factor loading.For native grasslands, PC1 included two variables that had highly weighted factor loadings: the slope gradient and elevation.Pasture grasslands had no variables with high factor loadings because it had no significant environmental variables except soil particle size and soil organic matter.The highly weighted factor loadings in areas with sea buckthorn were slope gradient, aboveground biomass, and plant density.In areas planted with Caragana korshinskii, two variables had high factor loadings: average annual rainfall and elevation (Table S5).
The PCA identified two PCs each for farmland and David's peach; the corresponding cumulative variances were 73.93 % and 81.07 %, respectively.The PC1 for farmland included two variables that had high factor loadings: slope shape and slope position; whereas PC2 only included slope gradient.In areas planted with David's peach, crown width, vegetation height, and vegetation coverage contributed to the high factor loading of PC1; whereas basal diameter alone had a high factor loading for PC2.In areas planted with black locust, the PCA identified three PCs that accounted for 70.25 % of the variance (Table S5).PC1 had slope position, elevation, and litter biomass as parameters with high factor loadings.The parameters with high factor loadings for PC2 were slope gradient and soil bulk density, and vegetation coverage had a high factor loading for PC3 (Table S5).
The MDS of soil erodibility included six environmental variables for black locust, four for David's peach, three each for farmland and sea buckthorn, two each for native grasslands and Caragana korshinskii, one for apple orchards, and none for pasture grasslands (Tables S1, S2, and S3).In addition to soil organic matter and soil particle size, which were included in the K-value estimation equations, the dominant factors affecting soil erodibility for farmland were slope shape, slope gradient, and slope position.For apple orchards, the only dominant factor affecting soil erodibility (other than soil organic matter and soil particle size) was plant density.For areas with native grasslands, the dominant factors affecting soil erodibility were soil organic matter, soil particle size, slope gradient, and elevation.For areas with sea buckthorn, the dominant factors affecting soil erodibility were aboveground biomass, slope gradient, and plant density in addition to the two soil properties.The dominant factors affecting soil erodibility in areas with Caragana korshinskii were soil particle size, soil organic matter, average annual rainfall, and elevation.For areas with black locust, the dominant factors were slope gradient, slope position, elevation, litter biomass, soil bulk density, and vegetation coverage in addition to soil organic matter and soil particle size.The dominant factors affecting soil erodibility in areas with David's peach included soil organic matter, soil particle size, crown width, vegetation height, and vegetation coverage.In this study, we found that different models resulted in different estimates of soil erodibility (Table 2).Since the different estimation methods use different soil attributes as input parameters, the coefficient of variation of the same input parameters will differ.For example, the EPIC model focuses on the features of the soil particles and soil nutrients, whereas the NOMO model focuses not only on soil particle size and soil nutrient characteristics but also on the soil structural characteristics, such as soil structure code and soil permeability code.The existing soil erodibility estimation equations are used to calculate soil erodibility based on data of physicochemical soil properties, such as soil texture, soil structure, soil permeability, and soil organic matter content (Wischmeier et al., 1971(Wischmeier et al., , 1978;;Williams et al., 1990;Torri et al., 1997;Shirazi et al., 1988).Among these factors, the main physical soil property is soil particle composition, such as the contents of sand, silt, and clay, and the main chemical soil property is soil organic matter content (Wei et al., 2017).
Our results showed that the K values based on the Torri, NOMO, and Shirazi models were located close to each other in the Taylor diagrams (Fig. 2) and that these three models could therefore represent soil erodibility in the Ansai watershed.Based on previous studies, these models have been recommended as the optimal models for China's subtropical zone, China's purple hilly region, Northeast China, and China's Loess Plateau (Table 4).However, we suggest that the Torri and Shirazi models are the best models based on estimated K values derived from these models and actual (measured) soil erodibility data from the Ansai watershed (Zhang et al., 2001; Table S6).The estimated K values based on the Torri and Shirazi models were closer to the measured soil erodibility data among those of the three possible appropriate models (Tables 2 and S6).Our findings are supported by a study by Lin et al. (2017) showing that the estimated K values based on the Torri and Shirazi models were closer to the measured value than NOMO and M-NOMO models.

Environmental factors that influenced soil erodibility
Based on the definition of K factor by Wischmeier et al. (1971), soil erodibility is estimated from texture data, organic matter content, soil structure index, and the soil permeability index.While soil erodibility does not directly depend on environmental factors, soil properties such as soil particle size distribution and soil organic matter can be affected by environmental factors.Soil erodibility thus has indirect relationships with environmental factors, particularly vegetation type, which influences the generation of soil organic matter and the composition of soil particles.Soil erodibility had various correlations with the selected environmental variables, which affected the dominant factors influencing soil erodibility (Tables S1-S5 and 3).In native grasslands, soil erodibility had significant correlations with terrain factors (e.g., elevation, slope degree; Tables S1 and S4), and the dominant factors influencing soil erodibility were soil properties and topography.Terrain factors have close relationships with soil properties.With changes in elevation and slope, the physical and chemical properties of soil (e.g., soil permeability, soil bulk density, and soil nutrients) and soil surface conditions (e.g., roughness, litter layer) change, leading to changes in soil particle size composition and soil erodibility (Zhao et al., 2015).For example, Li et al. (2011) found that the silt content was higher than the sand content in low but not high elevations, and Liu et al. (2005) found that slope gradient was negatively correlated with soil nutrients (e.g., soil organic matter, available nitrogen).
For most artificially managed vegetation types (apple orchards and David's peach) and artificially restored vegetation types (e.g., sea buckthorn and black locust), soil erodibility had significant correlations with vegetation properties (Tables S1, S3, and S4).By affecting physicochemical soil properties and soil structure stability, vegetation properties affect soil erodibility.For example, the dominant factors influencing soil erodibility were plant density for apple orchards; aboveground biomass for sea buckthorn litter biomass; vegetation coverage for black locust; and crown width, vegetation height, basal diameter, and vegetation coverage for David's peach (Table S1).Human activities (e.g., pruning) affect vegetation recovery and land cover change.These changes may www.solid-earth.net/9/1507/2018/Solid Earth, 9, 1507-1516, 2018 then influence vegetation properties and thereby impact soil erodibility.

Conclusions
We evaluated soil erodibility in the Ansai watershed using five estimation models.The estimated K values differed among the different models and ranged between 0.009 and 0.092 t • hm 2 • h/(MJ • mm • hm 2 ).Based on Taylor diagrams and previous studies, we considered the Shirazi and Torri models as the optimal models for the Ansai watershed.Since soil erodibility is estimated by soil properties, it has indirect relationships with environmental factors, including elevation and slope degree and, to a lesser extent, human activities.By changing vegetation density, biomass, and cover, humans can indirectly affect soil erodibility.
Data availability.Please contact the corresponding author if you need relevant data.
Author contributions.WZ and HW mainly conducted article structure design, data processing, and article writing; LJ, SD, XZ, and YL mainly engaged in article modification and polishing.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Locations of the study area and the sampling points.

Figure 2 .
Figure 2. Taylor diagram used to compare estimated K values among models.

Table 1 .
Landscape and soil characteristics in the study area.: NG denotes native grassland, AO denotes apple orchard, FL denotes farmland, PG denotes pasture grassland, SB denotes sea buckthorn, CK denotes Caragana korshinskii, DP denotes David's peach, BL denotes black locust, Ele denotes elevation, SP denotes slope position, SA denotes slope aspect, SG denotes slope gradient, SS denotes slope shape, Cla denotes clay, Sil denotes silt, San denotes sand, OM denotes organic matter, SBD denotes soil bulk density, Por denotes porosity, AAR denotes average annual rainfall, VC denotes vegetation coverage, AB denotes aboveground biomass, VH denotes vegetation height, LB denotes litter biomass, PD denotes plant density, Cro denotes crown, BD denotes basal diameter, and BN denotes branch number. Annotation

Table 2 .
Statistics of soil erodibility in the Ansai watershed.: EPIC denotes the erosion-productivity impact model, NOMO denotes the nomograph equation, M-NOMO denotes the modified nomograph equation, Torri denotes the K-value estimation model established by Torri, Shirazi denotes the K-value estimation model established by Shirazi, SD denotes the standard deviation, and Cv denotes the coefficient of variation. Annotation

Table 3 .
Principal component analysis (PCA) of environmental attributes.
Annotation: SS denotes slope shape, SP denotes slope position, SG denotes slope gradient, PD denotes plant density, Ele denotes elevation, AB denotes aboveground biomass, AAR denotes average annual rainfall, LB denotes litter biomass, SBD denotes soil bulk density, VC denotes vegetation coverage, Cro denotes crown width, VH denotes vegetation height, and BD denotes basal diameter.

Table 4 .
Suggested soil erodibility estimation models in China.