We sincerely appreciate your insightful feedback on our manuscript. Please find our detailed, point-by-point responses to your comments below.

Major Comments

1. Standardize terminology to avoid ambiguity. (1) The manuscript defines post-retirement SOC stock change as carbon capture and storage (CCS), which is easily confused with engineered carbon capture and storage. Please use this term cautiously and, at first mention, clearly constrain it to the study-specific definition of ΔSOC; alternatively, consider a less ambiguous term such as carbon sequestration (or SOC stock change, ΔSOC). (2) The phrase "length of farmland retirement" for the time elapsed since retirement should be refined; consider using the more common "years since retirement" and apply it consistently.

Response: We very appreciate the review’s suggestion. ΔSOC is a very suitable substitution for CCS. “Years since retirement” is more clearly represent the meaning of our intention. So we made the change through the manuscript.

2. Introduction—focus on the case for necessity. The introduction adequately motivates the importance of soil restoration and summarizes ecological restoration programs worldwide, but the discussion of uncertainty in estimating impacts of farmland retirement on SOC and the specific necessity of this study remains limited. Please sharpen the study aims by clarifying, relative to existing datasets/studies, where the present dataset and methods add value and which key problems the study is designed to resolve.

Response: While the increasing trend of soil organic carbon (SOC) following farmland retirement has been established in previous studies, critical gaps remain in quantitatively capturing its spatiotemporal dynamics. Existing datasets on retired farmland distribution (Xu et al., 2018; Bai et al., 2024; Yang and Huang, 2021) are limited to a few discrete years, failing to provide continuous annual patterns. Similarly, while SOC measurements after retirement exist (Li et al., 2020; Yi et al., 2023), no generalized model has been developed to account for spatial heterogeneity, nor have fine-resolution SOC accumulation maps been produced. These limitations result in three key research gaps: (1) no fine-scale heterogeneity of the Loess Plateau, (2) no pattern of difference in SOC storage between retired and cultivated farmlands and (3) no year-by-year dynamics of SOC accumulation. Therefore, this study intends to obtain a year-by-year dynamics of retired farmland, construct models evaluating the yearly change in SOC, and generate high resolution SOC stock change by considering the heterogeneity of the Loess Plateau with assistant of high resolution remote sensing data. Here are the detailed changes as follows for objectives by sharpen the purpose and advance in datasets:

“hile 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.”

3. Section 2.2 (Identifying Retired Farmlands): classification/validation details. Retired farmland is mapped via SVM supervised classification, but classification and validation details are insufficient, which undermines the credibility of subsequent CCS (ΔSOC) estimation. Please add the classification scheme, sizes of training/validation sets, per-class accuracies, and the full confusion matrix.

Response：Land use classification across the Loess Plateau was performed using ENVI remote sensing software and ArcGIS Pro, supported by high-resolution imagery such as Google Earth. Seven land cover categories were defined according to existing classification systems, local environmental characteristics, and research objectives: farmland, forestland, grassland, shrubland, water body, building land, and bare land. The analysis spanned the period from 1999 to 2021. For each land use classification scheme, 1,050 regions of interest (ROIs) were used for training, and 315 independent ROIs (30% of the training set) were used for validation. Classification accuracy, measured in overall accuracy or kappa coefficient ranges across years, was as follows: farmland (0.63–0.86), forestland (0.75–0.97), grassland (0.79–0.94), shrubland (0.83–0.97), water body (0.89–0.97), building land (0.84–0.98), and bare land (0.76–0.94). So we added the following explanation in section 2.2:

“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 6,930 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 to 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).”

4. Section 2.3 (Field Sampling and SOC Measurements): sampling and screening rules. The text states that initial sampling points were laid out at 5 km intervals and that unsuitable points were removed. Please clarify (1) the design of paired sites (how pairs were planned and located) and (2) the criteria for removing sites and the number removed.

Response：Paired sampling points were designed following the principle of spatial proximity to ensure homogeneous coverage across the study area, encompassing varied ecosystems, climatic zones, and years since farmland retirement. A systematic grid sampling approach was applied to the farmland retirement distribution map to generate retired farmland points. For each retired farmland point, the nearest long-term cultivated farmland site was identified to form spatially paired sampling units. Potential sampling locations were pre-screened using high-resolution imagery to minimize human disturbance, excluding areas near roads, villages, or ditches. During field surveys, land use status was verified, and inaccessible or unsuitable points were excluded. A total of 133 sample pairs were removed throughout the process from design to implementation. We added the information in section 2:

