Evaluation of Gully Erosion Susceptibility Using a Maximum Entropy Model in the Upper Mkhomazi River Basin in South Africa

Abstract

Soil erosion is one of the most challenging environmental issues in the world, causing unsustainable soil loss every year. In South Africa, several episodes of gully erosion have been documented and clearly linked to the presence of Quaternary colluvial deposits on the Drakensberg Mountain footslopes. The aim of this study was to identify and assess the triggering factors of gully erosion in the upper Mkhomazi River basin in KwaZulu-Natal, South Africa. We compiled a gully inventory map and applied remote sensing techniques as well as field surveys to validate the gully inventory. The gullies were subdivided into slope gullies and fluvial gullies. We derived susceptibility maps based on the spatial distribution of gully types to assess the most important driving factors. A stochastic modeling approach (MaxEnt) was applied, and the results showed two susceptibility maps within the spatial distribution of the gully erosion probability. To validate the MaxEnt model results, a subset of the existing inventory map was used. Additionally, by using areas with high susceptibilities, we were able to delineate previously unmapped colluvial deposits in the region. This predictive mapping tool can be applied to provide a theoretical basis for highlighting erosion-sensitive substrates to reduce the risk of expanding gully erosion.

1. Introduction

Soil erosion is one of the most severe environmental problems worldwide [1,2]. Soils subjected to a series of degradation processes, such as compaction, nutrient loss and a loss in water storage capacity, result in soil erosion that leads to the loss of fertile arable land [3].

Gully erosion is a significant environmental problem in arid and semi-arid regions [4,5,6,7,8], in Mediterranean countries [9,10,11,12,13,14,15,16] and in a wide range of climatic and environmental situations [17,18,19,20]. Over the past few decades, numerous studies have addressed the identification, mapping [9,21] and modeling [4,16,17,22,23,24] of gully erosion. Several parts of eastern South Africa are affected by a range of water erosion processes [1]. The hinterland of Kwa-Zulu Natal (KZN) province is characterized by wide areas of sheet (rill–interrill) erosion and has been affected by multiple episodes of dendritic gully erosion [25]. A gully is defined as a channel with steep walls caused by the removal of soil by concentrated turbulent flow of water after heavy rains [26], with dimensions that preclude remediation by tillage operations. Gully erosion represents a geo-environmental problem, causing severe soil loss, increasing the connectivity in the landscape, therefore transporting large quantities of sediment into the drainage systems [16,21,27], and are considered to be an indicator of desertification [28,29]. Furthermore, gullies can also cause damage to roads, buildings and infrastructures [16].

Given the implication of gully erosion, in order to have proper soil management, preserve the soil and mitigate gully erosion processes, especially in those areas where agriculture is the main source of income, it is essential to examine the many factors that drive gully formation [30,31,32,33,34].

Especially in South Africa, several elements contribute to the development of gully erosion, including soil type, bedrock lithology and structure, precipitation, slope angle, vegetation and land use [1,23,35,36,37,38,39]. Moreover, human activity and climate change can increase this phenomenon in many regions [16,40,41]. In order to understand the causes of gully formation as well as successfully apply land use planning strategies, the influence of the soil erosion factors need to be identified and quantified.

Therefore, the purpose of this study is to (1) differentiate the gullies of the study area based on their morphology and evolutionary development stage, (2) assess the gully triggering factors and derive susceptibility maps as well as (3) map the distribution of the colluvial deposits in the upper Mkhomazi River basin, where most of the gullies are located.

Gully formations are known locally as “dongas”, which have been studied from a geological and geomorphological point of view [25,33,42,43,44]. Moreover, gully classifications based on the specific morphological, morphogenetic and evolutionary characteristics have been applied [45,46,47]. Based on the morphological classification of Rowntree [45], derived from our own field observations, we classify the gullies of the upper Mkhomazi River basin into two distinct morphological types: (1) type A gullies, represented by V-shaped gullies where the erosion processes are still active, and (2) type B gullies, represented by more mature landforms and fluvial channels characterized by a U-shape that has reached the erosive base level (i.e., the bedrock). The two gully types may occur independently but can also be connected due to the retrogressive erosion of type B gullies or downslope evolution of type A gullies.

In the study area, both gully types incise the Quaternary colluvial deposits [25,42,48,49] that are widespread along the foot slopes of the Drakensberg mountain range [25]. A colluvial deposit is described as poorly sorted sediment characterized by a wide range of different grain sizes that are transported by unconfined overland flow or in shallow ephemeral channels. These sediments are often eroded by concentrated surface runoff and subsurface soil piping processes on the upper and mid slopes and re-deposited on lower hillslopes, where either transport capacity decreases or overland flow is dissipated [48].

The focus of many previous regional investigations has been the Thukela River basin in central KwaZulu-Natal, whereas in the Mkhomazi River catchment study area, a comprehensive description of colluvial deposits has not been undertaken so far, and the geological map [37] of the area represents only a small and incomplete part of the distribution of the colluvial deposits [43].

In an associated research project, Bosino et al. (2020) [43] recently presented a more accurate distribution of the Masotcheni Formation in the upper Mkhomazi River basin. However, this map illustrates an area characterized by thick colluvial deposits of the Masotcheni Formation that contain buried paleosols. Anyway, this map does not cover the entire area especially colluvia cropping out in non-accessible remote areas are not reported.

The Masotcheni Formation represents successive cycles of colluvial deposition, soil formation and truncation by gully erosion during the late Quaternary [25,50,51].

The goal of this project was to distinguish the two gully morphologies described above though a morphological and morphometric characterization process by using GIS and remote sensing analysis.

Over the years, different statistical approaches have been developed worldwide to assess specific erosion processes and, in particular, gully erosion [4,9,11,21,30,32,52,53,54,55,56]. GIS- and machine learning-based models have been applied to select the gully triggering factors, derive susceptibility maps, implement land management decisions and establish future strategies. Through machine learning methods, indeed, it is possible to evaluate the role of different factors and their interactions, offering a great potential, and they have been increasingly used in recent years [16].

In this study, a pre-existing inventory map of gully erosion was employed for the evaluation and validation of a stochastic approach based on mechanical statistics (MaxEnt) [57] to assess the driving factors for the two types of gully erosion and to derive their spatial susceptibility distribution. Finally, we delineated the colluvial deposits starting from the gully susceptibility map in order to refine the distribution of colluvium in different terrain settings across the study area.

2. Materials and Methods

2.1. Study Area Description, Landform Inventory and Classification

The study area is located in the province of KwaZulu-Natal in South Africa, specifically in the upper part of the Mkhomazi River basin (Figure 1). The area covers about 500 km² and is dominated by the Drakensberg Mountain foothills. The elevation ranges between 1051 m a.s.l. and 2269 m a.s.l. (Figure 2a). The major drainage system of the area is the Mkhomazi River with two main tributaries: the Mkhomazana and Lotheni Rivers (Figure 1).

The region is characterized by a subtropical highland climate (Cwb), according to the Köppen climate classification [58]. The average rainfall amounts to 920 mm/year based on the 1970–2005 precipitation time series [59]. The most common land use is stock grazing of the grassland that covers approximately 66% of the study area, followed by thickets or shrubland covering 7% and commercial afforestation (Pinus patula and Eucalypts) covering 6% (Figure 2b) [60].

Figure 2

Characteristics of the study area: (a) elevation range (derived from the TanDEM-X) and (b) land use classes [60].

[Figure omitted. See PDF]

