Carbon Sequestration

We analyse output from 18 DGVMs that are part of a recent model intercomparison project, TRENDYv10 and the Global Carbon Budget 2021 (GCB2021)⁸ . The models included in the analysis here are CABLE-POP, CLASSIC, CLASSIC-N, CLM5.0, DLEM, IBIS, ISAM, ISBA-CTRIP, JSBACH, LPJ-GUESS, LPJ, LPX-Bern, OCN, ORCHIDEE, ORCHIDEEv3, SDGVM, VISIT, and YIBs (see ref. 8 for model descriptions and setup). One model (JULES-ES-1.1) from TRENDYv10 is not included due to incomplete data. Note, CLASSIC-N^{68 (link)} and ORCHIDEE^{69 (link)} were not part of GCB2021 due to the inclusion of alternate model versions (CLASSIC and ORCHIDEEv3) but are included in TRENDYv10 and this study.
The models are forced with a merged monthly Climate Research Unit (CRU)^{70 (link)} and 6-hourly Japanese 55-year Reanalysis (JRA-55)^{71 (link)} data set. The models are also forced with atmospheric CO₂⁷² , gridded nitrogen deposition⁷³ and nitrogen fertiliser⁷⁴ . DGVMs use the HYDE (v3.3) land-use change data set^{75 (link),76} , which provides annual pasture and cropland areas at a global scale, and includes improvements in the spatial distribution of agricultural regions^{77 (link)}. Several models (CABLE-POP, CLM5.0, JSBACH, LPJ-GUESS, LPJ, and VISIT) also use harmonised land-use change data (LUH2-GCB2021⁸ ), which provides information on sub-grid-scale land-use transitions.
To isolate the response of the land to each driver (CO₂, climate, LULCC), each model performs four simulations: S0 (fixed pre-industrial atmospheric CO₂ and land-use, recycled 1901–1920 climate), S1 (transient atmospheric CO₂, recycled 1901–1920 climate, and fixed pre-industrial land-use), S2 (transient atmospheric CO₂, transient climate, and fixed pre-industrial land-use), and S3 (transient atmospheric CO₂, transient climate, and transient industrial land-use). Nitrogen deposition varies temporally in simulations S1–S3. Therefore, the transient CO₂ (+N deposition) effect on the ‘natural’ land sink is calculated by S1−S0, S2−S1+S0 is the climate effect on the ‘natural’ land sink, S3−S2 is the LULCC effect, and S3 is the net effect. There exists an artefact in the HYDE3.3 data causing a large land-use transition and emission peak around 1960. To correct this, we replace the LULCC estimates for 1959–1961 with the average of 1958 and 1962 in each DGVM.
Net Biome Productivity (NBP) from the S3 simulation represents the net land sink and can therefore be compared to the ‘observed’ Global Carbon Budget land constraint, which is calculated in the GCB as fossil fuel emissions—atmospheric carbon growth rate—ocean carbon sink (see Fig. 1 in the main text and Table 5 in ref. 8 and data taken from 10.18160/gcp-2021). Atmospheric CO₂ concentration measurements began in 1959⁷⁸ , and so this independent constraint on the land sink covers the period 1959–2020, which defines the study period used throughout our analysis.

To model the regrowth of secondary forests we applied a space-for-time substitution method. Instead of tracking the associated Aboveground Carbon (AGC) regrowth over time, the regrowth was estimated by considering the available ages of the standing secondary forest area in 2017 and the associated AGC at the same time. Here we explain the methods used to determine secondary forest AGC using the ESA-CCI Aboveground Biomass (AGB) product (100-m) for the year 2017²³ (see Supplementary Notes 1 and 2). All analysis was carried out in the original product units (AGB) but expressed as AGC by assuming a 2:1 ratio of biomass to carbon²⁴ .
The ESA-CCI AGB product was only released in late 2019 and was in its early phases of development at the time of use. However, given that its spatial resolution was high enough to separate areas of only secondary forest and its recent acquisition warranted its use for this research. Only areas of secondary forest greater than 9000 m² were considered for further analysis, an area approximately equal to 1 pixel of the ESA-CCI product. Despite limiting the study to these larger secondary forest polygons, we were still left with just under 2.5 million polygons of secondary forest to analyse. The secondary forest map was laid over the AGC data and the modal AGC was extracted for each secondary forest polygon using the “zonal_stats” function available in the “rasterstats” module for the programming language “Python” (v3.6). We then aggregated the AGC values by the age of secondary forest and used the median AGC value for each age in further analysis. We applied a bias correction to the median AGC values, subtracting the smallest median value from all values to shift the data to begin at or near 0 Mg C ha⁻¹ AGC for a 1-year-old secondary forest.
Following this, we used six remote sensing products of driving variables widely accepted to influence regrowth of forests. The data products included four environmental drivers (1–4) and two anthropogenic disturbance drivers (5–6): (1) Mean annual downward shortwave radiation (for the period 1985–2017)^{26 (link)}, (2) Mean annual precipitation (for the period 1985–2017)²⁷ , (3) the mean Maximum Cumulative Water Deficit (MCWD) (for the period 1985–2017)^{65 (link),66} , (4) Soil Cation Concentration^{30 (link)}, (5) Annual burned areas (between 2001 and 2017)³¹ and (6) Number of times a secondary forest area was deforested between 1987 and 2017 (repeated deforestations) (this study). These products all have different spatial resolutions (Supplementary Table 1) and so had to be resampled to the size of secondary forest pixels (30-m spatial resolution) using the “resample” package in the Geographic Information System programme, ArcMap10.6. We calculated the key zonal statistics of these variables such as the mean value of the driver affecting a specific area of secondary forest, again using the “zonal_stats” function in Python.
The drivers were then grouped according to numerical limits, such as the 25, 50 and 75th percentiles. We then modelled the AGC for the age of secondary forest under these groupings using the commonly used Chapman-Richard model for regrowth^{67 (link)}: $Y_t=A{\left(1-e^{-kt}\right)}^c\pm \epsilon ;A,k\mathrm{and}c>0$
where Y_t refers to the AGC at age t; A is the AGC asymptote or the AGC of the old-growth forest; k is a growth rate coefficient of Y as a function of age; c is a coefficient that determines the shape of the growth curve; and ε is an error term. We applied the “nls” function available in the “nlstools” package for the statistical software R (v4.0.2)^{68 ,69} . We assumed that after a given amount of time, the AGC could return to levels equivalent to old-growth forests, and reach a pre-calculated asymptote. As such, we extracted the median, bias-corrected AGC value of old-growth forests under each variable condition from the ESA-CCI AGC product to represent the value of the asymptote (Supplementary Fig. 6 and Supplementary Tables 8 and 9). From this, we could also determine if and when the modelled AGC of secondary forest regrowth would reach those equivalent to old-growth forest levels. Forcing the models to “fit” to an expected value for the asymptote value naturally increases the error of our model, partly due to heterogeneity in old-growth forest values within each variable condition.