“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. 1-b), 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.”

5. Sample size and power by stratum. The design is “3 profiles per site × 3 depths per profile × 3 samples per depth,” yielding 135 paired profiles (cropland and adjacent retired land as a pair). Because modeling is stratified by multiple factors (ecosystem type × climate zone), please report the effective sample size (n) in each stratum, the train/validation split, and a brief power assessment to support model robustness.

Response：Based on temperature, precipitation, and ecosystem characteristics, retired farmlands on the Loess Plateau were systematically stratified into seven conversion types: forestland in the SH and SA zones; shrubland in the WT-SH, WT-SA, and MT zones; and grassland in the WT and MT zones. The corresponding numbers of valid soil sample pairs for each stratum were 32 (SH for forestland), 21 (SA for forestland), 22 (WT-SH for shrubland), 13 (WT-SA for shrubland), 20 (MT for shrubland), 16 (WT for grassland), and 11 (MT for grassland), respectively. We added the sample size for each stratum in Table 1. The training/validation was conducted by the leave one out-cross validation. To support model robustness, we added the following changes in the main text (section 2.4):

“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’.”

6. Section 2.4 (Model Development and CCS Mapping). (1) Results indicate stratified mapping primarily via multiple linear regression (MLR). Please state explicitly which model(s) are used for final mapping and why they were chosen. (2) Specify the temporal scale of variables such as temperature and precipitation (multi-year means, growing-season metrics, or WorldClim bioclimatic variables). (3) List the final candidate predictors after correlation/collinearity screening and confirm that the VIF threshold is 10. (4) If terms like “optimal/best model” are used, provide the selection basis (multi-model/parameter comparisons and cross-validation metrics), or replace with “final selected model(s)” for rigor and transparency.

Response：The study employed a layered multiple linear regression model as the final mapping approach, selected based on its scientific rationality, statistical appropriateness, interpretability, and operational feasibility. This method aligns with current theoretical understanding of soil organic carbon dynamics, corresponds to the sample size characteristics across different layers, supports ecological mechanism analysis and policy application, and is readily implementable and validateble. The final models used for mapping were presented in Table 1 with detailed formula. Those model were chosen due to their good performance after validation and evaluation.

“The final selected models of ΔSOC in different ecosystem types were shown in Table 1 based on the results of evaluation and validation.”

All variables were aligned with the spatiotemporal context of the GFGP period. Temporal resolution was defined as mean interannual variation, and spatial resolution was uniformly resampled to 30-meter grids. Following correlation and collinearity screening, the final candidate predictors included: years since retirement, latitude, longitude, elevation, soil bulk density, and bioclimatic variables BIO1 to BIO19. All the bioclimatic variables for each grid cell of retired farmlands were the average of all the years since retirement. All selected predictors satisfied a variance inflation factor (VIF) threshold of 10. We added the following specification in section 3.3.2:

“For every grid cell of retired farmlands, the bioclimatic factors were calculated as the average of the years since retirement.”

“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.”

In line with the reviewer’s suggestion, the term “final selected model” has been used throughout the manuscript instead of “optimal/best model” to enhance rigor and transparency.

7. Section 3.1, Fig. 2—readability and area accounting. Increase font sizes and color contrast; consider trimming the number of years shown. In addition, please present side-by-side series (or a stacked chart) for annual retirement, reclamation, and final retained retirement areas to clarify that “cumulative annual retirement area > final retirement area” is mainly due to reclamation, avoiding misinterpretation as classification error.

Response：In response to reviewer comments, we have streamlined the visual presentation by relocating selected figures is Figure 2 to the supplementary information. The main text now retains only figures a, f, h, o, u, and v. Additionally, we have enhanced Figure 3 by incorporating a new visualization that displays annual land conversion from farming alongside recultivated land area. This addition presents cumulative yearly totals relative to previous years, providing a clear temporal perspective on land use dynamics throughout the study period. These modifications improve both the clarity and analytical value of our graphical presentation.

