Spatial assessment of soil degradation under intensive agricultural pressure in the central loess plateau

1 Introduction

Soil degradation (SD) has emerged as one of the most pervasive threats to soil resources in farmland and semi-dryland agroecosystems, where fragile soil structures are exposed to climatic fluctuations and intensive land-use practices (Khosravi Aqdam et al., 2023; Ma et al., 2024). As global food production increasingly depends on the sustained functioning of agricultural soils, understanding the mechanisms that drive SD has become a central challenge in environmental science and land management (Hossain et al., 2020). SD processes; such as erosion, nutrient depletion, salinization, structural breakdown, and soil organic matter (SOC) decline; are often interwoven with topographic complexity, hydrological variability, and anthropogenic disturbances (Hossain et al., 2020). These interacting forces generate spatially heterogeneous degradation patterns that cannot be captured adequately through traditional field-based assessments alone (Bajocco et al., 2018). Consequently, modern SD studies emphasize the need for integrated frameworks that link soil quality metrics, environmental covariates, and advanced modeling techniques to characterize degradation dynamics with higher precision (De Caires et al., 2025).

The Loess Plateau of China represents one of the most sensitive landscapes worldwide in terms of soil vulnerability (Ye et al., 2023). Its deep Quaternary loess deposits, inherently weak aggregate stability, and high erodibility create conditions where even modest disturbances can initiate severe degradation cascades. In the central part of the plateau, monsoon-driven precipitation, steep terrain discontinuities, and long-term agricultural intensification collectively shape soil transformation pathways. While numerous studies have examined erosion and soil loss in this region, there is still limited understanding of how physicochemical soil attributes; such as pH, SAR, CEC, organic carbon content, and nutrient dynamics; interact with environmental gradients to influence the progression of SD. The high spatial variability of soil texture, irrigation regimes, and cultivation practices further complicates the assessment, underscoring the need for data-driven approaches capable of revealing underlying patterns that traditional methods may overlook (Qiang-guo, 2001; Zhao et al., 2013).

Soil quality indicators (SQIs), particularly when structured through the minimum data set (MDS) approach, offer an effective means of reducing data redundancy and isolating the most meaningful parameters governing soil functionality (Bravo-Medina et al., 2021). When combined with multivariate statistical analyses such as principal component analysis (PCA), MDS allows for a systematic identification of key physicochemical and biological factors that regulate soil performance. Despite the strengths of these techniques, their predictive power remains limited unless they are integrated with environmental variables that capture terrain configuration, climate conditions, and spectral characteristics. Remote sensing imagery, especially multi-band and index-based products from platforms such as Landsat, has proven essential for characterizing land surface conditions, vegetation dynamics, and moisture gradients. Likewise, terrain derivatives derived from high-resolution DEMs illuminate geomorphological controls on soil redistribution and hydrological flow pathways. Yet, only a few studies have jointly incorporated these diverse data sources into a unified SD modeling framework for loess-based agricultural landscapes (Ghosh et al., 2025; Poudel et al., 2025).

With the advancement of machine learning, Random Forest (RF) models have become particularly valuable for SD prediction. Their ensemble architecture, ability to handle nonlinear and multicollinear relationships, resistance to noise, and strong generalization capability make them especially suitable for environmental systems where variable interactions are complex and multidimensional (Faramarzi et al., 2022; Wang and Sun, 2024). Furthermore, RF not only generates accurate predictions but also provides quantitative estimates of variable importance, offering insight into which environmental drivers exert the strongest influence on degradation patterns. Despite this, major gaps persist regarding the spatial uncertainty of such predictions. Without explicit uncertainty mapping, decision-makers risk misinterpreting model outputs, especially in heterogeneous terrains where prediction intervals can vary widely across space (Dadgar and Faramarzi, 2024; dos Santos et al., 2025; Khosravi Aqdam et al., 2022).

Given these knowledge gaps, there is an urgent need for a comprehensive, high-resolution, and uncertainty-aware SD assessment that integrates soil quality metrics with environmental predictors and advanced machine-learning methods. Accordingly, this study aims to develop an integrated SD assessment framework for the central Loess Plateau by: (1) identifying the most influential soil quality indicators through MDS and PCA; (2) quantifying the relationships between SD and environmental variables using correlation analysis; (3) modeling SD through an optimized RF approach; and (4) generating spatially explicit SD. By combining field-derived SQIs with multi-source environmental data and robust modeling, this research provides a refined understanding of degradation dynamics and supports the development of targeted, region-specific soil conservation strategies.

2 Materials and methods

2.1 The study area

