Maps with 1 km resolution reveal increases in above- and belowground forest biomass carbon pools in China over the past 20 years

To quantify the ecological consequences of recent nationwide restoration efforts in China, spatially explicit information on forest biomass carbon stock changes over the past 20 years is critical. However, long-term biomass tracking at the national scale remains challenging as it requires continuous and high-resolution monitoring. Here, we characterize the changes in the above- and belowground biomass carbon (AGBC and BGBC) of forests in China between 2002 and 2021 at 1 km spatial resolution by integrating multiple types of remote sensing observations with intensive field measurements through regression and machine learning approaches. On average, 8.6 ± 0.6 and 2.2 ± 0.1 PgC were stored in above- and belowground live forests in China. Over the last 20 years, the total forest biomass carbon pool in China has increased at a rate of 114.5 ± 16.3 TgC yr⁻¹ (approximately 1.1 % yr⁻¹). The most pronounced forest biomass carbon stock gains occurred in central to southern China, including the southern Loess Plateau, Qinling mountains, southwestern karsts and southeastern forests. While the combined use of multi-source remote sensing data provides a powerful tool to assess the forest biomass carbon changes, future research is also needed to explore the drivers of the observed woody biomass trends and to evaluate the degree to which biomass gains will translate into biodiverse, healthy ecosystems that are sustainable. Annual forest above- and belowground biomass maps for China are now available at https://doi.org/10.6084/m9.figshare.21931161.v1 (Chen, 2023).

1 Introduction

Forest biomass carbon stock contributes to over 90 % of the global vegetation biomass carbon pool (Ma et al., 2021). As a net outcome of carbon gains from photosynthesis and carbon losses from respiration, mortality and disturbances, forest biomass carbon stock (approximately 50 % of biomass) is a critical indicator of ecosystem function and ecosystem services, such as carbon sequestration, wood production and resource allocation (Kumar and Mutanga, 2017). Accurate forest biomass carbon stock monitoring over space and time is thus essential for assessing ecosystem management strategies and mitigation policies (Kumar and Mutanga, 2017).

In recent decades, remote sensing tools have been integral in our efforts to map aboveground biomass (AGB) or carbon stock (AGBC). By combining satellite imagery (e.g. MODIS) and airborne lidar signals, forest cover and canopy height can be mapped across large spatial scales (Hu et al., 2016; Saatchi et al., 2011; Su et al., 2016; Tong et al., 2020; Xu et al., 2021). Apart from optical images and lidar signals, microwaves can provide more detailed insights into subcanopy forest structure and AGBC due to their ability to penetrate the canopy. Active microwave techniques, i.e. Synthetic Aperture Radar (SAR) backscatters, facilitate high-resolution (e.g. 100 m) AGB mapping, but the temporal coverage is limited (Cartus et al., 2012; Bouvet et al., 2018). Conversely, vegetation optical depth (VOD) retrieved from multiple passive microwave sensors can be used to produce long-term continuous AGB maps (Frappart et al., 2020; Liu et al., 2011, 2015), yet at a coarse spatial resolution (e.g. 0.25^∘). Because different remote sensing techniques have their advantages and pitfalls, combining these techniques and complementing them with direct ground measurements is integral to maximizing the accuracy and precision of biomass carbon estimations across space and time.

Another source of uncertainty in vegetation biomass carbon stocks is the extent of biomass that is stored belowground as roots. While AGBC mapping is facilitated by a suite of emerging remote sensing techniques, investigating the spatiotemporal variation in the belowground biomass carbon pool (BGBC) remains challenging despite the large contribution of roots to total carbon storage (Huang et al., 2021; Ma et al., 2021). To map BGBC, the commonly used approach is to combine aboveground biomass information with vegetation-type-specific ratios of belowground biomass (BGB) to AGB (i.e. root–shoot ratio or RSR; Xu et al., 2021; Saatchi et al., 2011). Because field studies indicate a near-linear relationship between log-transformed BGB and AGB (Enquist Brian and Niklas Karl, 2002), BGB variations at large scales have often been approximated using this relationship (Spawn et al., 2020). To capture the complex relationship between BGB and biotic or abiotic variables (e.g. stand age, heat and water availability), machine learning algorithms have been applied to map BGB (Huang et al., 2021) and root mass fractions (Ma et al., 2021) globally. However, the reference plots were unevenly distributed across the world and limited in developing countries, leading to some uncertainties in BGB and BGBC estimation within those regions (Huang et al., 2021).