8. Section 3.2—SOC change statistics without uncertainty. In addition to reporting means, please provide confidence intervals or standard errors and briefly describe the uncertainty estimation method.

Response：For analysis the sample points, uncertainty estimation is calculated based on the variability of ΔSOC sample data, specifically through standard error calculation, confidence interval calculation, and hypothesis testing. The standard error is the sample standard deviation divided by the square root of the sample size. Given the large sample size (n=135), a normal distribution approximation is used, with the 95% confidence interval calculated as the mean±1.96×standard error. This method assumes strong representativeness and independent identically distributed samples, disregarding factors such as model error or spatial autocorrelation. Mean: 2.86 g C kg⁻¹, Standard Error: 1.17 g C kg⁻¹, 95% Confidence Interval: [0.56, 5.15] g C kg⁻¹.

“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.”

9. Figure 4—units/terminology, outliers, n by stratum, and consistent zoning. (1) The caption reads “SOC stocks,” but the unit is g C kg⁻¹ (concentration). Please standardize terminology and units: if the focus is on stocks, convert using 0–30 cm thickness and bulk density to report kg C m⁻². (2) There are apparent outliers; please describe outlier detection and handling (data cleaning and/or robust estimation) to limit their influence on modeling and conclusions. (3) Report sample sizes (n) for each climate × ecosystem stratum to support assessment of model reliability. (4) Zoning standards are not fully consistent across sub-panels (forest: SH/SA; shrub: WT-SH/WT-SA/MT-SA; grassland: WT/MT). If a three-class climate grouping is adopted, consider unifying to WT-SH, WT-SA, and MT-SA.

Response：1) We changed the unit of “Soil organic carbon stocks,” to kgC/m². 2) To maintain complete sample integrity, all data points were retained and utilized in the analysis. Robust statistical methods—including the use of medians and interquartile ranges for descriptive statistics, as well as robust regression techniques—were applied to minimize the influence of outliers on statistical inference and modeling. 3) The sample sizes (n) for each climate × ecosystem stratum are as follows: 32, 21, 22, 13, 20, 16, and 11, respectively (Table 1). Since each sampling site includes paired measurements from both cultivated and retired farmlands, the sample sizes for these two categories are identical. These sample sizes support the reliability of the model and ensure statistical adequacy within each stratified group. 4) The model is designed to assess changes in soil organic carbon during the Grain for Green Program, which are strongly influenced by local climatic conditions. The zoning criteria were intentionally differentiated—forestland (SH/SA), shrubland (WT-SH/WT-SA/MT-SA), and grassland (WT/MT)—to enhance within-group homogeneity and between-group heterogeneity, while maintaining sufficient sample sizes for robust statistical analysis. This stratified approach improves the scientific rigor, reliability, and practical applicability of the model outputs. In contrast, a unified three-category climate grouping (e.g., WT-SH, WT-SA, MT-SA) would not adequately support the learning of clear and stable patterns within groups.

“To facilitate the ΔSOC estimation by area, we converted the SOC stock to area based content by soil bulk density.”

10. Figure 5—consistent zoning for response curves. The climate-zone groupings are not fully aligned across sub-figures. Please indicate whether this is due to sample size or model-fit constraints, and make the zoning standard as consistent as possible across figures and text.

Response：The apparent inconsistencies in climatic zone delineation reflect deliberate adaptations to model-fit requirements, balancing ecological accuracy with statistical feasibility. As described in the previous response, these adjustments were made to develop a more reliable, parsimonious and predictive statistical framework. Therefore, based on the analysis, we obtained 7 combinations for different ecosystems as shown in Figure 5: forestlands in the SH and SA zones, shrublands in the WT-SH, WT-SA and MT-SA zones, grasslands in the WT and MT zones. It should be emphasized that the study maintains methodological coherence by applying a consistent set of refined zoning criteria across all components—from analytical procedures to result interpretation and graphical visualization.

11. Line 226 wording (“seven regression equations were the best representative”). If no formal multi-model comparison was conducted, please avoid “optimal/best representative.” Use “final selected regression equations” and provide the selection basis (e.g., cross-validation metrics).

