25-year, quarterly land change maps of China's Loess Plateau reveal long-term and substantial water-induced soil erosion mitigation

Unsustainable human activities have driven global ecological degradation. In China, decades of restoration policies have been implemented to reverse this trend in severely degraded regions with catastrophic soil erosion, transforming them into landscapes of ecological recovery. However, the evolution of soil erosion in these regions remains poorly quantified due to the absence of high-resolution, long-term, and high-frequency monitoring data. Here, to address this gap and provide a reliable spatiotemporal benchmark dataset, we conducted the first 10–30 m quarterly wall-to-wall land change mapping for China's flagship ecological restoration site: the Loess Plateau, based on the developed cross-temporal consistency-constraint deep learning framework. The dataset was generated using over 10 TB of Sentinel and Landsat imagery and documents land-cover dynamics across 100 quarterly time steps from 2000 to 2024, showing an overall accuracy of 81.44 % based on 40 000 annotated samples and 79.8 % for third-party validation sources. The resulting maps record pronounced land-cover dynamics, including forest expansion (+13 131 km²), cropland expansion (+28 095 km²), and bare land reduction (−65 029 km²) over the past decades. Furthermore, the produced dataset was combined with environmental factors to measure the 25-year water-induced soil erosion, where comparison with government survey data shows strong consistency, with a mean absolute error of 4.50 %. The dataset further illustrates that long-term ecological interventions have substantially reduced erosion intensity in the region by 30 % over the past 25 years, from 13.34 to 9.35 t (hm² a)⁻¹. Based on this benchmark, the long-term, fine-grained soil erosion becomes possible to estimate. The data-driven analysis indicates that current erosion is most severe in the central and southwestern Loess Plateau, and scenario modeling based on multiple factors suggests that optimized vegetation distribution – including grassland expansion and cropland-to-forest conversion – could potentially reduce future erosion intensity to 6.42 t (hm² a)⁻¹. This dataset provides a comprehensive benchmark for erosion mitigation in the Loess Plateau and its underlying drivers, providing critical insights for sustainable land management, ecological restoration, and policy development both in China and across fragile ecosystems worldwide. The land-cover maps and soil erosion maps are available at https://doi.org/10.57760/sciencedb.33656 (Cheng et al., 2026).

1 Introduction

Soil is the bedrock of human sustenance, providing more than 98 % of global food supply and hosting nearly 25 % of Earth’s biodiversity (Kraamwinkel et al., 2022). Meanwhile, Soil ecological functions result from the combined effects of aboveground vegetation and land-cover / land-use dynamics, both of which are highly sensitive to unsustainable human activities (Amin et al., 2025). As a consequence of their vulnerable environmental resilience, soil systems worldwide are experiencing unprecedented degradation (Amundson et al., 2015), a trend that threatens not only agricultural productivity but also the continued functioning of key ecosystem services (IPCC, 2019). Notably, existing evidence worldwide has identified soil erosion as the leading cause of global land degradation (UNCCD, 2022), a process that is further intensified by human-induced land-cover changes (Borrelli et al., 2017). Specifically, the spatial patterns of land-cover change, especially in regions susceptible to soil erosion, would strongly reflect the severity and distribution of soil erosion (Nut et al., 2021; Zhu et al., 2022). Based on the context, intensive monitoring of large-scale land-cover change and soil erosion at high temporal and spatial resolution can effectively reveal their underlying connection, filling the knowledge gap in fine-scale and long-term monitoring of soil erosion.

Over the past several decades, China has experienced severe soil erosion and associated environmental challenges (Wen and Zhen, 2020), particularly in the Loess Plateau, which has undergone sustained ecological decline characterized by widespread soil erosion, deforestation, and land degradation (Jiang et al., 2019). In particular, unsustainable agricultural practices applied to highly erodible soils have undermined agricultural productivity and disrupted local ecosystems in the region (Ge et al., 2020). Accordingly, to counteract these adverse environmental conditions, a suite of nationally mandated soil and water conservation measures has been implemented since the 1980s, including land-cover restructuring and large-scale vegetation restoration (Liu et al., 2008). As a result, these interventions have yielded substantial ecological improvements. Based on statistical records, vegetation coverage in the Loess Plateau has increased from 31.6 % in 1999 to 63.6 % in 2019 (Xie et al., 2015). This ecological transformation has fundamentally reshaped the sediment regime of the Yellow River, which drains much of the plateau, contributing to an approximately 90 % reduction in its annual sediment load (Wang et al., 2015). In general, the governance strategies developed in the Loess Plateau over the past decades offer insights of broad relevance for soil erosion mitigation in arid and semi-arid regions worldwide (Feng et al., 2013; Yang et al., 2024). Nevertheless, systematic evaluation of erosion mitigation and its ecological consequences remains hindered by the persistent lack of long-term, large-scale satellite monitoring data products (Pettorelli et al., 2014).

Time-series satellite imagery constitutes a vital resource for monitoring large-scale land-cover dynamics because of its extensive spatial coverage and frequent revisit capability (Li et al., 2026). Thus, leveraging remote-sensing technologies to conduct a spatiotemporal modeling of land-cover patterns (Huang et al., 2025; Tang et al., 2025) has become increasingly important for soil-erosion research, such as disaster-induced erosion (Di et al., 2010) and agricultural-driven erosion (Rafiei-Sardooi et al., 2025). However, limitations in the spatial coverage and temporal continuity of existing land-cover datasets constrain their utility for large-scale long-term monitoring of soil erosion dynamics (Kodl et al., 2024; Zweifel et al., 2019). Specifically, current map products are typically characterized by a trade-off between spatial resolution and temporal continuity. Products offering high spatial resolution usually lack long-term consistency, whereas long-term time-series products often operate at coarse spatial resolutions. In general, current approaches fall into two categories:

1. High-spatial-resolution mapping. This category of approaches focuses on land-cover mapping conducted at fine spatial resolutions, often finer than 10 m (Gong et al., 2019; Zhang et al., 2021; Zanaga et al., 2022) and in some cases reaching approximately 1 m (Li et al., 2023, 2024). Such datasets excel at capturing fine-scale landscape patterns and support detailed analyses, especially in urban environments (Balhas et al., 2025). However, they typically provide only a single or a few temporal snapshots, which limits their ability to represent long-term land-surface dynamics.

2. High-temporal-resolution mapping. In contrast, high-temporal-resolution datasets allow systematic analyses of land-cover dynamics and enable interpretation of temporal trends. These datasets feature revisit cycles ranging from several days to months and have been adopted in hydrological modeling (Bohn and Vivoni, 2019), vegetation monitoring (Liu et al., 2024), and agricultural management (Zhang et al., 2022). However, since they are commonly derived from satellite platforms such as MODIS and AVHRR, with spatial resolutions between hundreds and thousands of meters (Friedl et al., 2010), they lack the capacity to depict fine-scale landscape heterogeneity and thus fall short for detailed soil erosion assessments.

Monitoring soil erosion dynamics in the Loess Plateau hinges on the availability of high-resolution and dense temporal mapping data. Specifically, high spatial resolution is indispensable for resolving the heterogeneous landscape of the Loess Plateau, where erosion processes are modulated by fine-scale features such as fragmented croplands, narrow gullies, and localized vegetation restoration zones (Yang et al., 2019). Meanwhile, dense temporal observations are required to capture rapid intra-annual transitions driven by precipitation seasonality, vegetation phenology, and cropping cycles (Sharma and Ehlers, 2023), with quarterly granularity being particularly valuable because most annual soil loss occurs within a narrow seasonal window. Consequently, integrating fine spatial detail with quarterly temporal coverage enhances the ability to delineate erosion patterns and interpret their spatiotemporal evolution, thereby supporting long-term monitoring and process-based understanding. Nevertheless, conducting quarterly high-resolution land-cover mapping is fundamentally limited by the difficulty of obtaining reliable annotations. In practice, producing multi-temporal labels typically requires expert interpretation (Hermosilla et al., 2022), while high-resolution mapping demands detailed labeling to describe inter-pixel relationships accurately (Li et al., 2024). Therefore, making large-scale, time-series annotation is prohibitively time-consuming and labor-intensive. To address these limitations, leveraging existing high-resolution map products as supervisory signals and extending them along the temporal dimension is an effective strategy, alleviating the need for extensive manual labeling.

Predicting long-term quarterly land-cover maps from existing high-spatial-resolution products inherently falls within a semi-supervised learning paradigm, where only a limited set of labeled samples is available relative to the abundance of multi-temporal observations. In this context, a central difficulty arises from the pronounced appearance variability observed in images of the same location across different times. Although only a small proportion of the Earth’s land surface undergoes actual land-cover change each year (Zhu et al., 2017), multi-temporal satellite imagery still exhibits substantial fluctuations arising from illumination differences, atmospheric conditions, and sensor-related inconsistencies (Lin et al., 2023). In turn, this variability obscures true land-cover signals and potentially misleads models into focusing on transient, non-semantic differences. Consequently, developing methods capable of learning stable and temporally consistent representations, while remaining sensitive to meaningful land-cover transitions, is essential for capturing the dynamic interaction between soil erosion and its potential divers.

In this work, to address the inherent trade-off between spatial resolution and temporal continuity in large-scale, long-term mapping, we presented a 25-year quarterly land-cover dataset for the Loess Plateau (LP-QLC10) through the combination of over 10 TB of multi-source satellite imagery and a proposed Consistency-Constrained Time-Series (CCTS) framework that enables high-resolution land-cover information to be extended along the temporal dimension (Fig. 1). The resulting maps provide continuous coverage from 2000 to 2024 across 648 700 km² of the Loess Plateau, with a three-month interval, resulting in 100 phases in total. Meanwhile, spatial resolution is maintained at 10 m from 2016 onward and 30 m for earlier years, consistent with the native characteristics of the contributing sensors. Furthermore, to fill the knowledge gap of soil erosion monitoring in China's Loess Plateau, we integrated multi-source data – including precipitation, soil properties, the Digital Elevation Model (DEM), slope length and steepness factors, and Fractional Vegetation Cover (FVC) – to derive environmental and biophysical factors of erosion. Then, in accordance with the standards issued by the Ministry of Ecology and Environment of China, the spatiotemporal variations in water-induced soil erosion were quantified from 2000 to 2024. Finally, the resulting spatiotemporal dataset enables systematic investigation of erosion patterns, temporal evolution, and interactions among environmental factors, thereby providing a unified data basis for addressing three knowledge gaps: (1) What are the dominant land-cover change patterns over the past 25 years in the Loess Plateau? (2) Where and when did the erosion occur? (3) How did the erosion dynamically interact with environmental factors?

Figure 1 Workflow of the 25-year land-cover mapping and soil-erosion estimation. A hybrid model combining convolutional neural networks and Vision Transformers (CNN-ViT) was employed to generate quarterly land-cover maps. Subsequently, five environmental variables were integrated to estimate soil erosion from 2000 to 2024. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google.

2 Study area and data

2.1 Study area

The Loess Plateau, as shown in Fig. 2a, is located in north-central China between latitudes 32–41° N and longitudes 107–114° E, covering more than 648 700 km². As the largest loess accumulation region in both China and the world, it is distinguished by thick loess deposits and pervasive soil erosion. Elevations generally range from 500 to 3000 m, with mountains, hills, and plateaus constituting the dominant landforms (GPRC, 2010). The soils of the region originate from wind-blown loess that has accumulated since the Quaternary period, resulting in a loose texture that is highly susceptible to water erosion (Qiang-Guo, 2001). Besides, mean annual precipitation ranges from 150 to 800 mm, with 55 %–78 % of rainfall concentrated between June and September (Liu et al., 2011). Consequently, water-induced erosion constitutes the principal erosion type on the Loess Plateau, accounting for approximately 86 % of the total erosion area (Gao et al., 2016). The combined influence of climate, topography, soil properties, and long-term human activities has made the Loess Plateau one of the most severely eroded regions in the world.

