Path coefficient analysis, a different approach to identify soil quality indicators

ABSTRACT Soil quality indicators related to water erosion reduction can assist with correct soil management. The objective of this research was to identify some variables that could be used as soil quality indicators with the aid of path coefficient analysis in order to reduce water erosion. The research was carried out in the field between May 2011 and April 2013 in southern Brazil on an Inceptisol. The following treatments were studied under simulated rainfall conditions: 1) no-tilled, cultivated and covered by cultural residue of ryegrass (Lolium multiflorum) (HCR); 2) no-tilled, cultivated and covered by crop residue of vetch (Vicia sativa) (HCV); 3) cultivated and scarified soil containing ryegrass roots (HRR); 4) cultivated and scarified soil containing vetch roots (HRV); and 5) bare and chiselled soil (BHR). Eight simulated rainfalls were applied in each treatment. Flow velocity, soil and water losses as well as variables or soil attributes influenced by management were quantified. Path coefficient analysis verified that the coverage, surface roughness, water infiltration rate and total organic carbon have the greatest direct or indirect relationships with soil and water losses or runoff velocity. These variables were indicative of soil quality, particularly its resistance to water erosion. In a rough soil, the total organic carbon, root mass and macroporosity of the soil are more important as indicators for soil resistance to water erosion.


Introduction
Soil management can influence soil quality, that affects the magnitude of soil erosion. More recently, several studies have applied multivariate statistical methods to select the most relevant indicators. These are often based on assumed but not assessed connections between indicators and soil functions (Tesfahunegn, 2014;Askari & Holden, 2015;Obade & Lal, 2016).
However, it is important to note that explicit evaluation of soil quality with respect to specific soil threats has rarely been done (Bünemann et al., 2018).
Related to soil erosion resistance, Clearwater et al. (2016) developed the Tillage Erosion Risk Indicator (TillERI). Bramorski et al. (2012) found that the soil surface variables that most influenced its resistance to water erosion were the soil cover and surface roughness. Among the subsurface attributes, Martínez-Trinidad et al. (2012) highlighted the macroporosity, total porosity, pore size distribution and geometric mean diameter -GMD. Tian et al. (2018) using path analysis found that the soil bulk density was a factor that affected sediment yields directly and indirectly in shelter forests.
The objective of this study was to identify the attributes of the soil that resulted in soil quality indicators with a view to reducing water erosion. For this, path coefficient analysis was used to statistically interpret the results of water and soil losses and runoff velocities as well as surface and subsurface attributes in different soil management systems.

