The study area is located in Changshun County of Qiannan Buyi Miao Autonomous Prefecture (106°13′6″E−106°38′48″E, 25°38′48″N–26°17′30″N), situated within a representative karst rocky desertification zone in China (Figure 1). The region experiences a subtropical monsoon climate, with an average annual precipitation of 1,260 mm, most of which falls between May and August. Elevation ranges from 700 m to 1748 m, resulting in substantial topographic variation. Within Guizhou Province, mountainous and hilly terrain accounts for approximately 93% of the total land area. Carbonate rocks are widely distributed throughout the study area, with dolomite and limestone being the predominant lithologies. The dominant soil type is calcareous soil. Rocky desertification is severe in the region, primarily occurring in areas underlain by dolomite and limestone, which together account for more than 70% of the total land area. Grassland resources are distributed in a patchy pattern, with low connectivity and limited spatial continuity, leading to highly fragmented grassland landscapes. The total grassland area in Guizhou Province is 189,400 ha, of which 112,000 ha are scattered grasslands, accounting for approximately 59 percent of the total grassland area.

The grasslands in the study area are generally distributed on gentle slopes, with more than 70% located on terrain with a slope of less than 15°. Dolomite grassland (DG) and limestone grassland (LG) are the dominant bedrock types in the region. Based on lithological characteristics, grasslands developed on DG and LG substrates were selected for this study. Field sampling was conducted from July to August 2023. All selected plots shared similar natural conditions. Within each lithological type, plots were established with varying levels of bedrock exposure rates (RER) to represent different degrees of rocky desertification. For both DG and LG, which are primarily affected by slight to moderate desertification, plots were established with RERs of 20%, 40%, and 55%, based on field investigations. Control plots with no exposed bedrock were also included. Each plot measured 10 m × 10 m, resulting in a total of eight plots. Sampling sites were selected under comparable environmental conditions, with major factors controlled and climatic and temporal variability minimized.

Plant roots in grassland ecosystems are primarily distributed within the 0–20 cm soil layer. Accordingly, soil samples were collected vertically from two depth intervals, 0–10 cm and 10–20 cm. Before sampling, surface litter and the dense root mat were carefully removed, and mineral soil was collected at the specified depths. In both DG and LG, the influence of bedrock exposure on surrounding soil was considered. A total of ten plots with exposed bedrock (50 cm × 25 cm) were established, and sampling was conducted along the rock–soil interface adjacent to the exposed rock. In contrast, two control plots without bedrock exposure (1 m × 1 m) were selected, and sampling was carried out in areas distant from rock surfaces.

Soil physical properties were determined using samples collected with a 100 cm³ ring knife. For erosion resistance testing, undisturbed soil samples were collected using a 500 cm³ ring knife. Prior to sampling, surface vegetation and litter were carefully removed. During collection, care was taken to ensure full contact between the ring knife and the soil, with no air entering the sample.

1.

Soil Aggregate

The water-stable aggregates in the soil were measured using the Le Bissonnais method (LB method), which includes three treatments: fast wetting (FW), slow wetting (SW), and wetting and shaking (WS).

Mean Weight Diameter (MWD): M W D = ∑ i = 1 n x i W i where xi is the mean diameter (mm) of soil aggregates in the ith size class (mm), and wi is the proportion of aggregates in that size class relative to the total dry mass of the soil sample (%).

Relative Dispersion Index (RSI): R S I = M W D S W ‐ M W D F W M W D S W × 100 %

Where MWD _SW is the mean weight diameter after SW (mm), and MWD _FW is the mean weight diameter after FW (mm).

Relative Mechanical Breakdown Index (RMI): R M I = M W D S W ‐ M W D W S M W D S W × 100 %

Where MWD _SW is the mean weight diameter after SW (mm), and MWD _WS is the mean weight diameter after WS (mm). 2.

Undisturbed soil scouring experiment

Undisturbed soil samples collected using 500 cm² ring cutters were directly installed in the flume for scouring tests. The undisturbed soil cores were carefully transferred and placed directly into the flume to preserve their natural structure during the scouring tests. Surface runoff erosion was simulated using a uniform top-down water supply system. Prior to the experiment, all soil samples were placed in a container with a constant water level and left to stand for 12 h. This preconditioning step ensured uniform initial moisture conditions and minimized the influence of moisture variability on erosion behavior.The scouring duration was set to 15 min, based on the typical duration of the most intense rainfall events recorded in the study area over the past 10 years. The water supply rate was maintained at 4 L/min, which corresponds to the maximum runoff rate observed under moderate-intensity rainfall conditions in the region. The slope of the flume was adjusted according to field monitoring data to closely replicate actual hillslope conditions, with a slope of 15°.