The study area is situated in the central part of Shaanxi Province, China, within the Loess Plateau, extending approximately from 106°40′to 110°40′E longitude and 34°00′to 35°50′N latitude, and covering nearly 30,000 km² (Figure 1). The region experiences a continental climate that ranges from semi-arid to semi-humid, with an annual precipitation of about 400–600 mm, most of which occurs during the summer monsoon season (China Meteorological Administration, 2015–2024). This pronounced seasonality in rainfall strongly shapes erosion dynamics and influences soil stability. Geologically, the area is dominated by deep, fine-textured, and highly erodible Quaternary aeolian loess deposits, which control both soil formation processes and patterns of land use. Agricultural land constitutes the predominant land-cover type, and the main cultivated crops include winter wheat, maize, rapeseed, soybeans, and potatoes—crops that are well adapted to the properties of loess soils and the prevailing rainfall regime.

2.2 Field studies and laboratory analyses

In this study, a total of 350 surface soil samples, collected from a depth of 0–30 cm, were systematically obtained from agricultural lands, taking into consideration key environmental and pedological factors such as crop type, elevation, slope, parent material, and soil classification (Khosravi Aqdam et al., 2021). To adequately capture the spatial heterogeneity of the study area, sampling was conducted using an irregular network design, ensuring that both uniform and variable terrain features were represented (Cochran, 1977). Upon returning to the laboratory, the samples were air-dried at room temperature and carefully sieved through a 2 mm mesh to remove coarse fragments and achieve a consistent sample matrix. Subsequent laboratory analyses were carried out using standardized methods to determine a range of physical and chemical soil properties (Table 1). Soil sampling was conducted under the framework and authorization of the Shaanxi Province Cultivated Land Quality Big Data Application Technology Research Project (No. 202533).

TABLE 1. List and description of the soil properties used in this study.

2.3 Soil quality indicators

In this study, soil quality assessment was conducted using the minimum data set (MDS) approach. To select the most appropriate soil properties of farming land for MDS calculation, principal component analysis (PCA) was performed on the standardized dataset. Components with eigenvalues ≥1.0 and explaining more than 5% of the total variance were retained. Within each retained principal component, variables with absolute factor loadings within 10% of the highest loading value, or exceeding an absolute threshold of 0.70, were designated as high-loading candidates for the MDS. When more than one variable met these criteria within the same component, Pearson correlation analysis was conducted among them. If two candidate variables were strongly correlated (∣r∣≥0.70), only the attribute with the highest factor loading was retained to eliminate redundancy. Conversely, candidate attributes with ∣r∣<0.70 were treated as independent functional indicators and both retained in the MDS. In addition, soil attributes exhibiting distinct, critical pedological and agronomic roles in semi-arid loess environments (e.g., calcium carbonate equivalent for soil buffering and structural stability) were preserved to avoid the omission of key soil functions. For properties in which higher values indicate improved soil conditions, a positive linear scoring function was applied as follows (Equation 1):

S L = X i − X min X max − X min (1)

Conversely, for properties in which lower values reflect better soil condition, a negative linear scoring function was used (Equation 2):

S L = X max − X i X max − X min (2)

In these equations, SL represents the dimensionless linear score (ranging from 0 to 1), Xi denotes the observed value of the soil attribute, and Xmin and Xmax correspond to the minimum and maximum values of that attribute within the dataset (Zeraatpisheh et al., 2020). Following the linear standardization of all soil variables, the overall soil quality was quantified using the Integrated Quality Index (IQI). This index aggregates the weighted contribution of each standardized property according to the following expression (Equation 3) (Biswas et al., 2017).here, Wi represents the weighting factor assigned to each soil attribute, Ni denotes the corresponding linear score, and n is the total number of soil attributes included in the MDS. The weighting factor (Wi) for each indicator was calculated based on the proportion of total variance explained by the respective PC from which the variable was extracted. When multiple non-redundant variables were retained from a single PC, the explained variance of that component was partitioned among them proportionally according to their communality values, and all resulting weights were subsequently normalized so that ∑Wi = 1.00. All calculations related to linear scoring and IQI computation were performed using Excel and SPSS to ensure consistency and analytical rigor in the evaluation process (Askari and Holden, 2014).

I Q I = ∑ i = 1 n W i N i (3)

2.4 Environmental parameters

Table 2 summarizes the environmental variables used in this study, along with their units and spatial resolutions (Demattê et al., 2020). To ensure consistent spatial modeling, all datasets were projected and resampled to a common 30-m grid. It is important to emphasize that this operation was performed solely to facilitate spatial overlay and analysis and does not generate genuine 30-m information for any variable. The native spatial resolution of each predictor remains unchanged, and the 30-m grid should therefore be regarded as the spatial framework for mapping rather than a representation of true 30-m resolution.

TABLE 2. Environmental variables used for spatial modeling.

