Evaluating the Efficiency of Different Regression, Decision Tree, and Bayesian Machine Learning Algorithms in Spatial Piping Erosion Susceptibility Using ALOS/PALSAR Data

Abstract

Piping erosion is one form of water erosion that leads to significant changes in the landscape and environmental degradation. In the present study, we evaluated piping erosion modeling in the Zarandieh watershed of Markazi province in Iran based on random forest (RF), support vector machine (SVM), and Bayesian generalized linear models (Bayesian GLM) machine learning algorithms. For this goal, due to the importance of various geo-environmental and soil properties in the evolution and creation of piping erosion, 18 variables were considered for modeling the piping erosion susceptibility in the Zarandieh watershed. A total of 152 points of piping erosion were recognized in the study area that were divided into training (70%) and validation (30%) for modeling. The area under curve (AUC) was used to assess the effeciency of the RF, SVM, and Bayesian GLM. Piping erosion susceptibility results indicated that all three RF, SVM, and Bayesian GLM models had high efficiency in the testing step, such as the AUC shown with values of 0.9 for RF, 0.88 for SVM, and 0.87 for Bayesian GLM. Altitude, pH, and bulk density were the variables that had the greatest influence on the piping erosion susceptibility in the Zarandieh watershed. This result indicates that geo-environmental and soil chemical variables are accountable for the expansion of piping erosion in the Zarandieh watershed.

1. Introduction

Poesen et al. [1] shed light on the importance of gully erosion for the proper land management and evolution of the geomorphology. They were first to discuss the relevance of modeling to properly understand the nature and fate of gully and piping erosion. Another key step in the understanding of the evolution of gully and piping erosion was the contribution of Vincent Chaplot [2] when he updated the importance of the topography and parent material (and soil types) to research gully and piping evolution. This is part of a dynamic experimental and modeling approach to badland, gully, and piping landscapes [3,4].

Soil is one of the essential natural resources. From the deforestation in the middle ages to the industrial revolution by anthropic factors, soil has been under several natural and man-made threats [5]. Agriculture, grazing [6], poor land management [7], and biological factors such as water and wind erosion frequently damage the soil [8,9]. Soil is a relevant component of the Earth system as it regulates the hydrological, biochemical, and life cycles, making it a definitive element in the achievement of the Sustainable Development Goals of the United Nations and the Land Degradation Neutrality challenge [10]. Piping erosion causes significant changes in the landscape and environmental degradation [11]. Piping is a subsurface drainage system that is often caused by dissolution but can also result from animal burrowing and nesting, and, if formed on a slope, with the collapse of the tunnel roof, gullies develop and eventually lead to badland formation and expansion [12,13]. Although piping erosion is not as widespread as agriculture erosion [10], it has recently become an environmental challenge and public issue in several communities worldwide [3].

Piping generally occurs in formations with low permeability and high dissolution minerals. With the absorption of water by the formation in the wet season as well as the presence of elements such as sodium that influence the hydrolysis of the clay, the clay expands and, during dry periods, numerous cracks develop and allow the percolation of the surface runoff. It should be noted, however, that tunnels are formed upon crack development [14]. Despite the simple physical structure of these erosional forms, their mechanism of formation is very complicated. Soil, physiographic, climatic, and biological factors all contribute to the development of piping erosion, and they have different degrees of importance in different areas and conditions [15].

Piping erosion results from the concentration of surface runoff and its uneven infiltration into the soil, the existence of a lower horizon with less permeability, and the creation of a hydraulic gradient [16]. Researchers have not been able to determine precisely which of these factors has the highest contribution to the development of piping erosion due to the formation pipes in a wide range of soil textures, different climates, fauna, different land use classes, and wide slopes [17,18,19,20]. Piping erosion has led to unusual changes in the landscape, ultimately leading to inadequate agricultural practices and the inefficient use of natural resources [12,21]. Therefore, to control this type of erosion, it is best to take protective measures before it begins. The formation and onset of piping erosion is influenced by environmental factors. To design and implement soil conservation decisions and strategies, it is important to identify areas that are susceptible to piping erosion and the likelihood of piping erosion to begin and to take appropriate precautions in each area [22].