Material and Methods
The research was conducted between May 2011 and April 2013 and was located in the municipality of Lages in the state of Santa Catarina in Southern Brazil (27° 47' S and 50° 18' W, at 923 m altitude). The climate is Cfb type according to the Köppen classification, with an average annual rainfall of 1,533 mm (Schick et al., 2014). The soil is an Inceptisol (Ramos et al., 2016).
The plots had an area of 38.5 m 2 (11 m in length parallel with the slope and 3.5 m in width). The plots were delimited with 0.20 m high galvanized sheets, embedded 0.10 m into the soil. At the lower part of the plots, the delimitation was done with a collecting trough connected to a PVC pipe used for concentrating the runoff flow for measurements.
The ryegrass and vetch crops were sown manually in the treatments, and fertilizer application with the formulation 7-30-15 (N-P 2 O 5 -K 2 O, 300 kg ha -1 ) was simultaneously carried out. Seed application rates were 60 kg ha -1 for ryegrass and 100 kg ha -1 for vetch. After reaching the blooming stage, the plants were harvested, and then the aerial parts were removed from the soil surface of the scarification treatments.
Five treatments were studied, with two field replicates, totalling 10 randomly arranged plots (completely randomized design). The following treatments were studied: i) notilled, cultivated and covered by ryegrass residue (Lolium multiflorum) (HCR); ii) no-tilled, cultivated and covered by vetch residue (Vicia sativa) (HCV); iii), soil cultivated with ryegrass and prepared with a scarification after removal of the aerial part of the crop residue (HRR); iv) soil cultivated with vetch and prepared with a scarification after removal of the aerial part of the crop residue (HRV); and v) bare and chiselled soil, with high roughness (BHR).
In each treatment, eight simulated rainfalls were applied for 90 min at a planned constant intensity of 65 mm h -1 , with the first being just after treatment installation (17/12/2011), and the last test was applied on 18/12/2012 (one year between the 1 st and the 8 th tests). The intervals were 24 days between the 1 st and 2 nd tests and 28, 32, 62, 99, 76 and 46 days between each of the other simulated rain tests (3 rd -8 th ).
A thrust-type rainfall simulator with rotating arms that simultaneously covered two plots was used to apply the water. To achieve the intensity planned, 15 open VEJEET 80/100 sprinklers were used, with pressure of 6.2 psi. The volumes and intensities of the rainfall applied were quantified after the end of 90 min of rainfall by using 20 rain gauges arranged around the plots.
Soil samples were collected from the 0.025 m layer immediately before each rainfall and divided into two groups: not deformed, which was used to determine soil bulk density, pore volume (macropores, micropores and total porosity) and gravimetric water content in the soil in the 0-0.10 and 0.10-0.20 m layers; and deformed for the determination of particle density, total organic C and geometric mean diameter (GMD). The undeformed soil samples were collected in metal rings with 0.025 m in height and 0.06 m in diameter with sharp edges.
The total pore volumes were calculated from the relationships between soil densities and particle densities, and the micropore volumes were calculated from water retention after saturation of the soil sample and submission to a 6-kPa tension in a sand tension table. The volumes of macropores were obtained by the difference between the total volumes of pores and micropores. Soil density was determined by the dry soil mass/volume ratio of the ring at 105 °C. Soil particle density was determined by the volumetric flask method (EMBRAPA, 2017).
The geometric mean diameter (GMD) was calculated according to the methodology described by Kemper & Chepil (1965), to express the effect of water on the stability of aggregates, during vertical stirring of soil samples using sieves with mesh sizes of 4.76; 2; 1 and 0.25 mm. Total organic C (TOC) was determined by oxidation with 1.25 mol L -1 K 2 Cr 2 O 7 in acidic media and titration with 0.25 mol FeSO 4 L -1 (Tedesco et al., 2016).
For the microrelief readings, a roughness meter (rugosimeter) with rods was used that had a photographic camera coupled for recording the rod heights. The readings were carried out at one point per plot, with the rugosimeter resting on stakes embedded in the ground at the same height in relation to the terrain level and in the same place for all of the measurements. The random roughness index (RR) was calculated by the method proposed by Kamphorst et al. (2000), without logarithmically transforming the data and without eliminating the extreme values.
Soil coverages by crop residues in the HCR and HCV treatments were determined by the line transect method. The dry mass of aerial parts of the crops was determined by collecting from a single random point in the plot in a 0.24 m 2 area, and later drying for 72 h at 70 °C.
Samples of the surface flow were collected every five minutes after the beginning of runoff, for later determinations of the instantaneous runoff rates and water and soil losses. The total mass of soil lost by erosion was obtained by integrating the instantaneous rates of soil loss and extrapolating to a hectare. The data for total water and soil loss were adjusted to the planned rainfall intensity of 65 mm h -1 . Soil loss data were also adjusted for terrain slope, following methods proposed by Wischmeier & Smith (1978), while standardizing them by the mean slope of the experimental plots, 0.134 m m -1 .
The rates of water infiltration into the soil were calculated by the differences between the rates of rainfall and water loss and calculating the average of the last three readings. The flow velocity was measured with methylene blue dye (2%) at 70 min after the beginning of simulated rainfall (Ramos et al., 2016).
The estimated initial of dry mass of roots was found by ratio of dry mass of root and dry mass of residues production in the same location, according to research done by Wolschick et al. (2016). The temporal dry mass of the root was estimated based on the relative temporal reduction of the root mass determined by Volk & Cogo (2008).
The results were compared using the SAS version 9.2 statistical package, with a mixed model procedure, and the Akaike information criterion was used to choose the variance and covariance matrices. The means were compared by Fischer LSD test at p ≤ 0.05.
The soil losses, water losses and flow velocities were related to the soil subsurface and surface attributes studied.
To determine the degrees of association between the attributes and to identify the possible influence of each attribute on erosion, first Pearson's correlation analysis was carried out and then path coefficient analysis was performed using the GENES program version 6.1. This type of analysis allows the clarification of direct and indirect effects of the simple correlation coefficients (Lima et al., 2014).
The path analysis was separated between the treatments with high soil coverage by crop residues (HCR and HCV) and those with high surface roughness (HRR, HRV and BHR). Thus, the consistency of the analysis was verified by a set of variables that had weak collinearity, following the classification of Montgomery & Pek (1981). Therefore, the selection of the data was performed after an analysis of multicollinearity performed by GENES, version 6.1. When collinearity was identified, the analysis was corrected using path coefficient analysis with a collinearity constant (k).