All variables were resampled to a 30-m mapping grid for spatial analysis; however, the native spatial resolution of the original data sources remained unchanged.

Normalized height refers to elevation values rescaled relative to the minimum and maximum elevation within the study area or a local moving window, thereby indicating the relative position of each pixel within the topographic range. In contrast, standardized height refers to elevation values transformed through statistical standardization based on the distribution of the data

The Synthetic Soil Image bands (SYS), derived from Landsat 8 Operational Land Imager (OLI) Collection 2, Level-2 surface reflectance images available through Google Earth Engine (GEE), were used as spectral predictors. These bands, Blue (Band 2), Green (Band 3), Red (Band 4), Near-Infrared (Band 5), Shortwave Infrared 1 (Band 6), and Shortwave Infrared 2 (Band 7), were obtained at a native spatial resolution of 30 m. To minimize the influence of vegetation phenology on soil spectral signatures, images acquired during the first 2 months of autumn (2014–2022), when soil surfaces are typically bare, were selected, and areas occupied by perennial crops, such as alfalfa, were masked. Missing data within the seasonal composites were filled using the ‘Close Gaps’ function and a Gaussian filter implemented in SAGA GIS version 7.2.0 (Hijmans et al., 2005).

Climatic variables, namely, annual mean temperature and annual precipitation, were extracted from the WorldClim Version 2 bioclimatic dataset (1970–2000 normals) at a native spatial resolution of 1 km. Topographic attributes were derived from the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) at a native resolution of 30 m. After hydrologically conditioning the DEM using the ‘Fill Sinks’ algorithm (Wang and Liu, 2006), terrain indices, including normalized height, standardized height, aspect, slope, mid-slope position, LS-factor, convergence index, plan curvature, profile curvature, relative slope position, topographic wetness index, and valley depth, were generated in SAGA GIS version 7.2.0 (Table 2).

Soil covariates, including clay content, sand content, bulk density, soil organic carbon, soil water content, and soil pH (H₂O), were obtained from the SoilGrids 2.0 database for the 0–30 cm depth interval. This global product provides predictions at a native spatial resolution of 250 m, and the 0–30 cm values were calculated as depth-weighted averages of the standard SoilGrids layers (Wang and Liu, 2006). The effective spatial resolution of these soil properties therefore remains 250 m, even though the layers were co-registered to the 30-m grid for integration with the other predictors (Table 2).

2.5 Soil degradation

Following the calculation of the SQI, soil samples were collected from ten sites located in grasslands and forest margins within the study area. These sites were selected as natural or minimally disturbed reference locations and were used to provide a practical baseline for evaluating the condition of agricultural soils. The reference sites were not intended to represent all natural soil conditions across the approximately 30,000 km² study area. Rather, they were used to characterize the available natural or relatively undisturbed soil conditions against which the quality of agricultural soils could be compared.

The reference samples were collected, prepared, and analyzed using the same procedures applied to the agricultural soil samples. This ensured that the comparison between agricultural and reference soils was based on consistent sampling and analytical methods. The mean SQI calculated from the ten reference sites was used as the reference SQ value (SQI _ref). SD at each agricultural sampling site was then estimated as the difference between the reference SQI and the SQI of the corresponding agricultural soil, as follows (Equation 4) (Ma et al., 2024).where, SQI _Ref represents the reference SQI under natural conditions, and SQI _i is the IQI-MDS at the ith sampling point. This approach allows the assessment of SD at each site, reflecting the deviation of current soil conditions from the natural reference state.

S D = S Q I R e f − S Q I i (4)

In this study, SD was interpreted as a measure of the relative decline in SQ compared with the available natural or minimally disturbed reference condition. It was not considered a direct measure of comprehensive SD across the entire study area. Given the limited number of reference sites, the results were interpreted with appropriate caution, and the reference-based SD values were used primarily to compare the relative spatial variation in soil degradation among the agricultural sites.

2.6 Modeling SD using RF model

In this study, SD was modeled using the RF algorithm, an ensemble machine-learning method that combines a large number of decision trees to generate stable and accurate predictions. The RF model was selected because of its ability to represent nonlinear relationships and complex interactions between soil degradation and environmental predictors. In addition, this algorithm is relatively robust to noise, can accommodate predictors with different scales and distributions, and generally has a lower risk of overfitting than individual decision-tree models. The response variable used in the modeling procedure was the SD value calculated from the IQI, whereas the predictor variables consisted of selected soil, climatic, topographic, spectral, and environmental covariates.