Figure 2 Illustration of the study area and using data. (a) The Loess Plateau and its six subregions. (b) The training data employed for land-cover classification. (c) Annotated samples used for validation across five temporal phases. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google, Basemap data: Esri 2026, USGS, World Hillshade | Powered by Esri.

Since the 1980s, the Chinese government has undertaken substantial efforts to control soil erosion and restore the ecosystem. Beginning in 1999, several large-scale ecological restoration initiatives were implemented to accelerate land rehabilitation, including the Grain for Green Program (GGP) and the Three-North Shelterbelt Forest Program (TSFP) (Wen and Zhen, 2020). Operating for nearly a quarter century, these programs have markedly improved ecological conditions across the Loess Plateau. According to government planning (GPRC, 2010), the Loess Plateau is divided into six comprehensive management areas: the irrigated agricultural region (IA), the sandy land and desert region (SL), the loess hilly and gully region (HG), the loess tableland and gully region (TG), the river valley plain region (RV), and the rock-soil mountainous region (RM).

2.2 Time-series images and multi-source data for mapping and validation

2.2.1 Satellite images (2000–2024) collected from Sentinel and Landsat missions

The 25-year imagery used for high-temporal-frequency land monitoring was compiled from multi-source satellite platforms, primarily Sentinel-2 (Drusch et al., 2012) and Landsat series missions (Chander et al., 2009; Roy et al., 2014). Sentinel-2 mission began launching satellites in 2015, comprising 13 spectral bands with spatial resolutions ranging from 10 to 60 m. For land-cover mapping, five representative bands were selected, including blue, green, red, near-infrared (NIR), and shortwave infrared-1 (SWIR1). These bands are vegetation-sensitive (Yu et al., 2019) and comparable to those of the Landsat series. Sentinel-2 images from 2016 to 2024 were accessed via the Google Earth Engine platform (Gorelick et al., 2017), where cloud removal, atmospheric correction, and median compositing were applied to generate one multispectral image every three months, yielding quarterly 10 m composites. To extend the dataset to the earlier period (2000–2015), imagery from Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) was incorporated, all acquired at 30 m spatial resolution. The same five spectral bands (blue, green, red, NIR, and SWIR1) were extracted, and the imagery were atmospherically corrected, cloud-masked, and converted to surface reflectance following the official guidelines (U.S. Geological Survey, 2020). In total, a dataset of quarterly composites from 2000 to 2024 was generated, comprising five-band imagery, with a spatial resolution of 30 m before 2016 and 10 m thereafter. The acquisition periods and the number of images used from multiple platforms are detailed in Fig. A2. Finally, all images were normalized by calculating band-wise mean and standard deviation to provide standardized inputs for the training and inference of the land-cover mapping framework.

2.2.2 Training and validation data collected from well-established datasets

To overcome the substantial effort required for large-scale and multi-temporal high-resolution annotation, we guide the multi-temporal learning process using verified single-temporal 10 m land-cover labels. Accordingly, we employ the European Space Agency (ESA) WorldCover v200 product (Zanaga et al., 2022), which was generated from Sentinel-1 and Sentinel-2 imagery for 2021 and reports a global overall accuracy of 77 %. The product provides 11 land-cover classes, 8 of which are present in our study area. Moreover, it has been well assessed by the land-cover community and was applied in downstream analyses, including global green-space studies (Chen et al., 2022), urban land estimation (Chakraborty et al., 2024), and cropland cultivation dynamics (Wuyun et al., 2025).

The accuracy of the predicted land-cover maps was evaluated through point-based validation using three datasets: one multi-temporal reference set constructed in this study and two publicly available validation sets. First, the multi-temporal reference set was constructed by random sampling and expert annotation across the Loess Plateau, as shown in Fig. 2c. Specifically, over 8000 points were randomly selected and visually interpreted at five time steps (2001, 2006, 2011, 2016, and 2021) by five remote-sensing experts, with cross-validation performed to ensure label reliability. This process yielded more than 40 000 manually annotated multi-temporal sample points, which were used to assess the accuracy of land-cover maps and to evaluate their temporal consistency.

Additionally, to provide an objective and fair evaluation, two additional public validation sets were employed. The SinoLC-1 validation set (Li et al., 2023) contains 106 852 manually annotated sample points distributed across China, most of which were labeled in 2020 into 12 land-cover categories, encompassing the eight land-cover categories relevant to our study area. Within the Loess Plateau, 7382 points are available. Furthermore, Geo-Wiki global land-cover reference set was also used (Fritz et al., 2017), which consists of 151 942 crowdsourced land-cover points collected worldwide and primarily labeled in 2012. The annotations include ten land-cover categories, which encompass the eight categories present in the Loess Plateau, where 1183 points are available for validation.

2.3 Environment and topographic data for soil erosion assessment

The intensity of water-induced soil erosion is governed by several environmental and biophysical factors, including rainfall erosivity, soil erodibility, topographic characteristics, cover management, and support practices. In this study, multi-source datasets were integrated to quantify these factors.

• Rainfall data. Daily rainfall data were obtained from the Climate Hazards Center Infrared Precipitation with Stations (CHIRPS) dataset (Funk et al., 2015). CHIRPS provides a global rainfall record from 1981 to the present by merging satellite-based thermal infrared observations with ground station measurements, resulting in a 0.05° gridded land-based time series. These rainfall data were used to derive rainfall erosivity.

• Soil composition data. Soil property layers were derived from OpenLandMap, which supplies global predictions of soil composition at six standard depths (0, 10, 30, 60, 100, and 200 cm) with a spatial resolution of 250 m (Hengl, 2018; Hengl and Wheeler, 2018). In this study, sand, silt, clay fractions, and organic carbon content in surface soils were extracted through the Google Earth Engine platform to calculate soil erodibility.

• Slope length and steepness (LS) factor. The combined LS factor, which quantifies the effect of topography characteristics on soil erosion, was represented by a 30 m resolution dataset developed for mainland China (Zhao et al., 2025).

• Elevation data. Elevation information was obtained from the Copernicus Global Digital Elevation Model (GLO-30), published by the ESA. This DEM provides global elevation coverage at 30 m resolution and was used to compute the support practice factor.

In addition, the China Soil and Water Conservation Bulletin, a official report compiled by the Ministry of Water Resources of China, was used to validate the accuracy of the derived erosion estimates. The bulletin provides authoritative annual statistics on soil erosion across the Loess Plateau and other regions of China since 2018, serving as an official benchmark for national soil conservation monitoring. Thus, data from 2018 onward were utilized in this study to quantitatively assess the consistency between our erosion estimates and officially reported values.

3 Methods

In semi-supervised learning, consistency regularization (Miyato et al., 2019) encourages the model to learn invariant but discriminative representations across varying input conditions. Meanwhile, enforcing such stability enables abundant unlabeled samples to serve as effective supervisory signals, thereby narrowing the performance gap between labeled and unlabeled data. Consequently, consistency-based methods have demonstrated strong efficacy across a range of image interpretation tasks, including image segmentation (Chen et al., 2024), object detection (Hua et al., 2023), and change detection (Yang et al., 2023). Furthermore, in the context of long-term land-cover mapping, this principle becomes particularly valuable, because stabilizing both feature embeddings and predictions not only mitigates overfitting to the single-temporal annotated subset but also promotes the learning of class-discriminative representations that generalize across imaging conditions and temporal phases.

Inspired by these advances, we integrate both inter-temporal and intra-temporal consistency constraints into time-series land-cover mapping. As illustrated in Fig. 3, the proposed Consistency-Constrained Time-Series (CCTS) framework is designed to enhance the temporal coherence and semantic stability of land-cover representations. Specifically, images acquired at different time phases over the same region are treated as style perturbations resulting from variations in observation conditions, upon which inter-temporal consistency constraint is imposed at the feature representation. Moreover, intra-temporal consistency is achieved by perturbing both the inputs and their latent representations, while prediction alignment is maintained through pseudo-label supervision, enabling the abundant unlabeled data to provide informative signals for model optimization. Together, by jointly optimizing these constraints, the model learns temporally robust and semantically consistent representations across multiple periods. Building upon this framework, we produced the LP-QLC10 dataset (Cheng et al., 2026) – the first quarterly land-cover record of the Loess Plateau – comprising 100 phases from 2000 to 2024. Furthermore, by integrating the LP-QLC10 dataset with multi-source environmental and topographic datasets, we derived the determinant factors of soil erosion. In accordance with the standards issued by the Ministry of Ecology and Environment of China, we produced quarterly estimates of water induced soil erosion across the Loess Plateau from 2000 to 2024. Finally, these results were analyzed to characterize the long-term spatiotemporal evolution of soil erosion.

Figure 3 Overall workflow of the long-term land-cover mapping and soil erosion assessment, comprising two main steps: (1) long-term land-cover mapping using the proposed CCTS framework, and (2) soil erosion estimation based on quarterly datasets. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google, Basemap data: Esri 2026, USGS, World Hillshade | Powered by Esri.

3.1 Quarterly time-series land-cover mapping

3.1.1 Data preprocessing and organizing

To construct reliable and class-balanced training tiles for time-series land-cover mapping, we first divided the Loess Plateau region into non-overlapping tiles of 256×256 pixels based on the land-cover products provided by the ESA WorldCover (Zanaga et al., 2022). According to the ESA's classification system, eight land-cover classes were identified in the study area: “Forest”, “Shrubland”, “Grassland”, “Cropland”, “Built-up”, “Bare land”, “Water”, and “Wetland”, with class definitions consistent with the ESA product. The detailed definition is explained in Table A2. To achieve class balance, we randomly sampled 1800 tiles for each class, individually. As shown in Fig. 2b, after filtering out tiles that overlapped with regions outside the study area, 11 013 training tiles (each with 256×256 pixels) were retained. Subsequently, the time-series Landsat and Sentinel imagery was then spatially aligned and cropped to match these tiles. Consequently, imagery acquired in 2021 was combined with ESA WorldCover tiles to construct labeled training pairs, and images from all other years were used as unlabeled training material.

3.1.2 Cross-temporal learning based on semi-supervised and consistency-constraint strategies

The preprocessed training data consists of limited labeled samples and a large amount of unlabeled data. Thus, exploiting the semantic relationships embedded within these multi-temporal data becomes essential for accurately inferring long-term land-cover maps. In our previous work (Cheng et al., 2024), we proposed a method that detects land-cover changes by constraining content-feature consistency within unchanged regions, thereby enabling the model to distinguish semantic variations associated with actual land-cover transitions. This approach has demonstrated potential for capturing cross-temporal semantic features in disaster assessment applications. In this work, as illustrated in Fig. 4, we extend consistency constraints to multi-temporal settings. Specifically, let ${\mathbf{D}}^{\mathrm{l}}=\mathit{\left\{}\right({\mathit{x}}_n^{\mathrm{l}},{\mathit{y}}_n\left)\mathit{\right\}}$ denote the labeled image-label pairs and ${\mathbf{D}}^{\mathrm{u}}=\mathit{\left\{}{\mathit{x}}_{n,t}^{\mathrm{u}}\mathit{\right\}}$ denote the unlabeled images, where n and t represent spatial and temporal indices, respectively. Meanwhile, to generate land-cover maps from satellite imagery, we employ the High-Resolution Transformer (HRFormer) (Yuan et al., 2021) as the segmentation model, owing to its strong performance in dense prediction tasks. In addition, to effectively incorporate both intra-temporal and inter-temporal consistencies, the training procedure was organized into three consecutive stages that jointly integrate labeled and unlabeled samples. The following subsections detail each stage of the training process.

