Modeling Soil Erodibility by Water (Rainfall/Irrigation) on Tillage and No-Tillage Plots of a Helianthus Field Utilizing Soil Analysis, Precision Agriculture, GIS, and Kriging Geostatistics ^†

Abstract

The aim of our study is the modeling at the field level of the soil erodibility (K factor) by water (rainfall and irrigation) on traditional tillage (CoTl) and no-tillage (NoTl) plots cultivated with Helianthus annuus utilizing plot observations, soil sampling laboratory analyses, GIS, precision agriculture (PA), and Kriging geostatistical modeling. A split-plot layout consisting of four handlings × three replicates of trial blocks (with a southeast facing 7.5% slope) was used. Grid template surface soil core (0.0–5.0 cm) samples were taken to characterize the textures (sandy, silty, clayey, very fine sandy, and gravelly), organic matter concentrations, and the soil’s microstructure and water permeability categories. One GPS satellite tracker system was utilized to define the sampled positions, and 40 soil cores were air-dried and sieved with a 2 mm sieve to identify the soil’s mechanical microtexture using the Bouyoucos methodology. The organic matter was extracted by chemical oxidation with 1 mol L⁻¹ K₂Cr₂O₇ and titration of the remaining reagent with 0.5 mol L⁻¹ FeSO₄. The soil microstructure and permeability categories were defined following the USDA classification system. The soil erodibility by water modeling of K (Mg·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹) was derived according to the Wischmeier nomographic method by incorporating it into a developed GIS geospatial model using Kriging geostatistics. The statistical results of the ANOVA test (p = 0.05) among the soil erodibility datasets showed significant differences between the two tillage systems, as well as between the four management treatments. Moreover, it was found that the no-tillage (NoTl) plots and the treatment of no tillage plus vegetative coverage were the best tillage and agricultural practices for hillslope farm fields and can be considered environmentally friendly farming methods to curb soil erodibility by water, reduce runoff hazard, and maintain the soil’s environment and its beneficial nutrients.

1. Introduction

The erosion of soil is the phenomenon of soil particles being separated and transported by water or wind [1]. Today, it is a major issue for agricultural growth and food safety at regional, country, and world levels [2,3]. Greece has a developed agricultural sector with a declining farmer population and heavy farming operations that have led to the increased erosion of soils. In an effort to study new ways to decrease soil erosion and preserve precious soil reserves, several erosion models have been successfully deployed and widely tested all over the world. Soil erosion and hazards are considered major problems of the environment in Greece. The soil’s erodibility (K) is a fundamental parameter in erosion forecasting methods such as the USLE (Universal Soil Loss Equation) [4] and the RUSLE (Revised USLE) [5,6]. The K factor is a complicated soil attribute, which is the ease with which the soil is degraded by waterdrop splashing during rain or irrigation (mainly by sprinklers or waterjets), water runoff, or their combination [3]. The capturing of erosion’s principal variable (K factor) in forecasting modeling has proved to be a difficult task [7]. To overcome this issue, implicit methods are used to assess the K factor and allow these studies to be carried out [8]. The aim of our study is the geospatial modeling at field level of the soil erodibility by waterdrops on traditional tillage (CoTl) and no-tillage (NoTl) Helianthus plots utilizing observations, soil laboratory analyses, precision agriculture, Kriging geostatistics, and GIS mapping under climate change in Greece.

2. Materials and Methods

2.1. Study Area and Site Description

The trial was carried out in the agricultural hilly erosion-prone area of the Gaiopolis University Campus of the University of Thessaly (Larissa, Central Greece). The region enjoys moderate continental climatic conditions with a hot arid summer and a gentle winter that is characterized as Csa (Koeppen climate classification) [2] and is further classified as a XERIC MOISTURE REGIME [9], with an average annual temperature and precipitation of 17.35 °C and 380.75 mm, respectively. The highest and lowest average monthly precipitation were pr(hi) = 113.40 mm (May) and pr(low) = 12.20 mm (November), respectively. The cumulative precipitation was 652.40 mm year⁻¹. A split-plot layout consisting of 4 handlings (treats) × 3 replicates of trial blocks (with a southeast facing 7.5% slope) was used. Helianthus annuus plants were seeded to facilitate plant coverage in a number of treatments: (a) the A-treatment was traditional tillage (CoTl) plus vegetative coverage (VCov), (b) the B-treatment was CoTl with no vegetative coverage (NoVCov), (c) the C-treatment was no tillage (NoTl) plus vegetative coverage (VCov), and (d) the D-treatment was no tillage (NoTl) with no vegetative coverage (NoVCov). The dimensions of the 12 trial plots were 6 m × 22.1 m downslope, with an overall plot area of 1591.2 m².