Forest-level variables were then calculated per square meter and re-scaled to account for the forest area considered (1 ha). This study focused on the Shannon index of the diversity of life-history traits (based on the relative density of trees belonging to every combination of LHT values: SLA/SWD/SS), the organic above-ground carbon sequestration of the trees (MgC/ha/y), the live biomass and deadwood in the forest (t/ha/y) and, finally, the average forest level of LHTs (SLA, SWD, SS). Similar to Pichancourt et al. (2014) (link), the model exemplified the projected state of the forest RI under 100-year climate scenarios where sub-tropical forests were parameterized based on the latitude around Brisbane in Queensland, Australia, and were subject to one of the predicted scenarios of hotter and drier climate change (CSIRO mk3.5 using CMIP3 model: see details in Pichancourt et al. (2014) (link): Table S1).

HLJP is the northernmost and easternmost provincial administrative region in China (43°26′ N~53°33′ N, 121°11′ E~135°05′ E), with 12 prefecture-level cities and the Great Khingan Mountains region, covering a total area of 473,000 km² (Figure 1). HLJP has a cold-temperate continental monsoon climate with an average annual temperature between −5 °C and 5 °C, precipitation between 400 and 650 mm, and 83% to 94% of the total annual precipitation falling during the growing season. The vegetation accumulates a great deal of organic matter that is not suitable for decomposition, and the unique conditions form the black soil humus layer [27 ]. HLJP Black soil resources are extensively spread over cultivated land, forest land, grassland, wetland, and lake ecosystems, and black soil has a high organic matter content and a tremendous carbon sequestration capacity [28 (link),29 (link)]. According to the White Paper on Northeast Black lands (2020) of the Chinese Academy of Sciences, black soils are at a high risk of degradation, including from frequent soil erosion on sloping arable land between 3 and 15 degrees and a progressive loss in ecosystem carbon reserves. Regional risks are exacerbated by factors including strong winds, tillage techniques, snowmelt runoff, and freeze–thaw cycles [30 (link),31 (link),32 (link)]. In 2021, the HLJP government and the Development and Reform Commission proposed the Regulations on the Protection and Utilization of Blackland in Heilongjiang Province and the 14th Five-Year Plan, respectively, to promote the conservation and intensive utilization of black land, construct a new pattern of spatial development and national protection, and enhance the service function and stability of the ecosystem.
Land-sat-TM/ETM images from 2000 and 2010 and Landsat 8 images from 2020 were interpreted to produce land use maps with a spatial resolution of 1 km that met the accuracy requirements of the study. According to the China Land Use Status Classification Standard (GB/T21010-2017), land use is divided into (1) cultivated land, (2) forests, (3) grasslands, (4) waters, (5) construction land, and (6) unutilized land. The Digital Elevation Model (DEM) data with a resolution of 90 × 90 m were obtained using the most recent SRTM V4.1 data by resampling, and the slope data were derived using the DEM data. The Chinese Academy of Sciences (http://www.resdc.cn accessed on 6 December 2022) was consulted for spatially interpolated annual mean temperature and precipitation, DEM data, and land use data. Black soil distribution, national roads, provincial roads and highways vector data were obtained from the National Earth System Science Data Center (http://www.geodata.cn accessed on 6 December 2022). The Heilongjiang Statistical Yearbook (http://tjj.hlj.gov.cn accessed on 6 December 2022) was consulted for its population and GDP figures. Additionally, all data were projected to the Krasovsky_1940_Albers coordinate system to maintain spatial consistency during data processing.