Figure 4 Training workflow of the CCTS framework. A three-stage training process is designed to incorporate inter-temporal and intra-temporal consistency. Imagery © 2024 European Space Agency, Map data © 2026 Google.

Inter-temporal consistency in feature levels

In this stage, we promote inter-temporal consistency by regulating the behavior of the encoded feature representations. Multi-temporal land-cover mapping is characterized by a large proportion of unchanged areas. Therefore, images acquired at different times but covering the same location should share similar semantic features. To exploit this property, all training samples, including both labeled and unlabeled images, are input to the segmentation encoder to obtain high-level feature embeddings. Subsequently, these embeddings are organized according to their spatial positions: those originating from the same location but at different times are considered positive pairs, while embeddings from different locations are treated as negative pairs. To further enhance discrimination, we draw positive pairs closer in the representation space while separating negative pairs. Accordingly, we apply contrastive learning and optimize the encoder using the InfoNCE loss (He et al., 2020), which is formulated as:

$$\begin{array}{}\text{(1)}& {\mathcal{L}}_{\mathrm{inter}}=-\log {\displaystyle \frac{\exp \left(\cos \right({\mathit{e}}_{n,t},{\mathit{e}}_{n,s})/\mathit{\tau })}{{\sum }_{(n,t)\ne (m,s)}\exp \left(\cos \right({\mathit{e}}_{n,t},{\mathit{e}}_{m,s})/\mathit{\tau })}},\end{array}$$

where the pair $({\mathit{e}}_{n,t},{\mathit{e}}_{n,s})$ denotes a positive feature embedding pair, whereas $({\mathit{e}}_{n,t},{\mathit{e}}_{m,s})$ denotes a negative pair. The operator cos(⋅) computes the cosine similarity, and τ is a hyper-parameter that modulates the sensitivity of the loss to similarity differences. This loss function encourages the encoder to align positive-paired embeddings closely in the representation space while driving negative pairs apart, thereby guiding the encoder to learn inter-temporal and semantically consistent representations.

Basic model training with full supervision

In this stage, the segmentation model is initialized with the pretrained encoder obtained from Stage 1, thereby equipping it with the capability of semantic feature extraction. Subsequently, the entire model is trained on the labeled dataset D^l. Following common practice in semantic segmentation, we adopt the cross-entropy loss as the supervision signal, which guides the model to perform dense land-cover prediction. The loss function ℒ_seg is formulated as:

$$\begin{array}{}\text{(2)}& {\mathcal{L}}_{\mathrm{seg}}=-{\displaystyle \frac{\mathrm{1}}{M}}\sum _{j=\mathrm{1}}^M\sum _{i=\mathrm{1}}^C{y}_{j,i}\log \left(p_{j,i}^l\right),\end{array}$$

where C denotes the number of land-cover classes (set to eight in this study), M is the number of pixels in the image, y_j,i is the ground-truth label of pixel j for class i, and $p_{j,i}^l$ is the predicted probability that pixel j belongs to class i. This pixel-wise supervision enables the model to learn class-discriminative representations for accurate land-cover prediction.

Intra-temporal self-refinement

The first two stages of training equip the model to extract robust semantic features and generate land-cover maps. However, variations in imaging conditions and sensor characteristics introduce disturbances that challenge the temporal consistency of predictions. To address this issue, the final stage adopts a intra-temporal self-refinement strategy. Specifically, labeled data are used to provide direct supervision and prevent performance degradation during self-refinement, while perturbation-based augmentation is applied to both the unlabeled images and the encoded features to increase data diversity, formulated as:

$$\begin{array}{}\text{(3)}& {\mathit{p}}^{\mathrm{ip}}=\mathrm{Dec}\left(\mathrm{Enc}\right({\mathcal{P}}_{\mathrm{i}}\left({\mathit{x}}^{\mathrm{u}}\right)),{\mathit{p}}^{\mathrm{fp}}=\mathrm{Dec}({\mathcal{P}}_{\mathrm{f}}\left(\mathrm{Enc}\right({\mathit{x}}^{\mathrm{u}}\left)\right),\end{array}$$

where 𝒫_i(⋅) and 𝒫_f(⋅) denote image-level and feature-level perturbations, Enc(⋅) and Dec(⋅) represent the encoder and decoder of the segmentation model, and p^ip and p^fp correspond to predictions under image and feature perturbations, respectively. Image perturbations include flipping, rotation, shifting, cropping, contrast adjustment, and Gaussian blur, whereas feature perturbations are implemented by randomly dropping half of the feature channels. Additionally, the non-perturbed prediction p^u serves as the pseudo-label. In turn, training is guided by a confidence-weighted cross-entropy loss, defined as:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathcal{L}}_{\mathrm{ic}}={\displaystyle \frac{\mathrm{1}}{M}}\sum _{j=\mathrm{1}}^M\mathbb{1}(p_j^{\mathrm{u}}\ge \mathit{\lambda })\mathrm{H}(p_j^{\mathrm{u}},p_j^{\mathrm{ip}}),\\ & {\mathcal{L}}_{\mathrm{fc}}={\displaystyle \frac{\mathrm{1}}{M}}\sum _{j=\mathrm{1}}^M\mathbb{1}(p_j^{\mathrm{u}}\ge \mathit{\lambda })\mathrm{H}(p_j^{\mathrm{u}},p_j^{\mathrm{fp}}),\end{array}\end{array}$$

where H(⋅) denotes the cross-entropy function, M is the number of pixels in the image, λ is the confidence threshold, and 𝟙(⋅) is the indicator function used to filter high-confidence pseudo-labels. ℒ_ic and ℒ_fc represent the consistency losses under image-level and feature-level perturbations, respectively. In total, the intra-temporal self-refinement loss is formulated as:

$$\begin{array}{}\text{(5)}& {\mathcal{L}}_{\mathrm{intra}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathcal{L}}_{\mathrm{seg}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathcal{L}}_{\mathrm{ic}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\mathcal{L}}_{\mathrm{fc}}.\end{array}$$

In this stage, multi-level perturbations improve the model’s robustness by forcing it to focus on intrinsic semantic information rather than cross-temporal variations. At the same time, the pseudo-label supervision mechanism allows the abundant unlabeled images to provide additional supervisory signals. Overall, the proposed three-stage training framework progressively strengthens the model’s capacity to exploit both labeled and unlabeled multi-temporal data. By optimizing the three stages, the model effectively leverages cross-temporal semantic information to produce reliable and temporally consistent land-cover maps.

3.1.3 Time-series mapping with the well-trained model

To produce the time-series LP-QLC10 dataset, a crop-predict-merge workflow was adopted to process the satellite data. Firstly, as illustrated in Fig. 3b, satellite images over the Loess Plateau were chronologically organized to provide quarterly coverage from 2000 to 2024, totaling about 10 TB of data across 100 phases. Subsequently, the images were cropped into non-overlapping tiles of 256×256 pixels, which were batch-processed by the three-stage trained CCTS to generate land-cover predictions. Then, the prediction tiles were mosaicked into seamless maps at 30 m resolution for 2000–2015 and 10 m for 2016–2024, ultimately producing 100 quarterly phases covering the Loess Plateau from 2000 to 2024. In addition, a post-processing temporal refinement was performed using a five-step sliding window, which replaces each pixel’s label with the most frequent class in its local temporal neighborhood, thereby enhancing the temporal consistency of the final land-cover maps. As a result, this workflow yielded the LP-QLC10 dataset, a large-scale, seamless, and temporally dense land-cover record, forming a robust foundation for long-term soil erosion analysis.

3.2 Long-term soil erosion assessment

3.2.1 Multi-criteria assessment factors

The time-series land-cover maps enable the investigation of long-term soil erosion dynamics in the Loess Plateau, thereby underscoring the broader research significance of the mapping. Accordingly, in compliance with the National Ecological Environment Standards of China (MEE China, 2021b), water-induced soil erosion of the Loess Plateau was estimated using the Revised Universal Soil Loss Equation (RUSLE) (Renard et al., 1991; Zhang et al., 2008) method:

$$\begin{array}{}\text{(6)}& A=R\times K\times L\times S\times C\times P,\end{array}$$

where A denotes the soil loss per unit area for a given calculation period (t (hm² a)⁻¹). In this study, all time-varying factors were estimated at quarterly intervals so as to capture intra-annual variability. Since soil erosion is cumulative over time, the annual soil erosion used for subsequent assessment was obtained by summing the four quarterly estimates within each year. This strategy preserves seasonal information during the calculation process while ensuring consistency with annual-scale evaluation of long-term erosion dynamics. The subsequent subsections provide detailed descriptions of all factors involved in the calculation, including the rainfall erosivity factor (R), soil erodibility factor (K), topographic factors (L and S), cover management factor (C), and support practice factor (P).

Rainfall erosivity factor (R)

The R factor reflects the potential of rainfall and surface runoff to cause soil erosion within a unit period. In this study, R was calculated following the method specified in the Technical Specification for Investigation and Assessment of National Ecological Status (MEE China, 2021a), an official guideline issued by the Ministry of Ecology and Environment of China. The calculation is given by:

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{‾}}{R}=\sum _{k=\mathrm{1}}^{\mathrm{24}}{\stackrel{\mathrm{‾}}{R}}_k,{\stackrel{\mathrm{‾}}{R}}_k={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\sum _{j=\mathrm{0}}^m(\mathit{\alpha }\cdot P_{i,j,k}^{\mathrm{1.7265}}),\end{array}$$

where $\stackrel{\mathrm{‾}}{R}$ denotes the multi-year average annual rainfall erosivity (MJ mm (hm² h a)⁻¹), ${\stackrel{\mathrm{‾}}{R}}_k$ represents the erosivity of the kth half-month ($k=\mathrm{1},\mathrm{\dots },\mathrm{24}$), and α is an empirical coefficient (0.3937 for the warm season and 0.3101 for the cold season). $P_{i,j,k}$ denotes the rainfall (mm) of the jth erosive rainfall event in the kth half-month of year i. Daily rainfall data were obtained from CHIRPS, as described in Sect. 2.3. Moreover, to support quarterly soil-erosion assessment, the R factor was computed independently for each of the four quarters.

Soil erodibility factor (K)

The K factor characterizes the susceptibility of soil particles to detachment and transport by water. In this study, it was computed in accordance with the national technical specifications (MEE China, 2021a) and is expressed as:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& K_{\mathrm{EPIC}}=\mathit{\{}\mathrm{0.2}+\mathrm{0.3}\exp [-\mathrm{0.0256}{m}_{\mathrm{s}}(\mathrm{1}-m_{\mathrm{silt}}/\mathrm{100})\left]\mathit{\right\}}\\ & \times [m_{\mathrm{silt}}/(m_{\mathrm{c}}+m_{\mathrm{silt}}){]}^{\mathrm{0.3}}\\ & \times \mathit{\{}\mathrm{1}-\mathrm{0.25}{m}_{\mathrm{o}}/[m_{\mathrm{o}}+\exp (\mathrm{3.72}-\mathrm{2.95}{m}_{\mathrm{o}})\left]\mathit{\right\}}\\ & \times \mathit{\{}\mathrm{1}-\mathrm{0.7}(\mathrm{1}-m_{\mathrm{s}}/\mathrm{100})/\mathit{\{}(\mathrm{1}-m_{\mathrm{s}}/\mathrm{100})\\ & +\exp [-\mathrm{5.51}+\mathrm{22.9}(\mathrm{1}-m_{\mathrm{s}}/\mathrm{100}\left)\right]\mathit{\left\}}\mathit{\right\}},\\ & K=(-\mathrm{0.01383}+\mathrm{0.51575}{K}_{\mathrm{EPIC}})\times \mathrm{0.1317},\end{array}\end{array}$$

where K_EPIC is the soil erodibility output from the Erosion-Productivity Impact Calculator (EPIC) model (Williams, 1990), and K is the revised soil erodibility factor (t hm² h (hm² MJ mm)⁻¹) calibrated for Chinese soils (Zhang et al., 2008). The variables m_s, m_silt, m_c, and m_o denote the proportions (%) of sand, silt, clay, and organic carbon in the surface soil, respectively. Soil composition data were obtained from OpenLandMap.

Topographic factor (LS)

The LS factor quantifies the combined influence of slope length and slope steepness on soil erosion. In this study, we adopted the 30 m resolution national slope length and steepness factor dataset for China (Zhao et al., 2025), from which the portion covering the Loess Plateau was extracted for subsequent calculation.

Cover management factor (C)

The C factor represents the influence of vegetation cover and land management on reducing soil erosion. To capture its seasonal variability, the C factor for each temporal phase was derived using the quarterly LP-QLC10 land-cover dataset. Specifically, According to the national technical specifications (MEE China, 2021a), fixed C values of 0, 0, 0.01, and 0.7 were assigned to water, wetland, built-up land, and bare land, respectively. Meanwhile, for other vegetation types, C values are determined according to the Fractional Vegetation Cover (FVC), which was calculated from the Normalized Difference Vegetation Index (NDVI) using a linear scaling approach between the bare-soil and full-vegetation thresholds (Carlson and Ripley, 1997), as shown in Eq. (9). Specifically, the C factors of forest, shrubland, and grassland are assigned according to Table 1, while cropland is estimated using Eq. (10):

$$\begin{array}{}\text{(9)}& {\displaystyle }c={\displaystyle \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\mathrm{bare}}}{{\mathrm{NDVI}}_{\mathrm{veg}}-{\mathrm{NDVI}}_{\mathrm{bare}}}},\text{(10)}& {\displaystyle }{C}_{\mathrm{crop}}=\mathrm{0.221}-\mathrm{0.595}\log c,\end{array}$$