China has been implementing national-scale afforestation and reforestation programmes since the late 1990s (Lu et al., 2018), promoting vegetation cover and carbon storage in the Loess Plateau and the southwestern karst regions, etc. (C. Chen et al., 2019; Niu et al., 2019; Tong et al., 2018). A spatial understanding of forest biomass trends can help evaluate the efficiency of ecological restoration programmes. High-quality, high-resolution and long-term continuous woody biomass monitoring in China has remained challenging (Zhang et al., 2019; Huang et al., 2019).

In this study, by integrating multi-source remote sensing data with large quantities of plot measurements, we produced 1 km resolution above- and belowground forest biomass carbon pool maps for China during the past 20 years (2002–2021). This dataset, which is available at https://doi.org/10.6084/m9.figshare.21931161.v1 (Chen, 2023) could provide new insights into forest carbon stock changes in China over the past 2 decades.

2 Materials and methods

To map above- and belowground forest biomass carbon stock in China during 2002–2021, we (1) calibrated a SAR-based high-resolution forest aboveground biomass map in China based on massive field measurements of AGBC during 2011–2015, (2) extended the AGBC time series to 2002–2021 by referring to the tree and short vegetation cover retrieved from optical remote sensing, (3) calibrated the AGBC time series in some specific areas using a long-term, integrated microwave-based VOD dataset and (4) mapped forestlands' BGBC through a random forest model developed based on the in situ records in the published literature. The basic procedure is shown in Fig. 1 and described below.

Figure 1 Workflow of forest biomass carbon pool monitoring in China during 2002–2021. AGBC and BGBC are the aboveground and belowground biomass carbon, VCF is the vegetation continuous field, LPDR VOD is the global land parameter data record vegetation optical depth, and CLCD is the China land cover dataset.

2.1 A benchmark map of forest aboveground biomass carbon (AGBC) in China

By combining multiple satellite observations of SAR backscatter, including the L-band ALOS PALSAR (Advanced Land Observing Satellite Phased Array L-band Synthetic Aperture Radar) and C-band Envisat ASAR (Advanced Synthetic Aperture Radar) around the year 2010, the first global high-resolution (100 m) forest AGB dataset, GlobBiomass 2010, was published through the European Space Agency (ESA)'s Data User Element (DUE) project (Santoro et al., 2021), whose relative root mean square error (RMSE) was below 30 % (Mialon et al., 2020). Apart from GlobBiomass 2010, another high-resolution (30 m) forest AGB for China was produced by relating the (lidar)-derived ICESat/GLAS (Ice, Cloud and land Elevation Satellite/Geoscience Laser Altimeter System) footprint AGB to various variables derived from Landsat optical images (Huang et al., 2019). Because the ICESat data in 2006 were applied as the training target of the random forest model, Huang et al. (2019)'s dataset refers to the AGB status in 2006. According to a recent validation study, GlobBiomass and Huang et al. (2019)'s AGB performed the best among all existing AGB datasets in China (Chang et al., 2021). Mean forest canopy heights and tree coverage are also good indicators of the spatial pattern of forest biomass. The high-resolution (30 m) forest canopy height map for China was developed by interpolating the ICESat-2 and Global Ecosystem Dynamics Investigation (GEDI) data in 2019 through a neural network (Liu et al., 2022), while the tree cover map at the same resolution was derived from cloud-free growing season composite Landsat 7 data in around 2010 (Hansen et al., 2013). We resampled GlobBiomass from 100 m resolution ($\mathrm{1}/\mathrm{1125}$^∘) to $\mathrm{1}/\mathrm{1200}$^∘ (approximately 90 m) and averaged Huang et al. (2019)'s AGB map, canopy height map and tree cover map to the same resolution.