Response: We made the change because those models still have space to improve. We replaced “best representative” by a more moderate expression in the main text (including Line 226): “final acceptable representative” / ”final selected model”.

12. Table 1—transparency. Please add per-stratum sample sizes (n), variable definitions/units and preprocessing (e.g., standardization), and VIF diagnostics to facilitate verification and reproducibility.

Response：We provided and added information on sample sizes (n) at each level, variable definitions, and preprocessing procedures, along with VIF diagnostic results.

13. Section 3.4—accounting conventions. (1) Clarify whether the mapping and totals represent potential increase or realized increase accounting for lag effects (i.e., years since retirement). (2) For pixels retired and later reclaimed, explain how they are handled in mapping and in area-integrated accounting. Please make the conventions explicit in Methods and keep them consistent in Results.

Response：1) The mapped distributions and total values represent the realized soil organic carbon increment attributable to the GFGP from 2000 to 2021, as derived from the models presented in Table 1. 2) Areas that were recultivated after initial retirement during the GFGP implementation period were excluded from the calculation of the final benefit in carbon stock increment, as these lands had less contribute to the total GFGP-derived carbon stock increment.

To prevent potential misinterpretation, we implemented the following revisions:

“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. 6-a). 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/m2 at 30 m raster level.”

14. Consistency between Table 2 and modeling strata. Table 2 summarizes carbon increases by MT-SA/WT-SA/WT-SH, whereas earlier modeling strata differ across forest/shrub/grassland. Please align climate-zone definitions between modeling and reporting (or provide a clear crosswalk).

Response：In response to sample size considerations and model significance requirements, we stratified climatic zones differently across vegetation types: SH and SA for forestlands, WT-SH/WT-SA/MT-SA for shrublands, and WT/MT for grasslands. As presented in Table 2, SOC increments are categorized into three aggregated climatic zones for each ecosystem type. To facilitate flexible analysis, the SOC increment for any specific climatic zone combination (as applied in Figure 4 and Table 1) can be derived through summation of corresponding values—for instance, ΔSOC for forestlands in SA equals the sum of MT-SA and WT-SA values. This refined presentation offers more detailed and adaptable information for assessing GFGP-derived carbon stock benefits. To prevent potential confusion, we have incorporated explanatory guidance in the manuscript.

“Significant variations in ΔSOC were observed across different ecosystem types (Fig. 6-b, Table 2). To provide detailed and vegetation-specific insights, Table 2 presents ΔSOC values for three climatic zone combinations associated with each vegetation type.”

15. Uncertainty quantification. The discussion of uncertainty is qualitative. Please add pixel-level uncertainty layers/ranges and confidence intervals for zonal/regional aggregates, or provide the methods and results in the Supplement.

Response: The main uncertainties in this study arise from land cover classification, which was used to identify retired farmlands, and from the model simulation of ΔSOC. For land cover classification, the uncertainty was evaluated using kappa coefficients, as described in Section 2.2. Based on validation with a large number of training and validating samples, the classification achieved a relatively high level of accuracy. So we added the following information:

“During the accuracy assessment phase, the classification performance over the study period consistently achieved kappa coefficients ranging from 0.76 to 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).”

Although it is difficult to determine the exact area of retired farmlands from previous studies on the Loess Plateau, we compared our results with reported trends in earlier research and found them to be generally consistent.

“During the study period, the Chinese central government implemented two phases of GFGP: the first from 1999 to 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).”

For the model simulation, the coefficients used are provided in Table 1. Based on the analysis of the seven multiple linear regression models for ΔSOC (y1–y7), the simulated ΔSOC values fall within the 95% prediction intervals corresponding to the years since retirement for each model as follows: y1: –33.62 to 53.05, y2: –9.30 to 24.57, y3: –50.40 to 55.30, y4: –8.53 to 6.86, y5: –12.68 to 40.21, y6: –23.60 to 22.04, and y7: –12.12 to 11.89.

To further clarify the uncertainties in this study and potential ways to address them, we have added a discussion in Section 4.4, Limitations and Uncertainties:

“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.”

16. Figure standards and language editing. Aim for “read at a glance”: define acronyms at first use, standardize axes and units, and increase font sizes and color separability. Professional English editing is recommended to improve terminological consistency and grammar.