Before model fitting, all SD observations and environmental predictor variables were integrated into a single georeferenced database. The resulting database contained 350 observations corresponding to the agricultural soil-sampling points. To reduce the influence of multicollinearity on model interpretation and prediction stability, Pearson’s correlation analysis was performed among the predictor variables. When two predictors exhibited a strong correlation, defined as an absolute Pearson correlation coefficient greater than 0.85, one of the variables was removed. The variable retained from each correlated pair was selected according to its ecological relevance, data quality, interpretability, and relative contribution to the representation of the underlying soil-degradation process. This procedure was applied before model training to ensure that the final predictor set contained informative variables with limited redundancy.

The complete dataset was randomly divided into two subsets using a 70:30 training–validation ratio. Accordingly, 245 observations (70% of the total dataset) were assigned to the training subset, whereas the remaining 105 observations (30% of the total dataset) were retained as an independent validation subset. The training subset was used to develop the RF model, optimize its hyperparameters, and identify the most influential predictors. The validation subset was not used during model fitting or parameter optimization and was subsequently used to evaluate the predictive performance of the final model. A fixed random seed was applied during the data-partitioning procedure to ensure that the training and validation subsets could be reproduced.

Hyperparameter optimization was conducted using a grid-search procedure combined with ten-fold cross-validation within the training subset. This internal cross-validation procedure was used to compare different combinations of RF parameters and to identify the model configuration that minimized the cross-validated RMSE and maximized the coefficient of determination (R ²). The principal parameters considered during model optimization included the number of trees and the number of candidate predictors evaluated at each node split. The final model was constructed using 500 regression trees (ntree = 500) and three randomly selected predictors at each split (mtry = 3). These settings provided the most suitable balance between predictive performance, model stability, and computational efficiency within the training dataset (R Core Team, 2020).

The relative importance of the environmental predictors was assessed using the mean decrease in accuracy (MDA) and the mean squared error (%IncMSE). The MDA measure evaluates the reduction in predictive accuracy that occurs when the values of a given predictor are randomly permuted, whereas the MDG measure reflects the contribution of a variable to reducing prediction error during the construction of the regression trees. These two measures were used to rank the predictors and to identify the environmental factors that contributed most strongly to the spatial variation of SD. After hyperparameter optimization, the final RF model was fitted using the selected predictor variables and the 245 observations in the training subset. The predictive performance of this final model was then assessed using the independent validation subset of 105 observations.

The model-validation procedure used in this study was therefore based on a conventional random 70:30 training–validation split. The 245 training samples were used for model development and internal ten-fold cross-validation, while the 105 validation samples were reserved for the final assessment of model performance. Accordingly, all reported validation statistics were calculated from the 30% validation subset and were not derived from spatial block cross-validation (Malone et al., 2017).

2.7 Estimating uncertainty in RF predictions

To evaluate the spatial uncertainty of the RF predictions, we quantified the spread of the ensemble predictions generated across the 500 individual regression trees. It should be clarified that these empirical boundaries represent the ensemble prediction uncertainty range rather than statistically calibrated prediction intervals. For each pixel, the 5th and 95th percentiles of the 500 individual-tree predictions were computed as the lower and upper bounds of the ensemble uncertainty range, respectively. The overall ensemble prediction uncertainty spread was characterized by the difference between the 95th and 5th percentiles (ensemble prediction interval width, W₉₀), while the pixel-level standard deviation of tree predictions was calculated to quantify local model variance (Equations 5 and 6):

W 90 = y ^ 95 − y ^ 5 (5)

Standard deviation = 1 B − 1 ∑ b = 1 B y ^ b − y ¯ 2 (6)

Where, y ^ 95 y ^ 5 and are the 95th and 5th percentiles of the ensemble predictions, B = 500 is the number of trees, y ^ b is the prediction of the bth tree, and y ¯ is the ensemble mean prediction.

To quantitatively evaluate the reliability of these uncertainty boundaries against independent observations, the prediction interval coverage probability (PICP) was determined on the validation dataset (Equation 7):where, n = 105 is the number of validation samples, and ε i = 1 if the observed IQI/SD value falls between the 5th and 95th ensemble percentiles (y ^ 5 . i ≤y _*i*≤ y ^ 95 . i), and ε i = 0 otherwise. The calculated PICP value was evaluated to verify whether the ensemble spread adequately captured the validation data variability (Malone et al., 2011).

P I C P = 1 N ∑ i = 1 n ε i × 100 % (7)

2.8 Model validation and performance assessment

The predictive performance of the RF model was evaluated using three statistical indicators, namely, the coefficient of determination (R ²), root mean square error (RMSE), and concordance correlation coefficient (ρ_c), calculated between the observed and predicted SD values in the validation dataset (Equations 8–10). These metrics were selected to provide complementary information on model performance, including explained variance, magnitude of prediction error, and the strength of agreement between model outputs and observations. The R ² was used to quantify the proportion of variance in the observed SD values explained by the model, whereas RMSE measured the average magnitude of prediction errors. In addition, the ρ_c was calculated to assess the strength of the linear relationship between predicted and observed values. Higher values of R ² and ρ_c, together with lower RMSE values, indicate better model performance. where, Oi and Pi denote the observed and predicted SD values for sample i, respectively, O i ¯ and P i ¯ are the mean observed and predicted values, and n is the total number of validation samples (Malone et al., 2017).