A reviewable, consistent ecosystem carbon stock inventory was conducted in China between 2011 and 2015 (Tang et al., 2018). We requested the AGB carbon stock (AGBC) data at more than 5000 30 × 30 m sized forest plots from the authors. Due to the scale mismatch between the maps of biomass, canopy height or tree cover and the field measurements, we dropped out the data within the $\mathrm{1}/\mathrm{1200}$^∘ resolution grids in which the standard deviation of tree cover was greater than 15 % (according to Chang et al., 2021), leaving 2444 homogeneous forest plots remaining (see Fig. 2 for the spatial distribution of these forest plots and Fig. S1a–b in the Supplement for the cumulative frequency curve and histogram of the AGBC records). The AGBC records in these forest plots were further multiplied by the mean fraction of forestland over 2011–2015 in the corresponding grid, which was computed from the annual 30 m resolution China land cover dataset (CLCD) (Yang and Huang, 2021). By comparison, GlobBiomass 2010 AGB matches the best with the grid-scale forest AGBC derived from plot measurements, with a correlation coefficient (CC) of 0.50, followed by tree cover (CC = 0.42), the product of canopy height and tree cover (CC = 0.38) and, finally, the canopy height (0.27) and Huang et al. (2019)'s AGB (0.25). Therefore, to obtain an improved benchmark map of forest AGBC in China for the period of 2011–2015, we chose the GlobBiomass 2010 dataset as our basis and calibrated it against the in situ observation-based, grid-scale forest AGBC. To build an equation for the calibration, we divided the grid-scale AGBC values into 16 equidistant subranges (0–15, 15–30, …, 225–240 tC ha⁻¹), calculated the median of grid-scale AGBC values that are within each subrange and then the median of GlobBiomass AGB values in the corresponding grids. According to previous studies, an exponential function would be suitable for calibrating the GlobBiomass map in a region such as China (Mialon et al., 2020). After the calibration, we averaged the benchmark AGBC map from $\mathrm{1}/\mathrm{1200}$ to $\mathrm{1}/\mathrm{120}$^∘ (approximately 1 km) to further reduce the uncertainties.

Figure 2 The spatial distribution of (1) 2444 homogeneous forest plots with aboveground ground biomass carbon stock measured between 2011 and 2015 and (2) 8182 forest plots with both above- and belowground biomass records collated in this study.

2.2 Temporally continuous forest AGBC mapping during 2002–2021

Because the benchmark AGBC was mapped based on SAR data, the spatial pattern accuracy is guaranteed, but the temporal coverage is limited to just a few years. Hence, to create a forest AGBC time series over the past 20 years, we integrated the benchmark AGBC with long-term continuous optical and passive microwave remote sensing data.