From a geological point of view, the KwaZulu-Natal hinterland region is characterized by the sedimentary rocks of the Karoo Supergroup succession. In the Drakensberg foothills region, the most extensive unit is the Beaufort Group, comprising Permo–Triassic mudstone and sandstones [61,62,63] of the lower Adelaide Subgroup and upper Tarkastad Subgroup [64]. The Adelaide Subgroup consists mainly of siltstones, mudstones and subordinate very fine- to medium-grained sandstone. The Tarkastad Subgroup comprises more abundant fine- to medium-grained sandstone units and brown-colored mudstones than the Adelaide Subgroup [61,62,63,64]. The Karoo sedimentary succession is intruded by Jurassic dolerite dykes and sills [44,65]. Widespread Masotcheni Formation deposits on the middle-to-lower hillslopes [25,49] of the study area comprise successions of colluvial sediments and buried paleosols [44,66]. The hillslope deposits are characterized by sediment derived from weathered claystone, siltstone and sandstone or reworked Quaternary cover sediments [44,49,50], and in many areas, these susceptible sediments are eroded by gullies [44].

In KwaZulu-Natal, water-driven soil erosion processes such as sheet erosion, extensive rill formation and gullying related to surface and subsurface runoff are the main soil erosion features formed, especially in the highly erodible colluvial deposits along the slopes [42]. Gully incision is caused by the removal of soil by concentrated water flow through channels (Figure 3), and the morphological difference between large gullies and ephemeral low-order stream channels is debatable [67]. However, gully erosion is very effective at deeply incising linear channels that feed sediment into the major drainage systems [67]. In this study, we used the gully and sheet erosion inventory map of Bosino et al. [43] to derive relevant input data utilized for the susceptibility assessment. The gully inventory was subdivided into two groups for training and testing the stochastic model. Gully erosion was classified based on field observation in gully type A and B. The first are represented by active V-shaped dendritic gullies incised into the colluvial slope deposits that had not, in most cases, reached the bedrock (Figure 3a,b). Conversely, type B gullies are characterized by U-shaped entrenched channels that reach the bedrock and valley’s bottom alluvial deposits (Figure 3c,d). The two gully types were classified using aerial photo interpretation based on Google Earth^TM satellite images (year 2016) and orthophotos (year 2009) provided by the Council for Geoscience of South Africa. The classification was validated during a field survey in 2017, 2018 and 2019. The resulting training database contained 122 gully erosion features: 85 gullies of type A and 37 locations of type B gullies.

2.2. Topographic Indices and Environmental Parameters

From the end of the nineteenth century and later for most of the twentieth century, soil erosion forms and features were only studied through extensive field investigations and using aerial photographs [68]. More recently, the development of GIS systems allowed the characterization and modeling of the land surface based on remote sensing image analysis techniques and digital terrain analysis, conducted using digital elevation models [69].