R 2 = 1 ‐ ∑ i = 1 n O i ‐ P i 2 ∑ i n O i ‐ O i ¯ 2 (8)

RMSE = 1 n ∑ i = 1 n P i ‐ O i 2 (9)

ρ c = ∑ i = 1 n O i ‐ O i ¯ P i ‐ P i ¯ ∑ i = 1 n O i ‐ O i ¯ 2 ∑ i = 1 n P i ‐ P i ¯ 2 (10)

3 Results

3.1 Data analysis

The descriptive statistics of the analyzed soil properties are summarized in the supplementary material (Table 3). According to the results, the particle-size distribution showed considerable variation across the agricultural lands, with sand ranging from 5% to 73%, silt from 10% to 65%, and clay from 7% to 47%. This level of variability resulted in the presence of eight distinct textural classes, including clay, silty clay, silty clay loam, clay loam, silt loam, loam, sandy clay loam, and sandy loam. Such heterogeneity in soil texture is commonly associated with differences in hydrological and physicochemical behavior, which in turn affect plant growth and crop performance (Vinhal-Freitas et al., 2017).

TABLE 3. Summary statistics of selected attributes in the agricultural lands.

As shown in Table 3, clay exhibited the highest coefficient of variation (CV = 42.80%) among the mineral particle fractions. This elevated variability is likely influenced by the widespread use of flood irrigation in the region, which can redistribute fine particles such as clay through erosion and sediment transport (Roy et al., 2022). The mean saturated moisture content (39.05%) and total porosity (49%) suggest moderate water-holding capacity and a favorable proportion of macropores. The bulk density averaged 1.24 g cm^-3, indicating that most soils fall within a suitable to slightly compact condition, consistent with loam- and clay-rich textures.

The average soil pH was 7.88, with CCE values averaging 16.84%, classifying the soils as mildly to moderately alkaline. The mean CEC (17.12 cmol kg^-1) places the soils within the moderate capacity class based on standard criteria (Hazelton and Murphy, 2025). Soil organic carbon (SOC), an important indicator of soil quality (Herrick and Wander, 2018), averaged 1.12%, positioning the soils within the moderate SOC range considering the thresholds suggested by Jat et al. (2018). The values recorded for EC (0.57 dS m^-1), SAR (2.75), and ESP (3.23%) collectively indicate that the soils are neither saline nor sodic. Regarding nutrient status, total nitrogen (TN) was very low (0.029%), phosphorus was in the moderate range (14.25 mg kg^-1), and potassium levels were relatively high (415.30 meq L^-1), consistent with typical fertility patterns in irrigated agricultural lands. Among all measured properties, phosphorus exhibited the highest variability (CV = 386.86%), while TN showed the lowest (CV = 3.45%). According to the classification of Wilding (1985), these values indicate high spatial variability for P and low variability for TN across the study area.

3.2 MDS analysis

The results of the MDS analysis revealed substantial differences in the calculated weights of the various soil parameters, clearly indicating the relative importance of each variable in determining soil quality (Table 4). Parameters such as pH (0.242), SAR (0.246), CEC (0.155), SOC (0.119), K (0.098), and Clay (0.078) received the greatest weights, underscoring the decisive role of soil alkalinity, sodicity, organic carbon stock, cation retention, and fine particle distribution in shaping overall soil condition. In contrast, attributes such as SBD (0.031) and CCE (0.030) exhibited lower but functionally complementary contributions to the final SQI index.

TABLE 4. Scores and computed weights of variables in TDS and MDS approaches.

In the MDS, approach, weights were calculated based on the variance explained by the corresponding principal components, normalized across selected attributes so that the total sum equals 1.000.

The PCA indicated that the first six principal components explained more than 85% of the total data variance (cumulative variance = 85.2%), demonstrating that a substantial portion of the information contained in the dataset can be summarized within a limited number of components. PC1 accounted for 30.1% of the total variance, heavily dominated by CEC (0.94), Na⁺ (0.90), and Clay (0.76). Although CEC and Clay originated from PC1, their correlation was moderate (r∣<0.70), prompting the retention of CEC as the primary chemical exchange metric and Clay as the essential physical textural indicator, while Na⁺ was omitted due to its high correlation with SAR and CEC. PC2 (20.1% of variance) was governed by SOC (0.94) and porosity (n, 0.92); SOC was chosen due to its direct role in aggregate stability and carbon storage. PC3 (14.4% of variance) exhibited the highest loadings on pH (0.98), SAR (0.98), and ESP (0.98); SAR and pH were retained as critical indicators of sodicity and alkalinity hazards, respectively, whereas ESP was excluded due to strict collinearity with SAR (r > 0.95). PC4, PC5, and PC6 captured secondary structural, nutrient, and buffering variations, leading to the selection of SBD (PC4), K (PC5), and CCE (PC6). The combined application of PCA, communality partitioning, and correlation filtering enabled the construction of a robust, non-redundant MDS that accurately represents the soil system (Table 4; Table 5).