where c is the FVC value. NDVI_veg and NDVI_bare correspond to the NDVI values of full-vegetation and bare-soil endmembers in the study area, respectively. NDVI was calculated from satellite imagery following its standard formulation (Tucker, 1979).

Table 1 C factor values under different levels of Fractional Vegetation Cover (FVC).

Support practice factor (P)

The P factor represents the ratio of soil loss under erosion-control practices to that under natural slope conditions, thereby indicating the effectiveness of support measures. In this study, we adopted the P values reported in national-scale soil erosion products (Yan et al., 2025) for the Loess Plateau. Accordingly, the assigned P values for different land-cover types and slope ranges are summarized in Table 2. Furthermore, slope values were computed from the GLO-30 DEM using the Horn method (Horn, 1981).

Table 2 P factor values assigned to different land-cover types and slopes.

3.2.2 Historical soil erosion map (2000–2024) derivation

Quarterly soil erosion estimates for the Loess Plateau were generated for the period 2000–2024. As outlined in Fig. 3d, the erosion modeling followed a structured multi-factor workflow. First, the time-invariant factors, including soil erodibility (K) and topographic influence (LS), were derived to characterize the inherent susceptibility of the regional soils and terrain to erosion. Second, rainfall erosivity (R) was calculated by aggregating half-monthly erosivity values within each quarter, generating four seasonal indices that capture intra-annual rainfall variability while suppressing short-term fluctuations. Third, for each of the 100 quarterly land-cover phases from 2000 to 2024, the cover-management factor (C) and support-practice factor (P) were assigned to represent seasonal vegetation conditions and conservation interventions, respectively. Subsequently, all factor layers were integrated into the RUSLE formulation (Eq. 6) to derive spatially explicit soil erosion intensity for each quarterly interval across the 25-year period. Finally, quarterly estimates were aggregated to produce annual erosion maps, which were used for quantitative evaluation and cross-validation against officially reported government statistics. By combining quarterly and annual perspectives, the framework enables a high-resolution characterization of both seasonal fluctuations and long-term trajectories of erosion over the Loess Plateau.

4 Results

To address the three scientific questions outlined in the introduction – (1) what are the dominant land-cover change patterns over the past 25 years, (2) where and when soil erosion occurred, and (3) how erosion interacted dynamically with environmental factors – we begin this section by rigorously validating the land-cover and soil-erosion mapping products. Following the evaluation, we present (1) the spatiotemporal patterns of land-cover transitions, and (2) quantify the 25-year dynamics of soil erosion. In total, these results establish a solid foundation for subsequently (3) disentangling the mechanisms that shape soil-erosion processes and their linkages with environmental drivers.

4.1 25-year quarterly land-cover maps of China's Loess Plateau

In this study, we present a land-cover dataset for the Loess Plateau covering 2000–2024 with a three-month temporal resolution, referred to as LP-QLC10. Specifically, the dataset comprises 100 temporal phases and was produced at a spatial resolution of 30 m for 2000–2015 and 10 m for 2016–2024, reflecting the native characteristics of the underlying satellite sensors. Furthermore, to provide an overview of its temporal patterns and classification scheme, a subset of the quarterly maps is illustrated in Fig. 5. Spatially, forest is mainly distributed in the southern and eastern portions of the Loess Plateau, whereas grassland dominates the central and western regions. Cropland occurs largely along the regional margins and is typically situated near water bodies and built-up areas. In contrast, the northwestern Loess Plateau is characterized by extensive bare land interspersed with patches of grassland. Temporally, bare land has exhibited a pronounced decline over the past 25 years, while cropland and built-up areas have expanded, particularly in the southern and eastern regions. These trends reflect intensified land use and increasing human activities. In addition, the LP-QLC10 dataset developed under the CCTS framework significantly outperforms the baseline model in terms of both spatial accuracy and temporal continuity. Quantitatively, this advantage leads to an average overall accuracy (OA) improvement of 2.88 %. Qualitatively, it demonstrates enhanced multi-temporal mapping stability for land-cover types, such as grassland, cropland, and forest, with details provided in Appendix A1.

Figure 5 Demonstration of the LP-QLC10 product, showing representative temporal phases: (a) the first quarter of 2000, (b) the second quarter of 2010, (c) the third quarter of 2017, and (d) the fourth quarter of 2024.

To demonstrate the performance of the CCTS-derived LP-QLC10 products and its capability to capture temporal dynamics, we illustrate two representative regions for detailed analysis. As shown in Fig. 6, these regions are situated in Yulin and Yan’an cities, Shaanxi Province. Moreover, land-cover maps from four temporal phases (2000, 2007, 2016, and 2024) are presented to visualize long-term changes in land-cover composition. In the first region, Yulin City lies within a transitional zone between the sandy lands and desert margins of the northern Loess Plateau, where extensive ecological restoration programs have been implemented to combat desertification, such as aerial seeding, afforestation, and grassland rehabilitation. Correspondingly, the LP-QLC10 maps clearly depict these transformations. As illustrated in Fig. 6a, grassland areas have expanded markedly owing to vegetation restoration projects, while the growth of cropland reduced the extent of bare land, together indicating a consistent decline in barren surfaces. In the second region, Yan’an City has been a key demonstration area for China’s Grain for Green Program since its initiation in 1999. Consequently, following more than two decades of continuous implementation, the region has undergone profound land-cover transitions driven by farmland-to-forest conversion. As shown in Fig. 6b, cropland areas have progressively decreased, whereas large tracts of forest have expanded into zones formerly dominated by grassland. Collectively, these findings highlight the ability of the LP-QLC10 dataset to faithfully capture long-term, human-induced ecological transformations across the Loess Plateau.

Figure 6 Representative examples of the LP-QLC10 land-cover maps: (a) Yulin City and (b) Yan’an City, illustrating land-cover evolution across four temporal phases (2000, 2007, 2016, and 2024). Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google.

In addition, a qualitative comparison was conducted to compare the performance of the LP-QLC10 dataset with other widely used products. In particular, GLC_FCS30D (Zhang et al., 2024), CLCDs (Yang and Huang, 2021), and YRCC_LPLC (Wang et al., 2025) as the time-series, large-scale land-cover products, were selected for visual comparison. In detail, GLC_FCS30D is a global 30 m land-cover dataset containing 10 primary categories and 35 subcategories, covering 1985–2022 with 26 temporal layers updated every five years before 2000 and annually thereafter. Meanwhile, CLCDs provide national-scale annual land-cover data for China at 30 m resolution from 1990 to 2021, comprising nine categories. Moreover, YRCC_LPLC was developed to characterize land-cover changes in the Loess Plateau induced by ecological restoration, providing annual land-cover maps at 30 m resolution from 1990 to 2022. All the datasets include the eight major land-cover types present in the Loess Plateau, making them suitable for direct comparison.

Representative results are shown in Fig. 7 to provide a comprehensive comparison. First, Fig. 7a illustrates an urban area. The Loess Plateau contains numerous small cities where built-up surfaces are tightly interwoven with cropland, grassland, and forest, creating a highly fragmented landscape. In this setting, LP-QLC10 accurately identified built-up areas and effectively distinguished adjacent forest and grassland patches. In contrast, the other three products misclassified forest as cropland and grassland. Second, Fig. 7b depicts an agricultural landscape. In the Loess Plateau, cropland is typically located near water sources and interspersed with large patches of grassland. In this context, LP-QLC10 successfully distinguishes cropland, water bodies, and grassland with high spatial fidelity. By contrast, GLC_FCS30D misclassifies substantial grassland areas as cropland or bare land. Third, Fig. 7c presents a mountainous region. LP-QLC10 delineated cropland boundaries accurately, whereas GLC_FCS30D mislabeled cropland as grassland. CLCDs and YRCC_LPLC exhibited similar performance and incorrectly classified grassland as cropland. Moreover, since the other three datasets were produced based on pixel-based classification methods such as the Random Forest algorithm (Breiman, 2001), they tend to produce fragmented and spatially inconsistent cropland patterns. In contrast, the object-based CCTS framework underlying LP-QLC10 benefits from long-range contextual modeling enabled by self-attention mechanisms (Yuan et al., 2021), yielding spatially coherent land-cover results. Finally, Fig. 7d illustrates a wetland environment. GLC_FCS30D successfully identified wetland areas. However, errors were observed along rivers and their banks. Due to the fine-grained wetland classification system in GLC_FCS30D, confusion occurred among wetland subclasses. Specifically, cropland was misclassified as marshes, while water-covered surfaces were misidentified as flooded flats (sub-classes of wetland), particularly in sediment-rich river channels. In contrast, CLCDs and YRCC_LPLC failed to detect wetlands. Overall, LP-QLC10 demonstrates high classification reliability across diverse land-cover types and landscape settings, highlighting its superior spatial detail and temporal robustness compared to existing products.

