**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

**Citation:** Filintas, A.; Gougoulias, N.; Hatzichristou, E. 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. *Environ. Sci. Proc.* **2023**, *25*, 54. https://doi.org/10.3390/ ECWS-7-14254

Academic Editor: Athanasios Loukas

Published: 16 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Correspondence: filintas@uth.gr

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<sup>−</sup>1. 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 m2.

#### *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−<sup>1</sup> K2Cr2O7 and titration of the remaining reagen<sup>t</sup> with 0.5 mol L−<sup>1</sup> FeSO4 [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−1·MJ−1·mm<sup>−</sup>1) was derived according to the Wischmeier nomographic method [4,12–14], by incorporating it into a developed GIS geospatial model using Kriging geostatistics. The *K* factor Equation (1) [4,12–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(S - 2) + 2.5(P - 3)\right)}{100}\right] \times 0.1317 \tag{1}$$

where *K* = soil erodibility of the USLE method (Mg·ha·h·ha−1·MJ−1·mm<sup>−</sup>1), *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–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–29] was applied to evaluate the tillage systems and the treatment effectiveness. The Levene statistics test for the homogeneity of variants [14–22] was employed to verify the hypothesis of variance equality for the soil *K* datasets. The LSD test [14–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–19,21–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–19,21–23,27,30–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–19,21–23,27] modeling results are depicted in Figure 1a–c.

**Figure 1.** (**a**) Modeling outcomes on a soil organic matter digital GIS map of the *Helianthus* plots. (**b**) Diagram of the *OrM* classes vs. the percentage of the *OrM* area. (**c**) Semivariogram of the model.

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−1·MJ−1·mm<sup>−</sup><sup>1</sup> (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.

**Figure 2.** (**a**) The soil erodibility modeling results on a digital GIS map of the Helianthus plots. (**b**) Diagram of erodibility categories vs. percent of *K* factor area. (**c**) Semivariogram of the model.

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%.