TABLE 5. The PC analysis results for soil quality factors.

The results of the MDS analysis revealed substantial differences in the calculated weights of the various soil parameters, clearly indicating the relative importance of each variable in determining soil quality Parameters such as SAR, pH, CEC, SOC, and Clay exhibited the highest weights, underscoring their dominant role in shaping overall soil quality. In contrast, CCE and SBD contributed less to the final SQIs. These findings suggest that chemical properties and indicators related to soil salinity and alkalinity exert the greatest influence, whereas soil physical attributes and exchange-related properties play a more complementary role. The PCA indicated that the first six principal components explained more than 85% of the total data variance (cumulative variance = 85.2%), demonstrating that a substantial portion of the information contained in the dataset can be summarized within a limited number of components.

PC1, which accounts for the largest proportion of variance (30.1%), primarily reflects cation exchange capacity and sodium-related soil properties, with CEC (0.94), Na+ (0.90), and Clay (0.76) exhibiting the highest factor loadings. PC2 (20.1% of the variance) was mainly associated with SOC and n, where the high loadings of SOC (0.94) and n (0.92) highlight the role of organic matter and porosity in soil quality. PC3 (14.4% of the variance) is influenced by pH and alkalinity-related indices such as SAR and ESP, emphasizing the role of salinity and acidity in soil quality. The remaining components (PC4–PC6) contribute less to the overall variance, yet the specific combinations of variables within these components reflect the inherent complexity of soil properties (Table 5). A comparison of the MDS and PCA outputs indicates that the key parameters influencing soil quality include SAR, pH, CEC, SOC, Clay, and K, while the remaining variables play complementary roles. The combined application of MDS and PCA allows for the identification of the most influential parameters.

3.3 Correlation of environment and SD

The analysis of Pearson correlation revealed that several environmental factors, are positively associated with the SD index. In other words, higher values of these variables correspond to higher SD scores. Conversely, factors including soil pH (in H₂O), soil crusting, moisture stress, mean annual precipitation, mean annual temperature, and SYS I bands exhibit a significant negative relationship with SD, indicating that increases in these parameters lead to a reduction in the SD value (Figure 2).

3.4 Relative importance of the environmental parameters