The spatial resolution of optical remote sensing is higher and is thus preferred in this study. By adopting the MODIS vegetation continuous field (VCF) data (MOD44B v061), which includes three ground cover components, we obtain the percent tree cover, percent non-tree vegetation (i.e. short vegetation) cover and percent non-vegetated area (DiMiceli et al., 2022). We first calculated the mean tree cover (hereinafter TC_mean) and short vegetation cover (hereinafter SVC_mean) during 2011–2015 and resampled them from 250 m to $\mathrm{1}/\mathrm{120}$^∘, which is the same resolution as the benchmark AGBC map for 2011–2015. Because the canopy heights of trees are usually similar within a small area, the regional AGBC per TC_mean can be assumed to be the same, which is referred to as the “homogeneous assumption” hereinafter. Accordingly, for each grid, we searched the TC_mean, SVC_mean and AGBC within a 3 × 3 window ($\mathrm{1}/\mathrm{40}$^∘ × 1/40^∘) and then regressed the AGBC values in nine grids against both TC_mean (the primary, or key, predictor of AGBC) and SVC_mean (assumed as a supplementary predictor) linearly. Specifically, when the regression coefficient of SVC_mean was negative or the fitting efficiency was low (R²<0.5; significance p value > 0.05), we excluded the supplementary predictor from the regression and only explored the linear relationship between TC_mean and AGBC. Afterwards, if the regression between TC_mean and AGBC was still invalid, we enlarged the searching window size to 5 × 5, then 7 × 7 and, finally, 9 × 9 until the regression and the coefficients became valid. Then, the grid annual AGBC from 2002 to 2021 can be estimated from the TC or both the TC and SVC in each year, following the regression results. If the regression failed even if the window size reached 9 × 9, we stopped expanding the searching window to avoid the homogeneous assumption being invalid. In those grids, following a previous study (Xu et al., 2021), we divided the estimated AGBC by the TC_mean during 2011–2015 and then multiplied the TC in each year to obtain the AGBC time series. The above method utilized spatial information to estimate the temporal variation, and can thus be referred to as the “space-for-time” method.

Long-term continuous microwave VOD can also reflect forest biomass changes, although the relationship was nonlinear (Jackson and Schmugge, 1991; O'Neill et al., 2021; Liu et al., 2015; Wigneron et al., 1995). We selected the global land parameter data record (LPDR) v3, 0.25^∘ resolution, VOD product, which was generated using similar calibrated, X-band brightness temperature retrieved from the Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E) and the Advanced Microwave Scanning Radiometer 2 (AMSR2; Du et al., 2017). As revealed by a recent evaluation study, LPDR VOD is better correlated with AGB than other long-term VOD products, especially in less-vegetated areas (Li et al., 2021). Because X-band VODs are still more sensitive to canopy cover than stem biomass, and there is a data gap between October 2010 and June 2011, while the plot investigations were all conducted in summer (Tang et al., 2018). We averaged the VOD data from mid-July (the 206th day) until the end of September (the 274th day) in each year to represent the annual AGB status. We also aggregated the benchmark AGBC map and the VCF data (TC_mean and SVC_mean) to 0.25^∘ resolution. After each round of searching, we applied the shape language modelling algorithm (D'Errico, 2022), to fit the nonlinear but monotonous relationship between AGBC and VOD values within the searching window, and then fitted the bivariate linear regression between AGBC and VCF. If the nonlinear regression between AGBC and VOD is valid, and the R² is superior to the regression between AGBC and VCF data, then the LPDR VOD data are expected to outperform VCF in predicting the interannual AGBC changes in the corresponding 0.25^∘ grid. Therefore, in these areas, we calibrated the VCF-derived high-resolution ($\mathrm{1}/\mathrm{120}$^∘) annual AGBC by incorporating the ratio between the VOD-derived 0.25^∘ AGBC and the aggregated VCF-derived AGBC in that year.

2.3 Forest belowground biomass carbon (BGBC) mapping during 2002–2021