Figure 7 Visual comparison of land-cover classification results over four representative regions: (a) urban areas, (b) agricultural zones, (c) mountainous regions, and (d) wetlands. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google.

4.2 Quantitative validation with multiple independent sources

To quantitatively evaluate the LP-QLC10 dataset, we conducted validation using (1) a human-annotated time-series reference dataset developed in this study; two publicly available validation point sets collected from (2) SinoLC-1 evaluation set (Li et al., 2023) – the first 1 m land-cover map of China, and (3) Geo-wiki – widely used global land-cover reference set (Fritz et al., 2017). For comparison, three widely used long-term land-cover products were also evaluated, including GLC_FCS30D, CLCDs, and Esri_GLC10 (Karra et al., 2021). Esri_GLC10 is a global 10 m dataset comprising 10 land-cover categories, available annually from 2017 to 2024. Meanwhile, three standard metrics were employed for evaluation, including Overall Accuracy (OA), Kappa coefficient, and Frequency Weighted Intersection over Union (FWIoU). OA and Kappa assess overall classification performance, whereas FWIoU gives more weight to infrequent classes, thereby offering a more balanced assessment of classification reliability.

Table 3 reports the results for 40 000 manual annotations developed in this study, obtained by labeling 8000 validation points at five-year intervals. Overall, LP-QLC10 consistently outperformed the comparison products across all three metrics. Specifically, its average OA (0.8144), Kappa (0.7258), and FWIoU (0.6454) were all substantially higher than those of CLCDs, GLC_FCS30D, and YRCC_LPLC, with LP-QLC10 ranking first in most evaluated years and showing particularly strong performance from 2001 to 2016. Although Esri_GLC10 surpassed LP-QLC10 in FWIoU for 2021, its applicability to long-term research is limited because the Esri product covers only the period 2017–2024. Besides, the transition of LP-QLC10 from 30 m imagery to 10 m imagery after 2016 resulted in a substantial improvement across all three metrics, highlighting the importance of higher spatial detail for accurately distinguishing fine-scale land-cover patterns in the Loess Plateau. Taken together, the long-term validation illustrates that LP-QLC10 provides both high classification accuracy and strong temporal stability.

Table 3 Quantitative results of the land-cover products on the time-series valuation points. The highest value in each column is marked in bold.

Table 4 Quantitative results of the land-cover products on the publicly available validation points. The highest value in each column is marked in bold.

To ensure a fairer and more objective evaluation, the comparision results of two publicly available validation datasets are reported in Table 4, where our LP-QLC10 product consistently outperformed the other products across both benchmarks. The SinoLC-1 validation set, derived from 1 m resolution Google Earth imagery, provided highly reliable reference labels. Therefore, despite the strong performance across all compared products, LP-QLC10 attained the highest OA, exceeding the second-best product by nearly 10 %. In contrast, the Geo-Wiki dataset yielded lower accuracies for all methods, because the quality of crowdsourced annotations may decline when contributors lack sufficient expertise (Moltchanova et al., 2022). Nevertheless, LP-QLC10 still achieved the highest scores, demonstrating its robustness across diverse validation scenarios.

The confusion matrix provide a detailed class-specific evaluation of model performance. Accordingly, Table 5 presents the confusion matrix of LP-QLC10 derived from the SinoLC-1 validation set. Based on this matrix, Overall Accuracy (OA) and Intersection over Union (IoU) were computed to assess overall performance, while User Accuracy (UA) and Producer Accuracy (PA) were derived to characterize commission and omission errors, respectively. In terms of PA, the forest class achieved the highest accuracy (0.944), followed by grassland, cropland, water, and bare land. In contrast, shrubland, bare land, and built-up areas exhibited relatively low accuracies. First, the low PA for shrubland is mainly attributable to its extremely limited spatial extent in the Loess Plateau, accounting for about 0.2 % of the total area (Zanaga et al., 2022), which results in severe sample imbalance. Meanwhile, the 10–30 m spatial resolution leads to substantial spectral confusion and reduces separability between shrubland and grassland. Second, the performance of built-up area is primarily constrained by image resolution. The imagery provide limited spatial detail to adequately characterize small buildings and narrow roads, thereby causing omission errors and blurred urban boundaries. Third, the low PA of bare land mainly arises from differences in class definitions between LP-QLC10 and SinoLC-1, as sparsely vegetated areas are classified as grassland in our dataset but are included in the “barren and sparse vegetation” category in SinoLC-1, as shown in Table A2, resulting in substantial confusion between bare land and grassland in the validation. When the mapping results were evaluated using the time-series validation set, which follows consistent class definitions, the PA for bare land reached 0.741 (Table A3). Moreover, the IoU results further supported these findings: vegetation types such as forest, grassland, and cropland exhibited high overlap scores, whereas shrubland and built-up areas were the most challenging categories for accurate delineation. Collectively, these superior and stable results demonstrate that the LP-QLC10 dataset provides the most accurate and temporally continuous land-cover maps of the Loess Plateau to date, thereby offering a robust foundation for subsequent analyses.

Table 5 Confusion matrix on SinoLC-1 validation set with per-class Producer Accuracy (PA), User Accuracy (UA), and Intersection over Union (IoU).

4.3 Spatiotemporal patterns of land-cover transitions

4.3.1 Temporal trends of land-cover transitions

The long-term land-cover transitions of the Loess Plateau reveal a pronounced shift toward increased vegetation cover and reduced surface degradation over the past 25 years. Overall, cropland and forest expanded by 22 % and 15 %, respectively, whereas bare land declined substantially by 51 %. These broad trends were largely shaped by vegetation restoration and land management interventions. In particular, 8.6 % (with estimates ranging 7.1 %–10.6 % when accounting for land-cover uncertainty) of the region was transformed from bare land to grassland, 3.5 % from grassland to forest, and 5.7 % from grassland to cropland. Collectively, these transitions illustrate a sustained regional greening process and a progressive replacement of degraded surfaces with more productive or ecologically stable land-cover types.

Detailed land-cover dynamics are illustrated in Fig. 8a and b. Over the past 25 years, the region has undergone systematic landscape restructuring, characterized by a coordinated expansion of forest, cropland, and built-up land alongside a substantial reduction in bare land. Forest, cropland, and urban areas all exhibited sustained upward trends, with cropland showing the largest net gain (28 095 km²), followed by built-up land (17 163 km²) and forest (13 131 km²). Conversely, bare land decreased by 65 029 km², indicating the substantial success of large-scale ecological restoration and land-management initiatives. Although grassland remained the dominant land-cover type, its area fluctuated dynamically in response to both ecological restoration and land-use conversion. Besides, water and wetland together account for less than 1 % of the Plateau’s area, underscoring the region’s pronounced scarcity of surface water resources. When aggregated across categories, vegetated surfaces (forest, shrubland, grassland, and cropland) collectively expanded by approximately 6 % of the Plateau’s total area, underscoring a sustained shift toward greener and more stabilized landscapes.

Figure 8 Temporal trends of land-cover composition from 2000 to 2024. (a) illustrates the proportional evolution of land-cover types, with the proportions in Q1 2000 and Q4 2024 annotated. (b) shows the change in coverage area of each class relative to the baseline year (2000), where points represent observed values and a fitted curve illustrates the long-term trend. (c)–(h) present region-specific variations across six subregions of the Loess Plateau.

Region-specific analyses in Fig. 8c–h further reveal a spatially heterogeneous yet policy-responsive pattern of land-cover transition across the Loess Plateau. In the northern and central restoration zones (IA, SL, and HG), bare land declined sharply, reflecting the sustained impact of desertification control and ecological rehabilitation initiatives. These reductions were accompanied by substantial gains in cropland and built-up land in IA, pronounced grassland expansion in SL, and concurrent increases in both forest and cropland in HG. Together, these shifts indicate large-scale conversion of previously degraded or sparsely vegetated surfaces into productive or revegetated landscapes. By contrast, TG exhibited a conversion of grassland to cropland, highlighting the pressure of agricultural intensification. In the southern subregions RV and RM, rapid urbanization emerged as the dominant driver of land-cover change, with built-up areas expanding considerably. This urban growth occurred predominantly at the expense of grassland and forest in RV, whereas in RM, most newly urbanized land originated from grassland, accompanied by localized forest recovery in recent years. Collectively, these regionally differentiated dynamics underscore the interplay between ecological restoration programs, agricultural demand, and urban development.

To strengthen ecological restoration, the Chinese government launched the second phase of the Three-North Shelterbelt Forest Program in 2001 and later implemented the National Soil and Water Conservation Plan (2015–2030), both of which designated the Loess Plateau as a priority region for vegetation recovery and erosion control. Accordingly, the land-cover transitions revealed by the Sankey diagrams (Fig. 9) clearly reflect the ecological impacts of these policies. From 2000 to 2014, 8.6 % (with estimates ranging 7.1 %–10.6 % when accounting for land-cover uncertainty) of the region shifted from bare land to grassland, consistent with the program’s early emphasis on large-scale grassland restoration prior to extensive afforestation. During the same period, forest and cropland expanded primarily through conversions from grassland, contributing 1.3 % and 5.7 % of the regional area, respectively. After 2014, 3.5 % of the area converted from grassland to forest, aligning with the intensified afforestation and slope-land conservation measures specified in the 2015 national plan. Meanwhile, dynamic exchanges persisted among grassland, cropland, and bare land, with grassland-to-cropland transitions (4 %) becoming dominant. Notably, sparsely vegetated bare land in semi-arid zones often exhibits spectral similarity to low-cover grassland, causing the model to temporarily classify such areas as grassland. Consequently, some reclamation activities may appear as grassland-to-cropland transitions. Overall, forest and cropland expanded, bare land declined continuously, and grassland remained relatively stable. These temporally coherent patterns indicate that sustained national ecological policies have substantially reshaped land-cover patterns and improved environmental conditions across the Loess Plateau.

Figure 9 Sankey diagram of land-cover transitions in the LP-QLC10 dataset. Data from Q3 2000, Q3 2014, and Q3 2024 were used to construct the diagrams.

4.3.2 Spatial distribution of land-cover change

In this section, we analyzed the spatial distribution of individual land-cover transitions and identified three distinct regional patterns across the Loess Plateau. Forest expansion occurred predominantly in the southeast, where relatively high precipitation supports woodland development. In contrast, the northwest exhibited extensive increases in grassland, closely corresponding to the widespread reduction in bare land. Meanwhile, cropland expansion was concentrated in the central Loess Plateau, especially within gully-dense landscapes. These representative patterns are mapped in Fig. 10, including: (a) Forest gain, (b) Grassland gain, (c) Cropland gain, and (d) Bare land loss.

Figure 10 Spatial distribution of typical land-cover changes from 2000 to 2024: (a) Forest gain, (b) Grassland gain, (c) Cropland gain, and (d) Bare land loss. Blue and yellow represent the timing of gains and losses, respectively, with darker tones indicating later occurrences of change. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google.

