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

: 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 ﬁeld surveys to validate the gully inventory. The gullies were subdivided into slope gullies and ﬂuvial 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.


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.

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 2 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].  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 loworder 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

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

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.

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 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 1 , f 2 , . . . , f n . Gibbs distributions are exponential distributions parameterized by a vector feature weights λ = (λ 1 , . . . , λ n ) [85] and described as where Z λ is a normalization constant. The MaxEnt q λ 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.

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 2 , while type B gullies covered an area of about 1 km 2 for a total of 1.7 km 2 .
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 Figures 5 and 6 for gully types A and B, respectively.
-Inf. 2021, 10, x FOR PEER REVIEW 9 of 22 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 Figures 5 and 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) 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     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).

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

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 Figures 4 and 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). 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).   [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 susceptibil-ity 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). 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).    Figure 11. Integration of the two susceptibility maps that potentially represent the spatial distribution of colluvial deposits.

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 high- Figure 11. Integration of the two susceptibility maps that potentially represent the spatial distribution of colluvial deposits.

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.