This study mapped belowground forest biomass carbon (BGBC) following the random forest (RF) model approach (Huang et al., 2021). To reveal the forest above- and belowground biomass allocation rules in China, this study collated both AGB and BGB records at 8729 forest plots throughout China, which were obtained using allometric equations or clear-cutting methods from published papers, including Luo (1996), Luo et al. (2014), Guo and Ren (2014) and Wang et al. (2014). Because forest stand age and tree species (forest type) information are also available at 8182 plots, while the climatic backgrounds are available from the WorldClim v2.1 dataset (Fick and Hijmans, 2017), the forest plot AGB, forest type (hereinafter FOR_T), stand age, mean annual temperature (MAT), temperature seasonality (standard deviation of monthly temperature × 100; abbreviated as Tsea), mean annual precipitation (MAP) and precipitation seasonality (coefficient of variation in monthly precipitation; Psea) were applied as predictors of the forest plot BGB. For simplicity, we distinguished all forests into the following five types: evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), evergreen needleleaf forest (ENF), deciduous needleleaf forest (DNF), and mixed forest (MF). Using the data records at these 8182 plots (see Fig. 2 for the locations of these forest plots and Fig. S1c–f for the cumulative frequency curves and histograms of the AGB and BGB data), we trained 10-fold RF models using MATLAB R2021a. The number of regression trees was set to 500.

Because the $\mathrm{1}/\mathrm{120}$^∘ resolution grids in which forest AGBC data were available are often mixed with forestland and some other land cover types, e.g. waterbodies, bare ground and croplands, we converted the annual grid-average AGBC into the AGBC per area forestland by incorporating the annual fraction of forestland computed from the CLCD at 30 m resolution. Considering the potential uncertainties in the forestland fraction and the inclusion of shrub or herbaceous plant AGB in the SAR-derived AGB, we only calculated the annual AGBC per area forestland in grids that were dominated by forestland (forestland fractions were consistently over 50 %). In these forestland grids, we simulated the forest BGBC per area forestland during 2002–2021 by inputting the estimated annual AGB (approximately 2 times the AGBC) per forestland, annual forest type map derived from the ESA Climate Change Initiative's land cover classification dataset (Li et al., 2018), forest stand age (Besnard et al., 2021) and climatic background variables into the RF model. Afterwards, we multiplied the simulation results in every forestland grid with the annual forestland fractions to obtain the forest BGB and BGBC (0.5 × BGB) time series. Finally, for grids with forests, but which are not dominated by forestlands, we sequentially searched for at least five valid RSR values (the ratio of forest BGBC to AGBC) nearby (Y. Chen et al., 2019) and then multiplied the annual forest AGBC in the grid with the median of nearby RSR values in each year to estimate the annual forest BGBC.

2.4 Evaluation and assessment

We compared the interannual trend of forest biomass carbon calculated in this study against that of existing global/regional long-term woody biomass datasets, including the well-received global long-term terrestrial biomass data between 1993–2012, which were developed mainly based on a long-term integrated VOD dataset (Liu et al., 2015), and an updated woody biomass dataset covering 2001–2019, whose long time series was derived from optical remote sensing data (i.e. MODIS VCF dataset; Xu et al., 2021).

To justify the random forest models for BGBC predictions, we drew partial dependence plots (PDPs) in MATLAB R2021a to show the marginal effect that one predictor has on the training target (e.g. BGB at forest plots; Hastie et al., 2009). Here, for each predictor, we excluded the extreme values (the lowest 1 % and the highest 1 %) before calculating the corresponding PDP to avoid roughly extending the PDP lines to data-scarce areas. We performed 10-fold RF training to derive the mean PDP values and the standard deviations.

3 Results and discussion

3.1 Evaluation of forest AGBC and BGBC estimation

First, according to Fig. 3a, an exponential function, y = 1.63 × x^0.73, can fit the relationship between the actual grid-scale forest AGBC over 2011–2015 (y) and the AGB values predicted by GlobBiomass 2010 (x). Hence, this function was applied to derive the benchmark map of forest AGBC across China.

Figure 3 Evaluation of the forest AGBC and BGBC mapping in this study. (a) The regression relationship between the grid-scale forest AGB carbon stock derived from plot measurements during 2011–2015 and the GlobBiomass AGB dataset for 2010. (b) The minimum searching window sizes of every $\mathrm{1}/\mathrm{120}$^∘ grid when the spatial variation in MODIS VCF was applied as the predictor of AGBC changes. (c) The spatial pattern of the relative performances of MODIS VCF and LPDR VOD data in predicting the variation in AGBC. (d) Comparison of the mean tree cover between the grids where VOD data were more suitable for predicting the variation in AGBC and the grids where VCF data were the better predictor.