2.2. Soil Sampling, Laboratory Analyses, and Classification

Grid template surface soil core (0.0–5.0 cm) samples were taken to characterize the textures (sandy (Sa), silty (Si), clayey (Cl), very fine sandy (vfSa), and gravelly (Gra)), organic matter (OrM) concentrations, and the soil microstructure plus water permeability categories. One GPS (Global Positioning System) satellite tracker system was utilized to define all the sampled positions, and 40 surface soil cores were air-dried and sieved with a 2 mm sieve to identify the soil’s mechanical microtexture using the Bouyoucos methodology [10,11]. The organic matter was extracted by chemical oxidation with 1 mol L⁻¹ K₂Cr₂O₇ and titration of the remaining reagent with 0.5 mol L⁻¹ FeSO₄ [11]. The soil microstructure (that is the assemblage of soil particulates and agglomerates via identifiable particles or granules) categories [9] and water permeability categories were defined following the USDA classification system [9,12]. The soil erodibility by water modeling of the K factor (Mg·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹) was derived according to the Wischmeier nomographic method [4,12,13,14], by incorporating it into a developed GIS geospatial model using Kriging geostatistics. The K factor Equation (1) [4,12,13,14] was derived for all soils consisting of less than 70% silt plus vfSa:

$$K=\left[\frac{\left(2.1\times {10}^{-4} \left(12-OM\right)M^{1.14}+3.25\left(S-2\right)+2.5\left(P-3\right)\right)}{100}\right]\times 0.1317$$

(1)

where K = soil erodibility of the USLE method (Mg·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹), M = product of the percentage of silt plus vfSa and the other soil components except clay (0.002 mm > clay, 0.05 mm > silt > 0.002 mm, and 0.1 mm > sand > 0.05 mm), OrM = soil organic material concentration (%), S = soil microstructure category, and P = water permeability category.

2.3. Statistical and Geostatistical Data Analysis, Soil Erodibility Modeling, and Methodology

Data analysis was performed using the IBM SPSS v.26 [15,16,17,18,19,20,21] statistics software package. The outputs showed the means of the soil samples analyses and field metrics. The ANOVA (analysis of variance) statistic test [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29] was applied to evaluate the tillage systems and the treatment effectiveness. The Levene statistics test for the homogeneity of variants [14,15,16,17,18,19,20,21,22] was employed to verify the hypothesis of variance equality for the soil K datasets. The LSD test [14,15,16,17,18,19,20,21,22] was used to separate means when significantly different outputs (p = 0.05) among treatments were obtained. In the present study, we used geostatistics (Kriging) and precision agriculture [14,16,17,18,19,21,22,23,27] for modeling and Geographical Information System (GIS) mapping of the soil’s textural class, organic matter concentration, and soil structure and permeability categories, respectively, as well as the soil erodibility. A GIS interface (ArcGIS © version 10.2) was employed to treat, store, and model the entry data variables in order to produce a soil erodibility digital map by using geospatial analytics and PA. The evaluation of the K factor requires residual errors’ analytics among the projected and observable values and the forecast identification range of overestimates and underestimates. For this purpose, we applied the statistic parameters reported earlier in other studies [14,15,16,17,18,19,21,22,23,27,30,31,32], such as the equations for the Mean Prediction Error (MPE), the Root Mean Square Error (RMSE), the Mean Standardized Prediction Error (MSPE), and the Root Mean Square Standardized Error (RMSSE). The soil erodibility modeling outputs of the plots were used in order to extract the K data for the validation procedure on the basis of the soil’s K trained and evaluation datasets.

3. Results and Discussion

Soil erodibility is a function of four parameters: the soil’s texture, structure, permeability, and the OrM concentration. The soil analysis outputs showed that sand with very fine sand had the ranges of 26.47–46.34% and 21.73–22.08%, respectively. The mean silt and clay contents were 19.91% and 20.22%, respectively. The soil’s organic matter [14,17,18,19,21,22,23,27] modeling results are depicted in Figure 1a–c.

