Rill Erosion and Soil Quality in Forest and Deforested Ecosystems with Di ﬀ erent Morphological Characteristics

: Rill detachment capacity is a key parameter in concentrated ﬂow erosion. Rill erosion generally turns into gully erosion with severe environmental impacts. Changes in land use and human activities can have heavy e ﬀ ects in rill formation, particularly in forests subject to deforestation; soil morphology plays a signiﬁcant role in these e ﬀ ects. However, literature reports few studies about rill detachment rates and their implications on soil quality in forest and deforested soils with di ﬀ erent morphological characteristics. To ﬁll these gaps, this study has evaluated the rill detachment capacity (D c ) and the main soil quality indicators in three areas (upper, middle and lower slope) of forest and deforested (for 10 years) hillslopes exposed to the north and south in Northern Iran. The variations of D c have been measured on soil samples under laboratory conditions through a ﬂume experiment at three slope gradients (12 to 19%) and ﬁve ﬂow rates (0.22 to 0.67 L m − 1 s − 1 ) with four replications. The large and signiﬁcant ( p < 0.05) di ﬀ erence (about 70%) detected for D c between forest and deforested hillslopes was associated to the higher organic matter content of forest areas; as a consequence, these areas also showed higher aggregate stability, porosity, root weight density, microbial respiration and available water. In the deforested hillslopes exposed to the south, the soil erodibility was higher by 12% compared to those exposed to the north. The di ﬀ erences in the monitored soil quality indicators were instead less noticeable and not always signiﬁcant ( p < 0.05). Conversely, D c did not signiﬁcantly change ( p < 0.05) among the upper, middle and lower hillslope areas investigated in this study. Simple but accurate models to predict the rill detachment capacity, erodibility and critical shear stress of soils from indicators of soil quality or the unit stream power using regression equations are suggested. Overall, the results can support land planners in prioritizing the actions for soil conservation in deforested hillslopes exposed to the south as well as in the extensive application of the proposed equations in erosion prediction models.


Introduction
The forest resources play a fundamental role in the survival of the entire planet. However, both forest plants and soils are subject to several degradation factors [1] due to natural and anthropogenic processes (e.g., climate change, erosion, deforestation, wildfire, urbanization) [2]. In particular, the Wang et al. [22] found that soil detachment capacity is significantly higher on lower positions (on average, by about 70% compared to upslope areas) [22]. Again, Wang et al. [22] and Zhang et al. [26] found that the rill detachment capacity increases with the slope by an exponential law. Moreover, very few studies have explored the relations between the rill detachment capacity and the changes in physico-chemical and biological properties. Li et al. [28] demonstrated that the variability of rill detachment capacity under different land uses in soils of the Loess Plateau (China) is positively related to the silt content (r > 0. 44), and inversely related (|r| < −0.24) to the sand content, cohesion, water stable aggregate, aggregate median diameter, organic matter and root density. However, although from this study it was evident that land use and soil types had significant effects on erodibility, the effects of soil morphology (e.g., aspect, slope, position) and forest activities (e.g., deforestation or lack of management) on soil detachment capacity as well as on soil quality have not been evaluated. Evidently, more research is needed to identify which morphological factors (e.g., hillslope aspects and positions) trigger rill erosion and soil quality degradation both in forest and deforested environments; moreover, reliable models to predict the soil detachment capacity from simple input parameters may support the actions of land managers to control erosion in forest and deforested ecosystems. To fill these gaps, this study aims at estimating the variability of rill detachment capacity and the associated changes in soil quality indicators in forest and deforested hillslopes with different morphological characteristics in Northern Iran. The forest hillslopes of these areas are a representative case study, because deforestation and inappropriate management have been very intense for many years [35,36]. Several factors, including urbanization, expansion of arable land and intensive practices, have reduced forest cover and soil quality, increasing soil erosion and degrading its quality [9,10,27,37]. Samples of soils were collected in forest and deforested hillslopes with different aspects (north and south) and positions (upper, middle and lower slope); the rill detachment capacity was measured on these samples using a laboratory flume at different water flow rates and soil slopes. Moreover, the main physico-chemical and biological properties of the soil samples were determined and related to the measured rill detachment capacity.
Overall, this study is useful to verify (i) whether the aspect and slope have significant effects on rill detachment capacity in forest and deforested hillslopes and (ii) which soil properties have significant effects on rill erodibility. Finally, simple models to predict rill detachment capacity from hydraulic parameters are proposed and validated for the experimental conditions.