In areas where the lithological and soil structure is composed of fine and discrete elements, the mechanism of the formation of piping erosion is such that this erosion causes washing and removal of fine particles by mechanical leaching and creates the conditions for tunnel erosion problems. Moreover, in areas with high solutes, water infiltration and solute dissolution accelerate the formation of piping erosion forms, or, in other words, dissolution plays a supportive role in creating these forms [23,24].

The climatic conditions of Iran and inadequate land management has led to accelerated land degradation, which has recently resulted in high soil erosion rates and sediment yield [25]. Iran shows active soil erosion rates due to the land uses [26], and climate and land management are the key factors [27]. Many authors focused on the evaluation of the soil erosion rates to quantify them [28]. The soil erosion studies in Iran focused on different factors and scales, from watershed scale soil erosion behavior [29] till the control that some factors such as texture exert [30], to gully erosion [31], soil detachment [32], or the use of the Universal Soil Loss Equation [33]. Although the research on soil erosion in Iran is a key environmental issue, piping erosion is a topic which was not addressed until recently [24] and literature about it is scarce. Although piping is relevant in Iran and other semiarid lands in the world [14,34,35], there is a need to deeply understand the spatial distribution, the origin and development, and the proper management in order to prevent it. Mapping piping and researching the control that different factors exert on piping development is an ongoing challenge for today’s scientists, as Bernatek-Jakiel and Poesen [14] highlighted in their article on the current state of the art in this field.

All the above-mentioned constraints in the scientific research of piping can be addressed with machine learning models, which are already widely used in predicting various phenomena such as gully, landslide, and land subsidence [31,36]. Machine learning models have the excellent ability to identify the occurrence behavior of phenomena in terms of the use of distributed estimation algorithms, data-driven nature, and high iteration of the modeling process. In numerous investigations, different machine learning models have been used to model zoning erosion phenomena, such as support vector machine (SVM) [37,38,39,40], random forest (RF) [31,41,42,43,44,45,46], naive Bayes [42,47], multivariate adaptive regression splines (MARS) [41,48,49], artificial neural network (ANN) [50], and logistic regression (LR) [51,52,53].

The scientific community is looking for proper and comprehensive methods to create susceptibility maps. For this, we need pilot areas (representative watersheds and slopes) on which to test the new strategies to forecast piping development. Various studies have been carried out on the use of a machine learning method in mapping sensitivity and risk phenomena [54]. In the study area, despite the extensive expansion of piping erosion and the susceptibility of the area to surface wash, no research has been conducted to determine the potential of piping erosion and the factors that cause piping formation and development. Therefore, in this study, we evaluate the factors affecting piping erosion and finally piping erosion susceptibility mapping. In this study, the performance of random forest (RF), support vector machine (SVM), and generalized linear Bayesian models (GLM Bayesian) in piping erosion hazard susceptibility using Advance Land Observatin Satellite/Phased Array type L-band Synthetic Aperture Radar (ALOS/PALSAR) digital elevation model (DEM) was evaluated.

2. Materials and Methods

2.1. Description of the Study Area