In this study, morphometric analysis was conducted on a high-resolution DEM with a 12-m cell size resolution, from the TanDEM-X satellite platform (provided by Deutsches Zentrum für Luft und Raumfahrt (DLR)) to derive the topographic indices [13] with SAGA GIS 7.9 [70]. The DEM was preprocessed in order to identify and fill the surface depressions using the SAGA GIS tool “Fill Sinks” [71]. In total, as shown in Table 1, 18 factors that affect gully erosion were derived from an extensive literature review [7,9,30,55,56,72,73] consisting of six basic morphometric parameters (slope, aspect, plan curvature, profile curvature and catchment area, as well as catchment slope). In addition, we applied the Topographic Position Index (TPI), which compares the elevation of each cell of the DEM to the mean elevation of a specified neighborhood around that cell [74], as well as the Vector Ruggedness Measure (VRM), which measures the roughness of the terrain surface [75]. Two other parameters based on the slope and specific catchment area were calculated: the Stream Power Index (SPI), describing linear soil erosion potential [76], and the Slope Length Factor (LS-factor), where the L-factor defines the accumulation of water and the S-factor that represent the slope steepness [77]. Here we use the 3D version of the LS-Factor, the Transport Capacity Index (TCI), substituting the slope length with the specific catchment area. Moreover, the parameters for solar radiation, such as direct insolation and diffuse insolation, were calculated [78], and the valley depth and Vertical Distance to Channel Network (VDCN) were also derived [70]. Two further environmental parameters are represented by the lithology, with 8 lithotypes derived from the 1:250,000 geological map [37], and the 2014 land cover classification data derived from the BGIS.SANBI.ORG website (http://bgis.sanbi.org/Projects/Detail/44, accessed on 9 October 2021). The Normalized Vegetation Index (NDVI) was calculated using the QGIS SCP plugin-in for a Sentinel-2 image from 23 June 2019. The NDVI yielded information on the distribution of the vegetation in the area. Finally, the erosion forms were transformed into a grid with a 12-m cell size, and the respective centroids were converted into a point dataset. In total, 17,065 points were identified. For each point, the values of the parameters described above were assigned.

2.3. Modeling

In this study, we used a maximum entropy approach (“MaxEnt”) to assess the two gully erosion types and to assess the extent of the Quaternary deposits. MaxEnt is a machine learning method based on statistical mechanics with a simple and precise mathematical formulation [57] that predicts species distributions using detailed climatic and environmental datasets [82]. We decided to utilize this model because it was applied worldwide for different types of erosion landforms, particularly rill–interrill, gully and badland erosion [7,68,72,73,83,84], and generally produced high-quality data.

MaxEnt expresses the susceptibility for each grid cell as a function of the environmental variables of that grid cell. In this way, it allows predicting the areas affected by gully erosion and assessing the susceptibility of the area.

The MaxEnt model generally works better than other statistical algorithms in terms of model performance (e.g., [13,37]). In addition, the presence-only dataset makes the model more robust to spatial error [68]. The MaxEnt model could estimate a variable’s probability distribution (π) on the X set of sites in the study area. Therefore, π assigned a non-negative value to each site x so that the values π(x) added up to one [85]. The probability distribution function of a target variable in cell x is given as

(1)$P \left(y=1 x\right)=\frac{P\left(x\right|y=1) P \left(y=1\right)}{P\left(x\right)}=\pi \left(x\right)P \left(y=1\right)\left|X\right|$

where |X| is the number of pixels or positions, P(y = 1) is the overall prevalence of species in the study area and π(x) is estimated by the MaxEnt algorithm and is equal to a Gibbs probability distribution derived from the set of features f₁, f₂, …, f_n. Gibbs distributions are exponential distributions parameterized by a vector feature weights λ = (λ₁, …, λ_n) [85] and described as

(2)$q_{\lambda }=\frac{\exp ({{\displaystyle \sum }}_{j=1}^n{\lambda }_{j }{f}_j\left(x\right))}{Z_{\lambda }}$

where $Z_{\lambda }$ is a normalization constant. The MaxEnt $q_{\lambda }$ model at a specific site x depends on the environmental variables at x. The environmental variables on which the model was formed are continuous predictive variables derived from the DEM analysis and the NDVI, as well as categorical variables (limited number of discrete values) such as lithology and land use.

High probability function values in a particular grid cell indicate that the grid cell is expected to have proper conditions for a specific type of erosion. The calculated model represents a probability distribution across all grid cells of the study area. For the MaxEnt modeling, only a part of the original dataset was utilized. In this case, we used the gullies which cropped out in the Mkhomazi River basin to train the model (N_train = 122), while gullies that had been mapped in the Mkhomazana River basin were used to validate the model (N_test = 36).

The model output provided information on the model performance as well as on the variables that were most important to explain the spatial distribution of gullies. Moreover, MaxEnt analysis yielded a susceptibility map for the two gully erosion types. Probabilities ranged between 0, meaning no susceptibility, and 1, representing a very high susceptibility for gully occurrence.

Model Validation

The MaxEnt output was represented by two susceptibility maps for gully type A and type B. In addition, we analyzed the performance of MaxEnt modeling and the relative importance of the 15 environmental parameters for the 2 types of gullies. The MaxEnt model performance was expressed though the Receiver Operating Characteristic (ROC) curve [86] or the Area Under the Curve (AUC) [87]. The ROC curve plots positive instances (sensitivity), which represent no omission error, and negative instances (specificity) [88]. The AUC = 1 value plot in the upper left corner represents the 100% sensitivity vs. 100% specificity representative of the high predictive performance of the model. Following the guidelines in [89], the prediction skill of the model was acceptable, excellent or outstanding if the AUC values exceeded 0.70, 0.80 or 0.90, respectively.

The original gully inventory dataset of Bosino et al. (2020) [43] was split by using the Mkhomazi catchment basin gully data and the Lotheni River catchment data that were used as training data. The Mkhomazana River catchment gullies were used for testing and validating the model.

3. Results

The initial results from this study are represented by the inventory map of the 122 gully erosion forms and features located in the study area: 85 features of gully type A and 37 features of gully type B. Type A gullies covered an area of about 0.7 km², while type B gullies covered an area of about 1 km² for a total of 1.7 km².

As shown in Figure 4, the following environmental layers played a significant role in the formation of the gully erosion: catchment area, NDVI, Vertical Distance to Channel Network, lithology, valley depth, TPI, land use and SPI. To understand the specific contributions of the variables, the individual variable response curves were derived, as illustrated in Figure 5 and Figure 6 for gully types A and B, respectively.

For each of the two gully types, the variable importance of the single independent variables was evaluated (Figure 4). The most important variables that describe type A gullies are represented by the (1) catchment area (33.1%), (2) Vertical Distance to Channel Network (12.3%), (3) lithology (11.7%), (4) valley depth (11.1%) and (5) NDVI (10.8%). Instead, the formation of type B gullies was mainly dominated by the (1) NDVI (18.8%), (2) catchment area (15.2%), (3) Topographic Position Index (11.6%), (4) valley depth (11.1%) and (5) Vertical Distance to Channel Network (10.7%).

Type A gullies develop in an elevation range between 1140 m and 1170 m, as well as above 1460 m, mainly where the valley depth is between 65 m and 75 m. The response curve for the catchment area, the most influencing parameter, showed low values. Additionally, they were characterized by values between 60 m and 90 m for the Vertical Distance to Channel Network. According to Figure 5, type A gullies were associated with siltstone, mudstone and sandstone lithologies belonging to the Adelaide Subgroup. Another parameter influencing type A gullies was the NDVI, characterized by low values between 0 and 0.1, which represented almost no plant cover. The Topographic Position Index, represented by positive values, also flagged the terrain context of gully erosion, characterized by significant elevation differences. The type B gullies (Figure 6) are mainly related to elevation ranges between 1340 m and 1360 m, with a valley depth higher than 20 m. Conversely, NDVI was represented by values higher than 0.3, representing denser plant cover. The catchment area was characterized by a value higher than 0.1, associated with large drainage areas. The Topographic Position Index had negative values that represented portions of the territory characterized by low altitudes (valleys). Gully type B also showed lower values for the Vertical Distance to Channel Network. Finally, the main lithologies influencing the formation of type B gullies were the colluvial deposits of the Masotcheni Formation and the sandstone and mudstone of the Tarkastad Subgroup.

The analysis of the response curves showed that gully type A and gully type B had opposite trends for the VDCN, showing high values for gully type A and low values for gully type B. The valley depth response curves showed low values for type A gullies and high values for type B gullies. The TPI values confirmed that gullies of type A generally were present at higher elevations, whereas gullies of type B were situated mainly at the base of the topographical profile and were represented by low NDVI values.

Additionally, the MaxEnt model was applied to (1) delineate the susceptibility of the landscape to the two gully types and (2) better identify uneroded areas underlain by colluvial deposits. The resulting susceptibility maps as well as the ROC curves are illustrated in Figure 7a,b. Details of the susceptibility values and the respective gully types are reported in Figure 8. The susceptibility values were classified into four susceptibility classes. Pixels with a susceptibility less than 0.60 were designated “non-susceptible”, the “low susceptibility” range was from 0.60 to 0.70, the “medium susceptibility” range was 0.70–0.80 and the “high susceptibility” range was from 0.80 to 1.

Following the work in [90], we summarized the number of pixels for each susceptibility class for the two erosional forms (Table 2). The model performances were assessed using the ROC curves. Following the work in [89], the type A gullies showed an AUC of 0.75, while the type B gullies had an AUC value of 0.70.

The two susceptibility maps (Figure 7a,b) show that 24% of the area was susceptible to type A gully development and 12% was susceptible to the development of type B gullies (Table 2).

4. Discussion

In the study area, type A gullies occupied a smaller portion than the gullies of type B, which were more extensive. Gullies of type A were located in higher elevation areas and along the hillslopes, while gullies of type B developed mainly along the flatter valley bottom areas where the incision exposed stratified colluvium or alluvium on low-level terraces. We demonstrated that the MaxEnt approach can differentiate between the two gully types with acceptable precision and the independent variables that drive the formation of gully erosion. Indeed, the ROC curves obtained for both types of gully erosion showed AUC values above 0.7, which indicate an acceptable level of model performance [89]. To reveal information about the driving factors that characterized the spatial distribution of the two gully types, we derived and evaluated the variable importance.

The outputs showed that the catchment area, NDVI, Vertical Distance to Channel Network, lithology, valley depth, TPI, land use and SPI had a significant impact on gully formation. Even if a standard methodology in the selection of the parameter was not established [91] because each study area is characterized by morphological peculiarities and data availability, which could not be exclusively included in a standard selection process, the above-mentioned factors have been selected and evaluated by several authors [7,32,35,54,55,56]. However, other works [34] include factors like aggregate stability and bulk density, which were not taken into consideration here.

From this study, we can assess that the MaxEnt model can distinguish between type A gullies and type B gullies, which often reach the bedrock and exhibit lateral erosion processes. This trend was confirmed by the individual response curves, in particular those for the catchment areas, VDCN and NDVI (see Figure 4 and Figure 5). The other environmental layers seemed to play a minor role in the model’s application, as shown in Figure 4 (e.g., slope), although their specific response curves may be useful for the interpretation of the spatial distribution of the gully erosion types and location. Indeed, as stated by several authors [23,33,35], the gullies developed on gentle slopes with gradients less than 10 degrees.

The type A gullies were dendritic erosional forms and showed a tendency for expansive growth, with low vegetation cover confirmed by NDVI values close to zero, while the type B gullies reached the bedrock, had larger catchment areas, and accumulate more runoff water. Therefore, persistent water flows at the gully bottom can be commonly observed, especially after larger rainfall events. These gullies are often characterized by a denser vegetation cover. The NDVI was detected as one of the most important variables in [92] as well. The valley depth showed lower values for type A gullies which had not yet completed the erosive cycle. Conversely, the type B gullies were characterized by higher valley depth values, since they were normally older and therefore more incised into the valley infill deposits. The importance of land cover and valley depth was also observed in [56].

In addition, the type A gullies showed positive TPI values, indicating more elevated areas, while the type B gullies were represented by the negative TPI values associated with the valley bottoms.

The catchment area also played a key role in the development of the two gully types. The low values of the response curves indicate that the type A gullies were characterized by small catchment areas or were situated close to interfluves. The relevance of the catchment area was also highlighted by the authors of [33], who linked the largest continuous gullies with high catchment area values. Additionally, the type A gullies were represented by high VDCN values, whereas the type B gullies were characterized by larger drainage areas and lower VDCN values, which expressed the vertical distance of the gullies in respect to the valley bottom. Thus, the type A gullies had higher erosional potential due to the generally steeper slope inclination and hence higher runoff velocities, but they displayed low-to-medium flow accumulation values. In contrast, the type B gullies were associated with larger catchment areas and higher runoff volumes but lower slopes and hence lower runoff velocities.

With respect to the lithology, type A gullies occurred predominantly in sheetwash colluvium mantles overlying the mudstones and sandstones of the Adelaide Subgroup, while type B gullies developed principally in the Tarkastad Subgroup, characterized by terrain influenced by the alternating sandstone and mudrock beds, which is evidence of the development of gullies in soils derived from these formations that was already found [33,93]). However, our detailed mapping showed that the published 1:250,000 geological map does not delineate the full extent of Masotcheni Formation colluvia in this area. The observed development of gullies in the colluvial deposits was widespread in South Africa [25,33,35,43,49,50] and in other study sites worldwide (e.g., [7,73,94,95]).

