Remapping carbon storage change in retired farmlands on the Loess Plateau in China from 2000–2021 in high spatiotemporal resolution

The soil organic carbon pool is a crucial component of carbon storage in terrestrial ecosystems, playing a key role in regulating the carbon cycle and mitigating atmospheric CO₂ concentration increases. To combat soil degradation and enhance soil organic carbon on the Loess Plateau, the Grain-for-Green Program (GFGP) has been implemented. Accurately quantifying change in soil organic carbon stock (ΔSOC) resulting from farmland retirement is essential for informing land use management. In this study, the spatial and temporal distribution of retired farmlands on the Loess Plateau was analyzed using Landsat imagery from 1999–2021. To assess the effects of the years since retirement, climate, soil properties, elevation, and other factors on ΔSOC, climate-zone-specific multivariable linear regression models were developed based on field-sampled soil data. These models were then used to map ΔSOC across the retired farmlands. Results indicated that a total of 39 065 km² of farmland was retired over the past two decades, with 45.61 % converted to grasslands, 29.75 % to shrublands, and 24.64 % to forestlands. The years since retirement showed a significant positive correlation with ΔSOC, and distinct models were developed for different climatic zones to achieve high-resolution (30 m) ΔSOC mapping. The total ΔSOC from retired farmland on the Loess Plateau was estimated at 21.77 Tg in carbon equivalent, with grasslands contributing 81.10 %, followed by forestlands (11.16 %) and shrublands (7.74 %). The associated datasets are freely available at 30 m resolution from 2000 to 2021 (Yang, 2025, https://doi.org/10.6084/m9.figshare.28785971).

1 Introduction

Soil organic carbon (SOC), as the largest terrestrial ecosystem carbon pool, plays a crucial role in regulating climate change (Mir et al., 2023). Global SOC was estimated at approximately 1400–1500 Pg C, about four times the organic carbon pool of terrestrial plants (Scharlemann et al., 2014). The high SOC is essential to support multiple ecological benefits, such as purifying water, increasing crop yields and maintaining primary productivity (Paustian et al., 2019). Currently, $\mathrm{1}/\mathrm{3}$ soil in the world is degraded, causing many socioeconomic (e.g., unemployment, poverty, immigration) and environmental (e.g., desertification, ecosystem degradation, biodiversity loss) issues (Ferreira et al., 2022; Ouyang et al., 2016). The large area of degraded soil also released more than 50 Pg carbon annually into the atmosphere which conflicts with the decarbonization target for mitigating global warming (Prăvălie et al., 2021). Therefore, restoring degraded soil is urgently needed for sustainable development and environment security.

Ecological restoration by nature alone is a lengthy process. Under the urgent need for restoring degraded soils and mitigating climate change, scientific management measures are necessary to accelerate the ecosystem restoration process (Lengefeld et al., 2020; Pape, 2022; Wang et al., 2021a). Many large-scale ecological restoration strategies around the world have showed encouraging ecological benefits. Brazil's Atlantic Forest Restoration Pact (AFRP) was established in 2009, and Argentina and Paraguay joined the impressive project in 2018, forming the Atlantic Forest Restoration Tri-national Network (Calmon et al., 2011). Hundreds of organizations have been actively involved in this decade-long efforts to protect and restore the forests, which recovered about 7000 km² forests and enhanced regional biodiversity (De Oliveira Faria and Magrini, 2016). Forests established by restoration in this project between 2010 and 2015 would have sequestered 1.75 Pg carbon if they were not re-cut (Piffer et al., 2022). The Development Project “Green Great Wall” in Africa was launched by the African Union in 2007, aiming at restoring savannahs, grasslands and farmlands across Africa to help biodiversity cope with climate change and desertification. The goals of the project are to restore 1 million km² forests in 2030 and sequester 250 Tg C (Graham, 2022; Macia et al., 2023). China has started ecological restoration practices and researches since the 1970s, and has implemented six national key ecological restoration projects (Cui et al., 2021). Among the projects, the Grain-for-Green Program (GFGP) is one of the most ambitious projects in the world with the highest investment and the largest implemented area (Xu et al., 2022). From 1999–2019, the GFGP implemented in 25 provinces and exceeded 0.343 million km² land area with 49 Tg sequestered carbon, indicating a significant potential of SOC stock by ecological restoration (Lu et al., 2018). Based on Deng et al.'s (2017) study, the total carbon stock in the GFGP implemented area was 682 Tg C in 2010, and projected to 1697 Tg C in 2020.

One of the primary area of the GFGP is the Loess Plateau, because the long-term indiscriminate cultivation and logging on the Loess Plateau has caused over 40 % of the total area (about 270 000 km²) in severe soil erosion and a significant loss of SOC (Shao et al., 2022). As the implementation of the GFGP, 96.1 Tg C was sequestered from 2000–2008 on the Loess Plateau (Feng et al., 2013; Xiao, 2014). Nonetheless, current estimations of SOC stock still have large uncertainties due to the technology and data limits (Zhang et al., 2022). On the Loess Plateau, the accumulation of SOC can be affected by many untested factors, such as ecosystem types and years since retirement. Moreover, most of the studies fail to differentiate the carbon sequestration between retired and currently cultivated farmlands, and caused an overestimate of SOC stock. Therefore, a more reliable estimation should be reached to quantify the benefit in SOC stock of the retired farmlands with the consideration of those issues (Deng et al., 2017; Sun et al., 2016).

While previous studies have confirmed the overall increasing trend of SOC changes following farmland retirement, significant uncertainties persist due to limited spatial resolution and insufficient temporal coverage. Existing datasets fail to provide the continuous spatiotemporal dynamics of retired farmland distribution on the Loess Plateau (Xu et al., 2018; Yang and Huang, 2021; Bai et al., 2024). Furthermore, existing SOC assessments (Li et al., 2020; Yi et al., 2023) lack the capacity to quantify fine-scale differences in SOC stock between retired and cultivated farmlands (ΔSOC). They also fail to capture the year-by-year dynamics of retired farmlands and SOC accumulation in high resolution by considering the heterogeneity of the Loess Plateau. To address these gaps, this study aims to: (1) reconstruct annual farmland retirement patterns (2000–2021) using multi-source remote sensing data; (2) develop a high-resolution ΔSOC model integrating terrain, climate and vegetation covariates based on the difference in SOC stock between retired and adjacent cultivated farmlands; and (3) generate 30 m resolution ΔSOC maps to quantify the impact of GFGP on carbon sequestration. Our spatially explicit approach provides unprecedented insights for optimizing ecological restoration strategies in heterogeneous landscapes.

2 Materials and Methods

2.1 Study Area

The Loess Plateau ($\mathrm{100}\mathit{°}{\mathrm{52}}^{\prime }-\mathrm{114}\mathit{°}{\mathrm{33}}^{\prime }$ E, $\mathrm{33}\mathit{°}{\mathrm{41}}^{\prime }-\mathrm{41}\mathit{°}{\mathrm{16}}^{\prime }$ N) is located in the north central part of China (Fig. 1a), in the middle reaches of the Yellow River, with a sensitive and fragile ecological environment, belonging to the warm temperate continental monsoon climate, characterized by dry and cold in spring and winter, warm and hot in summer and autumn. The average annual temperature is 3.6–14.3 °C. The average annual precipitation is 400–600 mm, of which is concentrated between July and September, and decreases from east to west and south to north (Zhou et al., 2016). The annual evaporation is 1400–2100 mm, with a trend of low in the south and east, high in the north and west. The elevation is 800–3000 m, and the original surface vegetation mostly is grassland, shrubland, deciduous broadleaf forest, and mixed broadleaf-conifer forest (Zhou et al., 2016). The total area of the Loess Plateau is 635 000 km², including Shanxi, Ningxia, Shaanxi, Gansu, Qinghai, Inner Mongolia, Henan provinces. The main terrain is hilly and gully, with soft loessial soil texture.

Figure 1 The map of the study area, (a) location, (b) soil sampling sites and climatic zones.

2.2 Identifying Retired Farmlands

To identify and confirm the spatial range of the annual retired farmlands on the Loess Plateau, Landsat remote sensing images (30 m resolution) from 1999–2021 were downloaded from the United States Geological Survey (USGS, https://EarthExplorer.usgs.gov, last access: 20 August 2022). The images with less cloud (lower than 10 %) in growing season (from May to September) were selected for further analysis. Those images were processed by the standard steps recommended by ArcGIS Pro 2.8 (Environmental Systems Research Institute, Inc., ESRI), including preprocessing, image classification and validation. To improve image readability, remote sensing images were first preprocessed in ENVI 5.3, including radiometric calibration, FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) atmospheric correction, gram-schmidt pan sharpening, seamless mosaic and subset data from ROIs (regions of interest). The image classification was then performed in ArcGIS Pro 2.8. In this study, we used the support vector machine (SVM) supervised classification method to classify the land cover types into the following seven categories: farmland, forestland, grassland, shrubland, water body, building land, and bare land. Training samples were selected through visual interpretation of high-resolution imageries and systematically managed using a training sample manager. A total of 23 100 ROI samples were used for model training, with an additional 6930 independent ROIs reserved for validation. During the accuracy assessment phase, the classification performance over the study period consistently achieved kappa coefficients ranging from 0.76–0.90 and overall accuracy values between 0.80 and 0.91. The average accuracies for different land cover types were as follows: farmland (0.71), forestland (0.87), grassland (0.86), shrubland (0.92), water body (0.97), building land (0.92), and bare land (0.87).

2.3 Field Sampling and SOC Measurements

To determine the ΔSOC in ecosystems established on retired farmlands, we implemented a systematic sampling design based on spatial proximity principles. Initial sample sites were systematically generated at 5 km intervals across the retired farmland distribution map (Fig. 1b), forming a comprehensive grid framework. For each retired farmland point, we identified the nearest long-term cultivated farmland counterpart to create a spatially paired sampling site. The sampling strategy incorporated stratification across different ecosystems, climatic zones, and years since retirement. To minimize human interference, we pre-screened all potential sites using ultra-high resolution imagery (0.5 m) to exclude areas near roads, villages, or irrigation ditches. Additional considerations included accessibility and sampling feasibility, leading to the exclusion of 133 site pairs from initial design to field implementation. Finally, 2430 soil samples from 135 sample sites were collected from fields. Nine soil samples (three 10 cm layers from top 30 cm soil in 3 sample points) were collected for every sample site, and nine soil samples from the nearest farmlands were also collected similarly. Each soil sample was individually bagged, labeled, and stored in cold storage for lab measurement. After drying and grinding through a sieve at 0.25 mm, SOC of each soil sample was measured by potassium dichromate external heating method. The difference in total SOC stock of the top 30 cm soil layer between retired farmlands and the nearest cultivated farmlands was defined as ΔSOC that contributed by the GFGP.

2.4 Model Development and ΔSOC Mapping on the Loess Plateau

ΔSOC is influenced by both natural environmental conditions and human activities, leading to variations across different climatic conditions of the Loess Plateau. Therefore, we developed different models based on the relationships between ΔSOC and variables such as years since retirement, geographic location, elevation, soil bulk density (BD), and 19 bioclimatic factors. Years since retirement were obtained from the annual spatial distribution data in retired farmlands on the Loess Plateau (Sect. 2.2). The data sources for climate information can be found in Sect. 2.5. The 19 bioclimatic factors were derived by following the formula in WorldClim (https://worldclim.org/data/index.html, last access: 12 January 2022). For every grid cell of retired farmlands, the bioclimatic factors were calculated as the average of the years since retirement. All the variables were extracted to the sample sites by the Kriging interpolation and prepared for model development.

Based on the factors introduced above, we combined ANOVA, single-factor regression, all subset regression and stepwise regression to select variables for multivariate linear models of ΔSOC. The steps included: data preprocessing, univariate analysis, multivariate analysis, model evaluation, and diagnostic checks. Finally, several key variables that co-occurred were selected. In consideration of wide climatic range on the Loess Plateau and different possible response of ΔSOC to the factors among climatic conditions (Zhang et al., 2018), we divided the Loess Plateau into different climatic zones for different ecosystem types (e.g., forestland, shrubland, grassland) based on climate regionalization in China–Climatic zones and climatic regions (GB/T 17297-1998) and climate data (Sect. 2.5). As Fig. 1b shows, we obtained middle temperature zone (MT, <8 °C) and warm temperate zone (WT, >8 °C) by the annual average temperature, and semi-arid zone (SA, <400 mm) and sub-humid zone (SH, >400 mm) by annual precipitation. In addition, five combined climatic zones were obtained: MT-SA (same as MT), WT-SA, WT-SH (same as SH), WT, and SA. A multivariate linear regression model was developed specifically for each ecosystem types in each climatic zone. Before regression analysis, diagnosis of multicollinearity is conducted, and the threshold is generally set at 10 to detecting correlations between the independent variables and identify those independent variables that were incorrectly included in the same regression model. The regression models were evaluated and validated by residual analysis, significance level (p-value), coefficient of determination (R²), root mean square error (RMSE) and mean absolute error (MAE), and the robustness of the models were validated by leave one out-cross-validation. Statistical power analysis indicates that the current stratified sampling design provides adequate power for detecting medium to large effects, though sensitivity for detecting small effects remains limited. Model robustness under this design is rated as “acceptable”.

Based on the results of model evaluation, the final selected models were used to estimate the overall ΔSOC of retired farmlands on the Loess Plateau. With the final selected multivariate linear regression models, the ΔSOC in the top 30 cm soil layer were mapped by raster calculation in different climatic zones and ecosystem types at 30 m resolution. And the total ΔSOC on the Loess Plateau contributed by the GFGP was obtained by summing up the ΔSOC in all the retired farmlands without recultivation within the study period.

2.5 Data Sources

The air temperature and precipitation data to calculate the 19 bioclimatic factors were from the China Meteorological Data Service Center (CMDC, http://www.geodata.cn, last access: 13 July 2023). Elevation data of every grid cell were from the Digital Elevation Model database (https://e4ftl01.cr.usgs.gov/MEASURES/, last access: 14 July 2023). Soil properties were retrieved from Harmonized World Soil Database (HWSD, https://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/, last access: 14 July 2023), and the boundary of the Loess Plateau was downloaded from the Resource and Environment Science and Data Center (https://www.resdc.cn/, last access: 15 January 2022). All the raster data were resampled to 30 m resolution.

3 Results

3.1 Distribution of Retired Farmlands

From 1999–2021, the final retired farmlands without recultivation on the Loess Plateau was 39 065 km² (Fig. 2a). The final retired farmlands were less than the area by summing up annually retired farmlands because of frequent recultivation (Fig. 3a). The annual area of retired farmlands has been fluctuating throughout the study period with no significant trend (Figs. 2a–f and 3b). The least amount of retired farmlands occurred in 2002 (28 003 km²; 4.41 % of the whole studied area), and the most was 78 653 km² in 2016 (12.39 % of the whole studied area). The retired farmlands were converted to different vegetation types, including forestlands, shrublands and grasslands. The ratios of different vegetation types in every year were in the ranges of 10.65 %–38.60 %, 14.63 %–47.70 % and 17.02 %–64.98 % for forestlands, shrublands and grasslands, respectively (Fig. 3b). Within the studied period in average, most of the retired farmlands were converted to grasslands (45.61 %) and shrublands (29.75 %).

Figure 2 Spatial distribution of annually retired farmlands on the Loess Plateau in (a) 2000, (b) 2005, (c) 2010, (d) 2015, (e) 2021, and (f) cumulative retired farmlands from 1999–2021.

Figure 3 (a) Cumulative retired farmlands and recultivated farmlands and (b) Annual area of different vegetation types from retired farmlands from 2000–2021.

The annual retired farmlands were unevenly distributed among different climatic zones (Fig. 2a–f). The annual retired farmlands in the other years can be found in Fig. S1a–p in the Supplement. For the final retired farmlands, the area in MT-SA, WT-SA and WT-SH were 20 299, 10 572 and 8194 km², respectively. In the MT-SA zone, the dominant ecosystem type from retired farmlands was grasslands which had 9705 km² (47.81 %), and followed by shrublands (5887 km², 29.00 %) and forestlands (4707 km², 23.19 %). In the WT-SA zone, grasslands were also the dominant ecosystem type which accounted for 4925 km² (46.59 %), and forestlands accounted the least (2384 km², 22.55 %). In the WT-SH zone, the percentages of different ecosystem types were 30.96 %, 30.16 % and 38.88 % for forestlands, shrublands and grasslands, respectively.

Among different years (Figs. 2a–f and S1), the highest areas for each ecosystem type were forestlands in the WT-SH zone in 2016 (12 846 km²), shrublands in the MT-SA zone in 2001 (15 441 km²), and grasslands in the MT-SA zone in 2007 (26 171 km²). The lowest areas were found in 2019 for forestlands in the WT-SA zone (813 km²), in 2013 for shrublands in the WT-SH zone (271 km²), and in 2013 for grasslands in the WT-SH zone (806 km²).

Among provinces, the retired farmlands in different years had significant differences (Table S1 in the Supplement), where Shanxi Province had the most in 2016 (30 912 km²) and Qinghai Province had the least in 2017 (438 km²). The final retired farmlands from 1999–2021 was the most in Inner Mongolia Province (8626 km²) and the least in Henan Province (739 km²). More forestlands could be found in warmer and wetter regions. The largest forestlands (15 073 km²) were found in Shanxi Province in 2016, while the least were found in Qinghai Province in 2016 (34 km²).

3.2 Analysis of Soil Samples

The results of soil samples showed that the SOC stock were 2.19–62.70 g C kg⁻¹ in retired farmlands, and 2.25–63.83 g C kg⁻¹ in adjacent cultivated farmlands. The average SOC were the highest in forestlands (4.84–62.70 g C kg⁻¹), followed by shrublands (2.62–54.72 g C kg⁻¹) and grasslands (2.19–21.83 g C kg⁻¹). The average ΔSOC of the all sample points was 2.86 g C kg⁻¹, with a standard error of 1.17 g C kg⁻¹, and a 95 % confidence interval of [0.56, 5.15] g C kg⁻¹. The findings indicated that the farmland retirement had significantly increased the SOC stock. To facilitate the ΔSOC estimation by area, we converted the SOC stock to area based content by soil bulk density. The highest value of ΔSOC after retirement was from forestlands in the SH zone (26.52 kg C m⁻²) and the lowest value was from sample in grasslands in the WT zone (0.91 kg C m⁻²). Forestlands and shrublands had significantly increased the SOC stock by 48.53 % and 20.34 %, respectively (p<0.05, Fig. 4a). Among different climatic zones (Fig. 4b), forestlands in the SA zone had the biggest increase (58.80 %), and followed by forestlands in the SH zone (44.53 %) and shrublands in the MT-SA zone (26.74 %).

Figure 4 SOC stocks in farmlands and retired farmlands (kg C m⁻²), (a) Comparison of SOC stocks on the Loess Plateau in farmlands retired to different ecosystem types (forestland, shrubland, grassland) with those in adjacent cultivated farmlands, and letters a and b are labeled to indicate significant differences in the ANOVA. (b) Comparison of different climatic zones are emphasized.
Note: 1-forestlands in the SH zone, 2-forestlands in the SA zone, 3-shrublands in the WT-SH zone, 4-shrublands in the WT-SA zone, 5-shrublands in MT-SA the zone, 6-grasslands in the WT zone, and 7-grasslands in the MT zone; Letters a, b and ab are labeled to indicate significant differences in the ANOVA, for same ecosystem in Fig. 1a and for same climatic zone combination Fig. 1b.

The ΔSOC of different ecosystem types in different climatic zones had significant relationship to the years since retirement (Fig. 5). The ΔSOC was negative in the first few years and significantly increased as the years since retirement increases, except forestlands in the SA zone and shrublands in the MT-SA zone. Most of the relationships indicated constant increase in ΔSOC except ΔSOC in grasslands in the MT zone which had a saturation point after 15 yrs of retirement.

Figure 5 Relationship between years since retirement and ΔSOC, (a) forestlands in the SH zone, (b) forestlands in the SA zone, (c) shrublands in the WT-SH zone, (d) shrublands in the WT-SA zone, (e) shrublands in the MT-SA zone, (f) grasslands in the WT zone, (g) grasslands in the MT zone.

3.3 Models of ΔSOC

Samples for different ecosystem types were divided by different combinations of climatic zones to find the final selected models by Backward Stepwise Regression. All variance inflation factor (VIF) diagnostic results were below the threshold of 10, including years since retirement, latitude, longitude, elevation, soil bulk density, and bioclimatic variables BIO1 to BIO19. The final selected models of ΔSOC in different ecosystem types were shown in Table 1 based on the results of evaluation and validation. In this table, t is the years since retirement, lat is latitude, ele is elevation, BD is soil bulk density, and BIO1–BIO19 are 19 bioclimatic factors, n is sample sizes at each level.

The analysis showed that seven regression equations were the final acceptable representative for the ΔSOC on the Loess Plateau when the study area was divided into SH and SA zones for forestlands, WT-SH, WT-SA and MT-SA zones for shrublands, and WT and MT zones for grasslands. The coefficients of determination (R²) ranged from 0.476–0.830 with p<0.05, The models with the highest R² were obtained for grasslands (0.830 in the WT zone and 0.790 in the MT zone), and the model with the lowest R² was for shrublands in the MT zone (0.476).

Table 1 Models of the ΔSOC in retired farmlands on the Loess Plateau.

3.4 Mapping ΔSOC

According to the regression models for ΔSOC and the distribution of retired farmlands, the ΔSOC in the retired farmlands on the Loess Plateau was quantified throughout the GFGP implementation period, excluding recultivated farmlands (Fig. 6a). The total benefit in ΔSOC on the Loess Plateau till 2021 was 21.77 Tg C with a range between −26.52 and 31.91 kg C m⁻² at 30 m raster level. Significant variations in ΔSOC were observed across different ecosystem types (Fig. 6b and Table 2). To provide detailed and vegetation-specific insights, Table 2 presents ΔSOC values for three climatic zone combinations associated with each vegetation type. Grasslands contributed the most ΔSOC increment (17.657 Tg C). Among the different climatic zones for grasslands, MT zone contributed the most (78.04 %, $-\mathrm{0.48}-\mathrm{3.04}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$), followed by WT zone (21.96 %, $-\mathrm{8.20}-\mathrm{31.91}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$). Forestlands contributed the second largest ΔSOC (2.429 Tg C) with 151.96 % from SH zone ($-\mathrm{26.52}-\mathrm{22.86}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$), and −51.96 % from SA zone ($-\mathrm{2.96}-\mathrm{8.67}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$). The shrublands only contributed 7.74 % of the total benefit in ΔSOC (1.685 Tg C) with 78.04 % from MT-SA zone ($-\mathrm{26.49}-\mathrm{30.57}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$), 45.07 % from WT-SA zone ($-\mathrm{4.00}-\mathrm{3.28}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$) and −23.11 % from WT-SH zone ($-\mathrm{4.60}-\mathrm{26.10}\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$).

Figure 6 Spatial distribution of the ΔSOC, (a) the distribution in the whole study area, and (b) raster level frequency of ΔSOC.

Table 2 The ΔSOC (positive and negative portion) in retired farmlands in different ecosystem types in different climatic zones (Tg C).

The potential ΔSOC by different provinces also changed significantly, but the potential ΔSOC in different ecosystem types by the same provinces were evenly changed (Table S3 in the Supplement). ΔSOC increased more in Shanxi and Shaanxi provinces, followed by Henan, Gansu, Inner Mongolia and Ningxia, and less in Qinghai province.

4 Discussion

4.1 Distribution of Retired Farmlands

In consideration of the topographic complexity and vegetational variation on the retired farmlands, a large-scale retrieve of retired farmland information from remote sensing images is challenging (Wei et al., 2021). For instance, farmlands and grasslands have similar spectrum characteristics in spring and summer seasons and can be easily confounded (Estel et al., 2015), which lead to inaccuracy in remote sensing image classification. The inaccuracy can be minimized by comparing with multi-source high-resolution remote sensing images (Yan et al., 2023). In this study, although different vegetation types were involved on the retired farmlands (e.g., forestland, shrubland and grassland), the accuracy in identifying retired farmlands could high to 90 % by combining visual interpretation of Landsat dataset, field observation, globeland30 database, and ultra-high resolution images from Google Earth.

Farmland retirement is the main land use change driver on the Loess Plateau. As classified in this study, retired farmlands on the Loess Plateau from 2000–2021 are unevenly distributed across different climatic zones, because of the significant hilly and gully terrain in the study area (Huang et al., 2007; Wen et al., 2015). We focused on forestlands, shrublands and grasslands from retired farmlands, and noticed that most forestlands were distributed in the SH zone due to higher precipitation than the SA zone. Grasslands were more distributed in the MT zone than in the WT zone, due to the temperature in the MT zone being more favorable for grasses than in the WT zone, and people may be more engaged in pastoral activities in the WT zone. Shrublands were more distributed in the MT-SA zone than in the WT-SH zone because the WT-SH zone is more suited to forest growth, thus having high percentage of tree cover and relatively low distribution of shrub. In this study, grasslands accounted for a large proportion in retired farmlands on the Loess Plateau, but the increase in forestlands were more significant.

The spatial–temporal patterns of farmland retirement varied significantly across years, primarily driven by policy orientation and farmers' participation willingness. During the study period, the Chinese central government implemented two phases of GFGP: the first from 1999–2013, and the second from 2014 onward. High rates of retirement were observed at the beginning of every phase due to promising subsides. High retirement rates were observed at the launch of each phase, largely due to attractive subsidy schemes. However, participation willingness declined afterward, as falling grain prices reduced the relative value of subsidies, leading some farmers to recultivate retired land (Xie et al., 2023). Additionally, population growth between 2000 and 2020 escalated local food demand, further motivating recultivation. Some abandoned farmland-induced misclassification also could introduce bios into the spatial analysis of retired farmlands. These dynamics are consistent with the findings of Wang et al. (2013), who reported a rapid decline in farmland area from 1999–2003 during the first GFGP phase, followed by a rebound due to recultivation and subsequent stabilization.

4.2 Model development for ΔSOC

Land use change due to GFGP can strongly affect SOC, and SOC tend to be lower in farmlands (Deng et al., 2014), which was proved in this study by comparing retired and adjacent cultivated farmlands. The increase in ΔSOC in retired farmlands shows a strong relationship with the years since retirement, although a slight decrease in SOC may occur during the early stages of land use change (Deng et al., 2017). During the study period, all vegetation types exhibited a consistent increasing trend in SOC after the initial few years. However, the accumulation tends to approach an upper limit as the ecosystem matures and stabilizes, as observed in grasslands that follow a logarithmic growth pattern. Some retired farmlands with decreasing SOC were found, which could be explained by interchange of recultivation and retirement (Qiu et al., 2018), but the deeper mechanism is still need to be explored. Moreover, the high SOC in adjacent farmlands due to good agricultural practice could also offset the benefit of ΔSOC from the GFGP (negative ΔSOC was mostly found in farmland with high SOC).

Based on the statistical analysis (Fig. 4), the range of the ΔSOC in grasslands was significantly smaller than that in forestlands and shrublands. This indicates the accumulation rate of SOC in grasslands was lower than that in forestlands and shrublands due to the low primary productive and the fine quality of grass litter for decomposition (Lukina et al., 2020), whereas woody litter contains more lignin and decomposes slowly (Xiao et al., 2022). Therefore, different models were developed according to vegetation types and climatic zones. Based on the models, the climatic factors had significant effect on ΔSOC besides the years since retirement. Among the climatic factors, the models showed that ΔSOC were more sensitive to precipitation-based bioclimatic factors (e.g., BIO12–BIO19). This is because most of the Loess Plateau is located in semi-arid and arid area with limited precipitation (Zhang et al., 2015). Moreover, increased precipitation and temperatures can enhance the decomposition of surface litter (Sharma and Sharma, 2022), and in turn reduce ΔSOC.

4.3 Benefits in ΔSOC on the Loess Plateau

Under climate change, ecological restoration is an urgent need to improve the healthiness of degraded ecosystems (Liu et al., 2023; Yang et al., 2023). As a major benefit from ecological restoration, ΔSOC brings a lot of interests due to SOC is the major carbon pool in the ecosystems. To illustrate ΔSOC from ecological restoration, only a comparison of restored and adjacent unrestored ecosystems should be persuasive (Francaviglia et al., 2019). Numbers of studies focusing on ΔSOC in retired farmlands has been conducted on the Loess Plateau, and found an increasing ΔSOC as a result of GFGP (Wang et al., 2021b), and the national SOC sequestration caused by retirement was estimated to be 14.46 Tg yr⁻¹ (Zhao et al., 2013). But they failed to make comparison with the adjacent farmlands. In this study, we analyzed the ΔSOC of retired farmlands and adjacent cultivated farmlands, and confirmed that the GFGP can provide significant amount of ΔSOC on the Loess Plateau, although negative ΔSOC was found in some areas.

Recently, studies have shown that SOC stocks in the GFGP region on the Loess Plateau increased by 20.18 Tg C between 1982 and 2017 (Li et al., 2022). The total ΔSOC (21.77 Tg C) of retired farmlands on the Loess Plateau estimated in this study was slightly higher than that value, which proved that the results of this study are reliable. The mechanisms driving ΔSOC vary across vegetation restoration types and climatic zones. While warmer and more humid regions generally exhibit higher carbon sequestration rates – owing to enhanced photosynthesis and plant growth under favorable temperature and precipitation regimes – these conditions also accelerate SOC turnover, potentially limiting long-term storage benefits compared to arid and semi-arid regions (Sierra et al., 2017). Therefore, selecting appropriate vegetation types is critical to prevent slow SOC accumulation and early saturation. Moreover, sustainable management practices – such as controlled grazing and systematic harvesting – are essential to maintain ecosystem health and maximize long-term soil carbon storage, thereby strengthening the role of retired farmlands in climate change mitigation.

4.4 Limitations and Uncertainties

Remote sensing images are widely used in studies of land use change because of their accuracy and timeliness. In this study, the use of Landsat dataset has practical feasibility to provide reliable distribution of retired farmlands. However, the Loess Plateau has a large spatial area, and has a fragmented and complex topography, which increases the difficulty of land use classification. Therefore, the 30 m resolution images can result in misclassification (e.g., abandoned farmlands vs. retired farmlands), although we obtained acceptable accuracy (80 %–91 %). Recently, the availability of ultra-high resolution images (sub-meter resolution) allows a more accurate classification, but lacks of long period records.

In this study, the direct comparison of retired farmlands and adjacent cultivated farmlands reflected a more persuasive ΔSOC. The multivariate linear regression models that developed for estimating ΔSOC can effectively reduce estimation errors by accounting for the spatial heterogeneity of the Loess Plateau. Increasing the number of sample points would further enhance model flexibility, allowing the incorporation of additional factors – such as slope, elevation, and soil properties – to stratify the study area into more representative subzones. Furthermore, establishing permanent observation points to monitor both retired and adjacent cultivated farmlands would provide reliable pairwise comparisons essential for robust model calibration. To more accurately project the future soil carbon sequestration potential of retired farmlands, the integration of process-based ecosystem models could be a more reliable approach, such as DLEM (Dynamic Land Ecosystem Model, (Tian et al., 2003)), LPJ–GUESS (Lund Potsdam Jena General Ecosystem Simulator, (Smith et al., 2001)), and CENTURY (Parton et al., 1987).

5 Data availability

The associated datasets “The 30 m resolution distribution of retired farmlands and their carbon sequestration on the Loess Plateau in China from 2000 to 2021” are available at Figshare (https://doi.org/10.6084/m9.figshare.28785971, Yang, 2025), including distribution of retired farmlands from 2000–2021, years since retirement, and high resolution ΔSOC from the retired farmlands.

6 Conclusions

Farmland retirement is an effective strategy to restore degraded ecosystems and increase carbon storage on the Loess Plateau. In this study, we found the total area of retired farmlands on the Loess Plateau during the study period was 39 065 km². The dominant ecosystem type was grasslands, followed by shrublands and forestlands. The area of retired farmlands showed significant interannual changes without a specific trend, and the retired farmlands varied in different climatic zones. Area of retired farmlands in the MT-SA zone were significantly higher than WT-SA zone and WT-SH zone. Based on soil samples, we found that ΔSOC increased with the years since retirement, and developed seven regression models for ΔSOC by years since retirement, temperature, precipitation, soil bulk density, latitude and longitude, and ecosystem types. According to the models, the total benefits in ΔSOC from retired farmlands on the Loess Plateau were estimated to be 21.77 Tg C, with the variation ranged from −26.52 to 31.91 kg C m⁻² at grid cell level. The most ΔSOC were contributed by retired farmlands in the MT-SA zone (15.120 Tg C), followed by WT-SA zone (6.440 Tg C) and WT-SH zone (0.210 Tg C). Therefore, Long-term implementation of GFGP brought significant impacts on increasing soil carbon sinks on the Loess Plateau, which contributed significantly in mitigating climate changes and promoting sustainability in the studied area.