During the experiment, runoff–sediment mixture samples were collected at 1 min intervals during the first 3 minutes. Thereafter, mixed samples were collected every 2 minutes using a 10L container until the end of the test. All collected samples were used to determine instantaneous sediment concentration and total sediment yield, which served as indicators for evaluating soil erosion resistance under different treatment conditions.

Anti-scourability coefficient (AS): A S = S r × T D v

Where AS is the anti-scourability coefficient (L/g), Sr is the water supply rate (L/min), T is the scouring duration (min), and Dv is the dry weight of the eroded sediment (g). 3.

Root Content

After the undisturbed soil scouring experiment, soil samples were washed using a 100-mesh sieve to remove soil particles and debris. The root systems were then gently blotted to remove excess moisture. All root washing, sieving, and measurement procedures were conducted following a standardized protocol to minimize operational variability. Root morphological parameters, including total length, average diameter, surface area, and volume, were measured using the WinRHIZO root analysis system.

Standard statistical methods were used to analyze the experimental data. One-way analysis of variance (ANOVA) was performed to evaluate differences among treatments (p < 0.05). Principal component analysis and Pearson correlation analysis were conducted using the SPSS statistical software (version 19.0). All figures were generated using Origin 2021.

Soil structural stability declined significantly as the RER increased, with distinct responses observed under different wetting treatments (Figure 2). The 0–10 cm layer showed the greatest sensitivity to changes in RER. Aggregate size distribution was relatively uniform at 0% RER. In contrast, the proportion of aggregates larger than 2 mm was substantially reduced when RER reached 55%, while the proportion of smaller aggregates in the 0.25–0.054 mm range increased. These shifts indicate a substantial breakdown of soil structure under high bedrock exposure conditions.

As degradation intensified, the stability of soil aggregates declined with RER increasing from 0% to 55% (Figure 3). In the 0–10 cm soil layer, MWD decreased by approximately 40%–60% under 55% RER compared to 0% RER. Under the fast wetting (FW) treatment, MWD dropped from 1.2 mm to 0.5 mm, representing a reduction of more than 58%, indicating that rapid wetting significantly accelerated aggregate breakdown. Although the SW treatment also led to a reduction in MWD, the decline was less pronounced, suggesting that slower wetting helps preserve the structural resilience of the soil.

Analysis of the RSI and RMI showed that increasing bedrock exposure led to a substantial rise in both indicators (Figure 4). At 55% RER, RSI was 1.5–2 times higher than that at 0% RER, indicating an enhanced tendency for aggregate disintegration upon wetting. RMI also increased significantly under the FW treatment, with values at 55% RER reaching approximately 2.2 times those observed under 0% RER. These results suggest that degradation intensifies the risk of mechanical breakdown in soil aggregates. Compared with DG, LG exhibited smaller increases in both RSI and RMI, suggesting relatively greater structural resilience. Average aggregate length decreased from 1,000 mm to 500 mm as RER increased from 0% to 55% (Figure 5). In addition, aggregate diameter, surface area, and volume all declined significantly with increasing exposure, further indicating the deterioration of soil structure under more severe degradation.

The sediment concentration in overland flow increased markedly with rising RER in both DG and LG (Figure 6). During the first minute, sediment yield was relatively low under 0% RER, reaching approximately 0.9 g/min in DG and 0.6 g/min in LG. Under 55% RER, these values increased to around 1.8 g/min and 1.5 g/min, respectively more than double the initial levels. Peak sediment yield were generally observed within the first 5 min of the scouring process. These results indicate that an increase in exposed bedrock area significantly enhances surface runoff and erosive force, intensifying soil particle detachment and transport. The presence of exposed rock appears to amplify slope runoff convergence and hydrodynamic disturbance, leading to a more rapid and severe erosion response.