Regarding gully forms and features, the susceptibility maps show that 7% of the territory was affected by medium susceptibilities and 7% had high values of susceptibility for type A gullies, while for type B gullies, 4% had medium susceptibility and 3% had high values of susceptibility. These values are not very high, but much of the territory is mountainous, while the susceptible areas are mainly concentrated on the lower gradient, with colluvium-mantled slopes characterized by arable land and pastures where the erosion has a direct impact on the livelihood of the local population. Gullies are often directly connected to higher order drainage systems with negative implications on water quality due to the sediment load washed into the reaches of the major channel (Figure 9a,b).

The trained model was validated with gullies of the Mkhomazana River basin, using the inventory proposed by Bosino et al. (2020) [43]. The results showed good correspondence with the higher values of susceptibility and the presence of gullies (Figure 10). The results indicate that the model provided appropriate predictions of susceptible gully areas prone to future gully erosion development. Proper land management activities should focus on these areas in order to reduce soil degradation and loss of soil.

Many of the gullies of the study area are associated with colluvial deposits [42,49], which are more erodible than the other lithotypes present in the area. This is also highlighted by the gully susceptibility maps. Even in small, non-typical slope areas, the maps delineate previously unmapped colluvial deposits. Consequently, the high susceptibility zones of the study areas indicate the existence of colluvial deposits that may not be exposed by gully erosion. Indeed, the two susceptibility maps illustrate the vast spatial distribution of potentially sensitive areas for gully erosion related to colluvial slope and valley fillings. In particular, the highest susceptibilities were registered along the valleys and in the depositional areas of the slope systems, while low susceptibilities represent the very high elevation and well-vegetated areas. The susceptibility maps show that 14% of the terrain had medium-to-high values for the type A gullies and 7% had such values for the type B gullies, which presumably indicate the areas where colluvial deposits are present. By joining the two susceptibility maps, we created a more accurate map illustrating a more detailed distribution of colluvial deposits (Figure 11).

5. Conclusions

During the past century, widespread erosion has developed in the colluvial deposits of the Drakensberg foot hills, leading to a loss of agricultural productivity, loss of biodiversity, release of soil organic carbon and damage to infrastructures [66]. This study highlights that a large portion of the upper Mkhomazi catchment study area is susceptible to gully erosion, and consequently, protective measures and new management strategies need to be adopted to prevent and reduce potential future soil loss. We applied the MaxEnt model approach, which is based on previously mapped gullies, in the study area [43]. The inventory was split into two types of gullies that generally show morphological and maturity differences. The gullies of type A generally incise the colluvial slope deposits, and the type B gullies are mainly concentrated in the valley’s bottom areas, characterized by colluvial or alluvial terrace deposits with interbedded paleosols. The latter often expose the more resistant bedrock substrates and are characterized by lateral gully side wall erosion. As independent variables, we derived several topographic indices from a TanDEM-X elevation model, as well as information about lithology, land use and vegetation.

The accuracy of the MaxEnt model results was evaluated using the ROC curves. The model generated acceptable performance values (AUC values > 0.7), showing that a differentiation between the two gully types was feasible and generally showed good accuracy.

The most important variables for gully erosion are the catchment area, NDVI, Vertical Distance to Channel Network, lithology, valley depth and TPI.

The slope facets most susceptible to type A gullies are along the steep and gentle upper hillslopes, where vegetation is sparse. The most susceptible areas for type B gullies are along the river valleys, where there is a higher surface runoff. This study indicates that 14% of the area is characterized by type A gullies and 7% by type B gullies, and it shows medium-to-high susceptibility to gully erosion. Even if these values do not represent a very large erosion susceptible area, much of the impact occurs in areas that are suitable for livestock grazing or agricultural purposes.

Furthermore, we indirectly identify areas covered by colluvial deposits by using the susceptibility maps, confirming the strong relationship between colluvial deposits and concentrated surface and subsurface runoff landforms that was observed in the field [25,42,44,50]. The susceptibility maps indicate that colluvial deposits are more widespread along the main river valleys and the lower slope valley bottom positions than it is indicated in the published geological maps.

This study demonstrates that an integrative approach based on stochastic modeling, remote sensing and GIS analysis yields valuable results, allowing a precise prediction of susceptible areas for gully erosion and detecting the presence of the susceptible colluvial deposits. However, the preliminary map of colluvial slope deposits and valley fillings should be further validated and confirmed with additional fieldwork.

The continuous development of erosion landforms may lead to a severe desertification of the area and reducing or even destroying its agricultural potential. Moreover, the connection between the gullies and the drainage system has potentially severe environmental consequences related to water quality. Predicting the spatial distribution of gully erosion through multi-factorial susceptibility analysis may provide a theoretical background for erosion risk reduction and may provide a basis for an enhanced land and drainage network management.

Author Contributions

Conceptualization, Alice Bernini, Alberto Bosino, Greg A. Botha and Michael Maerker; methodology, Alice Bernini and Alberto Bosino; software, Alice Bernini and Alberto Bosino; validation, Alice Bernini and Alberto Bosino; formal analysis, Alice Bernini and Alberto Bosino; investigation, Alice Bernini, Alberto Bosino, Greg A. Botha and Michael Maerker; resources, Alice Bernini, Alberto Bosino, Greg A. Botha and Michael Maerker; data curation, Alice Bernini and Alberto Bosino; writing—original draft preparation, Alice Bernini and Alberto Bosino; writing—review and editing, Greg A. Botha, Alberto Bosino and Michael Maerker; visualization, Alice Bernini and Alberto Bosino; supervision, Michael Maerker; project administration, Michael Maerker; funding acquisition, Michael Maerker. All authors have read and agreed to the published version of the manuscript.

Funding

The ROCEEH project financed traveling support by6 Heidelberg Academy of Sciences and Humanities. This research was also conducted with financial support of the University of Pavia.

Data Availability Statement

Data available on request from the authors.

Acknowledgments

The authors are grateful to the Council for Geoscience office in Pietermaritzburg for logistical support. We acknowledge the DLR and the TDX Science Team for providing the Tandem-X dataset of the study area.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Study area in the upper Mkhomazi River catchment.

Figure 3. Gully erosion in the upper Mkhomazi River basin: (a,b) examples of type A gullies on hillslopes; (c,d) examples of type B gullies eroded into stratified colluvial sediments and buried paleosols at the base of the hillslope.

Figure 4. Relative importance of the environmental variables.

Figure 5. Variable importance for gully type A.

Figure 6. Variable importance for gully type B.

Figure 7. MaxEnt model susceptibilities for the two gully erosion types: (a) susceptibility maps, ROC curves and AUC values for gully type A; (b) susceptibility maps, ROC curves and AUC values for gully type B.

Figure 8. (a,c) Details of the susceptibility map. (b,d) gully type A and gully type B, respectively.