Its concentration classes ranged from 1.44% to 3.22% (Figure 1b), indicating the soil’s OrM had medium to high content. The soil’s organic matter geospatial analysis showed that 34.887% of the soil plots’ area had medium OrM content (1.44–2.00%), while the remaining 65.113% had high OrM content (2.00–3.22%). The modeling and statistical outputs revealed that the K factor over the measuring time span ranged from a min 0.025 to a max 0.043 Mg·ha·h·ha⁻¹·MJ⁻¹·mm⁻¹ (average K = 0.034, standard deviation s = 0.0062). The soil characteristics of the Helianthus plots were sampled, analyzed, and digitized in accordance with their GPS-located field positions using the WGS 1984 geographic coordinate system (CS) and stored in a geodatabase. The soil parameters, tillage, and treatment datasets were projected to the UTM Zone 34N CS (Greece’s zone). The outputs of the geospatial erodibility modeling are visualized in a digital GIS map of the field in Figure 2a–c. Furthermore, the outcomes of the erodibility categories in relation to the percentage of the K factor area are illustrated in Figure 2b. The validation of the geospatial soil erodibility modeling (Figure 2c) resulted in the following geostatistical outcomes: mean prediction error (MPE) = −0.000000924, root mean square error (RMSE) = 0.00598019, mean standardized prediction error (MSPE) = −0.00518898, and root mean square standard error (RMSSE) = 1.0498154. These results are highly acceptable considering that the MPE, RMSE, and MSPE scores should be close to zero for an optimized forecast, and the RMSSE scores should be close to unity, suggesting an accurate estimate of the forecast variability. The abovementioned results confirmed the reliability and accuracy of the generated soil erodibility digital GIS map for the trial hillslope field of Helianthus annuus. Furthermore, these outcomes have proven that the ordinary Kriging exponential model demonstrated a good performance and is regarded as highly appropriate for geospatial modeling and mapping of the K factor as well as other soil parameters (clay, sand, silt, OrM, very fine sand, etc.). The output of the ANOVA test (p = 0.05) among the soil erodibility datasets in relation to the tillage method showed that the two tillage systems [traditional (CoTl) and no tillage (NoTl)] differed significantly in certain ways; so, it was necessary to further investigate the pattern of their differences. Therefore, in order to validate the equality hypothesis of variance for the erodibility dataset, the Levene statistical test for homogeneity of variances was conducted.

The findings of the Levene statistics for the soil erodibility in the tillage systems and treatments showed the variations in homogeneity of the K factor between the tillage systems (CoTl and NoTl) and also between the treatments’ (A, B, C, and D) datasets were not significantly different, which means that the hypothesis of equality of variation was confirmed. The Levene hypothesis was found to be true; so, ANOVA and LSD (Least Significant Differences) statistics were performed to evaluate the treatment effects and the mean separation of treatments. The optimum tillage system in Central Greece for hillside plots of high erosion hazard with a downward slope ≥ 7.5% was proven to be the NoTl system. The ANOVA results (p = 0.05) revealed that the soil datasets of the erodibility treatments (A, B, C, and D) significantly differed (Sig. = 0.029). The optimum treatment for limiting soil erodibility (K-factor) and maintaining a healthy soil environment was judged to be treatment C [(NoTl-VCov) (no tillage plus vegetative coverage)] for hillside plots with a high erosion hazard with a downward slope of ≥ 7.5%.

4. Conclusions

The prediction errors’ outcome of the validation of the geospatial and geostatistical modeling for the GIS soil erodibility mapping proved the validity and accuracy of the generated K-factor GIS digital map of the Helianthus annuus tested plots. All these results have demonstrated that the ordinary Kriging exponential model had good performance and is regarded as very well suited for modeling soil erodibility and many other soil parameters (clay, sand, silt, OrM, very fine sand, etc.) and digital mapping. In consideration of the ANOVA test results of the tillage systems and the treatment effects on the soil erodibility, the optimum tillage system found was NoTl (no tillage) and the superior treatment was C [(NoTi-VCov) (no tillage plus vegetative coverage)] for hillside plots of high erosion hazard with a downward slope ≥ 7.5%. These can be considered as environmentally friendly farming methods to curb soil erodibility by water, reduce runoff hazard, and maintain the health of the soil’s environment and its beneficial nutrients.