Four watersheds with different spatial scales in the granitic hill region of southern China were selected as the study areas (Figure 2), including southern China (20°13′-34°38′N and 104°28′-120°40′E), Jiangxi province (24°29-30°05′N and 113°34′-118°29′E), Ganxian county (25°26′-26°17′N and 114°42′-115°22′E) in Jiangxi province, and Tiancun township (26°04′-26°14′N and 115°02′-115°16′E) in Ganxian county, Jiangxi province. The four sites were the focus-study areas of gully erosion (Wei et al., 2021b) due to the high density, various sizes, and different development stages of gullies (2005 National Survey by the Yangtze River Water Resources Commission of the Ministry of Water Resources of China), suggesting the four sites can provide an ideal place for analyzing the spatial scale effects of gully erosion.

There are 239125, 48058, 4138, and 590 gullies in southern China, Jiangxi province, Ganxian county, and Tiancun township, respectively. These areas have a subtropical monsoon climate, favoring the accumulation of weathering materials (Xia et al., 2019). The granite weathering crust thickness ranges from 5 to 70 m, providing a material basis for gully development (Liao et al., 2019). The primary lithology is granite (syenite), especially the syenite and monzogranite, and the granite in these areas was formed mainly in the Early Caledonian and the Yanshan periods, with weathered granitic acidic red soil as the primary soil type. Field investigation found that most gullies are located in low mountains or hills with an old geological stage (Wei et al., 2022). The natural vegetation is evergreen broadleaf forest, and gully initiation was triggered by large-scale forest destruction (Woo et al., 1997). Pinus massoniana, Eriachne pallescens, Phragmites australis, and Dicranopteris linearis are the main tree species in the gully areas. Detailed characteristics of these areas are given in Table 1.

The methodology in the four different spatial scales of watershed consists of five steps (Figure 3): (1) preparing the gully inventory map; (2) extracting the geo-environmental factors; (3) analyzing the multicollinearity between gully susceptibility and potential conditions; (4) assessing the driving factors of gully erosion at four different scales; (5) modeling and mapping gully erosion susceptibility.

The receiver operating characteristic (ROC) curve and accuracy were employed to validate the performance of the models. The ROC is a graph established by sensitivity and 1-specificity with different cut-off thresholds on other susceptibility classes, where a high area under the ROC curve (AUC) can reflect a high goodness of fit and model consistency (Zabihi et al., 2018). In addition, the accuracy, ranging from 0 to 1, can show the model’s overall performance (Bui et al., 2016). These procedures were run in the R 3.5.1 environment using the “pROC” and “caret” packages. The accuracy can be described by equation 8:

where TP, TN, FP, and FN represent true positive, true negative, false positive, and false negative, respectively.

In Figure 4, the relative importance of the driving factors was seen to vary significantly at the four investigated spatial scales. According to the classification of the potential factors, the average contribution values were ranked in the order of lithology and soil (32.65%) > topography (22.40%) > human activities (22.31%) > climate (11.32%) > vegetation (11.31%) across the four spatial scales. Among these factors, lithology and soil properties dominated gully erosion at the scales of southern China, province, county, and township, and their contribution to gully erosion at these four scales was 28.49%, 29.27%, 31.06%, and 41.78%, respectively, indicating an emphasized impact of lithology and soil properties on gully erosion as spatial scale decreased. Similarly, vegetation was more critical at the scale of county (14.57%) and township (13.85%) than at the scale of southern China (8.32%) and province (8.54%). In contrast, the contributions of climatic factors showed a downtrend and were most important in the scale of southern China (20.77%). Moreover, with decreasing spatial scale, the contribution of topography to gully erosion followed a trend of increasing first and then decreasing, with its peak at the county scale (24.07%). The total contribution of human activities showed an N-shaped variation, with the highest value in the province scale (25.17%) and the smallest value in the scale of southern China (19.93%).

The relative importance of a single factor suggested that precipitation (7.82%), land use and land cover (LULC) (6.16%), rainfall erosivity (10.15%), and elevation (DEM) (11.59%) were the most critical factors for triggering gully erosion at the individual scale of southern China, province, county, and township, respectively (Figure 4). In addition, elevation (7.33%), population density (6.69%), rainfall erosivity (5.46%), gross domestic product (GDP) (5.43%), hypsometric integral (HI) (5.37%), distance from road (5.31%), and soil moisture (5.09%) were the top seven leading factors at all the four scales based on their average contribution values. Notably, with decreasing basin scale, the relative importance increased for DEM and GDP, but decreased for precipitation and rainfall erosivity. Meanwhile, the relative importance exhibited a trend of increasing first and then decreasing for LULC, HI, and distance from road, a V-shaped variation pattern for population density, and an inverse N-shaped trend for rainfall erosivity and soil moisture. Collectively, with the decrease of basin scale, the impact weakens for climate factors, but enhances for soil properties and vegetation.