Figure 9. (a) Sediment transport in the Mkhomazi River during a flood’s peak. (b) Evidence of higher connectivity between gully systems and a secondary drainage network.

Figure 10. Gullies of the Mkhomazana River used for external model validation.

Figure 11. Integration of the two susceptibility maps that potentially represent the spatial distribution of colluvial deposits.

Table 1

Environmental predictors for maximum entropy modeling.

Type	Variable Range	Reference
Topographic indices
Elevation	1051–2269 m a.s.l.	[71]
Slope	0−64°	[79]
Aspect	0–360°	[79]
Profile curvature	−0.04–0.05	[79]
Plan curvature	−23.36–6.98	[79]
TPI	−33.39–36.54	[74]
VRM	0–0.43	[75]
Catchment slope	0–52.4°	[80]
Catchment area	144–558,355,264 m2	[80]
SPI	0–2,591,456	[76]
LS-factor (TCI)	−6.16−14.95	[77]
Direct insolation	0–6.60 kWh/m−2	[78]
Diffuse insolation	0.57–0.91 kWh/m−2	[78]
VDCN	−2–975 m a.s.l.	[70]
Valley depth	−66–275 m a.s.l.	[70]
Environmental data
Lithology	8 classes	[37]
Land use	35 classes	[60]
Remote sensing data
NDVI	−0.94–0.99	[81]

Table 2

Number of pixels for every susceptibility class for the two types of gully erosion.

Susceptibility Classes	Npixel Gully A	Npixel Gully B	Gully A (%)	Gully B (%)
No susceptibility	3,033,849	3,502,309	76	88
Low	403,579	204,983	10	5
Medium	272,187	149,545	7	4
High	286,605	139,383	7	3

References

