Millennial land carbon emissions in China offset by carbon sinks of the past four decades

Domain: chemistry
Published: 2026-02-26
Editorial outcome: major_revision

Abstract
China experienced substantial millennial-scale climate and anthropogenic land use changes, yet their combined impacts on land carbon dynamics remain largely unexamined. Here, we quantify spatiotemporal changes in terrestrial organic carbon over 851–2022 using a land surface model driven by reconstructed climate and land cover forcings. Simulated results show China’s pre-industrial millennial land carbon dynamics aligned with global carbon stock and atmospheric CO2 fluctuations, such as the ~284 ppm peak in the 12th century linked to land use during the Medieval Climate Anomaly-warmed Song Dynasty. Notably, China’s total land carbon emissions (13±0.5 PgC) accounted for 22% of global land carbon emissions during 1700–1900, with Northeast and Southwest China experiencing the largest historical land carbon losses from intensive deforestation. Nevertheless, the 17.0±1.7 PgC emissions during 851–1980 were fully offset by rapid carbon sinks over 1980–2022, driven by CO2 fertilization and large-scale afforestation. These findings provide insights into China’s historical carbon dynamics and drivers.

Introduction
Atmospheric CO2 exhibited multidecadal to centennial-scale variations over the pre-industrial millennium (850–1850 Common Era, CE), with a peak concentration of ~284 ppm during the 12th century1,2 (Fig. 1a). These variations were likely governed by changes in terrestrial plant and soil organic carbon stocks prior to the era when fossil fuel CO2 emissions became dominant3-5. However, the contribution of regional land carbon emissions and absorptions to these atmospheric CO2 variations and its underlying drivers and mechanisms at millennial scale remain unclear. This knowledge gap stems from uncertainties in reconstructions of regional climate, land use and land cover change (LULCC), and associated land carbon dynamics. In particular, emerging evidence links historical Chinese human activities to millennial-scale atmospheric CO2 changes6. For instance, it has been proposed that the rise in atmospheric CO2 during the 11th century coincided with accelerating urbanization, deforestation, and widespread agricultural expansion south of the Yangtze River during China’s Song Dynasty7,8. Conversely, the sharp CO2 decline in the 17th century is thought to be associated with the fall of the Ming Dynasty, which led to cropland abandonment and subsequent reforestation7. These findings highlight the critical need to quantify past land carbon changes in China, as regional dynamics may have played a key role in global CO2 variability. Beyond historical relevance, China’s modern land carbon dynamics are equally noteworthy9-11. Since the 1970s, the country has implemented large-scale afforestation and reforestation programs, which have doubled forest cover from ~117 million hectares (Mha) in 1980 to ~248 Mha in 20229,10,12,13. However, the recovery trajectory and current status of China’s land carbon pool within the broader context of millennial-scale changes remain poorly understood. It is therefore critical to systematically assess historical and modern dynamics in China’s land carbon stocks and their driving factors, thereby providing insights for global carbon cycle research and climate mitigation strategies.

With the advent of agriculture to satisfy the growing food demand of humanity, anthropogenic land use has become increasingly intensive over the past 12,000 years14. To capture this LULCC dating back to 850 CE, the Land-Use Harmonization (LUH) dataset—integrated with the History Database of the Global Environment (HYDE)—was developed15,16. This dataset has been used in dynamic global vegetation models (DGVMs) to quantify land use impacts on carbon emissions since 190117. However, the LUH dataset exhibits notable biases in Chinese-specific data, including overestimation of forest areas over the past three centuries18 and underestimation of forest expansion during the last four decades19,20. To address this gap, more reliable historical cropland and forest gridded data for China have been reconstructed based on Chinese historical documents21-23. Recent estimates showed that China’s cropland area increased from ~37 Mha in 1000 CE to ~129 Mha in 1999 CE, with a peak of ~135 Mha in 1980 CE24. Some of these datasets have been combined with bookkeeping models to explore the subsequent changes in land carbon stocks for the past 300 years25,26. Nevertheless, these studies are limited to the post-1700 period and overlook the coupled impacts of CO2, climate change, and LULCC on land carbon density and stocks. Notably, DGVM simulations constrained by field experiments suggest that CO2 fertilization alone explains most of the observed increase in the global land carbon sink since the 1960s27. With CO2 concentrations rising rapidly since the 17th century, this factor should have played a pivotal role in shaping historical land carbon sinks.

Parallel to CO2 and land use changes, global climate over the last millennium has been characterized by multidecadal temperature fluctuations and two major climate anomalies, the Medieval Climate Anomaly (MCA) and Little Ice Age (LIA)28,29. Besides, temperature changes have exhibited notable seasonal variations, with model simulations suggesting stronger autumn and winter cooling likely driven by orbital insolation changes30. In China, simulation ensembles from the Paleoclimate Modeling Intercomparison Project Phase 3 (PMIP3), last millennium ensemble runs with CESM (CESM-LME simulations), and proxy-based reconstructions consistently represent the cooling trend over the pre-industrial millennium31-34. Paleoclimate records and climate model simulations further reveal a dramatic weakening of the East Asian summer monsoon and reduced precipitation during the LIA relative to the MCA, particularly in northern China31,35. Climate models have demonstrated the ability to reasonably reproduce these millennial-scale climate variations in China, such as the colder and drier conditions of the LIA29,31,33.

Together, the advancement of improved climate and LULCC reconstructions for the past millennium provides a robust foundation for quantifying long-term variations in China’s land carbon stocks—addressing the limitations of previous studies by integrating coupled drivers across a full millennial timescale. In this study, we use the ORCHIDEE-MICT land surface model and design factorial simulations to disentangle the impact of CO2, climate and LULCC on carbon stock changes in China since 851 CE (Materials and Methods). To drive the model, climate forcing data is first generated by bias-correcting the multi-model mean of selected climate simulations, which best reproduces reconstructed millennial-scale temperature changes across China. For vegetation cover, we calibrate the latest LUH-GCB2023 dataset using the historical-document-based forest and cropland gridded data of China8,13,17,21,22. The land cover changes are then prescribed for the ORCHIDEE-MICT model. With the simulation results, we explore the spatio-temporal changes in organic land carbon, and its response and sensitivity to each driver. Our results provide benchmarking data of historical carbon dynamics in China and highlight the impact of millennial LULCC on land carbon stocks and global carbon balance.

Results
Millennial changes of climate and land cover
The transient simulations for the past millennium in PMIP3/4 and CESM-LME have accounted for major forcing agents, including changes in orbital parameters, solar irradiance, stratospheric aerosols of volcanic origin, LULCC, and greenhouse gas concentrations34,36. Three PMIP simulations and nine CESM-LME simulations were ultimately selected for their ability to reproduce the cooling trend during the pre-industrial millennium (851–1850) and the warming since the LIA (Supplementary Table 1; Supplementary Fig. 1). The multiple-model mean of the 12 simulations exhibits a cooling rate of -0.028 °C/century during 851–1850 for China (Fig. 1a), which falls within the range of -0.044 to -0.022 °C/century derived from proxy-based reconstructions32,33 and climate data assimilations37 (Supplementary Table 1; Fig. 1a). The temperature difference between the warmest century (1051–1150) and coldest century (1751–1850) (i.e., the defined MCA and LIA specific for China in this study, respectively) is comparable between the multiple-model mean of ~0.37°C and reconstructions of 0.17-0.50°C32,33 (Fig. 1d; Supplementary Fig. 2). Although the simulation could reproduce the long-term trend correctly, the centennial-timescale variations in temperature show some differences between model and data. For example, the cold anomalies during the 13th and 15th century in all simulations due to volcanic forcing36 are not seen in reconstructions (Fig. 1a; Supplementary Fig. 1). The bias in the LUH dataset, the LULCC forcing employed in these climate models, may also introduce uncertainties into climate simulations through land-atmosphere interactions. From the end of the LIA, China has been warming by ~1.6°C together with a CO2 increase of ~132 ppm (Fig. 1a). Ultimately, the multi-model mean was employed to drive our ORCHIDEE-MICT model following bias correction, given its ability to reproduce the long-term temperature trend.

Cropland area in the LUH-GCB2023 is generally consistent with the historical-document-based reconstructions for the period 1000–1700 CE, indicating an increase of 31 and 37 Mha in Yang, et al.8 and LUH-GCB202317, respectively (Fig. 1c). Meanwhile, the forest area decreased from 270 to 242 Mha with the strongest deforestation and agricultural development from 1000 to 1200 during the Song Dynasty (Fig. 1b). During 1700–1980, cropland area increased by 69.9 Mha while forest area decreased significantly by -126 Mha, with the largest changes in Northeast China and Southwest China according to the document-based reconstructions21,22. Although the LUH-GCB2023 could reproduce the similar direction and pattern of changes between 1700 and 1980, this dataset significantly underestimates the magnitude of cropland expansion and the forest shrinkage in most regions of China (Supplementary Fig. 3). After 1980, LUH-GCB2023 does not capture the forest expansion and its effect on increasing carbon stock over China38. Our well-benchmarked LULCC dataset, calibrated from historical-document-based reconstructions, should enhance the estimation of millennial-scale land carbon cycling in China by accounting for the substantial role of LULCC.