Study Area
Four hillslopes in the Saravan Forestland Park, Guilan province (37 • 08 04" N, 49 • 39 44" E), Northern Iran, were selected for this study (Figure 1a). The climate of this area is typically Mediterranean, Csa type according to the Köppen-Geiger classification [38]. The average annual rainfall and temperature are 1360 mm and 16.3 • C, respectively [39]. Precipitation is mainly concentrated in winter, while the hottest months are in summer ( Figure 2).
The four hillslopes, about 100-m long and 20-m wide, were covered by forest trees and plants with high density; two hillslopes (hereinafter indicated as "forested hillslopes") were not subject to any management operation, while two others were deforested about ten years ago to install high-voltage power towers ("deforested hillslopes") [31]. One of the deforested hillslopes was exposed to south ("DS") and one to north "(DN"); the forest hillslopes also had different aspects ("FS", south, and "FN", north) (Figure 1b). Each of the four hillslopes, with uniform profile (that is, with no convex or concave shape), was divided in three positions from the top to the bottom: upper, middle and lower slopes ("US", "MS" and "LS", respectively) ( Figure 1b). The hillslopes were selected at altitudes between 50 and 250 m above the mean sea level. The mean slope gradient was about 15%. The vegetation cover was approximately 60% in LS, and lower in MS (50%) and US (30%). Vegetation was absent in the deforested hillslopes. about 15%. The vegetation cover was approximately 60% in LS, and lower in MS (50%) and US (30%). Vegetation was absent in the deforested hillslopes.  The soil of the study area was classified as ultic Hapludalfs with silty clay loam texture in all hillslopes [40,41]. Soil quality was variable among the investigated hillslopes, because of different conditions (forest or deforested) and morphology (aspect and position).