Response: We carefully reviewed the manuscript and added definitions for all acronyms upon their first occurrence. In addition, we standardized terminology throughout the text, replacing terms such as “in-use farmland” with “cultivated farmland”, “climate zone” with “climatic zone”, and “SOC storage” with “SOC stock”.

“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).”

We have standardized the units across all figures, converting SOC units from g C kg⁻¹ to kg C m⁻² in Figures 4 and 5 for consistency. Font sizes and colors have been adjusted to improve readability, particularly in Figure 2. To enhance visual clarity, some non-essential figures from Figure 2 have been moved to the supplementary information.

In addition, the manuscript has been professionally edited by an English language editing service and thoroughly refined for clarity and readability.

We sincerely appreciate your insightful feedback on our manuscript. Please find our detailed, point-by-point responses to your comments below.

Methodology

1.  One of the most significant contributions of this study is its innovative methodology for accounting carbon benefits from ecological restoration. The authors should strengthen the review of commonly used approaches in the literature and explicitly highlight the method in this study as an objective.

Response: Although several studies have examined the dynamics of retired farmlands and associated SOC changes on the Loess Plateau, existing analyses are typically constrained by two key limitations: spatially, they often focus on localized areas, and temporally, they rely on data from only a few discrete years. As a result, a continuous, high-resolution, year-by-year distribution of retired farmland remains unavailable. Moreover, most SOC assessments have been conducted solely on retired farmlands without comparative analysis with adjacent cultivated farmlands, limiting the ability to accurately quantify the net SOC change attributable to ecological restoration. To address these gaps, this study introduces a spatially explicit and temporally continuous monitoring framework. We have accordingly refined both the literature review and research objectives to better highlight the methodological advances and scientific contributions of this work:

“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. Greater clarity is needed regarding the sampling process, particularly how and why certain sampling points were removed.

Response：Paired sampling sites were selected based on spatial proximity to ensure representative coverage across diverse ecosystems, climatic zones, and years since farmland retirement. To minimize human disturbance and ensure accessibility and safety, all potential sites were pre-screened using high-resolution imagery to exclude areas affected by roads, villages, or irrigation ditches. During field sampling, land use status was verified in situ, and sites with objective constraints (e.g., inaccessible terrain or road closures) were excluded from the final selection. In total, 133 sites were removed during the final field sampling process.

“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. 1-b), 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.”

3. The study develops models for seven regional combinations. To support their validity, the sample sizes for each combination should be reported.

Response：The sample sizes (n) for each model stratum have been added in Table 1 in the study. The corresponding numbers of valid soil sample pairs for each stratum were 32 (SH for forestland), 21 (SA for forestland), 22 (WT-SH for shrubland), 13 (WT-SA for shrubland), 20 (MT for shrubland), 16 (WT for grassland), and 11 (MT for grassland), respectively.

4.  Numerous variables were introduced into the SOC change models. While the use of machine learning for variable selection is commendable, more detail is needed regarding the criteria and procedures used for selecting important variables.

Response：The study employed a comprehensive analytical approach combining ANOVA, single-factor regression, all-subsets regression, and stepwise regression to identify appropriate predictor variables for the multivariate linear models of ΔSOC. The variable selection process followed a structured workflow consisting of five distinct stages: data preprocessing, univariate analysis, multivariate analysis, model evaluation, and diagnostic checks. This rigorous procedure led to the identification of several key variables that consistently demonstrated explanatory power across multiple methods.

Although we initially applied random forest to assess variable importance as an additional validation step, the results primarily reinforced our confidence in including "years since retirement" as a key predictor. However, the importance values generated by random forest were not used as formal selection criteria for the final model. To prevent potential confusion regarding our variable selection methodology, we have removed the random forest results from the revised manuscript. This ensures clarity and maintains methodological consistency throughout our analysis.

“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.”

5. The "Materials and Methods" section of the main text omits an introduction to the Random Forest algorithm, despite its results being presented in Figure 6b. This constitutes an oversight in the methodological description and requires revision. It is recommended that the authors add content introducing this algorithm and its application details to the "Methods" section.

Response：Based on the previous response, we removed the methods and results for random forests for clarity.