As shown in Fig. 10a, forest gain was mainly distributed in the southeastern Plateau. This spatial pattern is strongly controlled by climatic gradients, as forest ecosystems generally require more than 450 mm of annual rainfall to sustain stable growth (Wang et al., 2011). Given that rainfall gradually increases from the arid northwest toward the humid southeast, forest expansion is naturally concentrated in the better-watered southeastern region. In contrast, Fig. 10b illustrates that grassland expansion occurred most extensively in the northwestern Plateau, with additional patches emerging in the southeast. This trend aligns with the broad climatic tolerance and ecological resilience of grasslands, which enable them to establish under wide ranges of temperature and precipitation (Sala et al., 2000). Moreover, the newly established grassland areas in the northwest were located within the Maowusu Desert, the largest sandy region of the Plateau. This transformation reflects the effectiveness of national ecological restoration policies, with sustained interventions including sand-control forest belts and grass-grid stabilization under the Three-North Shelterbelt Forest Program contributing to the rehabilitation of the sandy surface. As for cropland, its expansion was primarily concentrated in the central part of the Loess Plateau, exhibiting a northeast–southwest distribution pattern. This spatial trend is closely associated with the regional strategy of filling gullies to create cropland, whereby gully bottoms and gentle slopes were converted into arable land to mitigate erosion and enhance agricultural productivity (Liu et al., 2013). Conversely, cropland that is unsuitable for cultivation can be restored to grassland or forest through the Grain for Green Program. Finally, the loss of bare land, as shown in Fig. 10d, is predominantly concentrated in the northwest, coinciding with the areas of grassland increase – most notably within the Maowusu Desert. This spatial correspondence confirms that the reduction of bare land is driven by the expansion of grassland in this region.

Overall, the LP-QLC10 dataset demonstrates strong potential for revealing long-term land-cover changes. Its spatiotemporal patterns reflect the combined effects of climatic gradients and ecological restoration policies, which together have driven substantial vegetation recovery across the Loess Plateau.

4.4 Soil erosion estimation results

Based on the LP-QLC10 dataset and the RUSLE model, as shown in Fig. 11, quarterly water-induced soil erosion maps for the Loess Plateau were generated for the period 2000–2024. Following the classification system defined in the Technical Specification for Investigation and Assessment of National Ecological Status (MEE China, 2021a), erosion intensity was categorized into six levels, ranging from slight to extreme. Among these, the slight level represents permissible conditions with negligible soil loss, whereas the light to extreme grades are considered eroded categories. The unit t (hm² q)⁻¹ denotes the amount of soil mass eroded per hectare of land per quarter. Annual soil erosion was obtained by summing the erosion values from all four quarters, expressed in t (hm² a)⁻¹, representing the soil mass eroded per hectare per year.

Figure 11 Illustration of the Loess Plateau soil erosion maps: (a)–(d) show annual maps for 2000, 2008, 2016, and 2024, respectively, whereas (e)–(h) depict the four quarterly maps of 2024. The soil erosion level legend is shown at the bottom.

The spatial distribution of erosion exhibits a clear association with both topography and land cover. As shown in Fig. 11a–d, erosion is predominantly concentrated in the central and western Loess Plateau, particularly within gully-dense regions such as HG and TG. In contrast, areas with dense vegetation cover, especially the extensive forest in the south, remain largely in a slight erosion state, demonstrating the strong protective effect of forest ecosystems. From a temporal perspective, moderate to severe erosion in the central region improved substantially during the early years (before 2008), whereas the severely eroded zones in the west continued to shrink in the following years. In addition, seasonal erosion patterns are strongly controlled by the monsoonal climate of the Loess Plateau. Precipitation is highly uneven throughout the year and is concentrated mainly in summer, while vegetation cover reaches its minimum in winter and peaks in summer. These joint variations regulate the intra-annual evolution of erosion intensity. As shown in Fig. 11e–h, erosion is essentially negligible in the first and fourth quarters, with only localized moderate erosion in the second quarter. In contrast, the third quarter experiences widespread moderate to strong erosion.

As illustrated in the above analysis, the pronounced seasonal variations in rainfall and vegetation cover across the Loess Plateau highlight the need for high-temporal-frequency observations to accurately characterize soil erosion processes. Accordingly, to investigate the influence of temporal resolution on erosion estimation, we designed three experimental groups that differ in both their temporal update frequency and the source of their input data:

• Annual-scale group (CSWED-based). This experiment adopted the Hydraulic Soil Erosion Dataset (CSWED) (Yan et al., 2025) as a representative annual-scale benchmark. CSWED provides 30 m soil erosion estimates for mainland China from 1990 to 2022, derived from annual rainfall, fractional vegetation cover (FVC), and land-cover inputs. Thus, similar to conventional annual erosion products, this group produces a single erosion map per year.

• Quarterly-scale group (LP-QLC10-based). This experiment utilized the proposed LP-QLC10 dataset in combination with quarterly rainfall and vegetation cover data. By updating all RUSLE factors every three months, this group achieves four times the temporal resolution of the annual approach, thereby capturing pronounced intra-annual variations in vegetation phenology, rainfall erosivity, and surface conditions.

• Mixed-resolution control group (CLCDs-based). To isolate the specific contribution of land-cover temporal frequency, this control experiment employed the same quarterly rainfall and vegetation cover inputs as the LP-QLC10-based group, while substituting quarterly land-cover updates with annual CLCDs products. Correspondingly, this configuration allows us to directly quantify how quarterly land-cover information enhances the accuracy of erosion estimates.

A quantitative comparison is presented in Table 6, with official statistics from the China Soil and Water Conservation Bulletin serving as the validation reference. The comparison between the CSWED and LP-QLC10-based results clearly illustrates the advantages of higher temporal resolution. Specifically, CSWED, which relies on annual inputs, substantially underestimated erosion intensity by more than 10 % relative to the reference data, primarily because it failed to represent the pronounced seasonal fluctuations in rainfall and vegetation that govern erosion processes. In contrast, the LP-QLC10-based estimates achieved a mean absolute error (MAE) of only 4.50 across all examined years. In addition, the CLCDs-based estimates exhibited larger deviations from the reference values than those produced using LP-QLC10. This discrepancy arises from the seasonal bias inherent in annual land-cover products, which are typically derived from imagery acquired during peak vegetation periods (e.g., spring or summer). Consequently, such imagery over-represents dense vegetation, leading to systematic underestimation of erosion intensity. Overall, these findings demonstrate that incorporating quarterly land-cover information substantially enhances both the temporal representativeness and accuracy of soil erosion estimation, particularly in regions with strong intra-annual variability such as the Loess Plateau.

Table 6 Quantitative comparison of the soil erosion maps derived from annual data and quarterly data. The reference data come from the China Soil and Water Conservation Bulletin, which has been published annually since 2018. The minimum error in each row is marked in bold.

4.5 Spatiotemporal dynamics of soil erosion

In this section, we examined soil erosion dynamics from both temporal and spatial perspectives and identified a clear long-term mitigation trend across the Loess Plateau. Approximately 37 % of the region exhibits a statistically evident decrease in erosion, corresponding to a 30 % reduction in mean erosion intensity and a 12 % shrinkage of erosion-affected area over the past 25 years. Moreover, seasonal contrasts remain pronounced, with the third quarter alone contributes 68 % of the annual erosion.

The 25-year erosion dynamics reveal a clear and persistent improvement in ecosystem stability. As shown in Fig. 12a, both mean annual erosion and the proportion of eroded land exhibited substantial declines, with mean erosion decreased by 29.91 % (from 13.34 to 9.35 t (hm² a)⁻¹), and eroded area shrinking from 34.73 % to 22.97%, representing an 11.76 % reduction. This sustained weakening of erosion intensity reflects the cumulative effects of vegetation restoration, land-management interventions, and landscape adjustments. Moreover, The seasonal pattern displayed in Fig. 12b underscores the dominance of summer rainfall in shaping erosion dynamics. Erosion was negligible in the first and fourth quarters, but increased sharply in the second and third quarters, with the latter peaked at 6.41 t (hm² q)⁻¹, accounting for 68 % of the annual total. Correspondingly, Fig. 12c shows that erosion declined most strongly in the third quarter, with a reduction of −3.53 t (hm² q)⁻¹. This long-term decrease is primarily driven by the continuous expansion of newly established vegetation, whose cover reaches its seasonal maximum during summer and thereby provides the strongest protection against intense rainfall events.

Figure 12 Illustration of the temporal trends in soil erosion. (a) presents the annual erosion trends from 2000 to 2024. (b) shows the intra-annual variation in erosion area and mean erosion in 2024. (c) details the quarterly change of the mean erosion from 2000 to 2024. (d) and (e) depict the overall and subregional changes in detailed erosion level across the Loess Plateau, respectively.

Besides, Fig. 12d demonstrate that soil erosion across the Loess Plateau has undergone substantial mitigation over the past 25 years, reflected by clear declines in both erosion extent and severity. Specifically, light erosion decreased from 133 434 km² (20.4 % of the Loess Plateau) to 93 821 km² (14.4 %), while moderate erosion declined from 59 362 km² (9.1 %) to 32 483 km² (5.0 %). Strong and above erosion dropped from 34 049 km² (5.3 %) to 23 736 km² (3.6 %). In contrast, slight erosion expanded significantly, indicating that increasingly large areas have transitioned into stable, permissible erosion conditions. Moreover, subregional patterns in Fig. 12e further reveal pronounced spatial heterogeneity in erosion mitigation. The loess hilly and gully region (HG) exhibits the most substantial improvements, with the largest reductions in both erosion area and mean erosion intensity. This trend is consistent with the substantial forest expansion, cropland reduction, and urban land conversion documented in Fig. 8e, all of which contributed to reduced bare land and enhanced soil and water conservation capacity. Meanwhile, the loess tableland and gully region (TG) and the sandy land and desert region (SL) show the next most notable declines. Although the remaining subregions also exhibit modest improvements, TG remains the most severely affected zone, with more than one-third of its surface still experiencing erosion and the highest mean erosion intensity on the Plateau. Overall, these patterns underscore the clear effectiveness of ecological restoration and land-cover restructuring in reducing soil erosion intensity, while also highlighting the persistent vulnerabilities in areas with complex topography and intensive human activities.

From a spatial distribution perspective, the Mann–Kendall (MK) trend test and Theil–Sen slope estimator provide a statistically robust method for assessing long-term erosion dynamics across the Loess Plateau. The MK test (Hamed and Ramachandra Rao, 1998) is a non-parametric method designed to detect the presence and significance of monotonic trends in time-series data. Complementarily, the Theil–Sen estimator (Sen, 1968) quantifies the magnitude of these trends by calculating the median slope across all pairwise observations. Together, these methods offer a reliable basis for diagnosing spatial variations in soil erosion dynamics. Accordingly, applying the MK test to annual soil erosion from 2000 to 2024 reveals a clear regional contrast in erosion dynamics. As shown in Fig. 13, areas with no significant trend appear in white, while increasing and decreasing trends are displayed in red and blue, respectively. Overall, 36.8 % of the Loess Plateau exhibit a decreasing trend in soil erosion, primarily concentrated in the central and southwestern Loess Plateau within the HG and TG subregions. Most decreasing areas show mild declines of −1 to 0 t (hm² a²)⁻¹, collectively covering over 2.16 × 10⁵ km². Furthermore, several zones in the southern Plateau display more substantial reductions (< −1 t (hm² a²)⁻¹). In contrast, only 7.9 % of the plateau exhibit an increasing trend, mainly distributed along the northern and southern margins, typically with slopes of 0–1 t (hm² a²)⁻¹. A few localized patches in the southwest exhibit stronger increases (> 1 t (hm² a²)⁻¹), particularly in the TG region, where the conversion of grassland to cropland has reduced vegetation protection. Overall, the spatial pattern of MK trends underscores the dominance of long-term improvements in soil erosion across the Loess Plateau. The widespread decreasing trends reflect the cumulative effectiveness of large-scale ecological restoration initiatives, while the localized increasing trends highlight remaining vulnerability zones where targeted soil conservation measures are still needed.