Second, when using the spatial information of tree cover and short vegetation cover to estimate the temporal variation in AGBC in each grid, the spatial searching window was at its minimum of 3 × 3 in most (53 %) grids with forests. Across China, the temporal extension of AGBC in only 15 % of all grids with forest cannot be achieved, even when the searching window was enlarged to 9 × 9 (Fig. 3b).

Next, as shown in Fig. 3c and d, the grids in which the LPDR X-band VOD performed better than MODIS VCF in predicting the temporal change in forest AGBC are usually located in regions with low tree cover. These grids account for just 10.4 % of all grids with forests and may suffer from high uncertainty within the optical-based variation in tree cover. Therefore, microwave-based VOD is supposed to be more suitable for estimating the forest AGBC changes in these regions.

The RF model designed for forest plot BGB estimation (see Sect. 2.3) achieved a predictive R² of 0.89 ± 0.02, while the RMSE was 6.3 ± 0.5 t ha⁻¹. AGB explained 53 % of the variation in BGB among different plots. Long-term climate backgrounds, i.e. mean annual temperature, temperature seasonality, annual precipitation and precipitation seasonality accounted for 8 %, 6 %, 8 % and 7 %, respectively. Forest type and stand age also contributed 12 % and 8 % to the training efficiency, indicating that the effects of these factors are nonnegligible. The selection of predictors of BGB basically followed the existing knowledge (Huang et al., 2021), and the seasonality of temperature and precipitation made sense in the prediction (see Text S1 in the Supplement). On the other hand, although previous studies incorporated many edaphic factors as predictors of BGB (Huang et al., 2021), by comparing the training efficiencies with whether these edaphic factors are incorporated or not, we could justify the reasonability of our simplified set of predictors (Text S1).

According to the collected woody plot data, AGB is a key driver of BGB (Fig. 4). Yet, RSR changes among different forest growth stages, decreasing in general, as previously reported (Mokany et al., 2006). The overall negative impact of mean temperature on BGB or RSR agrees with the mechanism that higher heat promotes nutrient accessibility (Luo et al., 2012; Ma et al., 2021) and increases the turnover rates of roots at a higher magnitude than stems (Reich et al., 2014). The U-shaped relationship between precipitation and belowground biomass allocation follows the optimal biomass allocation theory because arid climates promote root extension, yet too heavy rainfall reduces nutrient availability through leaching and dilution effects (Luo et al., 2012). Other factors, including temperature seasonality, precipitation seasonality and forest type, have supplementary effects on biomass allocation (Fig. S2).

Figure 4 Influence of key factors on forest belowground biomass (BGB) and root–shoot ratio (RSR) in China. Subfigures (a–d) show partial influences of (a) AGB, (b) stand age, (c) MAT and (d) MAP on BGB and RSR values of all forest plots. The error bars represent the standard deviations of the 10-fold training. We did not draw the PDP for the impact of AGB on RSR, since the dividend of RSR calculation is AGB.

3.2 Forest biomass carbon pool, allocation and change in China

Between 2002 to 2021, the forest above- and belowground biomass carbon (AGBC and BGBC) pools in China were 8.6 ± 0.6 and 2.2 ± 0.1 PgC, respectively (Table 1). The mean RSR for all forests was 0.25, basically equal to the global average (Huang et al., 2021). Separated by forest type, evergreen conifer forests (ENFs) occupy the highest biomass carbon pool per unit area, mainly because ENFs are mainly located in southwestern China and are more mature and natural (Yu et al., 2020; Zhang et al., 2017). Deciduous forests (DBF and DNF) in northern China (see Fig. S3 for the distribution of different forest ecosystems) harbour less biomass carbon but higher BGBC (Fig. 5a), which can be attributed to the higher RSR values (Table 1).