Figures 

1. Figure 1 (Study area, sampling sites and climatic zones) requires significant revision to enhance clarity and information integration. Currently, Figure 1 separates the presentation of sampling site distribution (subplot b) and climatic zonation (subplot c). As both depict the same region and are crucial for understanding the study's specific climatic zone modeling approach, displaying them separately reduces the efficiency of the figure's expression. Subplot (a), as a macroscopic location map, provides limited and redundant information. The authors should consider merging subplots (b) and (c) into a single integrated main map, using the climatic zonation as a base layer and clearly overlaying the field sampling points. Concurrently, the original subplot (a) could be reduced to an embedded small map for macroscopic positioning only (with a transparent border indicating the main map area and simplified internal elements), and a unified main legend should be created, with all captions placed below the figure.

Response：Based on the reviewer’s suggestion, we simplified the figure 1-a and merged figure 1-b and c.

2.  The figure 2 is currently difficult to interpret. For improved readability, I recommend showing only a few representative years in the main text, while moving the remaining years to the supplementary material.

Response: We have moved most of the figures to the supplementary information, only left six typical years: a (1999-2000), f (2004-2005), k (2009-2010), o (2014-2015), u (2020-2021), v (1999-2021). The readability of the figures was largely improved, and we also mentioned the other year-by-year distribution of farmlands can be found in the supplementary information.

“The annual retired farmlands in the other years can be found in the supplementary material (Fig. S1 a-p).”

3.  The core datasets in Figure 2 (Spatial distribution of retired farmlands) and Figure 3 (Area of ecosystem conversion types from retired farmlands) suffer from scientific ambiguity and require revision. The captions and legends of Figure 2 and Figure 3 fail to clearly and unambiguously define what each annual layer represents regarding "newly added ecosystem types converted from retired farmlands each year" and the "cumulative conversion status" in Figure 2(v). It is recommended that the authors define all layer data in the main text methodology, and revise the captions and legends of Figure 2 and Figure 3 to be concise, accurate, and completely self-explanatory (the caption should explicitly include "Annual New Retired Farmland Conversion" or "Cumulative Distribution," and the legend should clearly label items such as "Converted Forestland" instead of just "Forestland," which can be misleading as it implies all existing forestland).

Response: In figure 2 and 3, each picture and bar indicates the new retired farmlands in every year. To avoid ambiguity, we changed the captions for the figure 2 and 3, and change the legend to “Converted to Forestland, Converted to Shrubland, Converted to Grassland”.

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

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

We also changed the main text accordingly.

“The final retired farmlands were less than the area by summing up annually retired farmlands because of frequent recultivation (Fig. 3-a). The annual area of retired farmlands has been fluctuating throughout the study period with no significant trend (Fig 2-a-u, Fig. 3-b).”

4.  Figure 4 (Comparison of soil organic carbon stocks) has critical omissions and ambiguities in the presentation and description of significance analysis, warranting revision. Figure 4 uses letter annotations to indicate potential significant differences, but the main text lacks a description of the statistical significance analysis methods used for these comparisons. In Figure 4a, the meaning of the letters and the precise scope of comparison are unclear; the X-axis label "1-7" in Figure 4b is difficult to understand, and the significance represented by its letter annotations is also unexplained. The authors should consider providing a detailed description of all significance testing and multiple comparison methods in the main text's methodology section. Concurrently, revise Figure 4's legend/caption to clearly explain the meaning of the letters and the scope of comparison, and correct Figure 4b's X-axis label to improve self-explanatory power.

Response：The letter labels indicating significant differences within groups have been added to the Figure 4 legend, derived from the ANOVA analysis. The X-axis labels “1-7” in Figure 4-b are detailed in the note. Corresponding modifications have also been made in the Methods section.

“Figure 4. SOC stocks in farmlands and retired farmlands (kgC/m2), (a) Comparison of SOC stocks on the Loess Plateau in farmlands retired to different ecosystem types (forestland, shrubland, grassland) with those in adjacent 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 within groups in the ANOVA, for same ecosystem in figure 1-a and for same climatic zone combination figure 1-b.”