Soil Sampling and Analysis
The soil samples were randomly collected in three positions (US, MS and LS) of the four hillslopes (DS, DN, FS and FN) with four sampling points for each position (total of 48 samples, 2 conditions × 2 aspects × 3 positions × 4 replicates) (Figure 1b). The 48 samples were transported to the laboratory and quickly analyzed to measure the following physico-chemical and biological properties (hereinafter indicated as "soil quality indicators"): (a) organic matter (OM), by Walkley-Black method [42]; (b) aggregate stability (AS) in water; (c) porosity (SP); (d) bulk density (BD) (b, c and d using the wet-sieving and oven-drying methods) [43]; (e) microbial respiration (SMR, as soil biological activity in terms of accumulated C-CO2 using Anderson [44] method; (f) CaCO3, according to Sparks [45], expressed as calcium carbonate equivalent (CCE). Moreover, the soil texture was determined using the hydrometer method [46]. Finally, the available water (AW) of soil samples was estimated as difference between the field capacity and permanent wilting point [47]. The values of soil quality indicators are reported on Table S1 (SD) of the Supplementary Materials.

The Experimental Device and Procedure
A hydraulic flume (3.5-m long and 0.2-m wide) with a rectangular cross section was used to measure the rill detachment capacity of soil samples collected at each hillslope, considering different values of flow rate (0.26, 0.35, 0.45, 0.56, and 0.67 L m −1 s −1 ) and soil slope gradient (12%, 16%, and 19%). Asadi et al. [48], Raei et al. [49] and Parhizkar et al. [27] have detailed the flume characteristics in previous studies, while the experimental procedure is reported in the study of Parhizkar et al. [27]. To summarise, the soil was collected from each hillslope using a steel ring with 0.1 m in diameter and The soil of the study area was classified as ultic Hapludalfs with silty clay loam texture in all hillslopes [40,41]. Soil quality was variable among the investigated hillslopes, because of different conditions (forest or deforested) and morphology (aspect and position).

Soil Sampling and Analysis
The soil samples were randomly collected in three positions (US, MS and LS) of the four hillslopes (DS, DN, FS and FN) with four sampling points for each position (total of 48 samples, 2 conditions × 2 aspects × 3 positions × 4 replicates) (Figure 1b). The 48 samples were transported to the laboratory and quickly analyzed to measure the following physico-chemical and biological properties (hereinafter indicated as "soil quality indicators"): (a) organic matter (OM), by Walkley-Black method [42]; (b) aggregate stability (AS) in water; (c) porosity (SP); (d) bulk density (BD) (b, c and d using the wet-sieving and oven-drying methods) [43]; (e) microbial respiration (SMR, as soil biological activity in terms of accumulated C-CO 2 using Anderson [44] method; (f) CaCO 3 , according to Sparks [45], expressed as calcium carbonate equivalent (CCE). Moreover, the soil texture was determined using the hydrometer method [46]. Finally, the available water (AW) of soil samples was estimated as difference between the field capacity and permanent wilting point [47]. The values of soil quality indicators are reported on Table S1 (SD) of the Supplementary Materials.

The Experimental Device and Procedure
A hydraulic flume (3.5-m long and 0.2-m wide) with a rectangular cross section was used to measure the rill detachment capacity of soil samples collected at each hillslope, considering different values of flow rate (0.26, 0.35, 0.45, 0.56, and 0.67 L m −1 s −1 ) and soil slope gradient (12%, 16%, and 19%). Asadi et al. [48], Raei et al. [49] and Parhizkar et al. [27] have detailed the flume characteristics in previous studies, while the experimental procedure is reported in the study of Parhizkar et al. [27]. To summarise, the soil was collected from each hillslope using a steel ring with 0.1 m in diameter and 0.05 m in height). The soil surface was wetted by light spraying and then inserted in a hole of the flume bed. Then, the water flow rate and bed slope in the flume were set at the desired values. The water discharge was measured five times per experiment, using a graduated plastic cylinder. The surface flow velocity was measured ten times using the fluorescent dye technique. More specifically, the time for the tracer to travel from the injection point to the observation point was visually estimated [50]. The average flow velocity was calculated reducing the water velocity by 0.6 for laminar flow, 0.7 for transitional flow and 0.8 for turbulent flow [51]. Water temperature was determined to compute water viscosity. The water depth was measured at three points (0.01 m from the left and right sides and in the middle) for two cross sections (located at 0.4 m and 1 m from the flume outlet) using a level probe with accuracy of 1 mm. Then, the mean value of these six measurements was the average flow depth. When the depth of the eroded soil in the steel ring reached 0.015 m, the experiment was stopped. After completing each experiment, the sample of wet soil was oven dried for 24 h at 105 • C to determine its dry weight. The flume experiments were carried out on all soil samples collected. Overall, 720 soil samples (2 conditions × 2 aspects × 3 positions × 4 sampling points × 5 water flow rates × 3 slope gradients) were tested.

Measurement of Hydraulic Parameters and Rill Detachment Capacity
The rill detachment capacity (D c , kg s −1 m −2 ) was estimated as the mean value of four replicates using Equation (1): where ∆M is the dry weight of detached soil (kg), A is the area of the soil sample (m 2 ) and ∆t is the experiment duration(s).
The rill detachment capacity of the hillslopes with different conditions, aspects and positions was calculated as the average among the values measured under the replicates of different water flows and profile slopes in the flume.
Regarding the modeling approach, D c can be predicted with accuracy by a power function of some important hydraulic parameters, such as the hydraulic shear stress, stream power and unit stream power [25,52]. Therefore, according to Foster [19] and Yang et al. [53] and assuming a rectangular cross section for rills, the hydraulic radius (R, m), mean velocity (V, m s −1 ), shear stress (τ, Pa) and unit stream power (ω, m s −1 ) were calculated as follows: where: h = flow width (0.2 m in this study) - The values of the calculated hydraulic parameters and D c are reported in Tables S2 (SD) and S3 (SD) of the Supplementary Data.
Parhizkar et al. [27] have demonstrated that, in forest areas, ω is the hydraulic parameter that is the best predictor of D c . Therefore, the measured values of D c (considered as the dependent variable) were regressed on ω (independent variables) using non-linear equations based on power function, to identify the most accurate predictive model for each hillslope condition and aspect. The model has therefore the following expression: where "a" and "b" are two coefficients that must be calibrated. The accuracy of these regression equations was evaluated using quantitative indexes, as the relative error (RE), coefficient of determination (r 2 ) and coefficient of efficiency (NSE) [54]. The optimal values of these indexes are one for r 2 and E as well as zero for RE, while the acceptance limits of r 2 and NSE are 0.50 and 0.35, respectively; the prediction capacity of a model is considered good when r 2 and are higher than 0.75 [55][56][57].
Rill erodibility (K r , s m −1 ) and critical shear stress (τ c , Pa) are key input parameters for erosion prediction using process-based erosion models [58], such as the Water Erosion Prediction Project (WEPP) [59]. These parameters, which reflect soil resistance to rill erosion, were calculated as the slope and intercept of the following regression equation interpolating τ and D c :

Statistical Analysis
First, a one-way analysis of variance (ANOVA) and Tukey's test (at p < 0.05) were used to evaluate the statistical significance of differences in rill detachment capacity and soil quality indicators among the studied hillslopes. The condition (forest or deforested), aspect (north or south) and position (upper, middle and lower) of hillslopes were considered as independent factors, while the rill detachment capacity and quality indicators were the dependent variables. The normality of sample distribution was checked using QQ-plots. To satisfy the assumptions of the statistical tests, the data were square root transformed whenever necessary.
Then, a Principal Component Regression (PCR) [60] was used as a predictive model of rill detachment capacity (dependent variable), using the soil quality indicators (independent variables) as input. PCR is a combination of Principal Component Analysis (PCA) and Multiple Linear Regression (MLR, often an Ordinary Least Squares, OLS, regression), where the Principal Component (PC) scores are used as predictor variables and linear combinations are constructed between predictor and response variables. PCR is advised when the factors are many and highly collinear [61]. As a matter of fact, many pairs of performance indicators adopted in our study were strongly correlated (as shown by the Pearson's correlation matrix) and therefore collinear. Since PCs are orthogonal, the multicollinearity problems of the original variables are removed [62]. PCR was implemented as follows: (i) a PCA was applied to D c and soil quality indicators using Pearson's matrix method [63]; (ii) an OLS regression was run on the first two PCs; (iii) the dependent variable (D c ) was calculated as a linear combination of the first three PCRs, explaining most of the variance of the original independent variables. All statistical analyses were performed using XLSTAT ® software (Addinsoft, Paris, France).

Spatial Variations of Rill Detachment Capacity and Associations with Soil Quality Indicators
ANOVA showed that D c was significantly different between forest and deforested hillslopes and north vs. south aspects, but not among upper, middle and lower positions. On average, D c of deforested hillslopes exposed to south or north was 0.040 ± 0.020 and 0.035 ± 0.019 kg m −2 s −1 , respectively. Forest hillslopes showed D c of 0.014 ± 0.007 (exposed to south) and 0.010 ± 0.006 kg m −2 s −1 (north aspect) ( Figure 3). Regarding the soil quality indicators, the differences between forest and deforested soils were always significant according to ANOVA. In more detail, forest soils were less compacted (mean BD of 1374 ± 59.6 kg m −3 ) compared to the deforested areas (BD of 1554 ± 61.7 kg m −3 ) and the AS was higher as well (1.47 ± 0.29% against 0.69 ± 0.18% of deforested areas). The deforested hillslopes had lower SP (41.4 ± 2.33% against 48.1 ± 2.25% of forest areas) and AW (11.0 ± 1.93% against 14.5 ± 2.15%); also the OM content and SMR were lower (1.94 ± 0.36% and 0.10 ± 0.03 g of CO 2 per kg of soil, respectively) compared to the forest hillslopes (3.23 ± 0.28% and 0.21 ± 0.02 g of CO 2 per kg of soil). RWD of the forest hillslopes was 0.60 ± 0.11% and zero in the deforested areas ( Table 1). The hillslope aspect influenced OM, RWD, AS, CCE, AW and SMR, but not the other soil quality indicators (texture, Resources 2020, 9, 129 8 of 24 BD and SP). More specifically, OM content and AS of soils were lower in the hillslopes exposed to south (2.40 ± 0.68% and 0.96 ± 0.38%, respectively) compared to north (2.77 ± 0.74% and 1.20 ± 0.50%). The latter hillslopes showed higher SMR (0.17 ± 0.05 against 0.14 ± 0.06 g of CO 2 per kg of soil) and RWD (0.11 ± 0.11 against 0.56 ± 0.63 kg m −3 ) compared to the soils exposed to south. Finally, the AW was higher in north (13.7 ± 2.51%) and lower in south hillslopes (11.8 ± 2.58%) ( Table 1).
Resources 2020, 9, x FOR PEER REVIEW 8 of 31 g of CO2 per kg of soil, respectively) compared to the forest hillslopes (3.23 ± 0.28% and 0.21 ± 0.02 g of CO2 per kg of soil). RWD of the forest hillslopes was 0.60 ± 0.11% and zero in the deforested areas ( Table 1). The hillslope aspect influenced OM, RWD, AS, CCE, AW and SMR, but not the other soil quality indicators (texture, BD and SP). More specifically, OM content and AS of soils were lower in the hillslopes exposed to south (2.40 ± 0.68% and 0.96 ± 0.38%, respectively) compared to north (2.77 ± 0.74% and 1.20 ± 0.50%). The latter hillslopes showed higher SMR (0.17 ± 0.05 against 0.14 ± 0.06 g of CO2 per kg of soil) and RWD (0.11 ± 0.11 against 0.56 ± 0.63 kg m −3 ) compared to the soils exposed to south. Finally, the AW was higher in north (13.7 ± 2.51%) and lower in south hillslopes (11.8 ± 2.58%) ( Table 1).   Notes: D = deforested hillslope; F = forest hillslope; S = south; N = north; SaC = sand content; SiC = silt content; ClC = clay content; OM = organic matter content; BD = bulk density; AS = aggregate stability; CCE = calcium carbonate equivalent; SMR = soil microbial respiration; RWD = root weight density; SP = soil porosity; AW = available water; different lowercase and capital letters indicate significant differences between D and F condition as well as N and S aspect, respectively, after Tukey's test (at p < 0.05) following ANOVA.
High coefficients of regression were detected between D c and many soil quality indicators. In more detail, the correlations were positive with SiC (r = 0.880), ClC (r = 0.787) and CCE (r = 0.760) and negative with OM (r = −0.860), RWD (−0.796), AS (−0.768), SaC (r = −0.840) and SMR (r = −0.759). The other correlations among D c and BD, SP and AW were lower (r < |0.662|), but always significant ( Table 2). Interesting correlations were also found between couples of soil quality indicators. For instance, OM was correlated to BD (r = −0.777), AS (r = 0.860) and SMR (r = 0.880); BD was inversely correlated to RWD (r = −0.823), on its turn strictly linked to SMR and SP (r = 0.899 and 0.823, respectively). The weakest correlations, although being always significant, were detected between AW and the other soil quality indicators (r = |0.678|) ( Table 2).
The application of PCR to the soil quality indicators provided a Principal Component (PCR1) that alone explains 81.75% of the total variance of the original variables; the variance explained by the other PCRs was much lower (6.82% for PCR2). All these variables had high factor loadings on PCR1 and, particularly, OM (loading of 0.939), RWD (0.969), SMR (0.948), SaC and ClC (0.980 and −0.967, respectively) and CCE (−0.941). The loadings of all physico-chemical properties of soils on PCR2 were much lower and never significant (Table 3 and Figure 4a).    When the scores of the soil samples collected under different aspects and conditions of hillslopes were plotted on the PCR1-PCR2 chart, two evident clusters were identified. A first cluster grouped all the forest soils, associated to positive values of PCR1, while the second cluster was related to samples of deforested soils, associated to negative PCR1. No evident clusters grouping soil samples collected in hillslopes with different aspects could be found (Figure 4b).

Modeling the Rill Detachment Capacity and Erodibility
PCR identified a reliable multiregression model (r 2 = 0.89) that predicts D c from the first two PCRs. This model has the following analytical structure: The scatter plot of Figure 5 shows that almost all the predicted values lie inside the 95%-confidence interval; moreover, the values of RE (3.43%), r 2 (0.83) and NSE (0.83) indicate good model accuracy.
Resources 2020, 9, x FOR PEER REVIEW 15 of 26

Modeling the Rill Detachment Capacity and Erodibility
PCR identified a reliable multiregression model (r 2 = 0.89) that predicts Dc from the first two PCRs. This model has the following analytical structure: The scatter plot of Figure 5 shows that almost all the predicted values lie inside the 95%-confidence interval; moreover, the values of RE (3.43%), r 2 (0.83) and NSE (0.83) indicate good model accuracy. High coefficients of determination were found correlating Dc with ω (r 2 > 0.938) ( Table 4 and Figure 6). Moreover, the regression models established between these variables were accurate, as visually shown by the scatter plot of Figure 7, where the measured-predicted Dc points are grouped close to the identity line. Also, the quantitative indexes, adopted to evaluate the model accuracy, were very satisfactory (RE < 1.4%, r 2 and NSE > 0.96).  High coefficients of determination were found correlating D c with ω (r 2 > 0.938) ( Table 4 and Figure 6). Moreover, the regression models established between these variables were accurate, as visually shown by the scatter plot of Figure 7, where the measured-predicted D c points are grouped close to the identity line. Also, the quantitative indexes, adopted to evaluate the model accuracy, were very satisfactory (RE < 1.4%, r 2 and NSE > 0.96). Table 4. Power regression equations between soil detachment capacity (Dc, kg m −2 s −1 ) and unit stream power (ω, m s −1 ) measured on soil samples collected in four hillslopes of Saravan Forestland Park (Northern Iran).

Hillslope Characteristic Power Regression Equation (5) r 2 Condition Aspect
Deforested South D c = 0.947ω 1  When Dc was regressed on τ using Equation (6), the differences in the slope (Kr) and intercept (τc) were higher between hillslopes with different conditions (forested vs. deforested) and lower between hillslopes with different aspects. These equations had coefficients of regression (r 2 ) of 0.86-0.88 (Table 5 and Figure 8). Deforested hillslopes exposed to the south and forest areas with north aspects showed the maximum (0.0093 s m −1 ) and minimum (0.0027 s m −1 ) values of Kr, respectively, corresponding to τc of 3.11 and 3.57 Pa (Table 5). Also, these regression models were reliable for Dc predictions, as shown by the limited scattering of the measured-predicted Dc around the identity line ( Figure 9).  When D c was regressed on τ using Equation (6), the differences in the slope (K r ) and intercept (τ c ) were higher between hillslopes with different conditions (forested vs. deforested) and lower between hillslopes with different aspects. These equations had coefficients of regression (r 2 ) of 0.86-0.88 (Table 5 and Figure 8). Deforested hillslopes exposed to the south and forest areas with north aspects showed the maximum (0.0093 s m −1 ) and minimum (0.0027 s m −1 ) values of K r , respectively, corresponding to τ c of 3.11 and 3.57 Pa (Table 5). Also, these regression models were reliable for D c predictions, as shown by the limited scattering of the measured-predicted D c around the identity line (Figure 9). Table 5. Linear regression equations between soil detachment capacity and critical shear stress measured on soil samples collected in four hillslopes of Saravan Forestland Park (Northern Iran).

Hillslope Characteristic
Linear Regression Equation (

Spatial Variations of Rill Detachment Capacity and Associations with Soil Quality Indicators
The soil detachment capacity was significantly lower (on average by 68%) in forest hillslopes compared to the deforested areas and this result was expected. As a matter of fact, the higher erodibility found in the deforested hillslopes is consistent with the results of several studies [64][65][66][67]. In the deforested hillslopes exposed to south the mean Dc was higher by 12% compared to the north, while the variability between the two aspects was slightly higher (14%), considering both forest and deforested hillslopes. Conversely, although a variability in Dc values was detected among areas with different positions (upper, middle and lower slope), the related changes were small and never significant. This should be due to the similar soil characteristics among the three hillslope positions, which were not influenced by external factors (such as climate and past erosion). This is in contrast with some other studies: for instance, Khormali et al. [8] found that the upward slopes were more prone to erosion, especially in the deforested hillslopes, where organic matter content is lower and bulk density is higher with consequent lower infiltrability and more runoff generation capacity. Li et al. [28,31] and Wang et al. [22], proposing accurate regression equations to estimate Dc from hydraulic parameters, showed that hillslope position with different slopes plays a significant influence on rill erodibility.  Table 5. Notes: D = deforested hillslope; F = forested hillslope; S = south; N = north.

Spatial Variations of Rill Detachment Capacity and Associations with Soil Quality Indicators
The soil detachment capacity was significantly lower (on average by 68%) in forest hillslopes compared to the deforested areas and this result was expected. As a matter of fact, the higher erodibility found in the deforested hillslopes is consistent with the results of several studies [64][65][66][67]. In the deforested hillslopes exposed to south the mean D c was higher by 12% compared to the north, while the variability between the two aspects was slightly higher (14%), considering both forest and deforested hillslopes. Conversely, although a variability in D c values was detected among areas with different positions (upper, middle and lower slope), the related changes were small and never significant. This should be due to the similar soil characteristics among the three hillslope positions, which were not influenced by external factors (such as climate and past erosion). This is in contrast with some other studies: for instance, Khormali et al. [8] found that the upward slopes were more prone to erosion, especially in the deforested hillslopes, where organic matter content is lower and bulk density is higher with consequent lower infiltrability and more runoff generation capacity. Li et al. [28,31] and Wang et al. [22], proposing accurate regression equations to estimate D c from hydraulic parameters, showed that hillslope position with different slopes plays a significant influence on rill erodibility.
Several studies have demonstrated that the variability of soil detachment capacity can be associated to changes in soil properties, such as texture, bulk density, cohesion and water stable aggregate [20,26,31]. A review proposed by Knapen et al. [20] evidenced that, among the soil properties, the texture, structural stability of aggregates and organic matter play the main role in the resistance of soil to concentrated flow erosion. Zhang et al. [26] stated that, when the texture of soils with different land uses is the same, the clay content, aggregate median diameter, soil strength, bulk density, in addition to plant root density, are responsible for the differences in D c among the land uses. According to Li et al. [31], the variability of D c under different land uses is positively related to silt content, and inversely related to sand content, cohesion, water stable aggregate, aggregate median diameter, organic matter and root density of plants. Also, in this study, the differences in D c among hillslope conditions and aspects may be ascribed to the variability of the soil quality indicators, which were often significantly different (mainly for the soil condition) among the samples. These relationships between D c and the soil characteristics are confirmed by the strong correlations revealed by Pearson's matrix. Generally, for the investigated hillslopes, the soil erodibility is directly associated with the particle fractions of the sampled soils, despite the same soil texture. D c increased ClC content, as shown by the positive correlations among these soil quality indicators. Although a higher content of clay particles enhances soil cohesion action, which should tie the finer particles and decrease the soil erodibility [26,68,69], the higher D c of clay-richer soils should be attributed to other factors. These results are in accordance with Li et al. [31] and Parhizkar et al. [27], Ciampalini et al. [68], Zhang et al. [26] and Vaezi et al. [70]. These authors stated that, in soils with higher silt and clay contents, particle detachability increases. Conversely, when the SaC increases, particle detachability decreases, as reported also by Parhizkar et al. [27]. However, the direct associations between soil erodibility and texture characteristics should be generalised with caution, since some studies suggested that soils with high clay [71] and silt [23] contents have low rill erodibility and the effects of sand particles on D c are ambiguous [31]. Compared to the deforested hillslopes, forest soil showed a higher content of OM (on average 67% more), which is considered as one of the most common soil quality indicators [72]. The negative correlation found between D c on one side and OM, RWD, AS, SP and SMR on the other side confirms that the forest soils, which have a noticeable vegetation cover, are less erodible compared to soils with lower or no vegetation.
It is well known that soils with noticeable vegetation cover show a high content of organic matter [27,[73][74][75]. This higher OM content is due to root presence of the vegetation in the forest areas [76], where plant roots bind soil particles at the soil surface [77,78]. The effectiveness of roots in reducing detachment to flowing water is closely related to its physically binding and chemically bonding effects [18,31,79]. Vegetation roots bind soil particles at or near the soil surface to protect soil from water erosion and enhance its stability [77]. This is a more general process due to the control that vegetation exerts on soil particle detachment [80].
Besides higher OM content and RWD of the forest areas, the larger cover of vegetation increases porosity and aggregate stability of soils, as shown in our study, where SP and AS of forest soils was more than 16% and 100% higher, respectively, compared to deforested hillslopes. The soil structure stability is generally recognised as having a positive effect in reducing D c [20,31]. The lower soil AS found in deforested hillslopes is in line with the results of Caravaca et al. [81] and Khormali et al. [8]. In more detail, Caravaca et al. [81] showed that the aggregate stability of cultivated soils is significantly lower than that of forest soils. In the study of Khormali et al. [8], deforestation resulted in a significant decrease in AS of the soil surface in all the investigated slope positions.
Moreover, many researchers, such as Shepherd et al. [82], stated that the disintegration of soil aggregates decreases the soil porosity and therefore increases its bulk density. According to McDonald et al. [83], a BD increase exposes the deforested areas to greater erosion rates, while Lemenih et al. [84] associates decreases in SP and increases in BD to a decline in OM content of soils and therefore to worsened soil quality. A negative correlation between D c and soil porosity was also found in this study. The latter soil property was higher in forest hillslopes, where AW is higher (by 32%), presumably due to the higher AS and OM content of soil. Organic matter is a cementing agent and can help to increase the water content of soils [85]. This is in accordance with other studies [86,87], which reported high water content in soils with high amount of OM. Moreover, according to Rawls et al. [86], water retention of soils with coarse texture is substantially more sensitive to the amount of organic matter as compared with fine-textured soils. Therefore, deforestation and inappropriate management operations can reduce organic matter and thus water retention of soil. The higher water content that soil can store is beneficial for plant growth [88,89].
The high OM content of forest hillslopes is also linked to a high SMR. Conversely, the reduction in OM following deforestation reduced the SMR by over 100%, and a reduced microbial activity in soils is related to lower levels of available organic carbon [8,90]. The high production of biomass in forest soils plays an important effect on its microbial population [91]. Kiani et al. [92] stated that, when the rate of fresh plant residues increases in forest areas, the soil microbial respiration rises, which strictly depends on the nutrient contents of soil (i.e., P, K, Ca and Mg) [93][94][95].
If the differences in D c between forest and deforested hillslopes can be explained by the changes in the main soil quality indicators monitored in this study, the reasons for the variability of D c between the two hillslope aspects are less evident. The higher D c of south-facing hillslopes compared to the areas exposed to north is in accordance with some studies [96,97]. More specifically, Cerdà [96], under simulated rainfalls, demonstrated that both runoff and erosion increase from north-facing to south-facing afforested slopes as the results of variation in soil cover and organic matter content. Also, Petersen et al. [97] reported that water infiltration is significantly lower and sediment load in runoff is higher in soil with south aspect compared to north. In our study, in hillslopes exposed to north the OM content and AS were higher compared to south-facing areas (by about 15% and 24%, respectively) and the same contrast was found for SMR (+28%), RWD (+12%) and AW (+16%). The lower OM in south-facing hillslopes could be due to the higher temperature to which the soils are exposed [98], which accelerate the oxidation of the organic compounds [99,100].
Also, the PCR basically confirms the relationships between D c and soil quality indicators. When OM, AS, RWD and SP decrease, the rill detachment capacity increases. The forest soils are associated to higher values of OM, SMR, RWD, AS and SP (having high positive loadings on PCR1) and the opposite situation happens for soil samples collected in deforested hillslopes (that is, lower OM, SMR, RWD, AS and SP). The latter soils are characterised by higher BD and CCE, which negatively weigh on PCR1. Moreover, forest soils show higher AW compared to deforested hillslopes. Also, Wang et al. [22] reported negative correlations among the soil detachment capacity and ClC, OM and RWD, and a positive correlation with SP.

Modeling the Soil Detachment Capacity and Rill Erodibility
It is important to develop reliable and effective models for estimating rill detachment capacity from some measurable parameters, including hydraulic, soil, and vegetation indicators [22]. Several researchers proposed equations to predict D c in different regions by using hydraulic and soil parameters [26,68]. The regression analysis (by PCR as well as power and linear models) carried out in this study is helpful to predict rill erodibility parameters using some important soil quality indicators or hydraulic characteristics of the overland flow.
The multiregression model (7) to predict rill detachment capacity from the first two PCRs showed that D c can be estimated with acceptable accuracy using a linear combination of the soil quality indicators (which are the original variables, to which the derivative variables of PCR are correlated [101]). This is confirmed by the fact that the predicted values of D c lie inside the 95%-confidence interval and, therefore, a reliable estimation of D c can be achieved when the main soil characteristics are known. This outcome is in line with findings of Zhang et al. [26], who indicated that D c in rill erosion could be well simulated by clay content, bulk density, aggregate median diameter and soil strength. Moreover, Wang et al. [22], proposed to predict D c using an exponential equation from BD, OM, root density, τ and cohesion in hillslopes with permanent gullies in southeast China, achieving R and NSE equal to 0.98. The first PCR of the multiregression model (7) proposed in this study can be interpreted as a combined index of soil quality, since it derives from some important indicators of soil quality (e.g., OM, RWD, AS, AW, SMR) and an increased PCR1 indicates an improved soil quality.
The rill detachment capacity also depends on the hydraulic characteristic of the overland flow [48]. The regression analysis carried out in this study confirmed that unit stream power is a good predictor of D c in the experimental conditions, and particularly for the deforested hillslopes, which are more sensitive to soil erosion [8,102]. This high accuracy may be due to the fact that ω simultaneously takes into account the flow velocity and soil slope, both influencing the soil detachability due to the overland flow [27]. Also, Zhang et al. [26] found that D c can be accurately predicted by hydraulic parameters, such as stream power, slope and runoff density. These authors achieved r 2 of 0.93 and NSE of 0.90 in predicting soil detachment rates from stream power in soils collected in woodland. However, the extrapolation of these equations to other soil conditions and types must be done with caution, as stated by Nearing et al. [52], who stated that stream power is sometimes not an appropriate parameter to uniquely describe detachment rates for different soils. For instance, the application of the similar equations proposed by Parhizkar et al. [27] for forestland and woodland and validated in soils of Northern Iran shows a model accuracy that is lower than in this study (RE between 12% and 42% as well as NSE between −0.25 and 0.66). However, those experimental conditions were different than in this study (moderate differences in soil texture and quality indicators, and different hillslope aspect and vegetation), and this may affect the model prediction accuracy. This does not mean that the models developed in these studies are useless, but they should be applied only in the experimental conditions and validated case by case in soils of different characteristics. Conversely, the proposed equations provide accurate evaluations of rill erodibility in soil subject to conservation or management practices, where the expected changes in soil quality indicators can not affect the reliability of model predictions.
The D c -ω curves showed that the slopes of the regression models in the forest hillslopes were noticeably higher compared to the deforested areas, indicating that D c increases more rapidly with unit stream power. Regarding the same soil condition (forest or deforested), the hillslopes exposed to the south showed higher D c compared to north aspect.
The rill detachment capacity of a soil due to the overland flow is strictly linked to the rill erodibility and critical shear stress. These parameters reflect soil resistance to runoff [59] and are sensitive input parameters in process-based erosion models [23,58,103], such as the WEPP model [104]. In our study, the high coefficients of determination in the linear regression fitting D c to τ show the accuracy of these models in predicting the parameters linked to the soil resistance to erosion (K r and τ c ). This accuracy is consistent with the findings of Zhang et al. [26], who stated that the regressions between D c and τ can be used to fit the observed data satisfactorily, with r 2 up to 0.9. The higher values of K r in the deforested hillslopes, compared to the forest areas, again confirm the essential role of forest in reducing soil erosion in vegetated areas [13].

Conclusions
This study has confirmed in the forest and deforested hillslopes of the case study the expected significant differences in the rill detachment capacity. The variability of this soil property between the areas exposed to the north and south was lower but equally significant. Conversely, the soil erodibility was quite similar and thus not significant among the different soil positions. The differences detected in rill detachment capacity have been associated to the higher organic matter content, aggregate stability, porosity, root weight density, microbial respiration and available water of the forest hillslopes in comparison with the deforested areas. All these soil quality indicators were also higher in the hillslopes exposed to the north compared to the soils with south aspect. However, these differences were less noticeable than the changes detected in soils with different conditions.
The simple models to predict rill erosion proposed in this study may be adopted in order to control and mitigate the risks of erosion and soil quality degradation. The detachment capacity, rill erodibility and critical shear stress of soils can be predicted in the experimental conditions from indicators of soil quality or unit stream power using reliable and accurate regression equations. However, further studies should test the reliability of the suggested equations in other environmental conditions.
Overall, this study has contributed to the understanding of the spatial variability of soil quality and rill erodibility in forest areas. The results can support land planners and watershed managers in prioritizing the actions for soil conservation as well as in the extensive application of erosion prediction models towards the conservation of natural resources.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2079-9276/9/11/129/s1. Table S1. Soil quality indicators of samples collected in the four hillslopes of Saravan Forestland Park (Northern Iran). Table S2. Hydraulic parameters of flume experiments on soil samples collected in the four hillslopes of Saravan Forestland Park (Northern Iran). Table S3