Results and Discussion
The ryegrass aerial part dry mass (RM) in the HCR provided 17 percentage points higher soil cover (SC) than that of vetch in the HCV, in the average research period (Table 1). This difference in coverage was a consequence of the greater mass and slower decomposition of the ryegrass. In addition, the ryegrass residue was composed of thin stalks that were lightweight and more numerous than those of the vetch were, and were therefore a more effective soil cover than the vetch stalks were.
The lowest random surface roughness (RR) value occurred in the uncultured and scarified soil (BHR, 9 mm); this was low in relation to the average of the scarified treatments containing roots (Table 1). This is justified by the low GMD (stability of aggregates in water) and because of the lack of cultivation, which reduced soil resistance to disaggregation in the BHR (Zangiski et al., 2018).
Total soil losses (Table 1) were lower in the HCR treatment (2.93 Mg ha -1 ). In the HCV, the soil losses (SL) were around 10 times greater than in the HCR, demonstrating the efficacy of ryegrass residue. Among the scarified and cultivated soils, the HRR treatment had a lower SL (33.57 Mg ha -1 ) than the HRV (71.62 Mg ha -1 ), even without the RR difference, which was explained by the better soil aggregation due to the thin grass roots as demonstrated by Baets & Poesen (2010). The highest SL values occurred in the HRV and BHR, with respective values of 71.62 and 72.27 Mg ha -1 , due to low and/or non-existent soil cover in these treatments as well as the low GMD, which agreed with Laufer et al. (2016).
The water losses (WL) were equivalent to 50% of the volume of rainfall applied in the HCR treatment, which were justified by the higher surface cover and better macroporosity. In the other treatments, WLs were statistically the same (Table  1). According to Bertol et al. (2015), the Inceptisols have a low water infiltration capacity. The results of soil cover (SC), random roughness (RR), soil loss (SL) and water loss (WL) are more detailed in Ramos et al. (2016).
The lowest runoff velocity (RV) occurred in the HCR treatment (Table 1), and the highest value occurred in the BHR (0.190 m s -1 ) (Table 1), as similar to results obtained by Ramos et al. (2016). The smaller RV in the HCR could be explained by the formation of physical barriers made from ryegrass residue, which increased the flow tortuosity, and agreed with results founded by Laufer et al. (2016) who compared the effects of strip tillage (ST), full-width reduced tillage (RT) and intensive HCR and HCV -Soil cropped and covered by ryegrass and vetch residue, respectively; HRR and HRV -Chiselled soil after ryegrass and vetch crop, respectively, removing aboveground residues and keeping only the root; BHR -Uncropped and scarified soil; RM -Residue mass; SC -Soil cover; RR -Random surface roughness; SL -Soil losses; WL -Water losses; RV -Runoff velocity; WC01 and WC12 -Soil water content in the layers 0-0.10 and 0.10-0.20 m; TOC -Total organic carbon; GMD -Geometric mean diameter; BD -Soil bulk density; Mi -Micropores; Ma -Macropores; CWI -Constant water infiltration rate. Averages followed by the same letter in the column do not differ by Fischer LSD test at p ≤ 0.05 Table 1. Mean values of variables after eight simulated rainfalls in each treatment and total soil losses, in an Inceptisol tillage (IT). The cultural residues over the soil (HCR and HCV treatments) reduced on average RV by 38% in relation to the roughness treatments (HRR and HRV).
The total organic C (TOC) was higher for the HCR treatment in the 0-0.025 m soil layer, with 2.89%, followed by HCV and HRR treatments (statistically non significant) than HRV and BHR (Table 1). The highest organic C content in the HCR was due to the non-ploughing of the soil, reflecting the positive effect of the cultural residue on the surface and of the roots under the soil surface (Laufer et al., 2016). Similar to the organic C, the highest GMD occurred in the HCR treatment (4.47 mm), demonstrating the most pronounced beneficial effect of the grass roots (ryegrass).
The mean BD of 1.25 g cm -3 was below the critical limit for this soil textural class. The lowest soil bulk density (BD), with 1.20 g cm -3 , occurred for the BHR treatment in the 0-0.025 m layer (Table 1), followed by the HRR, HCR and HRV treatments, with values of 1.23; 1.24; 1.26 g cm -3 , respectively. On the other hand, with a value of 1.30 g cm -3 , the BD was higher in HCV. Soil microporosity (Mi) varied among the treatments with the higher value in the BHR, and lower value in the HCR, which are in agreement with the results found by Andrade et al. (2010). For the macroporosity (Ma), the highest value (20%) occurred in the HCR treatment in the 0-0.025 m soil layer, followed by the values for HRR, HRV, BHR and HCV, with 18, 17, 16 and 15%, respectively (Table 1).
The constant rate of soil water infiltration (CWI) was higher for the HCR treatment, with a value of 16 mm h -1 and was different from the other treatments that had an average of 10.8 mm h -1 (Table 1), which was consistent with what had been observed by Andrade et al. (2010). Tang et al. (2019) also found a positive effect of the cover and root system on the CWI compared to bare soil at hill and gully sites.
The lowest water content in the 0-0.10 m soil layer (Table 1) occurred for the HRV treatment (0.21 kg kg -1 ). In the 0.10-0.20 m layer, the lowest water content occurred for the HRV with 0.27 kg kg -1 , and the highest occurred in the BHR with 0.34 kg kg -1 .
For treatments with high soil cover (SC), the variables explained 72% (R 2 ) of the variations in the SL, while for the treatments with high surface roughness (RR), explained 53%, with residual effects of 0.53 and 0.69, respectively (Table 2).
Among the variables that were used for the path analysis involving the high soil cover treatments, SC was the only variable that correlated directly with SL and which had a higher coefficient than the residual effect. The possibility of interpreting the effects of these attributes directly was discarded, as was argued by Lima et al. (2014).
Indirectly, in the treatments with high SC, the variable that most influenced SL was SC. SC had a high capacity to dissipate the kinetic energy of the raindrops and, in part, that of the runoff. However, its capacity depended on the amount and type of residue. In this research, cultural residues were shown to be more relevant than the natural resistance to disaggregation and transport of the soil.
For the treatments with high surface roughness (RR), no variable had a significant direct effect, which reflected the high residual effect, although it is worth mentioning that there was a direct effect of -0.57 on the organic C; in the opossite, RR WC01 and WC12 -Soil water content in the layers 0-0.10 and 0.10-0.20 m, respectively; TOC -Total organic carbon; GMD -Geometric mean diameter; BD -Soil bulk density; Mi -Micropores; Ma -Macropores; CWI -Constant water infiltration rate; SC -Soil cover; RR -Random surface roughness; ERM -Estimated root mass; r -Pearson's correlation Table 2. Path coefficient analysis between the basic soil loss (SL) variable with the explanatory variables studied, for the treatments with soil cover (HCR and HCV) and surface roughness (HRR, HRV and BHR) in an Inceptisol did not present a direct correlation with SL. In the treatments with high RR, the variables with strong indirect effects were organic C, Ma, RR and estimated root mass (ERM). This shows that, in soil with high roughness, the residual effect of soil management was important for the ability of the soil to resist hydric erosion (Karlen & Stott, 1994).
For the treatments with high SC, the variables explained 83% (R 2 ) of the WL variations, whereas for the treatments with high RR, the variables explained 40%, with respective residual effects of 0.41 and 0.77 (Table 3). Among the variables analysed for high SC treatments, the only characteristics that showed a significant direct correlation (values higher than the residual effect) were the CWI and the SC.
For treatments with high RR, no variable had a significant direct effect, although different roughness values could influence the WL, due to the high residual effect. However, it is worth mentioning there was a direct effect of -0.46 on the CWI for these treatments, corroborated with the treatments with high SC.
Indirectly, the variables that most influenced the WL in the treatments with high SC were Ma, CWI, SC and ERM. For the treatments with high RR, the variables that most indirectly influenced were organic C, Ma, CWI and ERM.
The water loss (WL) showed a stronger relationship with CWI and SC than with the other variables, which agreed with results found by Laufer et al. (2016). CWI and SC contributed directly and indirectly to the specific capacity of the soil to facilitate the influx of water into the soil. However, the direct relationship between WL and Ma resulted in coefficients of -0.25 and -0.19, respectively, for treatments with high SC and RR additionally was indirectly influenced by Ma more than the other variables. Thus, even though it was less important than CWI, Ma could be considered an indicator for the facilitation of water influx into the soil.
The path coefficient analysis related to the runoff velocity of flow (RV) are shown in Table 4. The path coefficient analysis explained 92% (R 2 ) of the variations in the RV for the treatments with high soil cover, and 78% (R 2 ) for treatments with high surface roughness, with residual effects of 0.28 and 0.47, respectively.
In the treatments with high SC, only the SC and the BD had a significant direct correlation with the RV, with coefficients of -0.91 and 0.28, respectively. Indirectly, the soil variables BD, SC and ERM had greater effects on the RV than the others ones. For the treatments with high RR, there was no significant direct effect due to its high residual effect, although it was possible to highlight the direct effect of -0.41 on the organic C and of -0.31 on the RR with the RV. The variables WC01, TOC, Ma and RR had the greatest indirect effects over the other variables in the correlation with the RV.
In general, the residual effect of the path coefficient analysis was high, especially in the treatments with high RR, eliminating the possibility of interpretation of the effects of these attributes (Lima et al., 2014;Tian et al., 2018).
For SL, path coefficient analysis showed the importance of the SC as an indicator of soil quality, while the RR did not have the same importance. Partly, this is explained because only linear relationships occur in correlation analyses and in the path coefficient analysis. In this sense, Ramos et al. (2016) observed that SL correlated exponentially with the RR, and a WC01 and WC12 -Soil water content in the layers 0-0.10 and 0.10-0.20 m, respectively; TOC -Total organic carbon; GMD -Geometric mean diameter; BD -Soil bulk density; Mi -Micropores; Ma -Macropores; CWI -Constant water infiltration rate; SC -Soil cover; RR -Random surface roughness; ERM -Estimated root mass; r -Pearson's correlation Table 3. Path coefficient analysis between the basic water loss (WL) variable with the explanatory variables studied, for the treatments with soil cover (HCR and HCV) and surface roughness (HRR, HRV and BHR), in an Inceptisol Table 4. Path coefficient analysis between the runoff velocity of flow (RV) variable with the explanatory variables studied, for the treatments with soil cover (HCR and HCV) and surface roughness (HRR, HRV and BHR), in an Inceptisol WC01 and WC12 -Soil water content in the layers 0-0.10 and 0.10-0.20 m, respectively; TOC -Total organic carbon; GMD -Geometric mean diameter; BD -Soil bulk density; Mi -Micropores; Ma -Macropores; CWI -Constant water infiltration rate; SC -Soil cover; RR -Random surface roughness; ERM -Estimated root mass; r -Pearson's correlation 2-mm variation in the RR index was reflected by SL variations above 10 Mg ha -1 , which agreed with Zangiski et al. (2018).
The random surface roughness (RR) had an indirect effect, showing a negative relationship with the organic C, GMD, Ma and ERM, whereas with BD the relation was positive. Tian et al. (2018) found that, according to path analysis, the soil bulk density was a direct and indirect factor that affected runoff and sediment yields in shelter forests.
The soil loss (SL) was strongly influenced by the organic C in the treatments with high RR, while also showing the indirect effects of TOC, GMD and ERM. This demonstrates that the previous use of the soil is important for soil quality even in conditions with high roughness. Aggregate stability has been one of the most used indicators to evaluate soil physical quality and its resistance to erosion (Martínez-Trinidad et al., 2012;Moncada et al., 2013;Bünemann et al., 2018). On the other hand, in a soil with high residue cover, the mass of residues and the type of crop were more important.
For the reduction of RV, SC was both directly and indirectly important. This shows that, in the presence of high SC, the other soil attributes were less important for RV reduction. The difference was mainly due to the effect of the management and type of residue, as also verified by Ramos et al. (2016).
The random surface roughness (RR) indirectly influenced the other variables, even though it had no direct effect on RV. This result disagreed with that observed by Bertol et al. (2010), where RR exerted a strong influence over RV.
The organic C indirectly influenced the other variables in the treatments with high RR, demonstrating that the previous soil cultivation influenced the RR and decreased the RV. Bünemann et al. (2018) observed that the total organic matter/carbon was the most usual chemical attribute studied in soil quality research. This residual effect affected GMD and the resistance of soil to groove formation and, consequently, RV. RV was moderately correlated with ERM. However, was not direct, due to the indirect effect of mainly organic C and RR. The high correlation of the RV with ERM, organic C and RR probably masked the indicator values of these variables.

Conclusions
1. Ryegrass residue maintained on the soil surface contributed to higher soil organic carbon and, in addition, a greater macropore volume, which resulted in greater soil water infiltration when compared to the other treatments.
2. The soil cover, water infiltration rate and organic carbon were closely related to soil and water losses as well as runoff velocity and were therefore indicators of soil quality related to water erosion resistance.
3. In a soil with high surface roughness, the organic carbon, geometric mean diameter of the aggregates, root mass and macroporosity of the soil were more important soil quality indicators for water erosion resistance than they were for a soil with high surface cover of cultural residues.