5. Figure 5 (Relationship between length of farmland retirement and CCS) contains redundant axis labels and confusing Y-axis units, requiring revision. The current Y-axis label "CCS/g C·kg⁻¹" in Figure 5 may lead to technical confusion, as CCS quantification units are typically area-based (e.g., kg C/m²). Furthermore, the repetition of identical X-axis and Y-axis labels across the seven subplots creates visual redundancy. It is suggested that the authors revise the Y-axis unit to be consistent with the definition and calculation method of CCS. Additionally, to enhance the professionalism and conciseness of the figure, a shared axis label layout is recommended.

Response：The vertical axis of the seven subplots in Figure 5 has been revised from “CCS/g C·kg⁻¹” to “ΔSOC/m²” to better fit the study’s objective of investigating changes/increases in soil organic carbon contents on the Loess Plateau.

6.  Figure 6 (Correlation matrix and variable importance) lacks self-explanatory power and information completeness, warranting revision. Figure 6a's correlation matrix has missing X-axis labels, and some variable names are displayed ambiguously; in Figure 6b's variable importance chart, all abbreviated variables (BD, BIO1-19, etc.) are not fully explained in the caption, affecting the figure's self-explanatory nature. The authors should consider completing the axis labels and variable display for Figure 6a, and append the full meanings of all abbreviated variables in the caption.

Response：We deleted figure 6 for clarity purpose.

Results

1. On the Loess Plateau, pairwise comparisons of SOC are rarely conducted. Although field sampling is challenging and published references are limited, the uncertainty associated with the final SOC estimates should be explicitly presented to strengthen model evaluation.

Response：The main uncertainties in this study arise from land cover classification, which was used to identify retired farmlands, and from the model simulation of ΔSOC.

For the model simulation, the coefficients used are provided in Table 1. Based on the analysis of the seven multiple linear regression models for ΔSOC (y1–y7), the simulated ΔSOC values fall within the 95% prediction intervals corresponding to the years since retirement for each model as follows: y1: –33.62 to 53.05, y2: –9.30 to 24.57, y3: –50.40 to 55.30, y4: –8.53 to 6.86, y5: –12.68 to 40.21, y6: –23.60 to 22.04, and y7: –12.12 to 11.89.

To further clarify the uncertainties in this study and potential ways to address them, we have added a discussion in Section 4.4, Limitations and Uncertainties:

“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.”

Discussion

1. The observed dynamics of farmland retirement and reclamation from 2000 to 2021 are striking. The authors should provide some discussion on possible underlying drivers of these trends.

Response: The reclamation of retired farmland is a commonly observed phenomenon, which can be attributed to several factors. First, while subsidies initially enhanced participation willingness following the launch of the GFGP, their effectiveness declined over time as falling grain prices diminished the real value of subsidy payments. Second, population growth escalated the local demand for food, further incentivizing the recultivation of land. Third, current remote sensing techniques possess a limited capacity to reliably distinguish between abandoned cropland and legitimately retired farmland. In light of these factors, we have strengthened the discussion in our manuscript to better address these recultivation trends. We also changed “reclamation” to “recultivation” for clarity.

“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 to 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 to 2003 during the first GFGP phase, followed by a rebound due to recultivation and subsequent stabilization.”

2. The discussion would benefit from practical recommendations on strategies to enhance SOC on the Loess Plateau, as well as reflections on potential ways to improve model accuracy beyond process-based approaches.

Response: Several strategies show potential for enhancing soil organic carbon (SOC), including reducing grazing intensity, improving subsidy mechanisms to sustain local participation, and selecting appropriate vegetation types coupled with sustainable management of retired farmland to prevent degradation. This study reveals that different vegetation types exhibit distinct SOC accumulation potentials across climatic zones (Figure 5), highlighting the importance of vegetation selection and adaptive management under climate change. Accordingly, we have revised Subsection 4.3 as follows:

“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.”

To enhance the accuracy of this study, increasing the number of sample plots is essential, as it would allow greater flexibility in model development—for instance, by enabling stratification of the study area according to slope, elevation, and soil properties. In addition, establishing permanent observation points to monitor both retired farmlands and adjacent cultivated farmlands would provide reliable pairwise comparisons essential for robust model calibration. Accordingly, we have incorporated the following changes in Section 4.4.

“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).”