Evaluation of our simulations over the past four decades
Four simulations were performed with the ORCHIDEE-MICT model, including an equilibrium simulation for the year of 851 (S0) and three transient simulations for 851–2022 driven by time-variant CO2 only (S1), time-variant CO2 and climate change (S2) and all forcings (S3), respectively (Materials and Methods; Supplementary Table 2). The results of the simulation S3 with all forcings presents multidecadal to centennial-timescale variations in carbon fluxes for the past millennium due to climate variability, and followed by rapid increases since 1850 because of rising CO2 (Fig. 2). The S3 run accurately captures the observed changes in carbon fluxes and stocks for China. Consistent with estimates from eddy covariance measurements and machine learning methods41,42, our simulated GPP increases from ~6.3 to ~7.4 PgC/yr during 1990-2022 (Fig. 2a). The NPP increase for China in our estimates agrees with the estimate from remote-sensing based GLASS product43, indicating a total increase by ~0.6 PgC/yr during the past four decades (Fig. 2b). According to a summary of China’s land carbon sink over the past few decades44, there was a divergence among different methods, indicating a sink of 0.21-0.33 PgC/yr from inventory methods, 0.12-0.26 PgC/yr from ecosystem process models and 0.17-0.35 PgC/yr from atmospheric inversion models. Our result shows a net biome production (NBP) of 0.21±0.06 PgC/yr for the past four decades in China, falling within the range of those estimates (Fig. 2c).

Our estimates of current vegetation carbon stocks are comparable with satellite-based mappings10,45 and carbon storage sums derived from field measurements46 (Fig. 2d). More importantly, our estimated rate of increase in China’s vegetation carbon stock over the past two decades (0.25 PgC/yr) corroborates a recent estimate of woody biomass accumulation (0.21 PgC/yr)10. As for the soil organic carbon (SOC) at a depth of 0–100 cm for the past four decades, our result shows an increase from 82.5 to 87.2 PgC during 1980–2022 (Fig. 2e). This increase is in good agreement with previous upscaling estimates from China’s National Soil Survey, which indicates a rapid increase before 2000 followed by a slowing growth rate thereafter47,48.

Historical dynamics of land carbon stocks
The simulated total land ecosystem carbon decreased by -2.4 PgC in China during 851–1700, 80% of which resulted from LULCC-induced vegetation carbon loss (Fig. 3d). During the deforestation period of 1700–1980, our result depicts a land carbon source of 14.6±1.7 PgC (Fig. 3d). This carbon source magnitude exceeds the estimated range of 4.50-9.54 PgC by Ge, et al.50 and that of 2.94~5.61PgC by Yang, et al.25, but comparable to the latest estimate of 6.2–12.4 PgC—derived from a bookkeeping model that accounts for vegetation and soil carbon losses26. Specifically, the LULCC-induced emission is 23.4 PgC, which were partially offset by the CO2 fertilization effect in our estimates (8.2±1.7 PgC) (Fig. 3d), highlighting the importance of accounting for the often-overlooked CO2 effect in historical land carbon change assessments. For the period 1980–2022, the cumulative carbon sink (16.0±3.9 PgC) in vegetation and soil in our S3 experiment (Fig. 3a, c) is larger than the modeled estimate of 8.9±0.8 PgC for the similar period (1980–2019) in China38. This discrepancy may partially arise from uncertainties in grassland cover changes13, as these ecosystems store substantial amounts of soil carbon51, even though forest and cropland cover dynamics are well constrained. In our estimates, the total land carbon stock increased by 17.5±4.3 PgC during 1980–2022, an amount equivalent to the historical land carbon emission (17.0±1.7 PgC) during 851–1980 (Fig. 3d). It means that the land biosphere carbon is near neutral in China over 851–2022.

From the MCA and the LIA, a period marked by pronounced changes in both climate and land cover while the net CO2 changes were negligible (Fig. 1), total land organic carbon declined by 3.8 PgC, primarily driven by LULCC (Fig. 2f). By comparison, climate changes (i.e., cooling and drying) alone reduced vegetation carbon but raised SOC, yielding a net carbon gain (0.02 PgC) (Supplementary Fig. 4; simulation S2 minus S1). Spatially, climate impacts on total land carbon closely matched SOC change pattern, indicating the dominating role of SOC responses in driving climate-induced land carbon changes (Supplementary Fig. 4). Temporally, our multi-model mean and previous studies consistently show notable seasonal variations with stronger cooling in autumn and winter from the MCA to the LIA over China (Supplementary Fig. 4), likely driven by orbital insolation changes30,31. Our simulations reveal that climate-induced heterotrophic respiration (RH) decreased significantly in autumn, primarily following the seasonality of temperature changes. Consequently, the reduction in RH fully compensated the impact of NPP decrease, resulting in SOC increase.

Spatially, Northeast China (-3 PgC) and Southwest China (-2.6 PgC) experienced the largest net land carbon losses between 851 and 2022, driven by substantial forest degradation, especially during 1700–1980 (Fig. 4; Supplementary Fig. 5). Typically, the arable lands in Northeast China increased from only 0.6% to 26.6% during 1661-199924. However, both regions have emerged as the strongest carbon sinks over the past four decades, sequestering 2.7 PgC and 3.9 PgC, respectively, benefiting from elevated CO2 levels and afforestation efforts (Figs. 1, 4; Supplementary Fig. 5). In contrast, North China and Qinghai-Tibetan Plateau have functioned as the largest net carbon sink during 851-2022, accumulating 2.2 PgC and 2 PgC, respectively, with CO2 fertilization being the dominant contributing factor (Fig. 4; Supplementary Fig. 5). Notably, recent studies indicate that the Tibetan Plateau currently serves as a substantial greenhouse gas sink, nearly offsetting its regional anthropogenic emissions52.

The atmospheric CO2 and δ13C records combined with a deconvolution approach suggested an organic land carbon decrease by ~29 PgC from 850 to 1900 CE, which was partly absorbed by Ocean (~6 PgC) and consequently contributed to an atmospheric CO2 increase by ~11 ppm3 (Fig. 5). This net land carbon source was the tradeoff between global LULCC emissions at a rate of ~60 PgC/kyr3,6,53 and peatland development by ~30 PgC/kyr54,55. Meanwhile, our simulated terrestrial carbon emissions for China amount to 15.5 PgC during 851–1900, accounts for ~26% of global LULCC emissions (~60 PgC), underscoring the important role of regional LULCC on global land carbon budget and further on atmospheric CO2 variability at millennial scale. In detail, a reduction of ~1.7 PgC in terrestrial carbon stocks occurred during the 11th century, within the Song Dynasty (Fig. 5). This reduction is attributed to intensive deforestation and agricultural expansion, and is likely driven by the favorable climate conditions associated with the warmth of the MCA (Fig. 1). In particular, this regional carbon stock reduction coincides with a global land carbon loss of ~14 PgC in the 11th century, contributing to the consequent peak atmospheric CO2 concentration of ~284 ppm over the pre-industrial millennium (850-1850)3. However, our results do not support the hypothesis that the decline in global atmospheric CO2 in the late 1500s and early 1600s was ascribed to cropland abandonment and reforestation of late Ming Dynasty7 (Figs. 1c, 5). Other processes, such as land use change and secondary succession after the large-scale depopulation of the Americas, organic carbon burying due to anthropogenic soil erosion, and peatland development have potential to explain this decline of atmospheric CO2 54,56,57. Subsequently, the rapid increase in atmospheric CO2 since 1700 results from emissions from land-use change and fossil fuel combustion17,26. China’s land carbon emissions (13 PgC) account for 22% of global land carbon emissions (58 PgC) during 1700–1900 (Fig. 5).

Sensitivity of simulated land carbon changes to drivers
The impacts of CO2, climate change and LULCC on vegetation and soil carbon are decoupled with our factorial simulations in order to understand the sensitivity of land carbon changes in China to drivers from the model (Supplementary Fig. 6). The sensitivity of vegetation and soil carbon in China to CO2 changes is 6.7±1.5 PgC/100ppm and 4.1±2.4 PgC/100ppm, respectively (Supplementary Fig. 7). As a result, simulated CO2 fertilization effect during 851-2022 led to a net carbon sink of 16.5±6.1 PgC, which totally compensated the carbon loss from climate change and LULCC (Fig. 3). Given the significance of the CO2 fertilization effect in this study, we conducted two sensitivity simulations by setting atmospheric CO2 concentrations at 380 ppm and 540 ppm, respectively. The results indicate a 20% increase in NPP, which is consistent with findings from free-air CO2 enrichment (FACE) experiments (Supplementary Fig. 8). These experiments reported a 22% increase in NPP and an 18% increase in dry matter accumulation (i.e., aboveground wood production) when CO2 levels were elevated from ~380 ppm to ~540 ppm58.