The forest biomass carbon stock in China increased at an average rate of 114.5 ± 16.3 TgC yr⁻¹ (p<0.01) during 2002–2021, and the annual biomass carbon gains were the greatest from 2014 to 2015, reaching 736 TgC (Fig. 5b). Changes in AGB and BGB accounted for 81.9 % and 18.1 %, respectively, of the forest carbon stock gains over the past 20 years.

Figure 5 Forest biomass allocation and biomass change in China during 2002–2021. (a) Aboveground biomass carbon (AGBC) and belowground biomass carbon (BGBC) density of different forest ecosystems in China. (b) The interannual changes in forest AGBC and BGBC in China. Total forest biomass carbon stock changes from the previous to the current year are represented by green columns.

Our estimates of the forest biomass carbon pool, forest RSR and the recent interannual trend of forest biomass carbon are generally consistent with previous estimates based on massive field investigations (Table 1).

Table 1 Agreement of the estimated various forest RSR and the trend of forest biomass carbon in China with existing studies.

3.3 Spatial pattern of the forest biomass carbon stock trend in China

The highest forest biomass carbon pools during 2002–2021 were observed in northeastern and southwestern China, especially South Tibet. Forest biomass carbon stocks were also high in the natural or semi-natural forests in the Qinling mountains, Hengduan Mountains, Hainan and Taiwan (Fig. 6a). Above- and belowground forest biomass allocation varies significantly among regions. RSR is highest in northeastern deciduous conifer forests and northern China's deciduous broadleaf forests but low in southern China (Fig. 6b). The strongest forest biomass carbon increases were found in central to southern China, including the Loess Plateau, Qinling mountains, southwestern karst region and southeastern forests. Slight declines in forest biomass carbon only occurred in some mature and natural forests, e.g. those in the Greater Khingan mountain, Hengduan Mountains and South Tibet (Fig. 6c). A total of 40.3 % of all forests in China showed significant biomass carbon stock gains over the past 20 years, whereas only 3.3 % of forests experienced significant biomass carbon losses (Fig. 6d).

Figure 6 Maps of forest biomass carbon pool, allocation and trend in China during 2002–2021. (a) Spatial pattern of the forest biomass carbon pool in China. (b) RSR of all forestland pixels. (c) Map of the forest biomass carbon stock trend from 2002 to 2021, with shaded areas representing statistically significant trends at the 95 % confidence level. (d) histogram and basic statistics of all forest biomass carbon stock trend.

4 Data availability

Annual forest above- and belowground biomass maps in China between 2002 and 2021 are now available at https://doi.org/10.6084/m9.figshare.21931161.v1 (Chen, 2023). This dataset is also available on PANGAEA (https://doi.org/10.1594/PANGAEA.955074, Chen et al., 2023) and will also be available on the National Tibetan Plateau/Third Pole Environment website. Other open datasets that made this research possible, and the related references, are attached in the Supplement (Text S2).

5 Discussion

5.1 Comparison of the estimated forest biomass carbon pool change in this study against the existing datasets

Although there is potential overestimation, the interannual variation in forest AGBC in China, according to Liu et al. (2015), and that of total biomass carbon, according to Xu et al. (2021), are both highly correlated with our results (R²=0.65 and 0.88). Liu et al. (2015) predicted a forest AGBC increase rate of 102.2 ± 35.8 Tg yr⁻¹ (p<0.01), which is slightly higher than our estimate of 80.8 ± 25.1 Tg yr⁻¹ during 2002–2012, while Xu et al. (2021) indicated a biomass carbon stock trend of 99.4 ± 23.2 Tg yr⁻¹ (p<0.01) from 2002 to 2019, which is slightly lower than the rate of 115.6 ± 20.2 Tg yr⁻¹ in this study (Fig. 7a–b). The spatial maps of the forest biomass carbon trends estimated by Xu et al. (2021) and Liu et al. (2015) were slightly patchy (Fig. 7c–d). Compared to this study, the two existing datasets (i.e. the datasets of Liu et al., 2015, and Xu et al., 2021) predicted higher biomass carbon stock trends in the Qinling mountains and the mature deciduous conifer forests in northeastern China. Meanwhile, they predicted lower carbon sinks in southern China (Fig. 7c–f), where reforestation and forest-management-induced short-term extensive carbon uptake (Tong et al., 2020) have been confirmed by atmospheric inversions (Wang et al., 2020; Yang et al., 2021).