The Zarandieh watershed is located in the northern area of the Markazi province in Iran. This watershed is one of the sub-catchments of the Shoor river that is located at 35°10′ to 35°38′ N latitude and 49°40′ to 51°04′ E longitude, respectively (Figure 1). The catchment area is around 3592 km², with an average elevation of 1974 m and an average slope of 3.2%. According to the climatic classification of the Zarandieh watershed, it is located in the arid and cold arid climates of the desert, with an average rainfall of 343 mm, with the highest rainfall occurring in winter. Land uses in the area include agricultural, bare land, rangeland, and residential land; most of the study area is rangeland. Most of the piping affected land is located in the eastern and center regions of the watershed, which are mainly agricultural and bare land. The soil of the Zarandieh watershed is deep, saline, and alkaline and, at the surface, saline. The most important limitations in this area are soil texture, salinity, alkalinity, and lack of drainage. The soils are generally deep and sometimes semi-shallow to calcareous and calcareous and articular (Administration of Natural Resources, Markazi Province https://markazi.frw.ir/00/Fa/default.aspx). The soils of the Zarandieh watershed are prone to erosion due to climatic conditions and the presence of highly saline soils, especially piping erosion due to the low slope of the region, and this erosion is one of the main problems in the region. Examples of piping erosion in the study area are shown in Figure 2.

2.2. Methods

The present study illustrates the usage of various machine learning methods in piping erosion susceptibility. Preparing the piping erosion susceptibility maps in the current research involved data collection, inventory map preparation, effective factors preparation, collinearity test, piping erosion susceptibility (PES) mapping using random forest (RF), support vector machine (SVM), and Bayesian generalized linear models (Bayesian GLM), factors importance analysis using the RF, and validation of the models (Figure 3). The piping erosion susceptibility maps were classified into five categories using the natural break method of Jenks [55].

2.3. Dataset Preparation for Spatial Modeling

Extensive field survey was carried out to determine the areas of piping erosion in the Zarandieh watershed in Markazi province. In this study, the location of each piping was determined using Garmin GPSMap (Garmin Ltd., Olathe, KS, USA) based on field vists. Moreover, a digital elevation model (DEM) with a spatial resolution of 12.5 m was downloaded from the Alaska Satellite Facility website (asf.alaska.edu) to investigate the location of the piping erosion susceptibility in the study area. Factors such as slope, aspect, distance from river, distance from road, topographic wetness index (TWI), topographic position index (TPI), rainfall, lithology, land use, drainage density, bulk density, pH, cation exchange capacity (CEC), silt, clay, and sand to map piping erosion susceptibility were used (Figure 3). Soil physical and chemical properties of the study area (pH, CEC, bulk density, sand, clay, and silt) were prepared from the SoilGrids site (https://soilgrids.org/) with a spatial accuracy of 250 m due to the large extent of the area and the lack of soil studies. Layers of topographic wetness index (TWI), topographic position index (TPI) were also prepared using SAGAGIS 3.2 software. Finally, 152 piping locations were recorded in the study area. In addition, 152 absence points were created using ArcGIS 10.5 software randomly (Figure 4). The piping points and the non-piping points were divided in two parts: 70% for training and 30% for validation.

2.4. Multi-Collinearity Analysis

Addressing the multicollinearity issue of the independent variables is quite necessary to understand their strength of predicting susceptibility and it can directly be depicted with the help of VIF [56]. VIF (variance inflation factor) is the second measure of multi-collinearity analysis, derived from the tolerance (reciprocal of VIF), which is the direct measure [56]. The consideration of VIF as less than 10 (for weaker models, it is less than 5) indicates no multi-collinearity problem and it is ready for further use in terms of prediction [36,50,56,57,58,59]. The equations of tolerance and VIF are as follows (Equations (1) and (2)):

Tolerance = 1 − R²J

(1)

VIF = 1/Tolerance

(2)

where R²J represents the regression coefficient of determination of explanatory variable J on all the other explanatory variables.

2.5. Machine Learning Method Used in Modeling the Piping Erosion

2.5.1. Random Forest (RF)

The random forest method is a relatively complex method among tree methods in which several decision trees are trained to increase model accuracy. The result is the prediction of a group of decision trees. In the random forest learning method, each decision tree is trained using a random sample selected from the training dataset [60]. Selecting a set of predictor variables used to segment the nodes is also done at random. In the random forest method, the mtry and ntree properties are determined for the covariates used in each subset and the number of trees used in the forest, respectively. The number of variables can range from one to the total number of covariates. The number of trees is selected by the user, usually from 500 to 1000 [60].

In the random forest algorithm, it is possible to select the most important and effective covariates. Estimation of the importance of a variable is done by queuing the values of each variable in the out-of-bag (OOB) samples, which are usually 1.3 observations, not used in the model, and by reclassifying the OOB samples using the row variable. OOB error change is used as a measure of variable importance. Variables whose deletion results in a relatively higher increase in OOB error are considered more important variables [43]. In the random forest method, each tree is formed as follows: (i) If N is the number of states in the training dataset, N state will randomly sample by pasting in the original data. This is an example of a work set for this tree; (ii) If there are M variables, consider m smaller than M such that at each node, m variables are randomly selected from M and the best separation of this m variable is used to separate the node. The value of m is assumed to be constant during forest construction, and (iii) each tree grows as large as possible. There is no pruning.

2.5.2. Support Vector Machine (SVM)

The structure of the support vector machine model was introduced by Vapnik [61] and it utilizes the principle of structural error minimization. The support vector machine models are divided into two main groups: support vector machine classification model and support vector machine regression model. Support vector machine classification models are used to solve data classification problems that are placed into different classes, and support vector regression models are used to solve forecasting problems. If the data are linearly separated, the equation is as follows (Equation (3)):

$$y=f\left(x\right)=\mathrm{sign}\left[{\displaystyle \sum }_{i=1}^n{y}_i{\alpha }_i\langle x_i,x_j\rangle +b\right]$$

(3)

If the data are not linearly separable, the samples can be pre-processed by transferring the samples to a non-linear space with high dimensionality, where the result is internal multiplication and proves that if a symmetric kernel of equation conditions Mercer application to the low-dimensional input space can be seen as the result of high-dimensional internal multiplication and greatly reduces computation [62]. In this case, Equation (3) is changed to Equation (4):

$$y=f\left(x\right)=\mathrm{sign}\left[{\displaystyle \sum }_{i=1}^n{y}_i{\alpha }_{i}K\left(x_i,x_j\right)+b\right]$$

(4)

The $K\left(x_i,x_j\right)$ function is a kernel function that produces internal multiplication to create machines with different types of nonlinear decision-making levels in the data space. Different kernels are used for the support vector machine regression model (

Table 1. Equations of different kernel functions.

Parameters	Equations	Kernel Function
-	$K\left(x_i,x_j\right)=x_i·x_j$	Linear kernel
$\gamma$ and $d$	$K\left(x_i,x_j\right)={\left(\gamma x_i·x_j+r\right)}^d$	Polynomial kernel
$\gamma$	$K\left(x_i,x_j\right)=\exp \left(-\gamma {\parallel x_i-x_j\parallel }^2\right)$	Radial basis function kernel

).

2.5.3. Bayesian Generalized Linear Models (Bayesian GLM)

Bayesian reasoning is a probability-based approach to inference. The basis of this approach is that there is a probability distribution for each quantity, which can be optimized by observing new data and arguing about its probability distribution [63,64]. Bayesian networks, also known as belief networks, belong to the family of probabilistic graphical models. These graphical structures are used to represent information in an uncertain field. Specifically, each node in the graph represents a random variable and the branches represent probabilistic dependencies between the variables. These conditional dependencies are often evaluated by specific statistical and probabilistic methods. Bayesian networks combine principles of graph theory, probability theory, computer science, and statistics. Bayesian networks provide efficient representation and computation of the probability distribution shared over a series of random variables. Besides this, Bayesian networks model the intensity of the relationship between variables quantitatively, allowing the conditional belief about them to be automatically updated by accessing new information [64,65].

The basis of Bayesian learning is Bayesian theory. This theory allows the calculation of the secondary probability based on the first probabilities:

$$P\left(h|D\right)=\frac{P\left(D|h\right)P\left(h\right)}{P\left(D\right)}$$

(5)

where $P\left(h|D\right)$ is a posterior probability, $P\left(D|h\right)$ is a likelihood, $P\left(h\right)$ is prior probability, and $P\left(D\right)$ is evidence. As can be seen, by increasing P(D), the value of P(h|D) decreases, because the higher probability of D being observed independently of h means that D has less evidence to support h.

2.6. Methods of Validation and Accuracy Assessment

In the present study, four statistical measures including sensitivity (TPR), specificity (TNR), precision (PPV), and negative predictive value (NPV) were used for evaluating the fitness of the used machine learning models. The positive and negative predictive values (PPV and NPV, respectively) are the proportions of positive and negative results in statistics and diagnostic tests that are true positive and true negative results, respectively. Greater values of these statistical tests indicate the better result of the models [43,66]. The model output was speculated by area under the curve (AUC) of receiver operating characteristic (ROC) [67].

$$\mathrm{TPR}=\frac{\mathrm{TP}}{(\mathrm{TP}+\mathrm{FN})}$$

(6)

$$\mathrm{TNR}=\frac{\mathrm{TN}}{(\mathrm{TN}+\mathrm{FP})}$$

(7)

$$\mathrm{PPV}=\frac{\mathrm{Number}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{positives}}{(\mathrm{Number}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{positives}+\mathrm{Number}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{false}\hspace{0.17em}\mathrm{positives})}$$

(8)

$$\mathrm{NPV}=\frac{\mathrm{Number}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{true}\hspace{0.17em}\mathrm{negatives}}{(\mathrm{Number}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{true}\hspace{0.17em}\mathrm{negatives}+\mathrm{Number}\hspace{0.17em}\mathrm{of}\hspace{0.17em}\mathrm{false}\hspace{0.17em}\mathrm{negatives})}$$

(9)

$$\mathrm{AUC}=\frac{{\displaystyle \sum \mathrm{TP}+{\displaystyle \sum \mathrm{TN}}}}{\mathrm{P}+\mathrm{N}}$$

(10)

Finally, for all modeling approaches and analyses, we used R 3.5.2 software. The randomForest, e1071, arm packages were used to run the RF, SVM, and Bayesian GLM models. The susceptibility values generated from RF, SVM, and Bayesian GLM models were taken to the GIS environment to prepare the PES maps.

3. Results

3.1. Multi-Collinearity Analysis

The variance inflation factor (VIF) values of the variables are shown in

Table 2. Multi-collinearity analysis to determine the linearity of the independent variables.

Row	Variables	VIF
1	Aspect	1.13
2	Altitude	2.52
3	Plan	1.73
4	Profile	1.62
5	Distance from river	1.28
6	Slope	1.79
7	TWI	1.40
8	Lithology	1.49
9	Land use	1.41
10	Rainfall	3.62
11	TPI	1.51
12	Silt	2.58
13	Sand	1.89
14	Clay	3.93
15	CEC	1.88
16	Bulk density	2.96
17	pH	2.57

. The VIF values of all the variables are less than 5, indicating their independence and readiness for use in the models. Among the seventeen variables, clay has the highest VIF (3.93) and aspect has the lowest (1.13).

3.2. Piping Erosion Susceptibility Modeling

The piping and non-piping points containing the values of piping erosion susceptibility factors (PESFs) are taken as the input datasets in R 3.5.2 software to construct the piping erosion susceptibility (PES) models. Then, the PES maps are classified into five ordinal categories, i.e., very low, low, moderate, high, and very high based on the natural break classification method (Figure 5). The PES map generated from the RF model shows that very high and high PES zones account for 11.21% and 11.93% of the watershed area, respectively, whereas moderate, low, and very low PES zones cover 14.97%, 18.95%, and 42.93% of the area, respectively (

Table 3. Piping erosion susceptibility classes’ area.

Susceptibility Class	GLM Bayesian		SVM		RF
	Area (Km2)	Area (%)	Area (Km2)	Area (%)	Area (Km2)	Area (%)
Very Low	196.67	5.48	904.50	25.18	1542.10	42.93
Low	510.74	14.22	1055.98	29.40	680.66	18.95
Moderate	876.75	24.41	694.27	19.33	537.86	14.97
High	1196.65	33.32	523.98	14.59	428.56	11.93
Very High	811.04	22.58	413.11	11.50	402.64	11.21

). The SVM model also shows quite a similar picture to RF as very high and high PES zones occupy 11.5% and 14.59% of the basin area whereas moderate, low, and very low PES zones cover 19.33%, 29.4%, and 25.18% of the area, respectively. The very low PES zone according to the SVM model is significantly lower than the result of the RF model. However, as per the Bayesian GLM output, the very high and high PES zones account for 22.58% and 33.32% of the area, respectively, which is more than half of the basin area. This model shows the decreasing percentage of areas with piping erosion probability classes like moderate, low, and very low PES zones, which cover 24.41%, 14.22%, and 5.48% of the area, respectively. All three models identified that the downstream eastern end of the watershed is very highly susceptible to piping erosion. Among the models, Bayesian GLM shows more areas of high to very high PES zones (Figure 4). On the other hand, the southwestern end of the basin falls under very low PES zones, although the SVM model shows some share of moderate PES zones there. The map generated from the RF model shows a gully-like pattern in the PES zones of the central and western areas of the watershed.

3.3. Validation of the Models

The piping erosion susceptibility (PES) maps by random forest (RF), SVM, and Bayesian GLM are justified by the receiver operating characteristics (ROC) curve using both the training and validation datasets (Figure 6 and Figure 7). As per the training datasets, the area under curve (AUC) values of RF, SVM, and Bayesian GLM are 0.98, 0.96, and 0.93, respectively (

Table 4. Predictive capability of piping models using train and test dataset.

Models	SVM		RF		GLM Bayesian
Evaluation Parameter	Test	Train	Test	Train	Test	Train
Sensitivity	0.918	0.891	0.89	0.95	0.91	0.80
Specificity	0.702	0.893	0.74	0.97	0.70	0.89
NPV	0.891	0.893	0.87	0.95	0.89	0.82
PPV	0.762	0.891	0.78	0.97	0.76	0.88
AUC	0.88	0.96	0.90	0.98	0.87	0.93

). The sensitivity (TPR) values of RF, SVM, and Bayesian GLM are 0.95, 0.891, and 0.80, respectively, which means that the areas detected as piping erosion susceptible zones by these models are very accurate. On the other hand, the specificity (TNR) values of RF, SVM, and Bayesian GLM models are 0.97, 0.893, and 0.89, respectively, which also indicates the high accuracy of the non-piping erosion susceptible zones detected by these models. The precision (PPV) values of these models are also high like RF (0.97), SVM (0.891), and Bayesian GLM (0.88), which shows the robustness of the models. Therefore, the parameters of validation for training datasets present excellent positive results in favor of these models.

Next, the validation datasets were used to examine the prediction capability of the three models. In the context of validation datasets, the AUC values of RF, SVM, and Bayesian GLM are 0.90, 0.88, and 0.87, respectively (Table 4). This suggests very good to excellent predictability of these three models. As per the ROC method, the sensitivity (TPR) values of RF, SVM, and Bayesian GLM are 0.89, 0.918, and 0.91, respectively, which means, as per the test datasets, these models also predicted correctly the piping erosion susceptible zones. The specificity (TNR) values of RF, SVM, and Bayesian GLM models are 0.74, 0.702, and 0.70, respectively, which suggests the decent prediction capability of the non-piping erosion susceptible zones also. The precision (PPV) values of these three models RF (0.78), SVM (0.762), and Bayesian GLM (0.76) indicate their decent ability of prediction. So, the results of the validation assessment suggest that the models used in the PES are very accurate and have high predictability.

Table 5. Variable importance analysis.

Variables	Importance
Aspect	1.59
Altitude	9.29
Plan curvature	0.51
Profile curvature	0.84
Distance from river	4.47
Slope	1.76
TWI	2.76
Lithology	0.78
Land use	1.29
Rain	0.18
TPI	1.00
Silt	1.46
Sand	1.65
Clay	0.90
CEC	3.99
Bulk density	6.81
pH	8.80

shows the importance of environmental factors in piping erosion modeling in the study area. Among these factors, altitude (9.29) has the highest importance for the piping susceptibility method, followed by pH, bulk density, and distance from the river (Table 5).

4. Discussion

4.1. Comparison of the Models

Piping erosion is a critical pedo-geomorphological process that contributes to land subsidence and gully formation in many parts of the world [68]. Piping erosion is more than a key topic in geomorphological research as it also affects agriculture production and triggers land degradation processes [12,15,69,70]. The subsurface dissolution of soil solutes can lead to the formation of underground channels, cavities which may trigger the sudden collapse of land or the creation of gullies and rills [71]. Therefore, the identification and proper prediction of the piping erosion susceptible zones in the watershed through mapping is crucial for a sustainable land management. Mapping can contribute to efforts to reduce the risk induced by sudden subsidence and collapses that can result in economic damages and casualties [72].

There are various methods and techniques which have been developed and executed all over the world for the prediction of piping erosion zones. The preparation of a rational and accurate susceptible map of piping erosion is a quite debatable matter for scientists. Recently, with the arrival of machine learning techniques, environmental modeling has become significantly accurate as these techniques can analyze the complex relationship between the predictors and responses like flood [73], gully erosion [74], landslide [75], and soil piping [76]. Many types of research are now being conducted to enhance the prediction performance of these models using various hybrid or ensemble algorithms. This work aims at applying machine learning models like RF, SVM, and Bayesian GLM for piping susceptible zoning of the watershed. Using sensitivity, specificity, PPV, NPV, and AUC values, the random forest (RF) model appeared as the best model in predicting the piping erosion susceptible zones. In the case of training datasets, the RF model performed best in all aspects. So, this model is best suited for the PES study and its success rate is excellent. However, in the case of sensitivity analysis of test datasets, the SVM model predicted the piping erosion zones more accurately than the other two.

The piping erosion can be addressed better with RF than SVM as per its design. RF generally uses the data as they are but SVM relies more on distance in a different linear or non-linear form, which might be the reason that RF is the best model here. Various non-linear kernels were used in SVM, which produced very good results in the end, and, in terms of predicting the future piping erosion zones, it was as the best model. Now, if we examine the results of Bayesian GLM, then we can see a similar picture to SVM. The validation values are nearly equal to SVM in all categories (Table 3) but the GLM model has produced more PES zones in high to very high categories than the other two. The sensitivity value of 0.80 and 0.91 in training and test datasets, respectively, suggests that the Bayesian GLM is also quite robust and suitable for PES assessment. One positive of using machine learning algorithms like RF, SVM, and Bayesian GLM is that they can automate the process of examining multiple databases to assemble the vital required information. These models can also identify the unique assumptions for decisions along with automatic large dataset analysis.

4.2. Variable Importance Analysis

The RF model, introduced by Breiman [77] by assimilating algorithms of decision tree [78,79], is useful for predicting and identifying independent variables’ contributions in the prediction process. Eighteen piping erosion susceptibility factors (aspect, altitude, plan curvature, profile, distance from the river, slope, TWI, lithology, land use, rainfall, TPI, silt, sand, clay, CEC, bulk density, and pH) were used in the modeling, evaluation of goodness-of-fit as well as performance of modeling processes. No variable showed a multi-collinearity problem, which indicates that all of them are independent factors and can be used in the modeling process. Among these factors, altitude (9.29) has the highest importance for the piping susceptibility method, followed by pH, bulk density, and distance from the river (Table 5). Lower elevation zones are more susceptible to piping erosion in comparison to the higher elevations, which is why we can see that very high PES zones exist in the lowest area of the watershed, in the eastern region. Bulk density has a positive relationship with the formation of pipes in the watershed, likely due to the formation of cracks on soils with high clay content, which are usually heavy soils. There is little known by other researchers about how bulk density behaves and impacts the origin and development of pipes [12,80] and our research shed light on this issue for the first time. Heavy textured soils can be seen now as a source of risk due to the development of piping via mainly the preferential flows developed. Researchers such as Jones [22] describe these as being responsible for quick catchment responses. Higher bulk density indicates higher compaction of the soil, which restricts the root growth to some extent, and this is also a factor to the development of bare land, gullies, rills, and badland areas, which is the landscape where piping is also very visible and abundant [23,81,82]. The influence of roots on piping development is relevant to understand whether plants contribute to piping development [83].

In the Zarandieh watershed, we found that the pH and CEC are also positive factors to develop pipes and this is because high rates of sodium are very common in piping areas. Other authors have shown that piping development is a physico-chemical process [84]. Our results demonstrate that the physico-chemical properties of the soil are more important to foresee the spatial distribution of piping erosion than the other parameters that were researched such as plan curvature, profile, aspect, or slope. The topographical characteristics of slopes and watersheds are not relevant to define the soil erosion losses and piping development as the soil chemical and physical properties are more relevant.

Another factor that was found relevant to foresee the risk of piping is the distance from the river. This is because the proximity to the river accelerates the process of piping erosion through the supply of sub-surface water. This process is well known in the formation of lakes and lagoons in semiarid land, where the formation of piping determines the drainage of the watersheds. A good example of this interaction between piping and groundwater is found in the paper of Castañeda et al. [85]. This has also been found as a source of risk, as was reported by Wang et al. [86], due to the formation of landslide dams as a consequence of piping erosion. Piping is a key landform to understand the rapid evolution of the landforms and our research confirms that few factors control the piping development. Masannat [87] already found in the Benson area in Arizona that overgrazing combined with the climatic conditions of long dry summers with intermittent short rain storms contributed to the initiation and development of piping erosion, and then piping triggered and developed gully erosion and badland morphologies. This is why Parker [88] already found that piping was a geomorphic agent in dry lands and that piping was a key actor in the semiarid landform development.

The study of piping erosion susceptibility in the study area shows that areas with agricultural land use are very sensitive to this type of erosion. Since the agricultural areas in the Zarandieh watershed have expanded to low slopes with soil which has a high salt content, agricultural operations in these areas have increased the infiltration of water and solute dissolution and led to the development of piping erosion. However, piping is a landform that, although widely present in arid and semiarid landscapes, is widespread around the world. It was found by Parker [89] on upland blanket peat, by Seppälä [90] and Carey and Woo [91] in the permafrost, and by Onda and Itakura [92] in rivers and temperate climatic conditions [93].

4.3. Implications and Soil Erosion Control

Our research demonstrates that it is possible to foresee where piping erosion will take place. However, to solve the problems caused by piping erosion, programs and strategies to reduce the soil losses should be developed. There is a need to control the high erosion rates found in piping areas around the world and strategies such as cover crops, mulches [94], catch crops [95], stone bonds [96], land abandonment [97], and reduced tillage should be applied on agricultural land to avoid the formation of soil erosion features such as pipes, rills, or gullies. It is on agricultural land that soil erosion problems are found around the world due to the abuse of the soil in vineyards and in persimmon, olive, citrus, and fruit tree plantations in semiarid land [9,98,99,100]. This is due to the low organic matter, which can solved with the use of organic farming strategies to restore the carbon content in the soil [101] and reduce the soil losses [102]. The above-mentioned strategies will contribute to improving the ecosystem services offered by the soil systems and to achieving sustainability challenges [94,103].

5. Conclusions

The Zarandieh watershed is very much susceptible to piping erosion and our research mapped the locations of it as well as the areas susceptible to land subsidence and gully erosion. Machine learning models (RF, SVM, and Bayesian GLM) have been used to assess the location of the piping areas. We considered 18 piping erosion conditioning factors from 304 locations. The models have shown that altitude, pH, bulk density, CEC, soil, and distance from the river are the key factors affecting piping erosion. We also noticed that the downstream area of the basin, along with fewer central areas, is highly to very highly susceptible to piping erosion. Therefore, from this study, it may be inferred that machine learning models can be used as a very accurate tool for managing piping erosion affected land.