The AS variation in grasslands was illustrated over time and across different levels of degradation (Figure 7). After 15 min, the soil resistance coefficient under 0% RER reached approximately 600 L/g for DG and 550 L/g for LG, indicating strong erosion resistance. In contrast, under 55% RER, the values decreased significantly to around 300 L/g and 250 L/g, respectively. Overall, the resistance coefficient remained consistently low under high RER conditions, with the most pronounced differences observed during the initial 1–5 min of erosion. This suggests that severe degradation reduces soil structural stability and shear strength, leading to the development of surface scour and concentrated flow paths. These results demonstrate that increasing RER not only intensifies sediment production on the surface but also significantly reduces the soil’s resistance to scouring. Therefore, RER is a key factor driving the acceleration of hillslope erosion. Overall, lithology consistently modulated the erosion response to increasing RER. Across the gradient, LG generally maintained higher scouring resistance and a less pronounced deterioration in key stability-related indicators than DG, suggesting a buffering effect of limestone-derived soils under comparable degradation levels. These patterns indicate that the magnitude of degradation-driven changes depends not only on RER but also on parent material.

Principal component analysis was conducted based on nine diagnostic indicators related to soil retention capacity. Three principal components were extracted, namely, F₁, F₂, and F₃, each with an eigenvalue greater than 1. The cumulative variance explained by these components reached 71.78%, indicating that they effectively represent the majority of the information contained in the original variables and provide a satisfactory dimensionality reduction outcome (Table 1). According to the component loading matrix, F₁ showed high loadings on RSI, soil resistance coefficient, root surface area, and root volume density, primarily reflecting soil structural stability and root coverage characteristics (Table 2). F₂ was mainly associated with the resistance coefficient and average root diameter, highlighting its role in representing erosion resistance. F₃ was dominated by runoff sediment yield, indicating its relevance to soil erosion output. These three principal component expressions collectively reflect the dominant mechanisms underlying soil retention capacity in grasslands across different dimensions. The expressions for the three principal components are as follows: F 1 = − 0.19 C 1 − 0.31 C 2 − 0.4 C 3 + 0.11 C 4 − 0.16 C 5 + 0.44 C 6 + 0.3 C 7 + 0.46 C 8 + 0.42 C 9 F 2 = 0.24 C 1 − 0.36 C 2 − 0.04 C 3 + 0 . 23 C 4 + 0.65 C 5 − 0.13 C 6 + 0.53 C 7 − 0.1 C 8 − 0.14 C 9 F 3 = − 0.57 C 1 + 0.2 C 2 + 0.14 C 3 + 0.76 C 4 + 0.14 C 5 − 0.06 C 6 − 0.09 C 7 − 0.05 C 8 + 0.07 C 9

As shown by the principal component scores and composite scores for different lithological plots was shown in Table 3, the LG outperformed the DG in both F₁ and F₃, indicating stronger performance in terms of structural stability and anti-scour resistance. Although DG showed a slightly higher score in F₂, its overall composite score was relatively low, suggesting that despite a localized advantage in erosion resistance, its overall soil retention capacity was weaker. These differences in functional performance are primarily attributed to the variation in dominant indicators within each principal component. In F₁ and F₃, the indicators with the highest weights were RSI, root length, and root surface area, all of which are strengths of LG. Therefore, the integrated analysis indicates that LG exhibits superior soil retention capacity compared to DG. The principal components F₁, F₂, and F₃ effectively capture and differentiate the soil retention performance of grasslands under varying levels of degradation.

Soil structural stability plays a critical regulatory role in hillslope erosion by influencing both particle cohesion and shear resistance (Bryan, 2000; Reubens et al., 2007). It directly affects the processes of runoff generation and sediment detachment. Numerous studies have demonstrated that MWD is a reliable indicator of soil resistance to hydraulic disturbance, with lower MWD values typically reflecting increased soil erodibility (Ding and Zhang, 2016; Gan et al., 2024; Le Bissonnais, 1996). In this study, MWD showed a consistent negative correlation with sediment yield across all time intervals from 1 to 15 min. This suggests that the weakening of soil structural scale not only intensifies the initial erosion response but also alters the progression of the scouring process over time. The breakdown of macroaggregates rapidly exposes finer internal particles, increasing the potential for sediment supply (Chang et al., 2006; Six et al., 2000).

Differences in wetting rate further amplified the structural effects. Under FW treatment, soil aggregates were subjected to rapid redistribution of internal and external pressures, making them more prone to disintegration and fragmentation, which resulted in a significant decrease in MWD. When soil is in a structurally fragile state, abrupt hydraulic disturbances can trigger irreversible structural collapse (Haeri et al., 2017; Li et al., 2022). In this study, the RMI increased significantly under the FW treatment, consistent with this theoretical expectation. These results indicate that the mode of water input is a key factor regulating both soil structural stability and erosion response. This process illustrates a clear mechanism: as aggregate structure weakens, the soil becomes more vulnerable to raindrop impact and overland flow. Fragile aggregates break down more easily, exposing fine particles that are readily detached and transported. The resulting increase in runoff and shear force further accelerates erosion. This interaction between structural fragility and hydraulic stress explains the rapid rise in sediment yield under degraded conditions.