As for the climatic impact, we observed delayed response of vegetation and soil carbon in our simulations. For example, the cold anomalies around 1274 and 1466 CE coincided with SOC peaks in 1288 and 1480 CE, respectively (Fig. 6). Consequently, cross-correlation analyses were performed to quantify lagged responses and sensitivities of carbon stocks to temperature and precipitation. Results show that vegetation carbon storage correlates positively with precipitation and forest area, with peak correlation coefficients at a lag of 12 yeas and 1 year, respectively (Fig. 6). In contrast, temperature increase and forest loss cause peak SOC reductions at longer lags of 15 and 4 years, respectively. Accounting for these lags, China’s SOC exhibits sensitivities of -0.36 PgC/°C to temperature and 0.07 PgC/Mha to forest area, while vegetation carbon shows sensitivities of 1.1 PgC/(100 mm/yr) to precipitation and 0.08 PgC/Mha to forest area (Supplementary Fig. 7). Spatially, vegetation carbon in all regions exhibits a positive correlation with precipitation, with a markedly higher sensitivity in northern China than in southern China (Supplementary Fig. 9). In contrast, temperature shows a negative correlation with both vegetation carbon and SOC across most regions, excluding two cold regions: Inner Mongolia and the Qinghai-Tibetan Plateau (QTP) (Supplementary Figs. 4, 9). These findings align with results from intensive field surveys and remote sensing studies of China over the past two decades51,59 and highlight the pronounced—but often delayed—negative impacts of warming and drying on carbon storage. Notably, such lagged responses by years to decades have been understudied due to limited long-term observations60,61, underlining the importance to pay more attention to ecosystem carbon turnover time and long-term land carbon dynamics62,63. Afforestation since 1980 has restored forest area as well as vegetation carbon stock back to the level of 1700, yet SOC still exhibits a ~3.4 PgC deficit in 2022 relative to 1700 (Fig. 3). The asymmetric response and sensitivity of SOC to forest loss and increase, i.e., the relative fast release from forest clearing and the slow uptake from afforestation64, explain the recent slow recovery of SOC (Supplementary Fig. 7).

In summary, our study reveals millennial land carbon changes for China by using the state-of-the-art land surface model driven by improved climate forcings and LULCC data. Between 851 and 1900, China’s land carbon emissions amounted to ~15.5 PgC, accounting for ~26% of global LULCC-driven emissions. This underscores the dominant role of anthropogenic land use changes in shaping terrestrial carbon stocks, as well as the influence of regional land carbon budgets on pre-industrial atmospheric CO2 variability. Nevertheless, current land carbon stocks in China have rebounded to their highest levels in the past millennium which is attributed to the rapid carbon sinks for the past four decades because of CO2 fertilization effects (47%) and large-scale afforestation initiatives (53%). We reveal the delayed and asymmetric response of land carbon stocks to climate and forest changes, which should be considered for predicting future land carbon sink potential and policy-making of land management strategies.