Figure 13 Spatial distribution of soil erosion trends across the Loess Plateau from 2000 to 2024. (a) Spatial patterns of the Theil–Sen slope and Mann–Kendall trend. Areas with increasing trends are shown in red, decreasing trends in blue, and areas without significant change are shown in white. (b) Statistical summary of the proportions of increasing, decreasing, and no trends.

5 Discussion

Based on the comprehensive validation and description of the land-cover and soil erosion mapping results presented in the Result section, we are able to address the last scientific question of this study: How did soil erosion interact with environmental factors? Subsequently, this section provides a detailed discussion of the mutual interactions among soil erosion, land-cover types, rainfall, and landform conditions, revealing the key processes that jointly shape the erosion dynamics across the Loess Plateau.

5.1 Dynamic interaction between soil erosion and environmental factors

In this section, we present a comprehensive assessment of how soil erosion interacts with major environmental factors. The analyses reveal several key mechanisms: the transformation from bare land to grassland produces the most erosion reduction in the Loess Plateau, while forests consistently exhibit the strongest soil conservation capacity. In contrast, cropland expansion tends to intensify erosion. Furthermore, pronounced quarterly variations in erosion are primarily regulated by the concentration of summer rainfall. Finally, erosion is preferentially amplified on sunny slopes, within mid-elevation belts, and on steeper terrain. Together, these findings illustrate the coupled influence of ecological, climatic, and topographic factors on shaping erosion dynamics across the Loess Plateau.

First, the interactions between vegetation dynamics and soil erosion exhibit a clear and quantitatively coherent pattern across the Loess Plateau. including forest, grassland, cropland, and bare land. As revealed by the land-cover transition pathways in Fig. 9, forest expansion originates primarily from grassland, grassland expansion from bare land, and cropland expansion mainly from grassland. This hierarchical conversion structure provides the foundation for interpreting the erosion responses associated with different land-cover gains and losses.

Figure 14a indicates that different transitions have distinct impacts on erosion intensity: grassland transitions to forest produced substantial soil-conserving effect (−14.86 t (hm² a)⁻¹). Similarly, stronger mitigation occurs when bare land is converted to grassland (−24.18 t (hm² a)⁻¹). In contrast, cropland expansion increased soil erosion by +15.44 t (hm² a)⁻¹, and bare land expansion results in the strongest erosion intensification (+41.75 t (hm² a)⁻¹). Collectively, these results indicate that erosion mitigation in the Loess Plateau is primarily driven by grassland and forest gains, whereas cropland and especially bare land expansion serve as the primary contributors to erosion escalation. Further evidence from the erosion in 2024 reinforces this pattern. As shown in Table 7, forest exhibits the highest soil-conservation capacity (< 1 t (hm² a)⁻¹), followed by grassland and cropland, whereas bare land experiences the most severe erosion. Besides, cropland exhibits pronounced seasonal sensitivity under the dominant single-season crops on the Loess Plateau. Dense summer crops enhances protection, while largely bare fields in winter and spring increase vulnerability to water-induced erosion. This seasonality is reflected in quarterly erosion patterns that most land-cover types show an approximately fourfold increase from Q2 to Q3, whereas cropland rises more moderately from 5.09 to 9.17 t (hm² q)⁻¹. Importantly, cropland reduction associated with the Grain for Green Program (GGP) has substantially contributed to erosion mitigation by converting retired cropland to forest or grassland, yielding an average erosion decrease of −36.64 t (hm² a)⁻¹ in the converted areas and highlighting the strong conservation benefits of this restoration pathway.

Figure 14 Soil erosion dynamics associated with the change of major vegetation types, based on comparative assessments from 2000 and 2024. (a) Erosion intensity change attributable to single-class gains and losses. (b) Total erosion reduction associated with the Three-North Shelter Forest Program (TSFP) and the Grain for Green Program (GGP).

Table 7 Erosion changes associated with land-cover gains and losses, derived from comparisons between 2000 and 2024, along with the current erosion levels for each land-cover category.

Furthermore, ecological restoration on the Loess Plateau has been primarily guided by two national programs: the Three North Shelter Forest Program (TSFP) and the GGP. Based on their objectives and implementation priorities, we attribute major land-cover transitions to these programs. Specifically, transitions from bare land to grassland, grassland to forest, and bare land to shrubland are mainly associated with the TSFP, whereas transitions from cropland to grassland and cropland to forest are mainly associated with the GGP. We note that these attributions may not capture all land-cover changes across the entire Loess Plateau, because additional environmental drivers and local human activities also contribute to the observed transitions. Nevertheless, this rule based attribution provides a practical approach for policy evaluation using the proposed LP-QLC10 dataset and offers a quantitative reference for assessing restoration effectiveness. As illustrated in Fig. 14b, the TSFP primarily promotes the conversion of bare land to grassland (4.74 × 10⁴ km²) and grassland to forest (2.43 × 10⁴ km²) across the Loess Plateau, with a smaller extent of bare land converting to shrubland. Among these transitions, conversion from bare land to grassland yields substantial erosion mitigation (−6.97 × 10⁷ t a⁻¹, ranging −7.66 to −6.50×10⁷ t a⁻¹ considering the land-cover uncertainty). In contrast, although the land cover changes associated with the GGP occur over a smaller total area, their erosion mitigation benefits are more pronounced. The cropland to grassland transition covers 2.51×10⁴ km² and reduces erosion by −11.54 ×10⁷ t a⁻¹. The cropland to forest transition covers only 0.55×10⁴ km², yet it decreases erosion by −2.37 ×10⁷ t a⁻¹. Overall, the estimated erosion reduction attributed to the GGP is approximately 1.5 times that attributed to the TSFP over 2000–2024.

Second, to assess the influence of soil erosion on land-cover changes, we mapped the cumulative soil erosion from 2000 to 2024 and analyzed the land-cover transitions within the top 10 % most eroded areas. As shown in Fig. 15a, the Loess Plateau experienced an average cumulative erosion of 265.02 t hm⁻². Meanwhile, the highest cumulative erosion concentrated in the southwest, particularly the loess tableland and gully regions. This spatial dominance arises primarily from the rugged terrain, where steep slopes and deeply incised gullies enhance runoff concentration and accelerate sediment transport. Furthermore, vegetation in these areas is largely composed of grassland and cropland, whose soil- and water-conservation capacities are markedly weaker than those of forest ecosystems, thereby further contributing to the elevated erosion intensity. In addition, Fig. 15b and c illustrate the detailed land-cover changes within the top 10 % erosion zones. Specifically, grassland experienced the dominate change, decreased by 7735 km², while cropland and bare land increased by 2768 and 5169 km², respectively. The Sankey diagram in Fig. 15c reveals that 10 354 km² of grassland degraded into bare land, and 5712 km² was converted into cropland, Although 4990 km² of bare land and 3644 km² of cropland were restored to grassland, the dominant process is still a net transformation from vegetated surfaces (especially grassland) toward more erosive land types. Collectively, these patterns indicate that grassland not only provide insufficient protection in areas experiencing severe erosion, but also undergoes progressive ecological degradation as soil loss intensifies. Meanwhile, this finding is supported by the visual example in Fig. 15a, where a slope initially covered by grassland in 2001 becomes dominated by bare surfaces by 2024. Consequently, intense erosion accelerates vegetation loss, which in turn further exacerbates erosion, forming a self-reinforcing degradation cycle. This finding underscores the urgent need for targeted ecological restoration strategies, such as promoting deep-rooted woody vegetation and implementing slope-stabilization practices, to interrupt this degradation cycle in extreme erosion zones.

Figure 15 Impacts of soil erosion on land-cover changes: (a) cumulative erosion from 2000 to 2024, (b) land-cover transitions in the top 10 % erosion areas, and (c) Sankey diagram of detailed transitions. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google, Basemap data: Esri 2026, USGS, World Hillshade | Powered by Esri.

Third, the quarterly analysis reveals that soil erosion in the Loess Plateau is fundamentally shaped by the pronounced seasonality of rainfall and its associated rainfall–runoff erosivity. As shown in Fig. 16, rainfall in 2024 was overwhelmingly concentrated in the second and third quarters, with the third quarter alone receiving 323 mm, more than 63 % of the annual total. This concentration of rainfall led to a substantial increase in rainfall-runoff erosivity. This concentration of rainfall greatly increased erosive energy, and the quarterly R factor reached 1548 MJ mm (hm² h q)⁻¹, which is about 4.2 times higher than in the second quarter. By contrast, the first and fourth quarters exhibited only minimal erosivity, with values remaining very low throughout both periods. Although vegetation cover is relatively high during summer and offers partial protection, it cannot fully offset the sharp rise in erosive power associated with intense seasonal rainfall. Consequently, soil erosion displays pronounced quarterly variability, with the second and particularly the third quarter dominating annual erosion, as illustrated in Fig. 12b. These findings highlight the essential role of seasonal rainfall patterns in governing erosion dynamics and emphasize the importance of conducting erosion assessments at a quarterly temporal resolution.

Figure 16 Quarterly rainfall in 2024 and the corresponding quarterly averaged R factor for the Loess Plateau.

Finally, terrain related controls on soil erosion intensity were assessed by examining how erosion varies with slope aspect, elevation, and slope gradient. In particular, Fig. 17a presents the distribution of erosion in 2024 across different slope aspects, expressed as the proportion of eroded area within each aspect relative to the total area of the Loess Plateau. Erosion displays a pronounced north–south contrast, with southern sunny slopes experiencing substantially higher erosion intensity than the shaded slopes to the northwest. This pattern reflects the microclimatic advantages of shady slopes, where reduced solar radiation and lower evaporative demand support denser forest and grassland cover, thereby enhancing surface protection and suppressing runoff driven erosion. Furthermore, elevation gradients modulate erosion intensity. After normalizing erosion area by land area within each elevation band, Fig. 17b shows that approximately 68 % of total erosion is concentrated between 800 and 2400 m, where grassland and cropland dominate and human disturbance is relatively high. In contrast, although high-elevation regions above 4000 m contribute only a small fraction of total eroded area, they exhibit a disproportionately large share of strong and above erosion due to steep slopes, sparse alpine vegetation, and highly efficient runoff pathways. Additionally, slope gradient imposes an strong control on erosion severity. As illustrated in Fig. 17c, gentle slopes are dominated by light erosion, whereas the proportion of strong erosion increases sharply with slope steepness. Notably, slopes steeper than 45 % are primarily covered by grassland rather than forest, and the limited canopy density, shallow rooting depth, and weak soil reinforcement capacity of grass vegetation are insufficient to withstand intense rainfall episodes, leading to accelerated soil loss. Collectively, these findings demonstrate that these terrain based patterns play a pivotal role in shaping the spatial heterogeneity of erosion intensity on the Loess Plateau, and underscore the need for topographically targeted soil and water conservation strategies in areas where steep slopes and sparse vegetation coincide.

Figure 17 Erosion distribution across landforms, including (a) erosion patterns across slope aspects, (b) erosion patterns along elevation gradients, and (c) erosion patterns across slope classes.

5.2 Potential conservation of the Loess Plateau

In this section, we discussed the potential conservation benefits achievable through a rational redistribution of land-cover types across the Loess Plateau. Under an optimized vegetation configuration, the potential minimum erosion is estimated to be 6.42 t (hm² a)⁻¹), representing a 31 % reduction compared to the erosion in 2024.