The area under receiver operating characteristic curves (AUC-ROC) was used to evaluate the performance of gully erosion susceptibility modeling, and gully susceptibility was further divided into four classes: very low, low, moderate, and high (Figures 5a-5l). In the three models, the RF model was shown to have the highest AUC value in predicting gully erosion susceptibility at each scale, fluctuating between 0.978 and 0.70, followed by the LMT model with its AUC averaging 0.925, 0.898, 0.740, and 0.687, respectively. The MLR model had the lowest mean AUC values (0.827-0.60), especially in the low and medium classes. In Figure 6, the total accuracy of models was also shown to decrease with decreasing spatial scale and data input. Moreover, the RF model was seen to have consistently higher accuracy values than the other two models, keeping pace with the AUC-ROC. Overall, in the three models, the RF model exhibited the highest accuracy in predicting gully erosion susceptibility at different spatial scales.

Gully erosion susceptibility was mapped using the RF model (Figures 7a-7d). The very low, low, moderate, and high groups accounted for 61.68%, 27.97%, 7.10%, and 3.25% in southern China, respectively. Notably, the south side (e.g., the center of Guangdong province, north of Jiangxi province, north and east of Guangxi, south of Fujian province, and center and southeast Hunan province) were more susceptible to gully erosion than the northern areas (e.g., Anhui province, Hubei province, and north of Hunan province). In Jiangxi province, gully erosion was concentrated in the north, (especially in the counties of Ganxian, Xingguo, Yudu, and Ningdu), with a susceptibility value of 33.08%, 21.84%, 19.08%, and 26.00% for the very low, low, moderate, and high groups, respectively. At the two smaller spatial scales, the very low, low, moderate, and high groups accounted for 21.89%, 15.71%, 29.56%, and 32.84% in Ganxian county, while 4.38%, 22.42%, 27.66%, and 45.54% in Tiancun township.

Scale, a hot topic in soil erosion, ecology, hydrology and other disciplines, is critical for measuring geographic phenomena (Vanmaercke et al., 2021). A scale effect occurs with the change of a given temporal or spatial scale, which is primarily limited by the complexity of erosion process, the heterogeneity of environments, and the size of monitored research areas (Peeters et al., 2008; Borrelli et al., 2020). The gully distribution with a spatial clustering pattern is non-uniform, suggesting that gully erosion has an apparent scale effect and is susceptible to environmental conditions (Wang et al., 2020). However, the standard guidelines for feature selection are controversial, especially regarding scale transformation (Bui et al., 2016). Additionally, the driving factors of gully erosion vary from case to case (Conforti et al., 2011). Therefore, determining the driving factors of gully erosion and examining their spatial scaling effect is extremely important (Conoscenti et al., 2013). In this study, we found lithology and soil properties as the characteristics most responsible for gully erosion (Figure 4). Previous studies have shown that lithology plays a crucial role in gully distribution, development size, and evolution process (Bernatek-Jakiel and Poesen, 2018). Especially the granite weathering crust thickness was considered as a material basis for gully initiation, which has been confirmed by similar erosion landforms worldwide, such as “Lavaka” in Madagascar (Cox et al., 2010), “Vocorocas” in Brazil (Bezerra et al., 2020), and “Calanchi” in Italy (Caraballo-Arias et al., 2014). Besides, soil physiochemical properties can determine gully occurrence and development by mainly affecting erodibility and hydrological function (Xia et al., 2019). Several field studies demonstrated a significant difference in soil properties along the granite soil profile due to varying weathering degree (Deng et al., 2020; Wei et al., 2021b). Once the soil layer has been eroded, the lower layers have low shear strength and are susceptible to collapsing (Tao et al., 2017).

Human activities can directly trigger gully formation (Borrelli et al., 2020). Our results showed that human activities are also a crucial determinant for gully development at all scales (Figure 4). For example, the LULC ranked first at the province scale, agreeing with the reports of Chen et al. (2019) and Zhou et al. (2021) that unreasonable land use type could influence soil quality, runoff yield, and soil loss rate. Previous studies also indicated that large-scale forest destruction could decrease vegetation cover, cause rapid topsoil loss, expose weathering crust, and trigger gully formation (Woo et al., 1997). Additionally, long-term agricultural activity on slope cropland could affect soil degradation and cause gully occurrence in southern China (Wang et al., 2017). Moreover, the distance from road could significantly affect gully formation, aligning with the conclusion of Arabameri et al. (2018). From the perspective of physical process, road building in hills and mountains can directly generate a steep slope and significantly increase the susceptibility to slope failure, especially during heavy rains (Nyssen et al., 2002). Furthermore, population density ranked top three for gully occurrence at one or more scales, consistent with the conclusions of other researchers (Yuan et al., 2019; Wei et al., 2022). However, Tarolli and Sofia (2016) mentioned that direct evidence for population impact on the gully initiation could hardly be obtained, indicating an uncertain relationship between gully erosion and human activities. Therefore, great attention should be paid to identifying the quantitative impact of human activities on gully erosion.