1. Mararakanye, N.; Sumner, P.D. Gully erosion: A comparison of contributing factors in two catchments in South Africa. Geomorphology; 2017; 288, pp. 99-110. [DOI: https://dx.doi.org/10.1016/j.geomorph.2017.03.029]

2. Lal, R. Soil degradation by erosion. L. Degrad. Dev.; 2001; 12, pp. 519-539. [DOI: https://dx.doi.org/10.1002/ldr.472]

3. Panagos, P.; Standardi, G.; Borrelli, P.; Lugato, E.; Montanarella, L.; Bosello, F. Cost of agricultural productivity loss due to soil erosion in the European Union: From direct cost evaluation approaches to the use of macroeconomic models. L. Degrad. Dev.; 2018; 29, pp. 471-484. [DOI: https://dx.doi.org/10.1002/ldr.2879]

4. Azareh, A.; Rahmati, O.; Rafiei-Sardooi, E.; Sankey, J.B.; Lee, S.; Shahabi, H.; Ahmad, B. Bin Modelling gully-erosion susceptibility in a semi-arid region, Iran: Investigation of applicability of certainty factor and maximum entropy models. Sci. Total Environ.; 2019; 655, pp. 684-696. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.11.235] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30476849]

5. Zweig, R.; Filin, S.; Avni, Y.; Sagy, A.; Mushkin, A. Land degradation and gully development in arid environments deduced by mezzo- and micro-scale 3-D quantification—The Negev Highlands as a case study. J. Arid Environ.; 2018; 153, pp. 52-65. [DOI: https://dx.doi.org/10.1016/j.jaridenv.2017.12.006]

6. Kariminejad, N.; Hosseinalizadeh, M.; Pourghasemi, H.R.; Bernatek-Jakiel, A.; Alinejad, M. GIS-based susceptibility assessment of the occurrence of gully headcuts and pipe collapses in a semi-arid environment: Golestan Province, NE Iran. L. Degrad. Dev.; 2019; 30, pp. 2211-2225. [DOI: https://dx.doi.org/10.1002/ldr.3397]

7. Zakerinejad, R.; Maerker, M. Prediction of gully erosion susceptibilities using detailed terrain analysis and Maximum Entropy Modeling: A case study in the Mazayejan plain, Southwest Iran. Geogr. Fis. E Din. Quat.; 2014; 37, pp. 67-76. [DOI: https://dx.doi.org/10.4461/GFDQ.2014.37.7]

8. Vanmaercke, M.; Poesen, J.; Maetens, W.; de Vente, J.; Verstraeten, G. Sediment yield as a desertification risk indicator. Sci. Total Environ.; 2011; 409, pp. 1715-1725. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2011.01.034] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21316738]

9. Conforti, M.; Aucelli, P.P.C.; Robustelli, G.; Scarciglia, F. Geomorphology and GIS analysis for mapping gully erosion susceptibility in the Turbolo stream catchment (Northern Calabria, Italy). Nat. Hazards; 2011; 56, pp. 881-898. [DOI: https://dx.doi.org/10.1007/s11069-010-9598-2]

10. Conforti, M.; Pascale, S.; Pepe, M.; Sdao, F.; Sole, A. Denudation processes and landforms map of the Camastra River catchment (Basilicata—South Italy). J. Maps; 2013; 9, pp. 444-455. [DOI: https://dx.doi.org/10.1080/17445647.2013.804797]

11. Conoscenti, C.; Angileri, S.; Cappadonia, C.; Rotigliano, E.; Agnesi, V.; Märker, M. Geomorphology Gully erosion susceptibility assessment by means of GIS-based logistic regression: A case of Sicily (Italy). Geomorphology; 2014; 204, pp. 399-411. [DOI: https://dx.doi.org/10.1016/j.geomorph.2013.08.021]

12. Conoscenti, C.; Agnesi, V.; Angileri, S.; Cappadonia, C.; Rotigliano, E.; Märker, M. A GIS-based approach for gully erosion susceptibility modelling: A test in Sicily, Italy. Environ. Earth Sci.; 2013; 70, pp. 1179-1195. [DOI: https://dx.doi.org/10.1007/s12665-012-2205-y]

13. Märker, M.; Pelacani, S.; Schröder, B. A functional entity approach to predict soil erosion processes in a small Plio-Pleistocene Mediterranean catchment in Northern Chianti, Italy. Geomorphology; 2011; 125, pp. 530-540. [DOI: https://dx.doi.org/10.1016/j.geomorph.2010.10.022]

14. Oostwoud Wijdenes, D.J.; Poesen, J.; Vandekerckhove, L.; Ghesquiere, M. Spatial distribution of gully head activity and sediment supply along an ephemeral channel in a Mediterranean environment. Catena; 2000; 39, pp. 147-167. [DOI: https://dx.doi.org/10.1016/S0341-8162(99)00092-2]

15. Martins, B.; Meira-Castro, A.; Nunes, A.; Lourenço, L. The development of gullies in a Mediterranean environment: The example of the Corgo gully (central Portugal). Energy Rep.; 2020; 6, pp. 794-799. [DOI: https://dx.doi.org/10.1016/j.egyr.2019.11.004]

16. Vanmaercke, M.; Panagos, P.; Vanwalleghem, T.; Hayas, A.; Foerster, S.; Borrelli, P.; Rossi, M.; Torri, D.; Casali, J.; Borselli, L. et al. Measuring, modelling and managing gully erosion at large scales: A state of the art. Earth-Sci. Rev.; 2021; 218, [DOI: https://dx.doi.org/10.1016/j.earscirev.2021.103637]

17. Sidorchuk, A.; Märker, M.; Moretti, S.; Rodolfi, G. Gully erosion modelling and landscape response in the Mbuluzi River catchment of Swaziland. Catena; 2003; 50, pp. 507-525. [DOI: https://dx.doi.org/10.1016/S0341-8162(02)00123-6]

18. Sidorchuk, A.Y.; Matveeva, T.A. Periglacial gully erosion on the east european plain and its recent analog at the Yamal Peninsula. Geogr. Environ. Sustain.; 2020; 13, pp. 183-194. [DOI: https://dx.doi.org/10.24057/2071-9388-2019-01]

19. Gong, C.; Lei, S.; Bian, Z.; Liu, Y.; Zhang, Z.; Cheng, W. Analysis of the development of an erosion gully in an open-pit coal mine dump during a winter freeze-thaw cycle by using low-cost UAVs. Remote Sens.; 2019; 11, 1356. [DOI: https://dx.doi.org/10.3390/rs11111356]

20. Maerker, M.; Schillaci, C.; Melis, R.T.; Kropáček, J.; Bosino, A.; Vilímek, V.; Hochschild, V.; Sommer, C.; Altamura, F.; Mussi, M. Geomorphological processes, forms and features in the surroundings of the Melka Kunture Palaeolithic site, Ethiopia. J. Maps; 2019; 15, [DOI: https://dx.doi.org/10.1080/17445647.2019.1669497]

21. Torkashvand, A.M. Investigation of some methodologies for gully erosion mapping. J. Appl. Sci.; 2008; 8, pp. 2435-2441. [DOI: https://dx.doi.org/10.3923/jas.2008.2435.2441]

22. Pourghasemi, H.R.; Sadhasivam, N.; Kariminejad, N.; Collins, A.L. Gully erosion spatial modelling: Role of machine learning algorithms in selection of the best controlling factors and modelling process. Geosci. Front.; 2020; 11, pp. 2207-2219. [DOI: https://dx.doi.org/10.1016/j.gsf.2020.03.005]

23. Flügel, W.A.; Märker, M.; Moretti, S.; Rodolfi, G.; Sidrochuk, A. Integrating geographical information systems, remote sensing, ground truthing and modelling approaches for regional erosion classification of semi-arid catchments in South Africa. Hydrol. Process.; 2003; 17, pp. 929-942. [DOI: https://dx.doi.org/10.1002/hyp.1171]

24. Märker, M.; Moretti, S.; Rodolfi, G. Assessment of water erosion processes and dynamics in semi-arid regions of Southern Africa (Kwazulu/Natal, RSA, and Swaziland) using the Erosion Response Units concept (ERU). Geogr. Fis. E Din. Quat.; 2001; 24, pp. 71-83.

25. Botha, G.A. The geology and paleopedology of late quaternary colluvial sedimensts in northern Kwazulu/Natal. Mem. Geol. Surv. South Afr.; 1996; 83, 165.

26. Poesen, J.W.; Vandaele, K.; Van Wesemael, B. Contribution of gully erosion to sediment production on cultivated lands and rangelands. IAHS-AISH Publ.; 1996; 236, pp. 251-266.

27. Ionita, I.; Fullen, M.A.; Zglobicki, W.; Poesen, J. Gully erosion as a natural and human-induced hazard. Nat. Hazard; 2015; 79, pp. 1-5. [DOI: https://dx.doi.org/10.1007/s11069-015-1935-z]

28. Shruthi, R.B.V.; Kerle, N.; Jetten, V. Object-based gully feature extraction using high spatial resolution imagery. Geomorphology; 2011; 134, pp. 260-268. [DOI: https://dx.doi.org/10.1016/j.geomorph.2011.07.003]

29. Sugathan, N.; Biju, V.; Renuka, G. Influence of soil moisture content on surface albedo and soil thermal parameters at a tropical station. J. Earth Syst. Sci.; 2014; 123, pp. 1115-1128. [DOI: https://dx.doi.org/10.1007/s12040-014-0452-x]

30. Saha, S.; Sarkar, R.; Thapa, G.; Roy, J. Modeling gully erosion susceptibility in Phuentsholing, Bhutan using deep learning and basic machine learning algorithms. Environ. Earth Sci.; 2021; pp. 1-22. [DOI: https://dx.doi.org/10.1007/s12665-021-09599-2]

31. Debanshi, S.; Pal, S. Assessing Gully Erosion Susceptibility in Mayurakshi River Basin of Eastern India; Springer: Dordrecht, The Netherlands, 2020; Volume 22, ISBN 0123456789

32. Garosi, Y.; Sheklabadi, M.; Pourghasemi, H.R.; Besalatpour, A.A.; Conoscenti, C.; Van Oost, K. Comparison of differences in resolution and sources of controlling factors for gully erosion susceptibility mapping. Geoderma; 2018; 330, pp. 65-78. [DOI: https://dx.doi.org/10.1016/j.geoderma.2018.05.027]

33. Le Roux, J.J.; Sumner, P.D. Factors controlling gully development: Comparing continuous and discontinuous gullies. L. Degrad. Dev.; 2012; 23, pp. 440-449. [DOI: https://dx.doi.org/10.1002/ldr.1083]

34. Kariminejad, N.; Hosseinalizadeh, M.; Pourghasemi, H.R.; Bernatek-Jakiel, A.; Campetella, G.; Ownegh, M. Evaluation of factors affecting gully headcut location using summary statistics and the maximum entropy model: Golestan Province, NE Iran. Sci. Total Environ.; 2019; 677, pp. 281-298. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.04.306]

35. Kakembo, V.; Xanga, W.W.; Rowntree, K. Topographic thresholds in gully development on the hillslopes of communal areas in Ngqushwa Local Municipality, Eastern Cape, South Africa. Geomorphology; 2009; 110, pp. 188-194. [DOI: https://dx.doi.org/10.1016/j.geomorph.2009.04.006]

36. Valentin, C.; Poesen, J.; Li, Y. Gully erosion: Impacts, factors and control. Catena; 2005; 63, pp. 132-153. [DOI: https://dx.doi.org/10.1016/j.catena.2005.06.001]

37. Geological Survey 1:250.000 Geological Series Map. 2928 Drakensberg 1981. Available online: https://maps.geoscience.org.za/portal/home/item.html?id=a01d123ed8f14390941b48a1b5164bcd (accessed on 8 October 2021).

38. du Plessis, C.; van Zijl, G.; Van Tol, J.; Manyevere, A. Machine learning digital soil mapping to inform gully erosion mitigation measures in the Eastern Cape, South Africa. Geoderma; 2020; 368, 114287. [DOI: https://dx.doi.org/10.1016/j.geoderma.2020.114287]

39. Chaplot, V. Impact of terrain attributes, parent material and soil types on gully erosion. Geomorphology; 2013; 186, pp. 1-11. [DOI: https://dx.doi.org/10.1016/j.geomorph.2012.10.031]

40. Li, Z.; Fang, H. Impacts of climate change on water erosion: A review. Earth-Sci. Rev.; 2016; 163, pp. 94-117. [DOI: https://dx.doi.org/10.1016/j.earscirev.2016.10.004]

41. Yang, D.; Kanae, S.; Oki, T.; Koike, T.; Musiake, K. Global potential soil erosion with reference to land use and climate changes. Hydrol. Process.; 2003; 17, pp. 2913-2928. [DOI: https://dx.doi.org/10.1002/hyp.1441]

42. Rienks, S.M.; Botha, G.A.; Hughes, J.C. Some physical and chemical properties of sediments exposed in a gully (donga) in northern KwaZulu-Natal, South Africa and their relationship to the erodibility of the colluvial layers. Catena; 2000; 39, pp. 11-31. [DOI: https://dx.doi.org/10.1016/S0341-8162(99)00082-X]

43. Bosino, A.; Bernini, A.; Botha, G.A.; Bonacina, G.; Pellegrini, L.; Omran, A.; Hochschild, V.; Sommer, C.; Maerker, M. Geomorphology of the upper Mkhomazi River basin, KwaZulu-Natal, South Africa, with emphasis on late Pleistocene colluvial deposits. J. Maps; 2020; [DOI: https://dx.doi.org/10.1080/17445647.2020.1790435]

44. Lyons, R.; Tooth, S.; Duller, G.A. Chronology and controls of donga (gully) formation in the upper Blood River catchment, KwaZulu-Natal, South Africa: Evidence for a climatic driver of erosion. Holocene; 2013; 23, pp. 1875-1887. [DOI: https://dx.doi.org/10.1177/0959683613508157]

45. Rowntree, K.M. Morphological characteristics of gully networks and their relationship to host materials, Baringo District, Kenya. GeoJournal; 1991; 23, pp. 19-27. [DOI: https://dx.doi.org/10.1007/BF00204405]

46. Soms, J. Development and morphology of gullies in the river Daugava Valley, South—Eastern Latvia. Landf. Anal.; 2011; 17, pp. 183-192.

47. Deng, Q.; Qin, F.; Zhang, B.; Wang, H.; Luo, M.; Shu, C.; Liu, H.; Liu, G. Characterizing the morphology of gully cross-sections based on PCA: A case of Yuanmou Dry-Hot Valley. Geomorphology; 2015; 228, pp. 703-713. [DOI: https://dx.doi.org/10.1016/j.geomorph.2014.10.032]

48. Watson, A.; Price Williams, D.; Goudie, A.S. The paleoenvironmental interpretation of colluvial sediments and paleosols of the late Pleistocene hypotermal in South Africa. Paleogeography Paleoclimatology Paleoecol.; 1984; 45, pp. 225-249. [DOI: https://dx.doi.org/10.1016/0031-0182(84)90008-7]

49. Botha, G.A.; Fedoroff, N. Palaeosols in Late Quaternary colluvium, northern KwaZulu-Natal, South Africa. J. Afr. Earth Sci.; 1995; 21, pp. 291-311. [DOI: https://dx.doi.org/10.1016/0899-5362(95)00072-2]

50. Temme, A.J.A.M.; Baartman, J.E.M.; Botha, G.A.; Veldkamp, A.; Jongmans, A.G.; Wallinga, J. Climate controls on late Pleistocene landscape evolution of the Okhombe valley, KwaZulu-Natal, South Africa. Geomorphology; 2008; 99, pp. 280-295. [DOI: https://dx.doi.org/10.1016/j.geomorph.2007.11.006]

51. Botha, G.A. Cenozoic stratigraphy of South Africa: Current challenges and future possibilities. South Afr. J. Geol.; 2021; 124, pp. 1-26. [DOI: https://dx.doi.org/10.25131/sajg.124.0054]

52. Pourghasemi, H.R.; Yousefi, S.; Kornejady, A.; Cerdà, A. Performance assessment of individual and ensemble data-mining techniques for gully erosion modeling. Sci. Total Environ.; 2017; 609, pp. 764-775. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.07.198]

53. Rahmati, O.; Haghizadeh, A.; Pourghasemi, H.R.; Noormohamadi, F. Gully erosion susceptibility mapping: The role of GIS-based bivariate statistical models and their comparison. Nat. Hazards; 2016; 82, pp. 1231-1258. [DOI: https://dx.doi.org/10.1007/s11069-016-2239-7]

54. Rahmati, O.; Tahmasebipour, N.; Haghizadeh, A.; Pourghasemi, H.R.; Feizizadeh, B. Evaluation of different machine learning models for predicting and mapping the susceptibility of gully erosion. Geomorphology; 2017; 298, pp. 118-137. [DOI: https://dx.doi.org/10.1016/j.geomorph.2017.09.006]

55. Arabameri, A.; Chen, W.; Loche, M.; Zhao, X.; Li, Y.; Lombardo, L.; Cerda, A.; Pradhan, B.; Bui, D.T. Comparison of machine learning models for gully erosion susceptibility mapping. Geosci. Front.; 2020; 11, pp. 1609-1620. [DOI: https://dx.doi.org/10.1016/j.gsf.2019.11.009]

56. Arabameri, A.; Nalivan, O.A.; Saha, S.; Roy, J.; Pradhan, B.; Tiefenbacher, J.P.; Ngo, P.T.T. Novel ensemble approaches of machine learning techniques in modeling the gully erosion susceptibility. Remote Sens.; 2020; 12, [DOI: https://dx.doi.org/10.3390/rs12111890]

57. Phillips, S.J.; Robert, P.; Andersonb, C.; Schapired, R.E. Maximum entropy modeling of species geographic distributions. Ecol. Modell.; 2006; 6, pp. 231-259. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2005.03.026]

58. Kottek, M.; Grieser, J.; Beck, C.; Rudolf, B.; Rubel, F. World map of the Köppen-Geiger climate classification updated. Meteorol. Z.; 2006; 15, pp. 259-263. [DOI: https://dx.doi.org/10.1127/0941-2948/2006/0130]

59. Nel, W.; Sumner, P. Rainfall and Temperature Attributes on the Lesotho-Drakensberg Escarpment Edge, Southern Africa. Geogr. Ann. Ser. A Phys. Geogr.; 2016; 90, pp. 97-108. [DOI: https://dx.doi.org/10.1111/j.1468-0459.2008.00337.x]

60. BGIS—SANBI. Available online: https://bgis.sanbi.org/ (accessed on 27 August 2021).

61. Rutherford, A.B. The Sedimentology and Stratigraphy of the Beaufort Group of the Karoo Supergroup in the Vicinity of Thaba Nchu, Central Free State Province. 2009; Available online: https://wiredspace.wits.ac.za/handle/10539/7235 (accessed on 8 October 2021).

62. Turner, B. Revised stratigraphy of the Beaufort Group in the southern Karoo Basin. Geology; 1981; 24, pp. 87-98.

63. Rutherford, A.B.; Rubidge, B.S.; Hancox, P.J. Sedimentology and Palaeontology of the Beaufort Group in the Free State Province Supports A Reciprocal Foreland Basin Model for the Karoo Supergroup, South Africa. South Afr. J. Geol.; 2015; 118, pp. 355-372. [DOI: https://dx.doi.org/10.2113/gssajg.118.4.355]

64. Catuneanu, O.; Elango, H.N. Tectonic control on fluvial styles: The Balfour formation of the Karoo Basin, South Africa. Sediment. Geol.; 2001; 140, pp. 291-313. [DOI: https://dx.doi.org/10.1016/S0037-0738(00)00190-1]

65. Grab, S.; Knight, J. Landscapes and Landforms of the Africa; 1st ed. Grab, S.; Knight, J. Springer International Publishing: New York, NY, USA, 2015; ISBN 9783319035604

66. Botha, G.A.; Temme, A.J.A.M.; Singh, R.G. Colluvial deposits and slope stability. Quaternary Enviromental Change in Southern Africa: Physical ad Human Dimensions; Botha, G.A.; Temme, A.J.A.M.; Singh, R.G. Cambridge University Press: Cambridge, UK, 2016; pp. 137-152.

67. Poesen, J.; Nachtergaele, J.; Verstraeten, G.; Valentin, C. Gully erosion and environmental change: Importance and research needs. Catena; 2003; 50, pp. 91-133. [DOI: https://dx.doi.org/10.1016/S0341-8162(02)00143-1]

68. Bosino, A.; Giordani, P.; Quénéhervé, G.; Maerker, M. Assessment of calanchi and rill–interrill erosion susceptibilities using terrain analysis and geostochastics: A case study in the Oltrepo Pavese, Northern Apennines, Italy. Earth Surf. Process. Landf.; 2020; 45, pp. 3025-3041. [DOI: https://dx.doi.org/10.1002/esp.4949]

69. Schillaci, C.; Braun, A.; Kropacek, J. Terrain analysis and landform recognition. Br. Soc. Geomorphol.; 2015; 2, 18. [DOI: https://dx.doi.org/10.13140/RG.2.1.3895.2802]

70. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Böhner, J. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev.; 2015; 8, pp. 1991-2007. [DOI: https://dx.doi.org/10.5194/gmd-8-1991-2015]

71. Wang, L.; Liu, H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modelling. Int. J. Geogr. Inf. Sci.; 2006; 20, pp. 193-213. [DOI: https://dx.doi.org/10.1080/13658810500433453]

72. Javidan, N.; Kavian, A.; Pourghasemi, H.R.; Conoscenti, C.; Jafarian, Z.; Rodrigo-Comino, J. Evaluation of multi-hazard map produced using MaxEnt machine learning technique. Sci. Rep.; 2021; 11, pp. 1-20. [DOI: https://dx.doi.org/10.1038/s41598-021-85862-7]

73. Zakerinejad, R.; Omran, A.; Hochschild, V.; Maerker, M. Assessment of gully erosion in relation to lithology in the southwestern Zagros Mountains, Iran using ASTER data, GIS and stochastic modeling. Geogr. Fis. E Din. Quat.; 2018; 41, pp. 95-104. [DOI: https://dx.doi.org/10.4461/GFDQ.2018.41.15]

74. Weiss, A. Topographic position and landforms analysis. Proceedings of the Poster presentation, ESRI User Conference; San Diego, CA, USA, 9–13 July 2001; Volume 64, pp. 227-245.

75. Sappington, J.M.; Longshore, K.M.; Thompson, D.B. Qunatifying Landscape Ruggedness for Animal Habitat Analysis: A Case Study Using Bighorn Sheep in the Mojave Desert. J. Wildl. Manag.; 2007; 71, pp. 1419-1426. [DOI: https://dx.doi.org/10.2193/2005-723]

76. Milevski, I.; Dragicevic, S.; Kostadinov, S. Digital elevation model and satellite images an assessment of soil erosion potential in the Pcinja catchment. Glas. Srp. Geogr. Drus. Serb. Geogr. Soc.; 2007; 87, pp. 11-20. [DOI: https://dx.doi.org/10.2298/GSGD0702011M]

77. Panagos, P.; Borrelli, P.; Meusburger, K. A new European slope length and steepness factor (LS-factor) for modeling soil erosion by water. Geosciences; 2015; 5, pp. 117-126. [DOI: https://dx.doi.org/10.3390/geosciences5020117]

78. Hofierka, J.; Súri, M. The solar radiation model for Open source GIS: Implementation and applications. Proceedings of the Open Source Free Software GIS—GRASS Users Conference; Trento, Italy, 11–13 September 2002; pp. 11-13.

79. Florinsky, I.V. Computation of the third-order partial derivatives from a digital elevation model. Int. J. Geogr. Inf. Sci.; 2009; 23, pp. 213-231. [DOI: https://dx.doi.org/10.1080/13658810802527499]

80. Sørensen, R.; Zinko, U.; Seibert, J. On the calculation of the topographic wetness index: Evaluation of different methods based on field observations. Hydrol. Earth Syst. Sci.; 2006; 10, pp. 101-112. [DOI: https://dx.doi.org/10.5194/hess-10-101-2006]

81. The Copernicus Open Access Hub. Available online: https://scihub.copernicus.eu/ (accessed on 8 October 2021).

82. Medley, K.A. Niche shifts during the global invasion of the Asian tiger mosquito, Aedes albopictus Skuse (Culicidae), revealed by reciprocal distribution models. Glob. Ecol. Biogeogr.; 2010; 19, pp. 122-133. [DOI: https://dx.doi.org/10.1111/j.1466-8238.2009.00497.x]

83. Cama, M.; Schillaci, C.; Kroacek, J.; Hochschild, V.; Märker, M. A Probabilistic Assessment of Soil Erosion Susceptibility in a Head Catchment of the Jemma Basin, Ethiopian Highlands. Geosciences; 2020; 10, 248. [DOI: https://dx.doi.org/10.3390/geosciences10070248]

84. Maerker, M.; Bosino, A.; Scopesi, C.; Giordani, P.; Firpo, M. Geoderma Assessment of calanchi and rill-interrill erosion susceptibility in northern Liguria, Italy: A case study using a probabilistic modelling framework. Geoderma; 2020; 371, 114367. [DOI: https://dx.doi.org/10.1016/j.geoderma.2020.114367]

85. Phillips, S.J.; Dudík, M. Modeling of species distributions with Maxent: New extensions and a comprehensive evaluation. Ecography; 2008; 31, pp. 161-175. [DOI: https://dx.doi.org/10.1111/j.0906-7590.2008.5203.x]

86. Hanley, J.A.; McNeil, B.J. A method of comparing the areas under receiver operating characteristic curves derived from the same cases. Radiology; 1983; 148, pp. 839-843. [DOI: https://dx.doi.org/10.1148/radiology.148.3.6878708] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/6878708]

87. Daskivich, T.; Luu, M.; Noah, B.; Fuller, G.; Anger, J.; Spiegel, B. Differences in online consumer ratings of health care providers across medical, surgical, and allied health specialties: Observational study of 212,933 providers. J. Med. Internet Res.; 2018; 20, pp. 29-36. [DOI: https://dx.doi.org/10.2196/jmir.9160]

88. Phillips, S.B.; Aneja, V.P.; Kang, D.; Arya, S.P. Modelling and analysis of the atmospheric nitrogen deposition in North Carolina. Int. J. Glob. Environ. Issues; 2006; 6, pp. 231-252. [DOI: https://dx.doi.org/10.1504/IJGENVI.2006.010156]

89. Hosmer, D.W.; Lemeshow, S.; Sturdivant, R.X. Applied Logistic Regression. Wiley Ser. Probab. Stat.; 2013; 479. [DOI: https://dx.doi.org/10.2307/2532419]

90. Bianchini, S.; Del Soldato, M.; Solari, L.; Nolesini, T.; Pratesi, F.; Moretti, S. Badland susceptibility assessment in Volterra municipality (Tuscany, Italy) by means of GIS and statistical analysis. Environ. Earth Sci.; 2016; 75, pp. 1-14. [DOI: https://dx.doi.org/10.1007/s12665-016-5586-5]

91. Conoscenti, C.; Maggio, C.; Rotigliano, E. Soil erosion susceptibility assessment and validation using a geostatistical multivariate approach: A test in Southern Sicily. Nat. Hazards; 2008; 46, pp. 287-305. [DOI: https://dx.doi.org/10.1007/s11069-007-9188-0]

92. Mokarram, M.; Zarei, A.R. Determining prone areas to gully erosion and the impact of land use change on it by using multiple-criteria decision-making algorithm in arid and semi-arid regions. Geoderma; 2021; 403, 115379. [DOI: https://dx.doi.org/10.1016/j.geoderma.2021.115379]

93. Weaver, A.V.B. The distribution of soil erosion as a function of slope aspect and parent material in Ciskei, Southern Africa. GeoJournal; 1991; 23, pp. 29-34. [DOI: https://dx.doi.org/10.1007/BF00204406]

94. Jacobson, R.B. Genesis and Distribution of Colluvium in Buffalo Creek Area, Marion County, West Virginia. Transp. Res. Rec.; 1986; 1089, pp. 62-67.

95. Lin, J.; Zhu, G.; Wei, J.; Jiang, F.; Wang, M.K.; Huang, Y. Mulching effects on erosion from steep slopes and sediment particle size distributions of gully colluvial deposits. Catena; 2018; 160, pp. 57-67. [DOI: https://dx.doi.org/10.1016/j.catena.2017.09.003]