Figure 7 Comparison of the estimated forest biomass carbon pool change in this study with two existing datasets. (a) Comparison of the interannual variation in forest biomass carbon in this study with the estimate by Xu et al. (2021) during 2002–2019. (b) Comparison of the interannual variation in forest AGBC calculated in this study with the estimate by Liu et al. (2015) over 2002–2012. (c) Map of the interannual trend of forest biomass carbon stock in China during 2003–2019, according to Xu et al. (2021). (d) Map of the forest AGBC trend in China during 2003–2012, according to Liu et al. (2015). (e) Map of the estimated trend of forest biomass carbon stock over 2002–2019 in this study. (f) Map of the estimated forest AGBC trend over 2002–2012 in this study.

Finally, by comparing Fig. 7e and f, we could also notice that the hot spot of forest biomass carbon gains has moved from the Loess Plateau over the first decade of our study period (2002–2012) to southern China (e.g. Guangxi province) later. This change was probably due to the large-scale implementation of the “Grain for Green” project on the Loess Plateau (Liu et al., 2020; Wu et al., 2019) before 2012 and the massive plantation of fast-growing trees in southern China after 2010 (Tong et al., 2020).

5.2 Some uncertainties in the forest biomass carbon dataset and future prospects

During benchmark AGBC mapping, we converted the in situ AGBC data at forest plots into the grid-scale average AGBC by multiplying by the fraction of forestland during the time period of field investigation. Considering the overall high quality of the China's land use/cover datasets developed via the human–computer interactive interpretation of Landsat images (Liu et al., 2014; Yang and Huang, 2021), and that the producer's accuracy (PA) and user's accuracy (UA) for forestland classification in the CLCD dataset used in this study were 73 % and 85 %, respectively, the errors within the benchmark AGBC mapping induced by the scale conversion based on the forestland area fraction were generally limited.

The variation in climatic conditions in the short term may have subtle influences on that in the BGB, but explicit knowledge on this effect is lacking. Instead, woody vegetation BGB is much more driven by AGB (vegetation density), as indicated by the very strong relationship between BGB and AGB (R²≥0.85). Moreover, the long-term climatic background is expected to have a stronger influence on the RSR of perennial woody plants than the meteorological conditions in only a few years, since above- and belowground biomass allocation is the result of plants' long-term adjustment to the environment (Qi et al., 2019). Therefore, it is reasonable not to consider the influence of the specific climatic conditions in a year on the variation in BGB.

In the near future, P-band microwave sensors, which have higher penetrability into the canopy than L-band microwaves, will further improve AGB mapping. For example, Biomass (Biomass Monitoring Mission for Carbon Assessment), a fully polarimetric P-band SAR, is scheduled to be launched in 2024 (Le Toan et al., 2011). Therefore, in the future, the relationship between P-band microwave retrievals and biomass should be addressed, in addition to the calibration of historical AGB datasets (e.g. the long-term AGB dataset in this study) against the P-band SAR-based AGB benchmark map to extend the time series. In addition, an intercalibration between the AMSR-E-based VOD and the AMSR2-based VOD will further reduce the potential bias within the long-term integrated VOD datasets (Wang et al., 2021a, b). On the other hand, more in situ AGB and BGB measurements in larger plots are needed to further improve the estimation of belowground biomass allocation.