Methods
Model description
Factorial simulations were performed using the land surface model ORCHIDEE-MICT v8.8.0 to reproduce and attribute change in carbon fluxes and stocks in China over 851–2022. The ORCHIDEE-MICT (https://orchidee.ipsl.fr/about-orchidee/) is a process-based dynamic vegetation model that simulates carbon, water and energy exchanges between the atmosphere and biosphere65,66. The model defines vegetation types as 13 plant functional types (PFTs), including one for bare soil type, eight for trees, and one each for C3 natural grass, C4 natural grass, C3 crops, and C4 crops, based on morphology, leaf type, phenology, photosynthetic pathway and climatic regions. All PFTs go through the same suite of processes: photosynthesis, phenology, carbon allocation to plant biomass, carbon flow from living biomass to litter pools after senescence and/or mortality, from litter to soil carbon pools, and heterotrophic respiration. Each PFT has its specific parameter values.

ORCHIDEE-MICT describes carbon exchanges between the atmosphere and the different ecosystem carbon pools and includes the effects of land use change on terrestrial carbon cycle64,67. There are biomass pools (leaves, roots, sapwood, heartwood, etc.), litter pools (structural and metabolic), soil carbon pools (active, slow, and passive) and wood product pools (with turnover times of 10 years and 100 years). During deforestation, a certain fraction of above-ground woody biomass (i.e. sapwood and heartwood) is lost as instant CO2 flux into the atmosphere (i.e., due to on-site disturbance), and the remaining wood is entering into wood product pools. Since the forest clearing and turnover of wood product during the early land use is uncertain, this study focuses on total organic land carbon of biomass, litter and soil. Given parts of crop production are harvested and consumed by humans directly, the model assumes 45% of biomass turnover is harvested and exported outside the ecosystem. The model has integrated a prognostic fire module to simulate open grassland and forest fires arising from both natural and anthropogenic ignitions68. The ORCHIDEE-MICT model has been widely applied in millennial-scale studies54,69, and has been extensively evaluated against present observed carbon and hydrological fluxes at global and regional scales65,70. In this work, we further validated the performance of the model in reproducing the carbon fluxes and stocks in China for the past a few decades and the CO2 fertilization effect on NPP (Fig. 2; Supplementary Fig. 8).

Climate forcing data
Historical data of atmospheric CO2 concentrations, climate and land-cover maps are required to force the model. The observation-based historical CO2 concentrations for the Common Era were collected following the CMIP6/PMIP4 past1000 simulation protocol1. For the meteorological forcing, we used CRUJRA version 2.4 at 0.5° spatial resolution for the years 1901–2022, which was created in the context of TRENDY 2023 project (https://orchidas.lsce.ipsl.fr/dev/protocols/). For years between 851 and 1900, we considered outputs from the 12 PMIP3/4 past1000 simulations and 13 CESM-LME simulations34,36,71 (Supplementary Table 1). The simulated temperature was compared with proxy-based temperature reconstructions32,33, observation-based reanalysis39, and climate data assimilations37,72 (Supplementary Fig. 1). Three PMIP simulations (INM-CM4-8, CCSM4, CSIRO-Mk3L-1-2) and nine CESM-LME simulations were ultimately selected for their ability to reproduce the cooling trend during the pre-industrial millennium (851–1850) and the warming since 1800 (Supplementary Table 1; Supplementary Fig. 2). Consequently, we calculated the multiple-model mean at monthly temporal resolution over these 12 simulations for variables that are available (i.e., 2-m air temperature, precipitation, near surface specific humidity, surface longwave and shortwave incident radiation, wind speed and surface pressure). To further improve the temporal resolution to 6 hours and spatial resolution to 0.5°, and eliminate systematic error in outputs of climate simulations, these variables of multiple-model mean were bias-corrected73. We used either absolute change or relative change for calibration of each variable using the baseline of the pre-industrial (PI) climate, as follows:
Vinput,yr = VPI × (Vmodel,yr ÷ Vmodel,PI) (1)
Vinput,yr = VPI + (Vmodel,yr - Vmodel,PI) (2)
where the climate field V is calibrated input data for our land surface model; VPI is base climate of PI period using CRUJRA gridded climate dataset between 1901–1910; Vmodel,yr and Vmodel,PI are simulated climate from multiple-model mean for the past millennium and for 1901–1910, respectively. Equation (1) was used for calibration of flux variables (incident shortwave and longwave radiation); while for state variables (wind speed, temperature, humidity, pressure) as well as precipitation, Equation (2) was used for calibration in case of negative values.

Land cover mapping
Past land cover changes were estimated using provincial cropland and forest area data from historical documents in China in previous work21,22. Cropland and forest area allocation model were built based on the assumption that the lands with good natural conditions will be cultivated first, gradually followed by marginal lands with harsh natural environments. The potential maximum distribution extent of forest in the absence of human activities constrains the allocation extent of historical forest area. Given the expansion of cropland in the past several centuries, almost all historical cropland was located within the present cropland scope. Spatially explicit reconstructions of cropland and forest cover in China since 1700 were developed at a resolution of 10 km based on these documents and allocation models21,22. This dataset has previously been used to study the impact of LULCC on ecosystem carbon stocks with bookkeeping model25.

As the model distinguishes eight forest PFTs, and two cropland PFTs, we need to disaggregate historical-document-based land types into corresponding ORCHIDEE PFTs. Historical land-cover maps (i.e., LUH-GCB2023) at 0.25° spatial resolution for the period from 851 to 2022 were reconstructed under the program of Global Carbon Budget of 202317. The land-cover dataset was reconstructed by merging the European Space Agency (ESA) Climate Change Initiative (CCI) dataset (www.esa-landcover-cci.org/) with the LUH2v2 dataset (http://luh.umd.edu/). The merging steps are described at https://orchidas.lsce.ipsl.fr/dev/lccci/. Specifically, the land-use dataset is based on the updated HYDE 3.3 reconstructions (https://public.yoda.uu.nl/geo/UU01/67UHB4.html)16,17, the performance of which has not previously been evaluated over China. As land-cover maps for the past 300 years (i.e., 1700–2022) from the LUH-GCB2023 still has some uncertainties for China according to our comparison (Fig. 1; Supplementary Fig. 3), we calibrated it with the above-mentioned reconstructions of forest and cropland cover changes for China from historical documents and allocation models21,22. In practice, the percentage of forest and crop PFTs for each grid cell were adjusted to match reconstructions from Li, et al.21 and He, et al.22 while the PFT types and relative abundance were taken from the LUH-GCB2023. For years prior to 1700, we performed calibration using the vegetation map generated for 1700 and adopted land-cover change data from LUH-GCB2023 to ensure the consistency and coherence of our produced PFT data. In particular, the improved PFT dataset is generally consistent with the latest cropland reconstructions derived from historical documents24 (Fig. 1c). Finally, the grass PFTs were adjusted to match the reconstructed total grassland cover in China over the past millennium13.

Simulation protocol
In order to decouple the effect of rising levels of atmospheric CO2, climate change, and land cover change, four model experiments (S0-S3) were conducted (Supplementary Table 2). First, a spin-up simulation was conducted (S0). In this simulation, a full model was run for 100 years to reach steady state for biomass, followed by a 20,000-year spin-up of the soil carbon submodule to reach steady state for soil carbon pools. In this step, meteorological forcing data for 1960–1990 was looped and used to approximate the warmer-than-preindustrial climate of the early-middle Holocene while the land-cover map was set as the year of 851. Then, a full model was run for another 100 years with looping climate of 851–870 and static CO2 concentration of 280 ppm to reach a steady state for land carbon fluxes and stocks for the year of 851. In the equilibrium state, the mean NBP over the latest 100 years is less than 0.01 PgC/yr and the trend is less than 0.01 PgC/yr per century.

Three transient simulations from 851 to 2022 were performed on the basis of the spin-up simulation S0 (Supplementary Table 2). The S1 simulation was conducted with observed atmospheric CO2 concentration for all years while looping climate of 851–870 and static land cover in 851 were used. The S2 simulation was conducted with observed atmospheric CO2 concentration and climate changes while static land cover in 851 was used for all years. The S3 simulation used all time-varying forcings of CO2, climate change and LULCC. Consequently, the impacts of atmospheric CO2 variation, climate change, and land cover change were decoupled with S1-S0, S2-S1 and S3-S2, respectively.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript quantifies the spatiotemporal changes in vegetation and soil organic carbon from 851 to 2022, using a land surface model driven by newly reconstructed climate and land cover forcings. The authors report a millennial-scale decline in China’s organic land carbon by 5.7 PgC, primarily attributable to early land use activities. Carbon emissions of 12.6 PgC during 851-1980 has been offset by rapidly growing carbon sinks over the last four decades owing to the CO2 fertilization and large-scale afforestation. I find this study to be original and of high significance. The manuscript is generally well written and logically structured. However, some aspects of model simulations and presentation of results require clearer explanation. Therefore, I recommend publication after some revisions.\nGeneral comments:\nExperimental design: The authors selected the CCSM4 model based on its temperature trend for 851–1850 CE. However, relying on a single model as the driving fields may introduce uncertainties, particularly in the spatial patterns of temperature and precipitation. Could this choice potentially affect the robustness of the results? Given the availability of multiple simulations for the last millennium, such as the CESM-LME (Otto-Bliesner et al., 2016) with multiple ensemble members, have the authors considered or evaluated these datasets? The CESM-LME could better reduce the influence of internal variability, thereby providing a clearer representation of millennial-scale trends.\nRegarding the millennial-scale trend results: Since the period around 850–1000 CE overlaps with the warming phase of the Medieval Climate Anomaly and coincides with a rising trend in atmospheric CO2 concentration, this may have attenuated the overall calculated trend across the entire millennium. Given this, why did the authors choose not to calculate the trend specifically for the period 1000–1850 CE, or from 1000 CE onward?\nThe study should clarify which LUCC forcing dataset was used to drive the CCSM4 model. This appears to differ from the LUCC forcing applied in the land surface model simulations (i.e., the modified LUH-GCB2023). It may be worth discussing whether this difference could influence the results.\nMinor comments:\nL73-74: The time periods assigned to the Little Ice Age and the Medieval Climate Anomaly should be more precisely defined. For instance, the interval 1850–1900 already experienced significant anthropogenic influence, with temperatures beginning to rise, and thus may not strictly fall within the LIA. Additionally, it would be preferable to refer to the “Medieval Climate Anomaly” rather than the “Medieval Warm Period”.\nL84-86: While this paragraph focuses on climate model simulations of past-millennium climate change, the final statement unexpectedly shifts to LULCC. To enhance coherence, it would be helpful to introduce some context beforehand.\nOver the past millennium, temperature changes have exhibited notable seasonal variations. For instance, model simulations suggest stronger cooling in autumn and winter, likely driven by orbital insolation changes (Lucke et al., 2021). Could such seasonal differences have an impact on carbon dynamics?\nL91-93: Please include a proper citation for the LUH-GCB2023 dataset when it is first mentioned.\nL107-109: Assimilation datasets, such as PHYDA (Steiger et al., 2018), have also been evaluated by scholars for their performance in simulating East Asian temperatures over the past millennium (Sun et al., 2024). The authors may consider comparing these datasets to further strengthen their analysis.\nL377-379: Please specify which climate variables are used as driving inputs here, including relevant details such as their temporal resolution.\nL416-420: Please include a reference to Figure S4 here to facilitate reader comprehension.\nFigures 3 and 4: Please provide panel labels (e.g., a, b, c, etc.) and corresponding descriptions for each subfigure.\nI suggest including a table summarizing the details of each experiment. This would help readers better understand the experimental design.\nFigure S8: Are these results derived from the CCSM4 simulation?\nFigure S9: The data sources for panels a-d should be clearly specified.\nReferences:\nLucke, L. J., Schurer, A. P., Wilson, R., & Hegerl, G. C. (2021). Orbital forcing strongly influences seasonal temperature trends during the last millennium. Geophysical Research Letters, 48, e2020GL088776\nOtto-Bliesner B L, Brady E C, Fasullo J, Jahn A, Landrum L, Stevenson S, Rosenbloom N, Mai A, Strand G. 2016. Climate variability and change since 850 CE: An ensemble approach with the Community Earth System Model. Bull Am Meteorol Soc, 97: 735–754\nSteiger N J, Smerdon J E, Cook E R, Cook B I. 2018. A reconstruction of global hydroclimate and dynamical variables over the Common Era. Sci Data, 5: 180086\nSun W, Chen D, Lü G, Ning L, Gao C, Zhang R, Wang B, Liu J. 2024. Impacts of major volcanic eruptions over the past two millennia on both global and Chinese climates: A review. Science China Earth Sciences, 67(1): 61‒78", "reviewer_2": "The authors present simulations of the terrestrial carbon stock changes over the period 851 to 2022, arguing that anthropogenic changes in land use decreased the carbon stock, primarily in vegetation and to some extent in soils until the 1980s when massive reforestation began. This reforestation, but also the CO2 fertilization effect, is shown to balance the past carbon losses.\nThe aim is to evaluate long-term changes in land carbon stocks in China and the dynamic response to anthropogenic land use and climate change.\nThe simulated data is compared to measurements and other model data, and show general agreement with GPP, NPP and NBP for the later decades. A comparison of the simulation S3, which is driven by CO2 fertilization, dynamic climate, and land use change, shows good agreement also with vegetation carbon. For the soil C, the picture is less straightforward, but soil C is a vast stock and difficult to upscale and constrain. It is not surprising to find a lower agreement there.\nThe study offers an improved understanding of the past terrestrial C dynamics in China, based on recently published improved data sets as input data.\nTo study the past in this way help constrain the past atmospheric CO2 dynamics, and the study suggests that China stood for 45 percent of net global land emissions over the period 800-1850. This is relevant in order to understand the role of past land management, which may provide insights important for current land use change.\nHowever, a large proportion of the gains of land C during the later decades, proposed by this study, is a direct result of CO2 fertilization. As this is a modeling study, the result hinges on the model’s sensitivity to CO2 fertilization. This is not discussed.\nIn fact, I’m missing a wider discussion about the possibility of overestimating (or underestimating) the CO2 fertilization effect, which really drives the results in this study (for the relevant time period). This needs to be more transparent to the reader. It would not be amiss to also include some type of sensitivity experiment here to actually show what the error bars are. No potential errors to any of the numbers are presented.\nSome more notes on the CO2 fertilization:\nLine 70: You should include the role of CO2 fertilization here, since it’s the main driver of recent change.\nLine 174 …is mostly attributed to increase in biomass associated with forest expansion.\nLooking at the attributions for in figure S7, the main driver of these NBP changes is CO2 fertilization. If this is a misunderstanding on my part, then you will have to explain better how to understand the data.\nThe only C pool where LUC has a larger effect than CO2 is the vegetation C (according to Fig 3)\nLine 186 – 190. There are large losses of land C during the 1700-1980 in your S3 simulation, but comparable to recent estimates. However, without the strong CO2 fertilization here, the losses would be far greater than what others have reported. Hence, the CO2 fertilization is of great importance, which should be discussed.\nLine 308 - 309: Please provide how much of the effect (in percent) is from CO2 fertilization and afforestation respectively. For transparency.\nThe entire section Land carbon stocks influenced by climate and land cover changes is all about how the model responds to forcings. This section feels more like a by-product of the analysis, and it’s written in a way that could be misinterpreted by a reader. You have to clarify that it is the model performance you are looking into. It is the model that is sensitive to these forcings with certain time-lags, not China.\nOtherwise, the authors have used relevant methodologies and have done extensive work to find the best climate model for their studies. They have also worked to better reproduce the land cover for the period. The work could be extended by adding a sensitivity experiment related to the CO2 fertilization.\nIn general it seems possible to reproduce the work with the provided method description.\nFurther technical comments:\nLine 47: Besides, China implemented… This is not given a time stamp. When was this, is it the recent reforestation? When did it begin? Provide more details here, please.", "reviewer_3": "Comments to NCOMMS-25-56683_Chen et al. “Millennial land carbon emissions in China offset by carbon sinks of the past four decades”\nMajor comments:\nThe manuscript presents a well-executed study exploring and quantifying the spatiotemporal changes in the vegetation and soil organic carbon over 851-2022 in China using a land surface model. Overall, the subject matter is of interest to a wide audience, but I have major concerns regarding the presentation of the results. Some descriptions are often overly general and lack precise quantitative support from the data. Specifically, there are several results appear inconsistent with the data shown in the figures. In addition, some figures are cited imprecisely. I strongly recommend that the authors carefully revise this entire section, ensure that all descriptions are accurately and directly supported by the data, and double-check all figures (panels) reference for correctness and appropriateness. Poor presentations make it hard for readers to understand the scientific importance of this work. Results should have been presented in a clearly time sequence, and the relative contributions of different driving forces could have been presented and explained more clearly.\nSpecific comments:\nLine 43-45: Rephrase this sentence. Introducing “850 to 1900 Common Era (CE)” suggests this work focuses on this time period. In addition, variations of CO2 concentration from ~279 to ~295 ppm at millennial scale is too small to attract readers interests. Moreover, too many time periods with different lengths in the first two paragraphs confuse the readers.\nLine 43: It’s too abrupt to shift from global to China.\nLines 47-48: which time period?\nLine 51-52: Delete “and their dynamic responses…”. There is a semantic overlap between 'changes' and 'their response', the changes in carbon stocks are inherently driven by land use and climate change.\nLine 54: What does 12 ka mean?\nLine 69-71: redundant with Lines 50-52 and Lines 36-43\nLine 121: Delete the point after regions\nLine 127: Please be specific, which period (1xxx-1xxx) during the Song Dynasty?\nLine 133 Move “over China” to the end of this sentence.\nLine 134-136: This comparison makes no much sense.\nLine 145-146: The authors have defined these two periods as the warmest (1051-1150) and coldest (1751-1850) century in Line 110-111, but not Medieval Warm Period (900-1400) and Little Ice Age (1400-1900) (Line 73-74).\nLine 155: This section “Historical dynamics of land carbon fluxes and stocks”. “three transient simulations for 851-2022 using the ORCHIDEE-MICT model”. However, the results just show for the last four decades until Line 184.\nLine 168: Rewrite “the past a few decades”\nLine 197-202: These sentences make me confused. The carbon losses in Northeast and Southwest China occurred during 851-1980, but not during 1980-2022 as shown in Fig. S5.\nLine 201-202: The authors just said “carbon losses were driven by substantial forest degradation between 851 and 2022”, which was be contradictory to “afforestation efforts over the past four decades (1980-2022)”.\nLine 204: Qinghai-Tibetan Plateau acts as a carbon source during 851-1700 as shown in Fig. S5.\nLine 205: Qinghai-Tibetan Plateau has sequestered 2.11 PgC during 851-2022 as shown in Fig. S5, which is inconsistent with 2.2 PgC in Fig. S6.\nLines 217-218: “Taken together, the millennial carbon emissions in China account for ~45% of 217 net global land emissions or ~22% of global LULCC emissions”. It’s too hard to convince the readers\nLine 217-220: Fig. 4 did not show these results (45% and 22%).\nLine 267-268: Could you please specify what these distinct trajectories are? A more detailed description of the differences between vegetation and soil carbon trends would be helpful.\nLine 271-272: Which Figure have shown these results: “climate cooling and 271 drying alone reduced vegetation carbon (-0.11 PgC) but raised SOC (0.23 PgC)”, please cite each figure or panel in the parentheses that follows the result.\nLine 281: -29 gC/m2/°C in Fig. 5."}, "author_response": "Response to reviewers\n(Ms. No. NCOMMS-25-56683)\nWe thank the three reviewers for their valuable time and diligent efforts in conducting\na thorough review of our manuscript. Their insightful comments and constructive\nsuggestions have been instrumental in enhancing the quality and rigor of the paper.\nIn this revision, we have carefully addressed all the points raised by the reviewers.\nBelow are our detailed point-by-point responses to each comment.\nResponse to reviewer #1\nReviewer #1 (Remarks to the Author):\nThis manuscript quantifies the spatiotemporal changes in vegetation and soil organic\ncarbon from 851 to 2022, using a land surface model driven by newly reconstructed\nclimate and land cover forcings. The authors report a millennial-scale decline in\nChina’s organic land carbon by 5.7 PgC, primarily attributable to early land use\nactivities. Carbon emissions of 12.6 PgC during 851-1980 has been offset by rapidly\ngrowing carbon sinks over the last four decades owing to the CO2 fertilization and\nlarge-scale afforestation. I find this study to be original and of high significance. The\nmanuscript is generally well written and logically structured. However, some aspects\nof model simulations and presentation of results require clearer explanation. Therefore,\nI recommend publication after some revisions.\nGeneral comments:\nExperimental design: The authors selected the CCSM4 model based on its\ntemperature trend for 851–1850 CE. However, relying on a single model as the driving\nfields may introduce uncertainties, particularly in the spatial patterns of temperature\nand precipitation. Could this choice potentially affect the robustness of the results?\nGiven the availability of multiple simulations for the last millennium, such as the CESM-\nLME (Otto-Bliesner et al., 2016) with multiple ensemble members, have the authors\nconsidered or evaluated these datasets? The CESM-LME could better reduce the\ninfluence of internal variability, thereby providing a clearer representation of millennial-\nscale trends.\nReply: Thanks for the suggestion and recommendation. We agree that a single model\nmay introduce uncertainties regarding the spatial pattern of climate change. Therefore,\nwe first compared 12 PMIP simulations and 13 LME simulations with data\nassimilations from LMR and PHYDA, and proxy-based reconstructions. As this work\nfocused on the long-term variations of climate on carbon stocks, we selected three\nPMIP simulations and nine LME simulations that best reproduce the cooling trend\nduring the pre-industrial millennium (851-1850) and the warming since 1800\n(Supplementary Table 1 of the revised manuscript). The average climate of these\nsimulations was then used to calibrate the climate forcing as input for ORCHIDEE land\nsurface model. All simulations have been re-run. We concur that the multi-model mean\ncan mitigate the uncertainties associated with the last millennium’s climate data, as\nyou have suggested.\nRegarding the millennial-scale trend results: Since the period around 850–1000 CE\noverlaps with the warming phase of the Medieval Climate Anomaly and coincides with\na rising trend in atmospheric CO2 concentration, this may have attenuated the overall\ncalculated trend across the entire millennium. Given this, why did the authors choose\nnot to calculate the trend specifically for the period 1000–1850 CE, or from 1000 CE\nonward?\nReply: As noted in the IPCC Sixth AR6, MCA and LIA lack standardized definitions\nand exhibit regional variability at global scale. Based on our analysis of millennial-\nscale climate data for China, there is a consistent decreasing temperature trend over\nthe pre-industrial millennium, but with inconsistencies in the delineation of warm and\ncold periods (Fig. 1a; Supplementary Figs. 1, 2). Therefore, in this study, we define\nthe Medieval Climate Anomaly (MCA) and Little Ice Age (LIA) specific to China using\nthe multi-model mean of simulated climate data that force our vegetation model.\nFurthermore, the defined MCA (i.e., 1050–1150 CE) roughly aligns with the CO peak\n2\nobserved during the 12th century (Fig. 1a). We have discussed the impact of the\ncooling transition from the MCA to the LIA on terrestrial carbon dynamics\n(Supplementary Fig. 4). In addition, as different drivers dominate each period, we have\nattributed land carbon changes across three periods to distinct driving factors (Fig. 3).\nThe study should clarify which LUCC forcing dataset was used to drive the CCSM4\nmodel. This appears to differ from the LUCC forcing applied in the land surface model\nsimulations (i.e., the modified LUH-GCB2023). It may be worth discussing whether\nthis difference could influence the results.\nReply: We have now used multiple model ensemble rather than only relying on\nCCSM4 simulation as you suggested. The LME and PMIP simulations were using\nLUH1 and LUH2 land use forcing, respectively. The LUH datasets exhibit notable\nbiases in its Chinese-specific data, which may introduce some uncertainties to the\nsimulated millennial climate changes. We have added a discussion of the potential\ninfluence in Lines 134–136 of the revised manuscript: “The bias in the LUH dataset,\nthe LULCC forcing employed in these climate models, may also introduce\nuncertainties into climate simulations through land-atmosphere interactions.”\nMinor comments:\nL73-74: The time periods assigned to the Little Ice Age and the Medieval Climate\nAnomaly should be more precisely defined. For instance, the interval 1850–1900\nalready experienced significant anthropogenic influence, with temperatures beginning\nto rise, and thus may not strictly fall within the LIA. Additionally, it would be preferable\nto refer to the “Medieval Climate Anomaly” rather than the “Medieval Warm Period”.\nReply: We have replaced \"Medieval Warm Period\" with \"Medieval Climate Anomaly\"\nconsistently throughout the manuscript. The time periods assigned to the MCA and\nLIA here were derived from the IPCC AR5; this inclusion was intended to assist non-\nspecialist readers. However, as noted in the IPCC Sixth AR6, MCA and LIA lack\nstandardized definitions and exhibit regional variability at global scale. Consequently,\nwe have removed these cited time periods to avoid controversy and instead provided\ndata-derived timeframes for both climate events in Fig. 1a and Supplementary Fig. 2.\nL84-86: While this paragraph focuses on climate model simulations of past-millennium\nclimate change, the final statement unexpectedly shifts to LULCC. To enhance\ncoherence, it would be helpful to introduce some context beforehand.\nReply: The second and third paragraphs introduce advances in LULCC and climate\nresearch, respectively. We have relocated this summary sentence to the last\nparagraph of the introduction: “Together, the advancement of improved climate and\nLULCC reconstructions for the past millennium provides a robust foundation for\nquantifying long-term variations in China’s land carbon stocks—addressing the\nlimitations of previous studies by integrating coupled drivers across a full millennial\ntimescale.”\nOver the past millennium, temperature changes have exhibited notable seasonal\nvariations. For instance, model simulations suggest stronger cooling in autumn and\nwinter, likely driven by orbital insolation changes (Lucke et al., 2021). Could such\nseasonal differences have an impact on carbon dynamics?\nReply: Many thanks for the very interesting suggestion. Indeed, our multi-model mean\nand previous studies consistently show notable seasonal variations with stronger\ncooling in autumn and winter from the MCA to the LIA over China (Fig. S4, also shown\nhere). We have added a paragraph to present and discuss the climate and carbon\nchanges from the MCA to the LIA in the section “Historical dynamics of land carbon\nstocks”. The main findings are as follows: “Climate-induced heterotrophic respiration\n(RH) decreased significantly in autumn, primarily following the seasonality of\ntemperature changes. Consequently, the reduction in SOC fully compensated the\nimpact of NPP decrease.”\nSupplementary Fig. 4. of the revised manuscript. Climate changes and climate-\ninduced land carbon changes from MCA to LIA in China.\nL91-93: Please include a proper citation for the LUH-GCB2023 dataset when it is first\nmentioned.\nReply: We have added the citation of original LUH as well as the recent LUH-\nGCB2023 dataset when it is first mentioned.\nL107-109: Assimilation datasets, such as PHYDA (Steiger et al., 2018), have also\nbeen evaluated by scholars for their performance in simulating East Asian\ntemperatures over the past millennium (Sun et al., 2024). The authors may consider\ncomparing these datasets to further strengthen their analysis.\nReply: Thanks for the suggestion. We have added a comparison of the PHYDA\nproduct with other datasets (Supplementary Table 1; Supplementary Fig. 1). Although\nPHYDA data could show reasonable temperature variability and response to volcanic\neruption events, its long-term variations are not consistent with LMR and proxy-based\nreconstructions. Specifically, it does not correctly reproduce the cooling trend during\nthe pre-industrial millennial and the warming since 1800 in China (Supplementary\nTable 1; Supplementary Fig. 1).\nL377-379: Please specify which climate variables are used as driving inputs here,\nincluding relevant details such as their temporal resolution.\nReply: We have now specified the variables used (i.e., 2-m air temperature,\nprecipitation, near surface specific humidity, surface longwave and shortwave incident\nradiation, wind speed and surface pressure) and their temporal resolution (i.e., 6 hours\nafter bias correction).\nL416-420: Please include a reference to Figure S4 here to facilitate reader\ncomprehension.\nReply: We appreciate your suggestion – the figure has now been referenced here.\nFigures 3 and 4: Please provide panel labels (e.g., a, b, c, etc.) and corresponding\ndescriptions for each subfigure.\nReply: The figures and their captions have been revised and improved accordingly.\nI suggest including a table summarizing the details of each experiment. This would\nhelp readers better understand the experimental design.\nReply: Thanks for the suggestion. We have added a Supplementary Table 2 to\nsummarize the details of each experiment. This will indeed enhance the clarity of the\nexperimental design for readers.\nSupplementary Table 2. Summary of the experiments designed in this study.\nClimate\nSimulation CO Vegetation Purpose\n2\nforcing\nS0 280 ppm 851-870 851 Spin-up of the model. The start for S1-S3.\n1172-year transient simulation considering\nS1 851-2022 851-870 851\nCO changes. Start from S0.\n2\n1172-year transient simulation considering\nS2 851-2022 851-2022 851\nCO and climate changes. Start from S0.\n2\n1172-year transient simulation considering all\nS3 851-2022 851-2022 851-2022\nfactors. Start from S0.\nTwo 322-year transient simulations with\nadjustments made to the model parameters\n1701-\nS3_eCO2 1701-2022 1701-2022 of CO fertilization effect at potential\n2\n2022\nmaximum and minimum levels. Start from\n1700 of S3.\nDiagnose the model’s sensitivity to elevated\nCO and compare it with FACE data\n2\n380 ppm, 2001-\nS_eCO2 2001-2010 (Supplementary Fig. 8). Two equilibrium\n540 ppm 2010\nsimulations of the model at CO concentration\n2\nof 380 ppm and 540 ppm.\nFigure S8: Are these results derived from the CCSM4 simulation?\nReply: The revised figure (Supplementary Fig. 4) now presents climate changes\nbetween LIA and MCA from multiple-model mean (i.e., climate input of our ORCHIDEE\nmodel), and the climate-induced land carbon changes in China from our simulations\n(simulations S2 minus S1). The caption has been updated.\nFigure S9: The data sources for panels a-d should be clearly specified.\nReply: The caption has been updated to specify that the data are derived from our\nfactorial simulations (now Supplementary Fig. 7). Besides, we have added the\nsensitivity of simulated carbon stocks to elevated CO and its uncertainties.\nLucke, L. J., Schurer, A. P., Wilson, R., & Hegerl, G. C. (2021). Orbital forcing strongly\ninfluences seasonal temperature trends during the last millennium. Geophysical\nResearch Letters, 48, e2020GL088776\nOtto-Bliesner B L, Brady E C, Fasullo J, Jahn A, Landrum L, Stevenson S,\nRosenbloom N, Mai A, Strand G. 2016. Climate variability and change since 850 CE:\nAn ensemble approach with the Community Earth System Model. Bull Am Meteorol\nSoc, 97: 735–754\nSteiger N J, Smerdon J E, Cook E R, Cook B I. 2018. A reconstruction of global\nhydroclimate and dynamical variables over the Common Era. Sci Data, 5: 180086\nSun W, Chen D, Lü G, Ning L, Gao C, Zhang R, Wang B, Liu J. 2024. Impacts of major\nvolcanic eruptions over the past two millennia on both global and Chinese climates: A\nreview. Science China Earth Sciences, 67(1): 61‒78\nReply: Thank you for recommending these highly relevant references and datasets,\nwhich have been appropriately cited and compared in the manuscript.\nResponse to reviewer #2\nReviewer #2 (Remarks to the Author):\nThe authors present simulations of the terrestrial carbon stock changes over the period\n851 to 2022, arguing that anthropogenic changes in land use decreased the carbon\nstock, primarily in vegetation and to some extent in soils until the 1980s when massive\nreforestation began. This reforestation, but also the CO2 fertilization effect, is shown\nto balance the past carbon losses.\nThe aim is to evaluate long-term changes in land carbon stocks in China and the\ndynamic response to anthropogenic land use and climate change.\nThe simulated data is compared to measurements and other model data, and show\ngeneral agreement with GPP, NPP and NBP for the later decades. A comparison of\nthe simulation S3, which is driven by CO2 fertilization, dynamic climate, and land use\nchange, shows good agreement also with vegetation carbon. For the soil C, the picture\nis less straightforward, but soil C is a vast stock and difficult to upscale and constrain.\nIt is not surprising to find a lower agreement there.\nReply: In this revision, we have adopted updated climate and land cover forcing data.\nFor the land cover component, we referenced a recent publication on millennial-scale\ngrassland cover reconstructions for China (He et al., 2025), which indicates that\ngrassland area only experienced a slight decrease since 1900. In the previous version\nof our manuscript, however, grassland cover compensates fluctuations in forest cover,\nresulting in stable total vegetated area and soil organic carbon (SOC) stocks. By\nincorporating this improved land cover dataset, our simulated SOC now shows\nsubstantially better agreement with improved SOC estimates from recent upscaling\nstudies for the past four decades (Fig. 2e).\nHe, F., Yang, F., Wang, Y., 2025. Reconstructing forest and grassland cover changes\nin China over the past millennium. Science China Earth Sciences 68 (1), 94–110.\nThe study offers an improved understanding of the past terrestrial C dynamics in China,\nbased on recently published improved data sets as input data.\nTo study the past in this way help constrain the past atmospheric CO2 dynamics, and\nthe study suggests that China stood for 45 percent of net global land emissions over\nthe period 800-1850. This is relevant in order to understand the role of past land\nmanagement, which may provide insights important for current land use change.\nReply: Thank you for your recognition of our work.\nHowever, a large proportion of the gains of land C during the later decades, proposed\nby this study, is a direct result of CO2 fertilization. As this is a modeling study, the\nresult hinges on the model’s sensitivity to CO2 fertilization. This is not discussed.\nIn fact, I’m missing a wider discussion about the possibility of overestimating (or\nunderestimating) the CO2 fertilization effect, which really drives the results in this study\n(for the relevant time period). This needs to be more transparent to the reader. It would\nnot be amiss to also include some type of sensitivity experiment here to actually show\nwhat the error bars are. No potential errors to any of the numbers are presented.\nReply: Many thanks for your constructive suggestions. We fully agree that the\nimportance and uncertainties of the CO fertilization effect on land carbon dynamics\n2\nwarrant further discussion. First, we conducted two global sensitivity experiments with\natmospheric CO concentrations set at 380 ppm and 540 ppm to quantify the net\n2\nprimary productivity (NPP) response to elevated CO (i.e., the S_eCO2 simulation in\nFig. 8 of the revised manuscript), which is consistent with findings from previous free-\nair CO enrichment (FACE) experiments. Second, we added two sensitivity\n2\nexperiments focusing on China over the past three centuries, involving adjustments\nmade to the model parameters of CO fertilization effect at potential maximum and\n2\nminimum levels (i.e., the S3_eCO2 simulation in Supplementary Table 2). The\nassociated potential uncertainties have been incorporated into Fig. 3, and\nSupplementary Figs. 6 and 7. Finally, we revised the subheading from “Land carbon\nstocks influenced by climate and land cover changes” to “Sensitivity of simulated land\ncarbon changes to drivers” to explicitly include discussions on the importance and\nuncertainties of the CO fertilization effect.\nLine 70: You should include the role of CO2 fertilization here, since it’s the main driver\nof recent change.\nReply: The importance of CO fertilization on land carbon sinks has been highlighted\n2\nin the introduction part. “Notably, DGVM simulations constrained by field experiments\nsuggest that CO fertilization alone explains most of the observed increase in the\n2\nglobal land carbon sink since the 1960s. With CO concentrations rising rapidly since\n2\nthe 17th century, this factor should have played a pivotal role in shaping historical land\ncarbon sinks.”\nLine 174 …is mostly attributed to increase in biomass associated with forest expansion.\nLooking at the attributions for in figure S7, the main driver of these NBP changes is\nCO2 fertilization. If this is a misunderstanding on my part, then you will have to explain\nbetter how to understand the data.\nThe only C pool where LUC has a larger effect than CO2 is the vegetation C (according\nto Fig 3)\nReply: For the past four decades, LULCC and CO fertilization have equal importance\n2\nfor the increase in terrestrial carbon (Fig. 3d). We have removed this sentence here.\nDiscussions on the impacts of CO , climate change, and LULCC have been\n2\nconsolidated into the final section.\nLine 186 – 190. There are large losses of land C during the 1700-1980 in your S3\nsimulation, but comparable to recent estimates. However, without the strong CO2\nfertilization here, the losses would be far greater than what others have reported.\nHence, the CO2 fertilization is of great importance, which should be discussed.\nReply: The importance of CO fertilization for historical land carbon change\n2\nassessments is emphasized here. “Specifically, the LULCC-induced emission is 21.2\nPgC, which was partially offset by the CO fertilization effect (8.8±1.7 PgC) (Fig. 3d),\n2\nhighlighting the importance of accounting for the often-overlooked CO effect in\n2\nhistorical land carbon change assessments.”\nLine 308 - 309: Please provide how much of the effect (in percent) is from CO2\nfertilization and afforestation respectively. For transparency.\nReply: Thank you for your suggestion. We have added the percentage contributions\nof the CO fertilization effect (47%) and afforestation (53%) to this sentence.\nis all about how the model responds to forcings. This section feels more like a by-\nproduct of the analysis, and it’s written in a way that could be misinterpreted by a\nreader. You have to clarify that it is the model performance you are looking into. It is\nthe model that is sensitive to these forcings with certain time-lags, not China.\nReply: In the revised manuscript, the subheading has been revised as “Sensitivity of\nsimulated land carbon changes to drivers”. In this section, we clarify that we present\nthe model-derived sensitivity of China’s terrestrial carbon changes to drivers. Since all\nsimulations were re-run using improved climate and LULCC data, the corresponding\nanalyses were redone. Notably, the key result that the delayed response of land\ncarbon by several years to more than a decade persists in the new simulations. We\ntherefore believe these findings are worth reporting and would help us to understand\nthe simulated long-term land carbon dynamics in China.\nIn addition, we added a paragraph to discuss the CO fertilization effect in this section.\nhere). “Given the significance of the CO fertilization effect in this study, we conducted\n2\ntwo sensitivity simulations by setting atmospheric CO concentrations at 380 ppm and\n2\n540 ppm, respectively. The results indicate a 20% increase in NPP, which is consistent\nwith findings from free-air CO enrichment (FACE) experiments (Supplementary Fig.\n2\n8). These experiments reported a 22% increase in NPP and an 18% increase in dry\nmatter accumulation (i.e., aboveground wood production) when CO levels were\n2\nelevated from ~380 ppm to ~540 ppm (Norby, 2025).”\nNorby, R.J., 2025. Forest productivity response to elevated CO2 in free-air CO2\nenrichment experiments: the 23 percent solution, revisited. New Phytol 246 (5), 1952–\n1959.\nSupplementary Fig. 8. Comparison of simulated and observed NPP changes in\nresponse to elevated CO .\nwork to find the best climate model for their studies. They have also worked to better\nreproduce the land cover for the period. The work could be extended by adding a\nsensitivity experiment related to the CO2 fertilization.\nReply: We have incorporated discussions on the importance and uncertainties of the\nCO fertilization effect in the revised manuscript, as detailed above.\ndescription.\nFurther technical comments:\nLine 47: Besides, China implemented… This is not given a time stamp. When was this,\nis it the recent reforestation? When did it begin? Provide more details here, please.\nReply: These large-scale afforestation projects began in the 1970s. We have made it\nclear here.\nResponse to reviewer #3\nReviewer #3 (Remarks to the Author):\nComments to NCOMMS-25-56683_Chen et al. “Millennial land carbon emissions in\nChina offset by carbon sinks of the past four decades”\nMajor comments:\nThe manuscript presents a well-executed study exploring and quantifying the\nspatiotemporal changes in the vegetation and soil organic carbon over 851-2022 in\nChina using a land surface model. Overall, the subject matter is of interest to a wide\naudience, but I have major concerns regarding the presentation of the results. Some\ndescriptions are often overly general and lack precise quantitative support from the\ndata. Specifically, there are several results appear inconsistent with the data shown in\nthe figures. In addition, some figures are cited imprecisely. I strongly recommend that\nthe authors carefully revise this entire section, ensure that all descriptions are\naccurately and directly supported by the data, and double-check all figures (panels)\nreference for correctness and appropriateness. Poor presentations make it hard for\nreaders to understand the scientific importance of this work. Results should have been\npresented in a clearly time sequence, and the relative contributions of different driving\nforces could have been presented and explained more clearly.\nReply: Many thanks for your meticulous review and insightful suggestions. We have\nsubstantially enhanced the presentation of our results: the data and figures have been\nappropriately cited in the relevant text to facilitate reader understanding, and both the\ncontent and structure of the manuscript have been refined. Specifically, we\naggregated some text into a new section titled “Evaluation of our simulations over the\npast four decades”. In the section “Historical dynamics of land carbon stocks”, we now\nfirst present the temporal evolutions of land carbon stocks in a clear chronological\nsequence, followed by an analysis of their spatial patterns, and finally a comparison\nwith historical global carbon budget estimates. We believe the result section is more\ncoherent and logical in the revised manuscript.\nSpecific comments:\nLine 43-45: Rephrase this sentence. Introducing “850 to 1900 Common Era (CE)”\nsuggests this work focuses on this time period. In addition, variations of CO2\nconcentration from ~279 to ~295 ppm at millennial scale is too small to attract readers\ninterests. Moreover, too many time periods with different lengths in the first two\nparagraphs confuse the readers.\nReply: We have carefully rephrased these sentences to avoid using too many time\nperiods. The logic and fluency of this paragraph have been significantly enhanced.\nLine 43: It’s too abrupt to shift from global to China.\nReply: Thank you for your insightful comment. We now explicitly connect global\nmillennial CO variability (linked to terrestrial carbon dynamics) with emerging\n2\nevidence of historical Chinese human activities influencing such global changes. This\nintegration highlights China’s potential role in shaping pre-industrial atmospheric CO\n2\nfluctuations.\nLines 47-48: which time period?\nReply: These large-scale afforestation projects began in the 1970s. We have made it\nclear.\nLine 51-52: Delete “and their dynamic responses…”. There is a semantic overlap\nbetween 'changes' and 'their response', the changes in carbon stocks are inherently\ndriven by land use and climate change.\nReply: Thanks for the suggestion. We have rewritten this sentence: “It is therefore\ncritical to systematically assess historical and modern dynamics in China’s land\ncarbon stocks and their driving factors, thereby providing insights for global carbon\ncycle research and climate mitigation strategies.”\nLine 54: What does 12 ka mean?\nReply: It means 12000 years ago. We revised it.\nLine 69-71: redundant with Lines 50-52 and Lines 36-43\nReply: The sentence has been removed.\nLine 121: Delete the point after regions\nReply: Thanks for the careful review. We revised it.\nLine 127: Please be specific, which period (1xxx-1xxx) during the Song Dynasty?\nReply: The specific period is 1000–1200 CE. We have added it.\nLine 133 Move “over China” to the end of this sentence.\nReply: Corrected.\nLine 134-136: This comparison makes no much sense.\nReply: This sentence has been removed.\nLine 145-146: The authors have defined these two periods as the warmest (1051-1150)\nand coldest (1751-1850) century in Line 110-111, but not Medieval Warm Period (900-\n1400) and Little Ice Age (1400-1900) (Line 73-74).\nReply: As noted in the IPCC Sixth Assessment Report (AR6), the Medieval Climate\nAnomaly (MCA) and Little Ice Age (LIA) lack standardized definitions and exhibit\nregional variability at global scale. Therefore, we have removed their definitions from\nthe introduction. Instead, in this study, we define the warmest and coldest centuries of\nthe pre-industrial millennium specific to China as the MCA and LIA, respectively. We\nconsider this will facilitate readers’ better understanding of our results.\nLine 155: This section “Historical dynamics of land carbon fluxes and stocks”. “three\ntransient simulations for 851-2022 using the ORCHIDEE-MICT model”. However, the\nresults just show for the last four decades until Line 184.\nReply: We aggregated this text into a new section “Evaluation of our simulations over\nthe past four decades”. We believe the structure of the Result section is more coherent\nand logical in the revised manuscript.\nLine 168: Rewrite “the past a few decades”\nReply: Corrected as “the past few decades”.\nLine 197-202: These sentences make me confused. The carbon losses in Northeast\nand Southwest China occurred during 851-1980, but not during 1980-2022 as shown\nin Fig. S5.\nReply: We have reviewed these sentences and confirmed that they accurately reflect\nthe carbon loss during 1700–1980 and the strong carbon sinks over the past four\ndecades. These sentences have been improved to eliminate potential\nmisunderstanding.\nLine 201-202: The authors just said “carbon losses were driven by substantial forest\ndegradation between 851 and 2022”, which was be contradictory to “afforestation\nefforts over the past four decades (1980-2022)”.\nReply: We have revised these two sentences to eliminate potential misunderstanding.\n“Spatially, Northeast China (-3 PgC) and Southwest China (-2.6 PgC) experienced the\nlargest net land carbon losses between 851 and 2022, driven by substantial forest\ndegradation, especially during 1700–1980 (Fig. 4; Supplementary Fig. 5). However,\nboth regions have emerged as the strongest carbon sinks over the past four decades,\nsequestering 2.4 PgC and 3.2 PgC, respectively, benefiting from elevated CO levels\n2\nand afforestation efforts (Figs. 1, 4; Supplementary Fig. 5).”\nLine 204: Qinghai-Tibetan Plateau acts as a carbon source during 851-1700 as shown\nin Fig. S5.\nReply: We have revised the sentence to ensure clarity and precision.\nLine 205: Qinghai-Tibetan Plateau has sequestered 2.11 PgC during 851-2022 as\nshown in Fig. S5, which is inconsistent with 2.2 PgC in Fig. S6.\nReply: Many thanks for your very careful review. This minor inconsistency arises from\nthe summation and rounding procedures utilized in the calculations. We have revised\nthe relevant values accordingly.\nLines 217-218: “Taken together, the millennial carbon emissions in China account for\n~45% of 217 net global land emissions or ~22% of global LULCC emissions”. It’s too\nhard to convince the readers\nReply: During 851–1900, China’s land carbon decreased continuously while global\nland cover shows fluctuations because of LULCC-induced emission and peatland\ndevelopment-induced absorption. We have improved this sentence to make a better\ncomparison.\nLine 217-220: Fig. 4 did not show these results (45% and 22%).\nReply: We have rewritten this sentence. “Our simulated terrestrial carbon emissions\nfor China amount to 15.5 PgC during 851–1900, accounts for ~26% of global LULCC\nemissions (~60 PgC)”. In additions, two shadings have been added to the Fig. 5 to\nbetter illustrate the two periods we discussed.\nLine 267-268: Could you please specify what these distinct trajectories are? A more\ndetailed description of the differences between vegetation and soil carbon trends\nwould be helpful.\nReply: We have re-run all simulations using improved climate and land use/land cover\nchange (LULCC) data. The simulated soil organic carbon (SOC) dynamics now aligns\nwith changes in vegetation carbon. This sentence has been removed. This section has\nbeen reorganized to present the results of examining the sensitivity of simulated land\ncarbon changes to individual drivers.\nLine 271-272: Which Figure have shown these results: “climate cooling and 271 drying\nalone reduced vegetation carbon (-0.11 PgC) but raised SOC (0.23 PgC)”, please cite\neach figure or panel in the parentheses that follows the result.\nReply: We have revised the supplementary figure (Supplementary Fig. 4), which\nillustrates changes between the MCA and the LIA, and cited it in the relevant text.\nLine 281: -29 gC/m2/°C in Fig. 5.\nReply: Thanks for pointing it out. We have double-checked to ensure consistency\nbetween the numbers presented in the text and those in the figures."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the authors' revisions and thoughtful responses. The manuscript has been significantly improved and is now suitable for publication.", "reviewer_2": "I applaud the authors of this study for their dedication to develop their already ambitious analysis even further. I find that my notes on the importance of CO2 fertilization has been dealt with in the current version. The paper is much stronger now.\nI only have some minor comments (text edits).\nLine 226-229: reveal that climate induced heterotrophic respiration (RH) decreased significantly in autumn, primarily following the seasonality of temperature changes. Consequently, the reduction in SOC fully compensated the impact of NPP decrease.\nThe last sentence doesn’t make sense. I presume you mean that the reduction of SOC led to less RH which then compensated the impact from NPP decrease? I don’t find it reasonable to draw the line directly from reduction in SOC to compensate NPP decrease. Really minor detail, but important to me. It is of course implied, but still.\nLine 310: In contrast, temperature increase and forest lose cause peak SOC.\nTypo – lose – loss", "reviewer_3": "The authors have addressed all my concerns, I have no other comments."}, "author_response": ""}]

Decision letter