The potential capacity of a landscape to conserve soil and water defines the theoretical upper limit of erosion reduction that can be achieved through ecological management (Singh et al., 2025). Building on this concept, we delineated environmentally suitable zones for different restoration measures by integrating precipitation regimes, topographic conditions, and patterns of human land use. These suitability zones represent areas where ecological interventions can operate, such as afforestation, grassland restoration, or cropland retention. Thus, given the relative erosion-mitigation capacities of land-cover types quantified in Sect. 5.1, we adopt the restoration hierarchy of Forest–Grassland–Cropland as the guiding principle for land-cover optimization on the Loess Plateau. By reallocating land-cover types according to this hierarchy and their environmental suitability, the region can theoretically attain its optimal soil conservation state, thereby defining the lowest achievable erosion level under realistic ecological constraints.

• Forest-suitable area. Forest development requires adequate rainfall and moderate temperature conditions. Studies have shown that forestry-favorable zones on the Loess Plateau typically located in areas with mean annual precipitation exceeding 450 mm (Liu et al., 2021). Therefore, existing forest regions, together with grassland and bare land areas receiving more than 450 mm of precipitation and situated below 4000 m elevation (Peili et al., 2020), are identified as suitable for forest expansion.

• Cropland-suitable area. Although cropland provides limited soil and water conservation benefits, it remains indispensable for sustaining the livelihoods of more than 100 million residents of the Loess Plateau. To ensure agricultural security while avoiding ecologically vulnerable terrain, existing cropland on slopes < 15 % is classified as highly suitable, while those between 15 %–25 % are considered conditionally suitable when the silt proportion is below 45 % (Rossiter, 1996). In addition, bare land areas with precipitation > 350 mm, elevation < 3000 m, and slopes < 15 % are also considered suitable for cropland use (Zeng et al., 2018). Areas failing to meet these suitability thresholds are reassigned to forest or grassland categories to improve ecological stability.

• Grassland-suitable area. Grassland exhibits broad ecological tolerance and can thrive under relatively low temperature and moisture conditions. Previous research indicates that grassland generally requires annual precipitation above 250 mm to maintain stable growth (Kang et al., 2007). Consequently, existing grassland, as well as cropland and bare land that meet this precipitation threshold but are unsuitable for forest or cropland development, are classified as grassland suitable areas.

Figure 18 Optimized land-cover distributions based on the current land cover: (a) current land cover in Q4 2024, (b) optimized land-cover map, and (c) statistics of land-cover transitions.

On the basis of the land cover patterns in Q4 2024, an optimized land cover map was generated using the defined environmental suitability conditions. Compared with the current land cover patterns in Fig. 18a, the optimized map in Fig. 18b indicates substantial landscape reconfiguration. A detailed statistical summary is provided in Fig. 18c. The transitions can also be grouped by the two national programs: Bare to forest and bare to grassland transitions are consistent with the primary objectives of the TSFP, whereas cropland to forest and cropland to grassland transitions align with the GGP. In addition, bare to cropland serves as a compensatory adjustment to offset cropland losses in the optimized scenario. Among these transitions, bare to grassland under the TSFP accounts for the largest converted area (2.4 × 10⁴ km²), concentrated in the central Plateau, including the SL and TG regions. Meanwhile, under the GGP, cropland to forest represents a particularly challenging transition: it not only involves 7 × 10³ km², but also typically requires sustained management over years to decades to establish stable forest cover, underscoring the need for long-term planning and continued implementation.

Figure 19 Potential minimum soil erosion in the Loess Plateau. (a) Current soil erosion in 2024, (b) estimated minimum soil erosion under optimized land-cover conditions, and (c) statistical distribution of erosion intensity levels.

Furthermore, applying the RUSLE model to the optimized land cover scenario demonstrates the ecological benefits of these transitions. As shown in Fig. 19, compared with the current erosion pattern in 2024, the theoretical minimum erosion map exhibits a pronounced reduction in strong and severe erosion, especially in the central and southwestern Plateau where many grassland and cropland areas are replaced by forest. Moreover, in the northeastern Plateau, extensive bare land restoration to grassland alleviates many zones previously classified as weak or moderate erosion. Besides, the erosion level distributions from 2000, 2024, and the optimized scenario reveal a clear long term improvement in Fig. 19c, with the minimum erosion state pushing most light and moderate erosion areas toward slight erosion. In general, the mean soil erosion of the Loess Plateau in 2024 is 9.35 t (hm² a)⁻¹, whereas the optimized vegetation configuration lowers this value to 6.42 t (hm² a)⁻¹, representing a 31 % potential reduction. Achieving this minimum state requires substantial mitigation in high-intensity erosion zones, including a 43 % reduction in strong erosion and a 61 % reduction in severe erosion, while extreme erosion must be completely eliminated. These findings illustrate the magnitude of the improvement achievable through ecologically guided land cover management.

Overall, the results demonstrate that soil and water conservation initiatives in the Loess Plateau have achieved substantial progress over the past 25 years, leading to a marked reduction in erosion intensity and a significant expansion of vegetation cover. However, despite these notable improvements, further ecological restoration remains challenging. High-intensity erosion zones, particularly those classified as strong, severe, or extreme, continue to impede sustainable land management. Looking ahead, the SL and TG regions, which contain the largest concentrations of such erosion hotspots, should be prioritized for future remediation. Targeted ecological interventions in these high-risk areas will be essential for consolidating previous achievements and advancing the region toward long-term ecological stability.

6 Data availability

The LP-QLC10 land cover product and the corresponding soil erosion maps generated in this study are available at https://doi.org/10.57760/sciencedb.33656 (Cheng et al., 2026).

7 Conclusions

Long-time-series monitoring of land-cover and soil erosion dynamics in the Loess Plateau remains a major challenge due to the lack of high-frequency temporal observation data. By filling three knowledge gaps – (1) What are the long-term land-cover change patterns? (2) Where and when did the erosion occur? (3) How did the erosion dynamically interact with environmental factors? – this study conducted the first quarterly land-cover and soil erosion mapping for the Loess Plateau (LP-QLC10), covering the period from 2000 to 2024 in 100 time steps. Specifically, a semi-supervised learning framework based on consistency constraints (CCTS) was employed to learn semantically invariant representations from publicly available single-temporal labels, and then time-series land-cover maps were generated using Landsat and Sentinel multi-spectral imagery. Moreover, accuracy assessment was conducted in several independent sources, demonstrating accurate and spatiotemporal-consistent results. In general, the proposed LP-QLC10 dataset effectively captured the spatial–temporal dynamics of land-cover transitions across the Loess Plateau.

Building upon LP-QLC10, we further constructed a quarterly hydraulic soil erosion analysis using the RUSLE model. Owing to its high-frequency temporal attribute, the modeled erosion estimates exhibited strong agreement with the China Soil and Water Conservation Bulletin (2018–2024), yielding a mean absolute error of 4.50 %. The results reveal a pronounced decline in soil erosion intensity over the past 25 years: the mean erosion intensity decreased by 30 %, from 13.34 to 9.35 t (hm² a)⁻¹, and the areas experiencing light, moderate, and strong or above erosion declined by 30 %, 45 %, and 30 %, respectively. Meanwhile, erosion hotspots remained concentrated in the central and southwestern Loess Plateau, and the observed spatial and temporal variations were strongly linked to land-cover distribution, terrain characteristics, and precipitation patterns. Furthermore, scenario-based modeling indicated a potential minimum erosion intensity of 6.42 t (hm² a)⁻¹ under optimized and rational management practices, particularly through grassland expansion and cropland-to-forest conversion.

Overall, the LP-QLC10 dataset provides a fine-grained, large-scale, and long-term representation of land-cover evolution across the Loess Plateau, offering valuable support for regional environmental modeling and dynamic analysis. The associated soil erosion mapping and analysis provide a comprehensive characterization of erosion mitigation and its driving factors, serving as a crucial reference for sustainable land management, ecological restoration, and policy decision-making in fragile ecosystems.

Appendix A

A1 Evaluation of the long-term stability of CCTS framework

To evaluate the improvement of the long-term stability of the proposed CCTS framework, we conducted both quantitative and qualitative evaluations to assess the improvement of the CCTS in terms of accuracy and temporal stability for long-term land-cover mapping on the Loess Plateau. Since CCTS uses HRFormer (Yuan et al., 2021) as its backbone segmentation model, we trained HRFormer on 2021 satellite imagery using the ESA WorldCover v200 product (Zanaga et al., 2022) as labels, and then used the trained model to generate baseline land-cover maps from 2000 to 2024. The quantitative results, evaluated using the multi-temporal reference set constructed in this study, are presented in Table A1. The LP-QLC10 dataset generated by the CCTS framework consistently outperformed the baseline model across all evaluated years. Notably, the baseline achieved performance comparable to LP-QLC10 in 2021 because the training and evaluation data came from the same year and therefore had highly consistent feature distributions. However, this advantage did not generalize to other years. When applied to earlier periods, particularly those dominated by Landsat imagery, HRFormer was more sensitive to temporal differences in imaging conditions and sensor characteristics, resulting in lower performance in the 2001, 2006, and 2011 validations. In contrast, CCTS improved both accuracy and temporal stability by learning more consistent cross-year representations, achieving approximately 4 % higher OA than the baseline in these years.

The qualitative comparison further reinforced this performance gap. As illustrated in Fig. A1, variations in imaging conditions and sensor characteristics substantially compromised the stability of multi-temporal predictions generated by the baseline land-cover model. Figure A1a presents an area that experienced only limited land-cover change. Although the actual semantic changes were minor, the HRFormer results showed pronounced interannual fluctuations, particularly in vegetation-related categories such as grassland and cropland. In contrast, by incorporating consistency constraints during training, the CCTS framework produced temporally stable and spatially coherent predictions throughout the entire multi-year inference period. Furthermore, Fig. A1b shows an area characterized by the gradual expansion of built-up area and forest. For categories such as forest, which generally require many years to establish, and built-up land, which typically changes incrementally,the baseline model produced markedly inconsistent predictions across years in some areas due to multi-temporal spectral variability, clearly indicating misclassification. This instability not only reduced spatial prediction accuracy but also introduced cumulative errors into area statistics, thereby limiting the ability to reliably characterize long-term land-cover change trajectories. Such limitations are especially problematic for analyses spanning multiple decades. By contrast, the CCTS framework enabled more temporally consistent inference and more faithfully captured the expansion of urban land and the increase in forest area from 2011 to 2015.

Figure A1 Qualitative comparison with/without CCTS. HRFormer serves as the baseline model without using CCTS, while LP-QLC10 denotes the dataset generated using CCTS. To facilitate the identification of actual land-cover categories in remote sensing imagery, all visualizations are based on composite images from the third quarter. Imagery © 2000–2024 USGS and ESA, Map data © 2026 Google.

Table A1 Quantitative comparison with baseline model. The highest value in each column is marked in bold.

A2 Additional materials

Table A2 Land-cover types in LP-QLC10: value, definition, and correspondence with existing products. YRCC_LPLC and CLCDs share the same classification system.

Figure A2 The obtained image count from multiple platforms within each temporal period. The statistic denotes the number of quality-screened satellite image tiles in Google Earth Engine that intersect the Loess Plateau. A Landsat scene covers approximately 3.3×10⁴ km², whereas a Sentinel tile covers approximately 1.0×10⁴ km².

Table A3 Confusion matrix of the 2021 mapping results, evaluated using the time-series validation set, with per-class Producer Accuracy (PA) and User Accuracy (UA).