We also found that the relative importance of other driving factors varies at different spatial scales. For instance, climatic factors showed a downtrend with decreasing spatial scale. It is fascinating to point out that the contribution of precipitation ranked first in southern China and then decreased with declining spatial scales, which may be related to the weak spatial heterogeneity (Zhang et al., 2022). In contrast, with decreasing spatial scales, the importance of elevation rose to the top ranking at the township scale. Lin et al. (2015) confirmed that elevation as a critical factor in the watershed scale, providing sufficient energy for gully initiation and expansion. In addition, elevation change is more drastic at a relatively small scale than at a larger scale, leading to a stronger scouring velocity (Wei et al., 2022). However, the resolution of digital elevation images could directly affect the accuracy of capturing spatial variability in geomorphic processes (Millares et al., 2012). Analogously, vegetation factors, especially NDVI, presented an upward trend with decreasing basin scale. Erktan et al. (2016) reported that vegetation growth and distribution in a smaller spatial scale are more sensitive to climate, topography, and human activities (e.g., agricultural cultivation, deforestation, and urbanization), which can significantly affect the protective effect of vegetation on soil and water. Furthermore, the control effect of soil properties increases with decreasing spatial scale, further supporting the widespread agreement that soil properties are the intrinsic and most important indicators at the plot scale (Castillo and Gomez, 2016). For the increasing importance of soil properties on gully erosion with reducing spatial scale, one possible explanation is that vertical differentiation of soil properties is more likely to cause collapse instability at a relatively small scale than at a larger scale (Liao et al., 2022). Although their correlation with gully erosion could be affected by field sampling accuracy of soil properties (Wu et al., 2021), our results could provide a reference for selecting potential factors to gully erosion at different spatial scales.

Numerous susceptibility prediction techniques have been employed in gully susceptibility mapping (Pourghasemi et al., 2017). However, accurate comparison of susceptibility models at different spatial scales remains challenging (Zabihi et al., 2016), and cross-comparison of various categories of approaches is widely accepted as a proxy for accurate risk assessment (Gayen et al., 2019). Additionally, the receiver operating characteristic (ROC) and other indices (such as total accuracy, precision, sensitivity, or specificity) have been shown as an effective way to select the optimal model and validate gully erosion susceptibility maps (Bui et al., 2016). In the present study, the performance of three data-driven algorithms (MLR, LMT, and RF) was evaluated and compared, and their performance was shown to gradually decrease with the spatial scale, probably due to decreasing data input (de Vente and Poesen, 2005). Moreover, AUC and total accuracy evaluation results revealed that the RF model consistently outperformed the other two models at different spatial (Figures 5 and 6), with its optimal results in line with most of the literature (Pourghasemi and Rahmati, 2018; Lei et al., 2020; Jaafari et al., 2022), probably due to its advantage of no data loss or outlier intervention (Naghibi et al., 2015). Furthermore, the RF model can handle mass databases (continuous and categorical variables) and is more sensitive to dataset noise than other boosting approaches (Ließ et al., 2012).

A few studies have reported that some other methods (such as boosted regression tree, artificial neural network, and multi-criteria decision-making) were more suitable for gully susceptibility modeling (Elith et al., 2008; Garosi et al., 2018; Zhang et al., 2022), probably due to different conditioning factor inputs and algorithm structure (Arabameri et al., 2020). The present study also found that the LMT model ranked second, followed by MLR, confirming a higher accuracy of machine learning models than traditional statistical methods (Azareh et al., 2019). However, machine learning is hampered by the increasing complexity of models and the difficulty in exploring driving factors (Micheletti et al., 2013). The traditional statistical method, MLR, was simple and effective for more than two types of data (Wei et al., 2021a). Therefore, research needs, data categories, and user’s operation difficulty should also be considered in model selection.

Potential determinants selection at different scales should carefully consider which factors are crucial to the results of importance ranking and accuracy of gully susceptibility mapping (Vanmaercke et al., 2021). Based on our findings, we put forward several recommendations. Specifically, soil and lithology, human activities, and topography should be adopted at all scales. Microtopography factors (such as plane curvature, profile curvature, and topographic wetness index) extracted from high-resolution digital elevation models are essential at the plot scale. Climatic factors (such as precipitation, temperature, and rainfall erosivity) should be prioritized at the regional scale. Comparatively, vegetation factors had better be considered at the watershed scale, and besides vegetation coverage, vegetation species and morphology should also be considered.

In this study, the dominant drivers of gully erosion were found to differ across different spatial scales, so the gully control strategies should have distinct orientations and priorities at different scales. For instance, lithology and soil properties ranked first at all scales, indicating that soil reclamation might be essential for gully control. Additionally, given the profound effect of precipitation at the scale of southern China, it is suggested to establish an extreme rainfall monitoring system to forecast gully formation. Moreover, land use has a significant impact at the province scale, so appropriate land use management should be emphasized. Furthermore, elevation has a greater effect at the scales of county and township, suggesting that microtopography alteration (such as slope cutting, drainage ditches, check dams, and terracing) should be addressed at these scales.