The relative importance analysis of the RF model indicated that several SYS spectral indices played a prominent role in distinguishing different levels of soil degradation. Among these variables, SYS I band 2 and SYS I band 6 exhibited the highest importance scores, largely due to their strong sensitivity to variations in surface moisture, declines in vegetation cover, and soil physical instability. These characteristics substantially enhanced the model’s ability to discriminate between varying degrees of degradation. Following these, SYS I band 7 and the Calcareous Sediment Index ranked next in importance, reflecting their direct association with erosional processes and surface alterations induced by shifts in soil properties. Additionally, variables such as Mid-Slope Position, soil pH (H₂O), and several other SYS bands fell within the medium-importance category, mainly capturing geomorphological variations and chemical characteristics of the soil. In contrast, parameters including Soil Water Content, Moisture Stress, and Soil Crusting demonstrated the lowest relative importance. This suggests that these factors either possessed limited discriminatory capability for separating degradation levels in the study area or exhibited overlapping patterns with more dominant variables, thereby reducing their unique contribution to the model (

3.5 Modeling SD using the RF model

The RF model demonstrated strong performance in predicting SD across the study area. Calibration using 245 training samples (70% of the total N = 350) resulted in a high coefficient of determination (R ² = 0.83), indicating that the selected environmental predictors captured the majority of the variability in SD. The model also exhibited low prediction error (RMSE = 0.031) and a high concordance correlation coefficient (ρc = 0.82), confirming both the accuracy of predictions and the agreement between observed and predicted values. Validation on 105 independent testing samples (30% of the dataset) showed slightly lower but still robust metrics (R ² = 0.76, RMSE = 0.034, ρc = 0.74), demonstrating that the RF model maintains good generalization capability when applied to unseen data. The minor reduction in performance during validation may reflect inherent spatial heterogeneity in soil loss processes or minor overfitting during model training (Table 6).

TABLE 6. Calibration and validation indices for modeling SD by RF model.

3.6 Spatial distribution of SD and its uncertainty

To map SD across the study area, the RF model was applied by integrating training data with relevant environmental covariates, producing a continuous spatial representation of SD (Figure 3). The resulting map reveals substantial spatial heterogeneity in SD intensity throughout the region. SD values ranged from 0.003 to 0.540, with higher values (purple and blue areas) indicating more severe soil loss, while lower values (yellow and light-colored areas) correspond to relatively stable conditions with minimal SD. Spatial patterns indicate that the central and southwestern portions of the study area exhibit the highest SD, whereas eastern and northern regions predominantly fall into low-SD categories, suggesting that SD distribution is strongly shaped by land use patterns and anthropogenic pressures.

The ensemble prediction uncertainty range maps of SD (Figure 4) illustrate a spatially heterogeneous and complex distribution of model predictions under the lower (5th percentile, Figure 4a) and upper (95th percentile, Figure 4b) limits, reflecting the spatial variation of model ensemble uncertainty across the study area. Quantitatively, the calculated PICP for the validation dataset was 88.6%, confirming that the vast majority of observed validation samples were successfully captured within this ensemble uncertainty range.

In the lower limit (Figure 4), most parts of the region appear in yellow and light tones, indicating lower predicted SD values under this conservative threshold. Only certain portions of the central and southwestern areas exhibit purple and blue patches, which signify higher predicted values even at the lower limit. Under the upper limit (Figure 4b), this spatial pattern undergoes a notable shift, with a substantial increase in the extent and intensity of purple and blue tones; extensive zones in the central and southwestern regions fall within high predicted SD classes. This contrast between the lower and upper limits highlights the spatial distribution of ensemble prediction uncertainty, with areas showing the largest differences between the two limits indicating higher model sensitivity to input variations and environmental conditions.

4 Discussion

The integrated assessment of SD in the central Loess Plateau revealed clear evidence of strong spatial heterogeneity driven by the interplay between soil physicochemical attributes, terrain configuration, and spectral–environmental conditions. The substantial variability observed in soil properties (Table 3), including particle-size distribution, pH, SOC, and CEC, demonstrates the dynamic nature of pedogenesis and the pronounced influence of land-use practices in the region. The presence of distinct textural classes, combined with high coefficients of variation for parameters such as clay and phosphorus, highlights the intrinsic complexity of loess-derived soils, which is consistent with previous studies in highly erodible aeolian landscapes. These heterogeneous conditions form the foundation on which SD processes evolve and intensify (Cao et al., 2021; Dai et al., 2018; Piotrowska-Długosz et al., 2016).

The MDS and PCA analyses (Table 4 and 5) further confirmed that a limited subset of soil attributes, specifically pH, SAR, SOC, CEC, and Clay, dominates the overall soil quality variability. These variables are strongly associated with soil chemical stability, cation retention capacity, nutrient dynamics, and structural stability, indicating that SD in the region is primarily governed by processes related to alkalinity, sodicity, carbon depletion, and textural degradation. The significant loadings of SOC and soil porosity (n) in PC2 underline the central role of organic matter dynamics and aeration/compaction status in SD, especially in cultivated loess soils where residue removal, oxidation, and limited biomass return can rapidly reduce carbon stocks and impair soil structure. Likewise, the prominence of pH and SAR in PC3 supports the interpretation that chemical degradation, particularly alkalization and sodicity in irrigated fields, is a major contributor to declining soil functionality (Andrews and Carroll, 2001; Khosravi Aqdam et al., 2023; Qiu et al., 2022).

Correlation analysis (Figure 2) revealed strong and consistent relationships between SD and environmental variables, emphasizing the combined effects of topography, spectral indices, and climatic conditions on degradation dynamics (Wei et al., 2020). Positive correlations with variables representing vegetation stress, thermal conditions, and specific spectral bands suggest that SD intensifies in areas experiencing moisture deficits, disturbance of surface aggregates, and reduced biological cover. In contrast, negative associations with mean annual precipitation, soil moisture, and soil pH indicate that relatively humid microenvironments and chemically stable soils are more resilient to degradation. These relationships illustrate the multidimensional controls on SD, where both soil-intrinsic factors and external environmental gradients collectively shape degradation (Mariano et al., 2018).

The variable-importance results of the RF model (Figure 5) provide deeper insights into the drivers of SD. The strong contribution of SYS I spectral bands, especially Bands 2 and 6, highlights the diagnostic value of surface reflectance indicators in capturing early degradation signatures such as crusting, SOC loss, and vegetation decline (Ben-Dor, 2002). The importance of the Calcareous Sediment Index similarly indicates that carbonate-rich soils, which are widespread in the Loess Plateau, respond sensitively to erosional processes and management-induced disturbances. The relevance of mid-slope position and soil pH underscores the role of landscape position and chemical equilibrium in shaping soil vulnerability. Conversely, the low importance of variables such as soil crusting and moisture stress may reflect either limited spatial variability in these parameters or their indirect influence in the presence of more dominant predictors (Zhang et al., 2019; Zhao et al., 2020).

The strong performance of the RF model (Table 6) further validates the suitability of ensemble-based machine-learning approaches for SD prediction in complex agroecosystems. High R ² (0.83 for training and 0.76 for validation) indicate that the model effectively captured nonlinear interactions between environmental variables and soil quality indices. The small difference between training and validation accuracies demonstrates minimal overfitting, suggesting that the selected covariates provide stable predictive signals across diverse landscapes. These results reinforce earlier findings that RF models are particularly advantageous for soil and environmental applications due to their resistance to noise, ability to handle multicollinearity, and robustness in high-dimensional datasets (Khosravi Aqdam et al., 2022).

The spatial SD map (Figure 5) clearly illustrates the concentration of severe SD in the central and southwestern sections of the study area (indicated by purple and blue tones), where intensive cultivation, sloped terrain, and lower vegetation cover interact to accelerate soil quality decline. In contrast, the eastern and northern zones exhibit lower SD (yellow and light-colored areas), likely due to gentler slopes, improved soil structural conditions, and better water-retention capacity. These spatial patterns align well with the region’s geomorphology and known agricultural dynamics, demonstrating that SD hotspots coincide with areas where anthropogenic pressures overlap with inherent soil fragility (Roy et al., 2022; Wei et al., 2020).

The uncertainty maps (Figure 4) provide an additional layer of interpretability. The lower-limit map (Figure 4a) shows widespread stability (lower predicted values) across most of the region, represented by predominant yellow tones. However, the upper-limit map (Figure 4b) reveals a substantial expansion of high-predicted SD values, particularly in the central and southwestern areas, where purple and blue tones become more extensive. This shift suggests that these zones are particularly sensitive to environmental fluctuations or data-input variability. The pronounced differences between the two limits highlight the value of uncertainty analysis in guiding spatial decision-making and prioritizing regions for field validation (Bravo-Medina et al., 2021; Cao et al., 2021).

The findings of this study indicate that the increasing intensity of cultivation and intensive land-use practices in the Loess Plateau constitute the primary drivers of SD in the study area. Observed soil patterns, including declines in organic carbon, chemical instability reflected in elevated SAR and pH variations, and reduced cation exchange capacity, clearly demonstrate that repeated cropping cycles, heavy tillage, imbalanced input applications, and successive harvests without sufficient nutrient replenishment have exceeded the natural resilience of the soils. The prominence of chemical indicators in the MDS and PCA analyses, along with the high weight of soil-quality-related variables in the RF model, indicates that intensive agricultural pressure has pushed soils toward structural degradation and disruption of biogeochemical processes. The increased uncertainty in some areas further highlights that these zones are particularly sensitive to management practices and require closer monitoring. Consequently, the results emphasize the urgent need for sustainable soil management strategies, including reducing cultivation intensity, increasing organic amendments, improving irrigation management, and implementing low-impact cropping systems, to mitigate degradation and maintain the long-term resilience of the agricultural ecosystem in the region.

5 Conclusion

The results of this study demonstrated that integrating SQIs, spectral–environmental variables, and the RF model provides an efficient framework for monitoring and predicting soil degradation in the Loess Plateau. Specifically, soil properties such as pH, SOC, SAR, and CEC were identified as the most influential factors controlling soil quality and long-term sustainability, while MDS and PCA confirmed that a limited number of key variables could effectively explain most of the spatial variability in soil quality. The RF model showed high predictive accuracy and strong performance in delineating critical degradation zones, and the spatial distribution maps enabled the identification of the most vulnerable areas across the region. In addition, the prediction interval maps provided important insight into the stability of model outputs and highlighted locations where additional sampling and monitoring may be required. Accordingly, the proposed framework can provide a reliable basis for risk assessment, land restoration planning, and data-driven environmental decision-making in vulnerable areas. Nevertheless, future studies should incorporate long-term temporal monitoring, biological soil indicators, higher spatial and temporal resolution remote sensing data, and advanced deep-learning approaches with spatiotemporal architectures to improve predictive reliability and better capture degradation dynamics. Moreover, integrative frameworks that account for future climate scenarios, land-use change, and multi-criteria decision-making should be developed to support climate-resilient policy formulation and sustainable land management strategies.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Yong Q. Tian, Central Michigan University, United States

Reviewed by: Yu Bai, Chang’an university, China

Quanping Zhang, Chinese Academy of Sciences (CAS), China