The structural equation model further revealed the sustained regulatory role of soil structural stability in the erosion process (Figure 8). MWD had a significant negative direct effect on sediment yield across all time points, with the effect being more pronounced in DG soils. This suggests that in soil systems with inherently weaker structure, aggregate stability plays a more critical role in constraining erosion dynamics. Compared to vegetation cover, physical structural parameters are often more direct determinants of erosion intensity (Bronick and Lal, 2005; Li et al., 2025; Sadeghian et al., 2021). In addition, changes in AS closely mirrored the trends observed in MWD and RMI. As the RER increased, AS declined sharply, indicating a systematic reduction in the slope’s ability to resist hydraulic forces. A decrease in aggregate stability is known to significantly lower the critical shear stress of the slope surface, thereby accelerating the detachment of soil particles (Doerr et al., 2000; Yuan et al., 2022). The findings of this study provide direct empirical evidence for this mechanism within karst grassland systems.

The response of slope systems to degradation disturbances often exhibits pronounced nonlinearity. When key structural support elements are sufficiently weakened, soil erosion resistance may undergo an abrupt decline. Previous studies have shown that once vegetation degradation or structural disruption reaches a critical threshold, erosion processes can shift from a gradual progression to a rapidly accelerating phase (Zheng, 2006; Ravi et al., 2010). In this study, when the RER increased to 55%, MWD decreased by more than 50%, while RSI and RMI rose significantly. Simultaneously, sharp changes were observed in both AS and sediment yield, clearly indicating the presence of a threshold response.The structural equation model further clarified the internal pathway underlying the erosion resistance threshold triggered by degradation (Figure 9). RER did not directly determine sediment yield intensity. Instead, it exerted an indirect effect by weakening MWD, a key structural variable, thereby amplifying the erosive impact of runoff. This chain response, characterized by degradation-induced structural deterioration leading to intensified erosion, is consistent with concepts proposed in various disturbed surface systems. Structural breakdown acts as a central mediator driving the abrupt increase in erosion risk (Gyssels et al., 2005; Ye and Liu, 2020).

The threshold mechanism exhibited distinct differences under varying lithological conditions. In DG soils, structural variables had a stronger path effect on sediment yield, indicating that erosion resistance in these systems relies more heavily on the maintenance of aggregate stability. In contrast, LG soils retained a certain degree of structural resilience under the same level of degradation, with a relatively smaller decline in scouring resistance. This divergence aligns with previous findings that parent material conditions significantly influence erosion sensitivity. This higher structural stability in LG is likely associated with differences in soil structural development under limestone conditions, which enhance aggregate stability and erosion resistance during scouring processes. The environment in which soils form can markedly affect both the manner and the rate at which erosion thresholds are crossed (Lal, 1998; Poesen et al., 2003).

The results of the principal component analysis provided integrated evidence for identifying erosion resistance thresholds (Table1–3). Variables such as MWD, RSI, soil resistance coefficient, and root surface area contributed most strongly to the principal components, indicating that the threshold in erosion resistance is not driven by a single physical factor. Instead, it reflects the combined deterioration of both structural stability and biological support elements. When root systems and aggregate structure are simultaneously compromised, the slope system’s capacity to buffer erosive disturbances is significantly reduced (Reubens et al., 2007; De Baets et al., 2007). This pattern is consistent with the integrated diagnostic results of this study. From a temporal perspective, structural parameters continued to exert significant influence on sediment yield during the later stages of the scouring process. This indicates that once a threshold in erosion resistance is crossed, its effects are persistent rather than transient. Such a dynamic, where slow-changing variables govern rapid system responses, is a key reason why degraded ecosystems often struggle to recover on their own (Pan et al., 2022). Therefore, implementing intervention measures before structural thresholds are exceeded is essential for maintaining the stability of slope systems.

This study was based on plot-scale experiments in a typical karst region and may not fully capture variability across broader spatial or climatic gradients. Future studies should include longer-term and landscape-scale observations. Nonetheless, indicators such as MWD, AS, and RMI show promise for early identification of degradation and can inform erosion control efforts, especially where rock exposure exceeds critical thresholds.