Dynamics of China's forest carbon storage: the first 30 m annual aboveground biomass mapping from 1985 to 2023

Accurate estimation and monitoring of forest aboveground biomass (AGB) are essential for understanding carbon dynamics, managing forest resources, and guiding environmental policies. However, the spatial and temporal patterns, dynamics, and driving factors of forest AGB in China over recent decades remain insufficiently understood, hindering ecosystem analysis and forest management strategies. This study combines multi-source remote sensing data with residual neural networks (ResNets) to develop the first 30 m resolution annual China Forest AGB dataset (1985–2023) with uncertainty quantification. Validation results confirm the robustness of the ResNets model, achieving an R² of 0.91, RMSE of 16.49 Mg ha⁻¹, and Bias of 0.50 Mg ha⁻¹ against GEDI footprint AGBD, and an R² of 0.63, RMSE of 68.26 Mg ha⁻¹, and Bias of −19.87 Mg ha⁻¹ against independent multi-year ground survey data. The dataset reveals a notable increase in China's average forest aboveground biomass density (AGBD) from 95.74 ± 11.30 Mg ha⁻¹ in 1985 to 122.69 ± 13.94 Mg ha⁻¹ in 2023. During this period, total forest aboveground carbon (AGC) stock rose from 5.50 ± 0.23 PgC to 13.97 ± 0.87 PgC, establishing China's forests as a significant carbon sink over the past four decades, with a net carbon sink of 0.22 ± 0.01 PgC yr⁻¹, offsetting 11.5 %–14.9 % of China's annual fossil fuel and industrial carbon dioxide emissions. Forest growth contributed 65.1 % (5.75 PgC) of the total AGC increase, while forest expansion accounted for 34.9 % (3.09 PgC). This dataset provides critical information for forest carbon accounting in China and offers valuable insights for climate change mitigation, ecosystem conservation, and sustainable land management.

The dataset is available at https://doi.org/10.5281/zenodo.12620984 (Part I: 1985–1993, Cai et al., 2025a), https://doi.org/10.5281/zenodo.12637101 (Part II: 1994–2001, Cai et al., 2025b), https://doi.org/10.5281/zenodo.12655492 (Part III: 2002–2008, Cai et al., 2025c), https://doi.org/10.5281/zenodo.12658255 (Part IV: 2009–2015, Cai et al., 2025d), https://doi.org/10.5281/zenodo.12742210 (Part V: 2016–2021, Cai et al., 2025e), https://doi.org/10.5281/zenodo.12747329 (Part VI: 2022-2023, Cai et al., 2025f).

1 Introduction

Aboveground biomass (AGB) is a critical metric for quantifying carbon stocks in forest ecosystems, representing approximately 70 %–90 % of total forest biomass (Cairns et al., 1997). Variations in AGB over time directly shape the dynamics of the global carbon cycle. Empirical studies estimate that forest AGB accounts for approximately 15 % to 30 % of the terrestrial carbon pool (Davidson and Janssens, 2006) and has sequestered an estimated 106.9 PgC over the past three decades, corresponding to roughly 17 % of global fossil fuel emissions (Pan et al., 2024). Given its critical role in the global carbon cycle, AGB has been incorporated into the compliance framework of the Paris Agreement and recognized as a key monitoring indicator under the United Nations Sustainable Development Goals (SDGs) (Baccini et al., 2012; Xu et al., 2021).

China ranks fifth globally in total forest area and possesses the largest expanse of planted forests worldwide (FAO, 2020) (see Sect. S1 in the Supplement for details on the study area). Since the 1970s, the Chinese government has launched a series of large-scale forestry initiatives aimed at improving ecological quality and mitigating environmental degradation (Liu et al., 2008; Lu et al., 2018). These programs have driven extensive afforestation and reforestation, expanding forest cover from 12.98 % in the 1980s to 24.02 % in 2021 (Fig. S1) (Cai et al., 2024). Between 1990 and 2019, the carbon sink capacity of global temperate forests increased by 30.2 %, with China's reforestation programs contributing substantially to this rise (Pan et al., 2024). Nevertheless, the degree to which these afforestation initiatives have augmented AGB, especially regarding the spatiotemporal dynamics of carbon accumulation across China's forested landscapes, remains inadequately quantified.

Annual, spatially explicit assessments of forest AGB provide a direct and effective approach to evaluating carbon sink enhancement resulting from large-scale afforestation. While mature forests typically experience gradual annual biomass changes that may not be detectable using yearly remote sensing data, monitoring intervals of 3–5 years are generally sufficient to capture meaningful growth trends (Le Toan et al., 2011). In contrast, young, fast-growing forests or recently disturbed areas require more frequent (e.g., annual) observations due to rapid carbon turnover in dynamic carbon pools, rendering high-frequency monitoring essential for accurately capturing rapid disturbance events and recovery processes (Quegan et al., 2019).

Traditional forest inventory approaches, which rely on permanent sample plots and periodic surveys (typically conducted at five-year intervals), remain widely employed at the national scale for estimating forest biomass and assessing its temporal dynamics. While these surveys provide valuable statistics, they are labour-intensive, lack spatial and temporal continuity, and are unsuitable for producing wall-to-wall, spatiotemporally consistent biomass estimations. Process-based ecosystem models (e.g., Dynamic Global Vegetation Model, DGVM), which simulate biogeochemical processes such as photosynthesis, carbon allocation, and uptake, often incorporate biotic, climatic, and anthropogenic drivers (Sitch et al., 2008). However, these models are limited by uncertainties in input datasets (e.g., climate and soil data), coarse spatial resolutions, and assumptions of stand homogeneity, resulting in significant uncertainties in biomass estimates and an inability to capture spatial variability in forest biomass. In contrast, remote sensing–based models is a promising alternative for achieving spatially explicit and annual assessments of forest AGB (Lu et al., 2016).

Optical remote sensing, synthetic aperture radar (SAR), and light detection and ranging (LiDAR) are the primary remote sensing technologies employed for estimating forest AGB. Spaceborne LiDAR systems (e.g., GEDI) provide accurate footprint-level forest AGB estimates but are limited in spatial and temporal coverage, making them unsuitable for nationwide, long-term monitoring (Coops et al., 2021; H. Nguyen et al., 2019). While synthetic aperture radar (SAR) data enable spatially continuous monitoring, historical satellite missions like ALOS-2 and Sentinel-1 lack sufficient temporal depth to capture the multi-decade impacts of large-scale afforestation. Vegetation optical depth (VOD) derived from SAR has demonstrated strong correlations with AGB over long timescales (Fang et al., 2024; Moesinger et al., 2020) but suffers from coarse spatial resolution, limiting its effectiveness in capturing fine-scale AGB changes caused by deforestation, degradation, and recovery.

To date, large-scale and long-term forest biomass monitoring systems have primarily relied on multispectral satellite imagery (e.g., Landsat time series) combined with field inventory data or airborne LiDAR to estimate AGB (Coops et al., 2021). Three main approaches have been employed for historical AGB reconstruction: (1) Space-for-time substitution: Models are trained using multitemporal reference datasets to predict historical AGB from spectral data, assuming that contemporary relationships are temporally transferable (Harris et al., 2021; Matasci et al., 2018; Wulder et al., 2020). (2) Parameter-based reconstruction: Historical AGB is inferred from related forest parameters, such as tree cover or canopy height (Chen et al., 2023). This approach depends on the robustness of parameter–AGB relationships and the quality of input datasets. (3) Continuous time-series estimation using machine learning: Continuous AGB dynamics are derived by integrating Landsat time series with advanced algorithms that extract temporal and spectral features (e.g., CCDC) and predict biomass via machine learning models, such as random forests (Fu et al., 2024). Among the available approaches, space-for-time substitution provides a practical solution for reconstructing historical AGB at large scales. By leveraging multitemporal reference data, this method can generate temporally generalizable models without requiring continuous field observations, thus balancing prediction accuracy, temporal coverage, and computational feasibility. Compared with approaches relying on continuous Landsat time series and machine learning, space-for-time substitution effectively captures both abrupt and gradual forest changes while maintaining manageable computational demands (H. Nguyen et al., 2019). However, a key challenge in predicting AGB from Landsat spectra is spectral saturation, wherein reflectance values become insensitive to biomass variations beyond certain thresholds, particularly in dense forest conditions (Zhang et al., 2019). This phenomenon restricts the capacity of traditional regression-based models to accurately capture biomass dynamics in high-biomass regions. Previous studies indicate that deep learning approaches outperform traditional machine learning algorithms in terms of cross-validation stability and generalization (Dixon et al., 2025; He et al., 2016; Huy et al., 2022). In particular, the ResNet architecture – through its residual connections – facilitates the training of deeper networks and the extraction of complex, non-linear relationships between spectral–structural features and AGB (He et al., 2016). This capability offers a promising pathway to mitigating the saturation effects inherent to single-source data.

Another limitation to achieving high-quality, high-resolution, and long-term AGB monitoring in China is the absence of forest-specific thematic maps that reliably reflect national forest dynamics. Forest extents derived from general land-use/land-cover products frequently underrepresent afforestation and plantation dynamics (Cai et al., 2024), which leads to omissions in historical records and underestimation of China's carbon sink potential (Wang et al., 2020; Yu et al., 2022). There is also a pronounced shortage of large-scale, geographically representative biomass survey data that cover China's diverse forest types, soils, climatic zones, and biomass gradients – datasets that are essential for developing spatially and temporally transferable models. Although several static AGB maps for China have been produced using different approaches, substantial inconsistencies persist (see Table S1); reported national AGB/AGC totals vary widely because of differences in observation epochs, input data, and upscaling/conversion procedures (Santoro et al., 2021). A recurring problem is temporal and spatial misalignment between field inventories and remote sensing inputs, which reduces the ability of models to capture forest change reliably (Araza et al., 2022; Babcock et al., 2016; Le Toan et al., 2011).

Representative global and China time series AGB products (see Sect. S1) exhibit complementary strengths but also systematic weaknesses that limit their direct applicability for high-fidelity, temporally explicit mapping in China. For example, coarse global products and compilations (e.g., European Space Agency's Climate Change Initiative (ESA CCI) biomass dataset, and GEOCARBON/GlobBiomass) provide extensive spatial coverage and useful baselines but tend to smooth local detail because of coarse native resolutions and multi-year compositing; some pan-tropical maps also exhibit systematic positive or negative biases outside the wet tropics (Liu et al., 2015; Santoro et al., 2021). Meanwhile, emerging spaceborne lidar (e.g., GEDI, and ICESat-2) provides higher-quality vertical structure information that can reduce saturation and structural ambiguity, but footprint sparsity requires careful upscaling to wall-to-wall maps (Duncanson et al., 2022). In short, available global products trade off spatial detail, temporal continuity, and structural sensitivity – none simultaneously deliver the high spatial resolution, dense temporal coverage, and unbiased response across the full biomass gradient that are needed for national-scale, temporally consistent AGB time series in China.

In response to the urgent need for improved national biomass mapping and the lack of a comprehensive dataset to support continuous monitoring of AGB spatial dynamics, this study proposes a bottom-up spatial framework to quantify changes in China's forest biomass over the past four decades (1985–2023). This framework integrates a ResNet-based deep learning algorithm, leveraging 72 150 multi-temporal GEDI AGB training samples (2019–2021) and remote sensing observations to train the model and generate accurate AGB time series, thereby reducing uncertainties in inventory change estimates. The primary goals of this study are: (1) to generate the first 30 m resolution annual time series dataset of China's forest AGB from 1985 to 2023, referred to as the China Forest AGB Time Series Dataset (CFATD); and (2) to monitor and quantify long-term and short-cycle changes in China's terrestrial forest biomass, improving estimates of carbon sources (e.g., from deforestation) and carbon sinks (e.g., from forest regeneration and afforestation). By achieving these objectives, this study seeks to answer four key questions: (1) What are the spatiotemporal trends of China's AGB carbon sink? (2) How much carbon have forests sequestered under China's large-scale afforestation initiatives? (3) How much of this sequestration is directly attributable to afforestation? (4) To what extent can forests contribute to China's carbon neutrality goals, and how much of the nation's anthropogenic carbon emissions have past forest carbon sinks offset? The resulting dataset provides a robust, data-driven foundation for understanding forest biomass changes and their implications for environmental policy, climate change mitigation, and sustainable forest management.

2 Materials and methods

2.1 Remote sensing data

2.1.1 Landsat imagery

This study utilized Landsat satellite imagery from the Google Earth Engine (GEE) platform, covering the period from 1985 to 2023. The dataset includes images from Landsat 4/5 Thematic Mapper (TM), Landsat 8/9 Operational Land Imager (OLI), and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) prior to 2003. Due to the Scan Line Corrector (SLC)-off failure in mid-2003, post-2003 Landsat 7 data were excluded from the analysis, except for 2012 when no other Landsat sensors provided coverage. All imagery was obtained from the Landsat Collection 2, Level-2, Tier-1 surface reflectance products at 30 m spatial resolution, which have undergone sensor-specific radiometric calibration, atmospheric correction, and orthorectification (Wulder et al., 2022).

To ensure temporal consistency, we first selected images acquired during the main growing season (Day of Year 150–300) and applied the cross-sensor spectral harmonization procedure described by Roy et al. (2016) to align the spectral responses of TM, ETM+, and OLI bands. Subsequently, annual medoid composites were generated for each year using QA band filtering (Zhu et al., 2015) to remove cloud, shadow, and snow contamination. Remaining gaps in the annual composites were filled using linear interpolation or extrapolation along the temporal dimension. Each composite contained three visible bands (red, green, blue), one near-infrared (NIR) band, and two shortwave infrared (SWIR) bands. Additionally, to capture vegetation growth status and canopy structure, we calculated several spectral indices, including the Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Normalized Burn Ratio (NBR). These spectral bands and indices served as input features for subsequent modelling. A more detailed description of the processing workflow can be found in Cai et al. (2025g).

2.1.2 China annual tree cover dataset

The China Annual Tree Cover Dataset (CATCD) is the first long-term dataset capturing annual tree cover changes across China at a 30 m spatial resolution from 1985 to 2023 (Cai et al., 2024). Developed using Landsat time-series imagery and advanced random forest-based ensemble learning techniques, the dataset has been rigorously validated against multiple reference datasets, achieving correlation coefficients of 0.70 to 0.96 and RMSE values ranging from 5.6 % to 25.2 %. The dataset aligns closely with National Forest Inventory (NFI) trends, effectively capturing forest area dynamics driven by major afforestation and reforestation initiatives. This capability enables the development of an AGB dataset that reflects temporal forest changes with high fidelity. In accordance with the technical standards of the National Forestry and Grassland Administration (National Forestry and Grassland Administration, 2019) and forest definition in the previous study (Fang et al., 2001; Piao et al., 2009), a 20 % tree cover threshold was applied to delineate forested areas, which formed the foundational layer for AGBD estimation. Furthermore, leveraging the strong correlation between canopy cover and forest biomass (Chen et al., 2023), tree cover was also used as an explanatory variable in AGBD modelling.

2.1.3 Meteorological and topographic data

To enhance the accuracy of AGBD predictions by capturing the environmental and spatial factors influencing forest biomass, we integrated climate, topography, and geographic location as explanatory variables in our AGBD prediction model (Su et al., 2016; Zhang et al., 2019). Climate variables were sourced from WorldClim 2.1, which provides global climate data for the period 1970–2000 at high spatial resolutions (∼ 1000 m). This dataset includes variables such as temperature, precipitation, and other bioclimatic factors. We selected four climate variables from WorldClim 2.1 that are relevant to forest growth: bio1 (annual mean temperature), bio4 (temperature seasonality), bio12 (annual mean precipitation), and bio15 (precipitation seasonality). To match the spatial resolution of other datasets such as Landsat, we resampled the climate data to 30 m resolution using bilinear interpolation. Topographic factors were derived from the SRTM DEM V3, which provides 30 m resolution elevation data. Based on the DEM data, we calculated slope and aspect using the GEE platform. To capture large-scale spatial patterns and residual spatial autocorrelation beyond environmental predictors, we also included latitude and longitude as auxiliary spatial variables, following Tobler's First Law of Geography. This practice has proven effective in improving model generalization and reducing spatial bias in regions with sparse training data (Yang and Huang, 2021).

2.2 Reference data

2.2.1 GEDI AGBD product

The Global Ecosystem Dynamics Investigation (GEDI) instrument aboard the International Space Station (ISS) provides footprint-level estimates of aboveground biomass density (AGBD) at ∼ 25 m spatial resolution, with footprints spaced ∼ 60 m along track. These estimates are derived from statistical models calibrated with field measurements and simulated waveforms, stratified by region and plant functional type to improve predictive performance (Duncanson et al., 2022; Kellner et al., 2023). The dataset also provides per-footprint uncertainty information, which enhances its applicability in underrepresented regions. We obtained GEDI L4A footprint-level AGBD observations for 2019-2021 and applied standard quality-control filters to retain only reliable measurements, excluding those from degraded orbits or steep slopes (Shendryk, 2022). To remove residual outliers, we employed the GEDI Level-4B (L4B) 1 km gridded product as a consistency benchmark. The L4B product generation implements a rigorous multi-tier screening process that goes beyond standard quality flags. It applies stricter beam sensitivity thresholds (> 0.98 in tropics) to ensure canopy penetration, filters out observations with degraded geolocation precision, and incorporates an iterative local outlier detection algorithm to remove spurious retrievals caused by atmospheric scattering (e.g., clouds or fog). Leveraging this robust baseline, we performed a consistency check by identifying and excluding L4A footprints that deviated significantly (> 2 standard deviations) from the high-quality local L4B grid mean. Finally, to minimize spatial autocorrelation, we applied a systematic spatial thinning strategy. Within each 1 km grid cell, we prioritized the retention of the single observation that exhibited optimal beam sensitivity and minimal deviation from the grid-level ensemble mean, thereby ensuring the selection of the most representative footprint for each 1 km cell. This rigorous process resulted in a curated dataset of 72 150 high-confidence samples distributed across China's major ecoregions (Fig. S2).

2.2.2 Field survey data

To independently assess the accuracy of forest AGBD estimates derived from remote sensing, we compiled ground-based measurements from 4243 forest plots across the study area from published sources (Fig. S3) (Avitabile et al., 2016; Luo et al., 2014; Usoltsev, 2010; Zhang et al., 2019). Given the scale mismatch between field plots and satellite pixels, strictly controlling data quality was critical. We applied a rigorous screening protocol: (1) Temporal Alignment: Survey years were aligned with satellite observation epochs to minimize errors from forest growth or disturbance; (2) Homogeneity Assessment: Plots located in highly heterogeneous landscapes or forest edges (identified using high-resolution imagery) were excluded to reduce the scale effect; and (3) Statistical Outlier Removal: Inconsistent records and potential outliers (e.g., biomass values exceeding biological potential for specific climate zones) were removed. After this quality control, 2109 representative plots were retained. These plots, surveyed between 1978 and 2008, cover diverse forest types (primary, secondary, and plantations) and dominant species such as Betula platyphylla, Cunninghamia lanceolata, and Pinus koraiensis.

2.2.3 Provincial statistics from national forest inventories

The National Forest Inventories (NFIs) provide provincial-level forest area and total forest volume data but lack biomass-specific information for direct validation. Here, we utilized provincial forest AGC statistics derived from Fang et al. (2001), based on forest volume data from five NFIs (3rd NFI: 1984–1988 to 7th NFI: 2004–2008). These AGC estimates were obtained by converting the NFI-reported forest volumes into biomass using continuous biomass expansion factors (BEFs) (Fang et al., 2005). The reliability of the AGC estimates has been validated in previous studies (Fang et al., 2001).

2.3 Satellite-based datasets for temporal trend validation

To cross-check the temporal dynamics of AGBD, we employed two independent datasets: vegetation optical depth (VOD) and gross primary production (GPP). VOD primarily reflects canopy structural and water characteristics, whereas GPP quantifies ecosystem carbon uptake; both are strongly and positively correlated with aboveground biomass dynamics. VOD was derived from the VOD Climate Archive (VODCA v2, CXKu band), which provides harmonized microwave retrievals from 1987 to 2021 at 0.25° resolution (Zotta et al., 2024). To better capture biomass-related signals, we filtered daily VOD by the main growing season (day of year 150–300) and composited annual medians. For ecosystem productivity, we used the Global Sunlit and Shaded GPP dataset (1992–2020, 0.05° resolution), which estimates photosynthesis with a two-leaf light use efficiency model. Annual GPP was adopted to evaluate AGBD variations at interannual to decadal scales (Bi et al., 2022). Together, VOD and GPP provide complementary perspectives on canopy structure, water content, and carbon uptake, supporting the validation of long-term AGBD trends.

2.4 Forest AGBD mapping using ResNet

2.4.1 Model development

To improve computational efficiency and remove redundant information, we applied Recursive Feature Elimination (RFE) combined with five-fold cross-validation to identify the optimal set of features for forest AGBD mapping (Cai et al., 2020). The importance of the selected explanatory variables was further evaluated using the Random Forest method. Subsequently, these RFE-selected features were utilized as explanatory variables in developing the forest AGBD estimation model based on the Residual Neural Network (ResNet).

Figure 1 illustrates the architecture of the ResNet employed for AGBD estimation. Building on the original ResNet framework (He et al., 2016), all 2D convolutional layers were replaced with 1D counterparts to accommodate the structure of the per-pixel input feature. Specifically, the input was constructed as a 1D vector comprising spectral bands, indices, and auxiliary variables. The network comprises two residual blocks, each containing a 1D convolutional layer, batch normalization, and a Rectified Linear Unit (ReLU) activation. Each block incorporates a skip connection that adds the block's input directly to its output feature map. This residual connection stabilizes training, alleviates the vanishing gradient problem, and enables the model to learn complex, non-linear relationships between multitemporal spectral–structural features and biomass. By fusing shallow and deep feature representations, the residual mechanism also enhances the model's ability to capture subtle spectral variations in high-biomass regions, thereby helping to mitigate the saturation effects inherent to single-source optical data. Following each residual block, max-pooling layers reduce the dimensionality of intermediate feature maps by half, decreasing memory consumption, accelerating inference, and gradually enlarging the receptive field so that deeper layers encode broader contextual information. To further improve generalization, a dropout layer with a rate of 0.5 is applied before the final fully connected component, which consists of two linear layers separated by a ReLU activation.

Figure 1 Architecture of the ResNet used for AGBD estimation.

We randomly selected 80 % of the multi-temporal GEDI AGBD samples (n= 57 720) as the training set, while the remaining 20 % (n= 14 430) were reserved for multi-year validation to test the robustness of the model over time. The training data were aligned with the observation dates of Landsat satellites to develop a temporally consistent deep learning model, which was then applied to the entire time series of images to generate the forest biomass time series dataset. During model training, we adopted the Adam optimizer, an adaptive version of stochastic gradient descent. In each iteration of the optimizer, the loss was calculated based on random batches of 64 data samples (batch size = 64). Based on this loss, the partial derivatives (gradients) of all model parameters were determined by backpropagation, using the chain rule to compute numerical gradients. The model parameters were updated in small steps along the negative gradient direction to minimize the loss, with the step size controlled by a hyperparameter known as the learning rate, initialized to 0.001. The model was trained for 1,000 epochs, where an epoch refers to the entire training data set going through a complete process of forward propagation and back propagation.

2.4.2 Uncertainty estimation

In this study, we employed quantile regression to quantify the uncertainty of AGBD estimates (Chung and Neiswanger, 2021; Shendryk, 2022). Assuming a Gaussian error distribution, the 95 % prediction interval was approximately within ±2 standard deviations (SD), and the 68 % prediction interval was approximately within ±1 SD of the AGBD estimates. Therefore, we estimated the AGBD values at the 2.5th percentile and 97.5th percentile for each pixel using the ResNet model with quantile loss function (Eq. 1). The difference between the quantile loss estimates at these percentiles provided the 95 % prediction interval, allowing us to quantify the uncertainty (±1 SD) of AGBD estimates for each pixel. The uncertainty derived from the quantile regression can capture the variability in predictions associated with the model, reflecting modelling uncertainty derived from the spread of prediction intervals.

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{ll}{\mathrm{L}}_{\mathrm{q}}(y,\widehat{y})& =q(y-\widehat{y}{)}_{+}+(\mathrm{1}-q\left)\right(\widehat{y}-y{)}_{+}\\ & (\cdot {)}_{+}=max(\mathrm{0},\cdot )\end{array}\right.\end{array}$$

where L_q represents the quantile loss function, q is the number of quantiles of interest, y and $\widehat{y}$ signify the observation and prediction, respectively. In this study, q was set to 2.5 % and 97.5 % to obtain AGBD estimates at the 2.5th and 97.5th percentiles, respectively.

2.5 Accuracy assessment

To assess temporal robustness, we implemented a Leave-One-Year-Out (LOYO) validation using the multi-year GEDI AGBD record (2019–2021) (Laffitte et al., 2025). In each LOYO iteration (i.e., leaving one year out for validation), the model was trained on the remaining two years, and the process was repeated three times to evaluate temporal generalization. In addition, a 5-fold cross-validation was performed within the full dataset to further assess model stability and ensure consistent predictive performance. Reported accuracy metrics represent the mean performance across LOYO and 5-fold evaluations. For assessing the accuracy of the AGBD estimation model, we utilized validation metrics such as the coefficient of determination (R², Eq. 2), root mean square error (RMSE, Eq. 3), and bias (Eq. 4). These metrics quantified the agreement between predicted AGBD values and reference values, providing insights into prediction accuracy and systematic errors.

$$\begin{array}{}\text{(2)}& {\displaystyle }{R}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{(n-\mathrm{1})\sum _{i=\mathrm{1}}^n({\widehat{y}}_i-y_i{)}^{\mathrm{2}}}{(n-\mathrm{2})\sum _{i=\mathrm{1}}^n({\widehat{y}}_i-y_i{)}^{\mathrm{2}}}}\text{(3)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\sum _{i=\mathrm{1}}^n({\widehat{y}}_i-y_i{)}^{\mathrm{2}}}{n}}}\text{(4)}& {\displaystyle }\mathrm{Bias}={\displaystyle \frac{\sum _{i=\mathrm{1}}^n({\widehat{y}}_i-y_i)}{n}}\end{array}$$

where n is the number of reference samples, ${\widehat{y}}_i$ represents the estimated AGBD, and y_i denotes the reference AGBD.

2.6 Comparison of existing AGB products

To further validate the effectiveness of our models, we compared the results against several external AGB datasets, focusing on both spatial and temporal dynamics (see Sect. S2 for details). For spatial comparison, we used static AGB datasets from global (e.g., Zarin map (Zarin et al., 2016), Santoro map (Santoro et al., 2021)) and regional (e.g., Hu map (Hu et al., (2016)), Yang map (Yang et al., 2023)) sources to assess the spatial distribution and accuracy of AGBD estimates at different resolutions (from 30 to 1000 m). For dynamic comparison, we utilized annual and multi-year datasets (e.g., Hengeveld map (Hengeveld et al., 2015), Liu map (Liu et al., 2015), Chen map (Chen et al., 2023), ESA CCI) spanning different periods from 1950 to 2021. This allowed us to evaluate how well the models captured forest biomass changes over time across various temporal scales and spatial extents. By comparing static and dynamic AGBD estimates, we provide a comprehensive view of the performance and robustness of the developed models under different conditions.

2.7 Forest AGC stock and change analysis

The AGB (unit: Mg) for each pixel was calculated by multiplying the estimated AGBD by the forest area derived from the CATCD for the corresponding year (represented as the product of forest cover fraction and pixel area). Carbon content was then derived from the biomass using a conversion factor of 0.5, allowing for the calculation of forest AGC for each pixel. The net annual change in forest carbon was calculated by aggregating the pixel-level changes.

2.8 Quantifying the impact of forest changes on AGC stock

To consistently quantify how forest dynamics influence AGC stocks, we established clear operational definitions of forest change modes based on tree cover transitions. Annual forest cover data from the CATCD were used to determine forest status, with pixels having tree cover ≥ 20 % classified as forest. In this study, four distinct forest change modes were defined as follows:

• Forest growth. pixels that were forested in both 1985 and 2023 and exhibited an increase in tree cover.

• Tree cover loss. pixels that were forested in both 1985 and 2023 but showed a decrease in tree cover.

• Forest expansion. pixels that transitioned from non-forest in 1985 to forest in 2023.

• Forest loss. pixels that changed from forest in 1985 to non-forest in 2023.

For each pixel i, the change in aboveground carbon was calculated as:

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }{\mathrm{AGC}}_i={\mathrm{AGC}}_{i,\mathrm{2023}}-{\mathrm{AGC}}_{i,\mathrm{1985}}\end{array}$$

Positive ΔAGC values indicate carbon gain, while negative values denote carbon loss. To facilitate interpretation, the four modes were further grouped into two broader categories: (1) Tree cover change (TCC)-induced changes, including forest growth and tree cover loss, which occur without a change in land-cover class; and (2) Land-use and land-cover change (LULCC)-induced changes, including forest expansion and forest loss, which involve transitions between forest and non-forest classes. Finally, we aggregated pixel-level ΔAGC to quantify the relative contributions of each change mode and category to the national AGC balance between 1985 and 2023.

3 Results

3.1 Feature selection and importance

The RFE results indicated that the model's cross-validation score increased with the number of features, stabilizing around 10 features, with the highest performance at 17 features (Fig. 2a). Therefore, we selected 17 explanatory variables to train the ResNet, including five spectral bands, three vegetation indices, three topographic factors, three climate variables, geographic location, tree cover, and land cover. Feature importance analysis showed that among the 17 explanatory variables, tree cover contributed the most to the model, followed by slope and longitude. The near-infrared band and EVI also demonstrated high importance, while other spectral features generally ranked lower than climate and topographic variables (Fig. 2b).

Figure 2 Feature selection and importance analysis. (a) Cross-validation scores of the model with different numbers of features, measured by negative mean squared error. (b) Feature importance distribution based on Random Forest.

3.2 Accuracy assessment

3.2.1 Model performance

The ResNet model exhibited strong predictive performance for AGBD. In 5-fold cross-validation, the model achieved R²= 0.91, RMSE = 16.49 Mg ha⁻¹, and a minimal bias of 0.50 Mg ha⁻¹, demonstrating its overall accuracy. Temporal validation using the LOYO approach yielded slightly lower but still robust performance (R²= 0.85, RMSE = 21.20 Mg ha⁻¹, Bias = 0.66 Mg ha⁻¹; Fig. 3), highlighting the model's ability to generalize across years. The marginally lower performance observed in the LOYO validation, likely reflects the temporal distribution imbalance of training samples, particularly the limited representation of high-biomass conditions in specific years. Nevertheless, the model still delivered satisfactory predictions across all years. Moreover, the ResNet model outperformed other machine learning ensemble models (e.g., Random Forest, XGBoost, and LightGBM), particularly in mitigating the effects of spectral saturation in high-biomass forests (see Sect. S3 for details).

Figure 3 Model evaluation of AGBD predictions. (a) Overall performance based on 5-fold cross-validation. (b) Temporal generalization assessed using Leave-One-Year-Out (LOYO) validation for 2019–2021.

3.2.2 Independent accuracy evaluation using field survey data

The validation using historical field survey data confirms that the AGBD estimates exhibit reasonable accuracy across different time periods. For the period 1985–2008, the model achieved an R² of 0.63, an RMSE of 68.26 Mg ha⁻¹, and a bias of −19.87 Mg ha⁻¹ (Table 1, Fig. 4). Among ecoregions, the cold temperate coniferous forest showed the highest accuracy, with the lowest RMSE (25.73 Mg ha⁻¹). Similarly, the temperate coniferous and deciduous broadleaf mixed forest achieved low RMSE (30.74 Mg ha⁻¹) and minimal bias (4.36 Mg ha⁻¹). In contrast, the tropical monsoon rainforest and rainforest regions presented greater challenges, with higher RMSE (90.82 Mg ha⁻¹) and a significant negative bias (−68.57 Mg ha⁻¹). Other ecoregions, such as the temperate desert and subtropical evergreen broadleaf forest, displayed intermediate accuracy levels, with RMSEs of 47.82 and 67.83 Mg ha⁻¹ and biases of 29.54 and −36.48 Mg ha⁻¹, respectively.

Figure 4 Scatter plot comparing field survey data with ResNet-estimated AGBD. The black dashed line represents the 1:1 line, while the red solid line represents the fitted regression line.

Table 1 Independent accuracy evaluation using field survey data across different ecoregions.

The correlation analysis between provincial-level AGC estimates and NFI forest AGC statistics reveals strong correlations, with coefficients consistently exceeding 0.85 (p<0.05) across different assessment periods, ranging from 0.85 to 0.89 (Fig. 5). These results suggest that the AGC estimates effectively capture the temporal changes in national AGC, as reflected in the NFI's multi-period data. The consistent correlations observed across the various assessment periods demonstrate that the AGC estimates are reliable indicators of forest carbon stock dynamics over time.

Figure 5 Scatterplots comparing provincial-level aboveground biomass carbon (AGC) from National Forest Inventory (NFI) data and ResNet-based estimates of AGBD for 1984–2008. Panels (a)–(e) correspond to five-year intervals: 1984, 1989, 1994, 1999, and 2004. Each point represents a province, and the 1:1 line (dash line) indicates perfect agreement.

3.2.3 Residual analysis

The comparison between our AGBD estimates and historical ground surveys (Fig. 4) indicates a systematic underestimation in some regions, as reflected by regression slopes lower than one and a negative mean bias. Additional diagnostics (Table S3) show that this bias increases with observed biomass density, implying greater underestimation in high-biomass plots. This pattern is consistent with the well-known spectral saturation effect, whereby optical indices lose sensitivity in dense canopies, as well as with potential scale and methodological mismatches between satellite-based estimates and historical field plots.

The residual analysis stratified by geographic location, tree cover, and terrain further underscores these tendencies (Fig. 6). Residuals are wider in low latitude (< 25° N) and high latitude (> 40° N) regions compared to mid-latitudes. As slope and tree cover increase, residual spread also broadens, indicating that predictions tend to be less accurate in steep or densely forested areas. By contrast, no significant biases were observed across longitudes, aspects, or elevations, suggesting that errors are more strongly associated with biomass density and slope rather than with all terrain factors. These findings highlight slope- and biomass-related challenges that should be addressed in future AGBD mapping efforts.

Figure 6 Residual analysis for forest AGBD estimation. Panels (a)–(f) depict AGBD prediction residuals under varying latitudes (° N), longitudes (° E), aspects (°), slopes (°), elevations (m), and tree cover (%), respectively. Positive residuals indicate predicted values higher than reference values, while negative residuals signify underestimation of AGBD compared to GEDI references. The box plot displays the median, quartiles, as well as the 10th and 90th percentiles of the residuals. The analysis was performed based on the 20 % of 2019–2021 GEDI AGBD footprints and their corresponding estimates.

3.3 Evaluation of AGBD temporal dynamics

3.3.1 Cross-comparison with VOD and GPP

The cross-comparison between AGBD trends and independent satellite-based datasets (GPP and VOD) further supports the temporal reliability of our product. Pixel-level correlation analysis (Fig. 7a, b) shows that positive correlations dominate for both GPP and VOD, while negative correlations are relatively limited. For VOD, 37.7 % of pixels exhibit highly significant positive correlations (p<0.01) and 6.1 % are significant ($\mathrm{0.01}<p<\mathrm{0.05}$), with highly significant negative correlations accounting for only 5.6 %. Similarly, for GPP, 34.3 % of pixels show highly significant positive correlations and 8.3 % are significant, while negative correlations remain rare (3.2 % highly significant, 2.6 % significant). Non-significant correlations comprise roughly 25 %–31 % of pixels.

Temporal trajectories of national mean AGBD closely follow those of GPP (1992–2020) and VOD (1987–2021) (Fig. 7e–f). Scatterplots of annual national means (Fig. 7g–h) demonstrate strong positive linear relationships, with R²= 0.65 for both GPP and VOD, indicating substantial consistency in interannual patterns. Collectively, the predominance of significant positive correlations and the alignment of temporal trajectories corroborate that the estimated AGBD trends are ecologically plausible and effectively capture large-scale interannual to decadal variations in vegetation biomass.

Figure 7 Cross-comparison of temporal dynamics between estimated AGBD and satellite-based datasets (GPP and VOD). (a, b) Pixel-level correlation maps between AGBD and GPP/VOD trends. (c–d) Fractions of highly significant (p<0.01), significant ($\mathrm{0.01}\le p<\mathrm{0.05}$), and non-significant (p≥ 0.05) pixels. (e, f) Temporal trajectories of national mean AGBD and GPP/VOD. (g–h) Relationships between national mean AGBD and GPP (1992–2020) or VOD (1987–2021).

3.3.2 Temporal validation via before–after–control–impact (BACI) analysis

The BACI analysis provides additional evidence that the estimated AGBD trends are consistent with known forest dynamics and management interventions. Typical examples of forest growth and natural forest degradation detected by CFATD highlight its effectiveness in capturing temporal changes in forest biomass (Figs. 8–9). In the Three-North shelterbelt regions (Fig. 8), BACI confirms the establishment of plantations, showing a sigmoidal growth curve typical of young forests: AGBD increased at 5.60 Mg ha⁻¹ yr⁻¹ during the first phase (1994–1999), accelerated to 8.45 Mg ha⁻¹ yr⁻¹ in the second phase (2000–2007), and stabilized at −0.03 Mg ha⁻¹ yr⁻¹ in the final phase (2008–2023). In contrast, natural dark coniferous forests in the alpine valleys of Sichuan and Yunnan (Fig. 9) experienced continuous degradation. Despite stable or slightly expanded forest area, AGBD declined from 292 Mg ha⁻¹ in 1985 to 136 Mg ha⁻¹ in 2023, with a degradation rate of −3.15 Mg ha⁻¹ yr⁻¹. BACI results reveal that this decline is temporally consistent with disturbances and policy interventions, with no recovery observed in the impacted areas, while nearby control regions remain relatively stable.

Figure 8 Temporal increase in AGBD of newly established forests under the Three-North Shelterbelt Program from 1985 to 2023, as captured by CFATD. (Satellite imagery: © Google Earth).

Figure 9 Temporal decline in AGBD of degraded natural forests in southwest China from 1985 to 2023, as captured by CFATD. (Satellite imagery: © Google Earth).

When combined with the cross-comparison with satellite-based proxies (GPP and VOD, Fig. 7), these BACI-based case studies reinforce the reliability of AGBD temporal trends. Significant positive correlations with GPP and VOD, consistent national mean trajectories, and strong linear relationships (R² = 0.65) across years collectively demonstrate that the estimated AGBD trends capture ecologically plausible interannual to decadal variations, responding appropriately to both natural disturbances and anthropogenic interventions.

3.4 Spatial pattern of forest AGB

The 30 m resolution mapping of China's forest AGBD reveals substantial spatial variability in biomass density (Fig. 10). In 2023, the national average forest AGBD was 122.69 Mg ha⁻¹, with a total forest AGB stock amounting to 13.97 PgC. Among the seven geographic regions, the Southwest (142 Mg ha⁻¹), Northeast (123 Mg ha⁻¹), and Northwest (123 Mg ha⁻¹) regions exhibited biomass densities above the national average. High-density areas with forest AGBD exceeding 300 Mg ha⁻¹ were predominantly located in these regions, including southern Tibet, the Qinling Mountains, and parts of Northeast China. Additionally, central Taiwan exhibited several hotspots with forest AGBD greater than 300 Mg ha⁻¹. These regions are notably rich in forest resources and represent significant areas of natural forest distribution in China. Conversely, the northern, central, eastern, and southern regions have historically been subjected to frequent human disturbances and high-intensity impacts, leading to younger forest stands and relatively lower biomass densities (Fig. 10e). Nonetheless, these regions present a substantial potential for carbon sequestration. The provinces of Yunnan, Sichuan, Heilongjiang, Guangxi, and Tibet are the top five in terms of forest biomass carbon stocks, collectively accounting for 41.5 % of the national total (Fig. 10f). This highlights the uneven spatial distribution of forest resources across China.

Figure 10 Spatial patterns of forest AGBD and regional summaries. (a–d) Spatial distribution of forest AGBD in China for 1985, 2000, 2015, and 2023, respectively. (e) Mean AGBD across different geographical regions. (f) Estimated forest AGC stock for different provinces.

In 2023, the tropical rainforest exhibited the highest AGBD (165.17 Mg ha⁻¹), followed by the temperate coniferous and deciduous broadleaf mixed forest (136.38 Mg ha⁻¹), subtropical evergreen broadleaf forest (123.88 Mg ha⁻¹), warm temperate deciduous broadleaf forest (111.08 Mg ha⁻¹), cold temperate coniferous forest (107.77 Mg ha⁻¹), temperate desert region (97.50 Mg ha⁻¹), and Tibetan plateau alpine region (72.98 Mg ha⁻¹).

The pixel-level analysis reveals that the average uncertainty in AGBD estimates across the study area is 13.94 Mg ha⁻¹, with individual pixel uncertainties ranging from 0 to 85 Mg ha⁻¹ (Fig. 11a). The vast majority of the study area (99.35 %) exhibits an uncertainty of less than 30 Mg ha⁻¹, while only a small fraction (0.65 %) has uncertainties of 30 Mg ha⁻¹ or greater (Fig. 11b). Notably, uncertainties are higher than the national average in the southwestern and southern regions, at 15.55 and 16.62 Mg ha⁻¹, respectively. In contrast, the northwestern and northern regions have relatively lower uncertainties, with estimates falling below 10 Mg ha⁻¹ (Fig. 11c). Overall, the uncertainty in forest AGBD across the study area is well-controlled, indicating a reliable dataset. However, the higher uncertainties in specific regions underscore the need for further refinement in these areas, particularly in complex terrains where topographic and vegetation heterogeneity may challenge the accuracy of AGBD estimates (Fig. 6d).

Figure 11 Uncertainty analysis of forest AGB in China for the year 2023. (a) Predictive uncertainty map (expressed as standard deviation) of forest AGBD across China. (b) Frequency distribution of predictive uncertainty, categorized into 10 Mg ha⁻¹ intervals. (c) Regional statistics of predictive uncertainty across different geographical areas.

3.5 Temporal trend and dynamics of forest AGB and AGC in China

The temporal analysis of China's average forest AGBD reveals a significant upward trend over the study period, with the average AGBD increasing from 95.74 ± 11.30 Mg ha⁻¹ in 1985 to 122.69 ± 13.94 Mg ha⁻¹ in 2023 (Fig. 12b). The average AGBD grew at an annual rate of 0.69 Mg ha⁻¹ yr⁻¹, representing a 28.1 % increase over the entire period. Most regions experienced an increase in AGBD during the study period, with the exception of certain areas in northeastern and northern China, where agricultural expansion and deforestation likely contributed to forest AGBD reduction (Fig. 12a). In contrast, significant increases in AGBD were observed in the Three-North Shelterbelt region, as well as in southwestern and southern China. In the Three-North region, extensive afforestation and the construction of shelter forests transformed non-forested areas into forested land, resulting in a substantial increase in AGBD. Similarly, large-scale afforestation efforts as part of desertification control in the southwestern and southern (e.g., Guangxi) regions, also contributed to a significant rise in AGBD. Among the seven geographical regions, the southern region experienced the fastest growth in AGBD (0.90 Mg ha⁻¹ yr⁻¹), followed by the northeast (0.81 Mg ha⁻¹ yr⁻¹), southwest (0.79 Mg ha⁻¹ yr⁻¹), north (0.73 Mg ha⁻¹ yr⁻¹), northwest (0.67 Mg ha⁻¹ yr⁻¹), central (0.41 Mg ha⁻¹ yr⁻¹), and eastern regions, which exhibited the slowest growth rate (0.26 Mg ha⁻¹ yr⁻¹) (Fig. 12b).

Figure 12 Long-term (1985–2023) trends in forest AGB and AGC stocks in China. (a) Pixel-level forest AGBD trend map derived from Sen's slope analysis (Sen, 1968). (b) and (c) display annual changes in forest AGBD from 1985 to 2023 across geographical regions and ecoregions, respectively. (d) and (e) illustrate annual changes in total forest AGC stock over the same period, also across geographical regions and ecoregions. AGC is calculated from biomass using a conversion factor of 0.5.

We observed increasing AGBD trends across most ecoregions, except for the Tibetan plateau alpine region, where forest AGBD declined from 80.07 Mg ha⁻¹ in 1985 to 72.98 Mg ha⁻¹ in 2023, with a trend of −0.12 Mg ha⁻¹ yr⁻¹. Alpine forests, being more sensitive to environmental changes, may have experienced degradation due to climate warming and extreme weather events. The tropical rainforest showed the most pronounced increase in AGBD (1.17 Mg ha⁻¹ yr⁻¹), followed by the temperate coniferous and deciduous broadleaf mixed forest (1.00 Mg ha⁻¹ yr⁻¹), subtropical evergreen broadleaf forest (0.67 Mg ha⁻¹ yr⁻¹), cold temperate coniferous forest (0.64 Mg ha⁻¹ yr⁻¹), temperate desert region (0.60 Mg ha⁻¹ yr⁻¹), and warm temperate deciduous broadleaf forest (0.48 Mg ha⁻¹ yr⁻¹).

The spatial distribution of China's forest aboveground carbon (AGC) stock in 2023 was uneven. The southwestern region held the largest share, with 4.84 PgC, accounting for over one-third of the national total (34.43 %). The northeast contributed 13.65 %, followed by the east (12.94 %), the south (11.32 %), the central region (9.80 %), the north (9.19 %), and the northwest (8.66 %) (Fig. 12d). Disparities in AGC stock were even more pronounced across different forest types. The subtropical evergreen broadleaf forest contained 8.82 PgC, making up over half of the national total (62.80 %), followed by the temperate coniferous and deciduous broadleaf mixed forest (10.32 %), warm temperate deciduous broadleaf forest (9.11 %), tropical rainforest (8.90 %), cold temperate coniferous forest (5.44 %), temperate desert region (2.91 %), and the Tibetan Plateau alpine region (0.55 %) (Fig. 12e).

Throughout the study period, China's forest AGC showed an upward trend, with total forest AGC stock expanding from 5.50 ± 0.23 PgC in 1985 to 13.97 ± 0.87 PgC in 2023 (Fig. 12d, e). Across all geographic and ecological regions, total forest AGC stock consistently increased from 1985 to 2023. The southwestern region experienced the fastest increase (0.09 PgC per year), while growth rates in other regions were relatively similar (0.02 PgC per year). Among ecoregions, the subtropical evergreen broadleaf forest showed the fastest growth (0.15 PgC per year), while the Tibetan Plateau alpine region exhibited the slowest growth (0.001 PgC per year). However, net carbon storage initially declined at a rate of −0.02 PgC per year during 1985–1989. This trend reversed after 1990, with a robust increase of 0.25 PgC per year. Overall, the total forest AGC stock increased at a rate of 0.22 ± 0.01 PgC per year, effectively doubling China's forest carbon reserves over the past 40 years, achieving a remarkable growth of 154.0 % (see Sect. S4 for discussion). The forest net AGC sink offset 11.5 %–14.9 % of fossil fuel and industrial emissions between 1985 and 2023.

3.6 Contribution of forest change to changes in AGC

This study highlights the dual role forests played in carbon sequestration and emissions in China during the period from 1985 to 2023, depending on the nature of forest changes. Forest growth, which involves increased tree cover in already forested areas, was a major contributor to carbon sequestration, particularly in regions south of the Qinling Mountains (Fig. 13a). The spatial changes in AGC induced by forest growth closely mirror the age distribution of China's forests. The carbon sequestration per unit area due to forest growth was significantly higher in southern regions compared to northern, northeastern, southwestern, and Taiwan, where natural forests dominate. These older forests with higher canopy cover had limited potential for further AGC increase, following the “law of constant final yield” (Adler et al., 2018), which explains why AGC growth in younger and less dense forests in central, southern, and eastern China outpaced that of other regions, reflecting the effects of long-term conservation efforts.

Forest expansion also contributed to additional carbon sequestration by converting non-forested land to forest, with significant gains observed in the southwest and western regions (Fig. 13b). This is closely related to reforestation efforts initiated since 1985. These areas, previously allocated for other land uses, significantly boosted carbon stocks through afforestation initiatives. Notably, we observed that carbon sequestration per unit area from forest expansion was considerably lower than that from forest growth, which can be explained by the fact that AGC accumulation rates follow a “slow-fast-slow” pattern as forests mature (see Fig. 8). Hence, the carbon sequestration in newly planted forests was naturally lower than that in pre-existing forests. However, the carbon sequestration due to forest expansion in the southwest was significantly higher than in other regions, aligning with recent studies showing that the carbon sink potential in southwest China may have been historically underestimated (Wang et al., 2020). Conversely, tree cover loss and forest-to-non-forest transitions were the primary contributors to carbon emissions, mainly concentrated in northern and northeastern regions (Fig. 13c, d). The spatial patterns of carbon emissions due to forest loss suggest that deforestation and land-use conversion were concentrated in the same regions experiencing high developmental pressure and agricultural expansion.

Figure 13 AGC Spatial changes induced by forest dynamics in China between 1985 and 2023. (a) and (b) represent the spatial distribution of carbon sequestration driven by forest growth and expansion, respectively. (c) and (d) illustrate the spatial distribution of carbon emissions resulting from tree cover loss and forest loss, respectively.

We observed that approximately 65.1 % (5.75 PgC) of the carbon sequestration was due to forest growth, while the remaining 34.9 % (3.09 PgC) resulted from forest expansion (Fig. 14). Although forest growth contributed more to the AGC stock increment, this finding suggests that afforestation activities since 1985 have significantly boosted AGC, accounting for nearly one-third of the total increase. These initiatives have involved planting new forests on previously non-forested land and restoring degraded forest areas, thus expanding China's forested landscapes and bolstering their carbon sequestration capacity. This substantial contribution underscores the positive impact of forest conservation policies and enhanced management practices, which have facilitated forest recovery, maturation, and increased carbon sequestration over time. In terms of carbon emissions, tree cover loss accounted for 47.77 % of the AGC reduction, while forest loss (forest to non-forest transitions) was responsible for 52.23 %. This indicates that complete deforestation had a more substantial impact on AGC loss than reductions in canopy density.

Figure 14 Quantification of AGC changes induced by forest dynamics in China between 1985 and 2023.

The contrasting spatial dynamics between carbon sequestration and emissions underscore the complexity of AGC changes over the past four decades. While large areas of China have acted as carbon sinks due to reforestation and forest regeneration, other regions have simultaneously contributed to significant carbon emissions through deforestation and land degradation. This dual impact highlights the critical importance of targeted forest management strategies aimed at enhancing carbon sequestration while minimizing deforestation and forest degradation, ensuring a balanced and sustainable approach to forest carbon dynamics in the future.

4 Discussion

4.1 Comparison with external AGBD datasets

4.1.1 Static and spatial comparison

In pixel-level comparisons of our AGBD map against datasets from Zarin et al. (2016), Hu et al. (2016), Su et al. (2016), Santoro et al. (2021), Yang et al. (2023), and ESA CCI (see Sect. S2 for more information about the external AGBD datasets), the average differences were −5.86, −45.23, 24.78, 47.91, 19.55, and 60.14 Mg ha⁻¹, respectively (Fig. 15). Our map aligns more closely with Zarin et al. (2016), Su et al. (2016), Yang et al. (2023), while larger discrepancies are observed with the Hu, Santoro, and ESA CCI datasets. Specifically, the percentage of pixels with absolute differences less than 50 and 25 Mg ha⁻¹, respectively, are as follows: Zarin et al. (49.5 %, 25.9 %), Hu et al. (50.7 %, 27.6 %), Su et al. (59.8 %, 33.9%), Santoro et al. (53.2 %, 26.4 %), Yang et al. (65.8 %, 39.0 %), and ESA CCI (38.1 %, 18.2 %).

Figure 15 Difference between the generated AGBD map and external AGBD datasets. (a)–(f) are pixel-level differences between the estimated forest AGBD from CFATD and the forest AGBD from Zarin et al. (2016), Hu et al. (2016), Su et al. (2016), Yang et al. (2023), and ESA CCI, respectively (the estimated AGBD minus external AGBD).

In southern, eastern, and southwestern China, our estimates are more than 50 Mg ha⁻¹ lower than those of Zarin, Hu, and Su. The coarser spatial resolution of these products likely leads to aggregation bias, potentially overestimating AGBD in fragmented forest landscapes. Temporal mismatches and intensive human activities in these regions exacerbate the differences. Conversely, our estimates in northern and northwestern regions are higher than Zarin and Su but align well with Hu, with minimal differences (< 25 Mg ha⁻¹). Compared to Santoro and ESA CCI, our map consistently reports higher AGBD across China, as their global focus and lack of region-specific samples likely lead to underestimations. The Santoro and ESA CCI maps, which report average AGBD of 57.05 and 64.88 Mg ha⁻¹ for 2010 and 2020, respectively, represent the lowest levels listed in Table S1. The Yang map shows strong spatial consistency with ours but underestimates AGBD in hotspots like the Qinling Mountains, southern Tibet, and Taiwan province due to limited field samples and spectral saturation.

Key differences among datasets stem from their capacity to capture spatial variability, manage spatial resolution, and comprehensively represent forested regions. Our map captures fine-scale patterns, such as higher AGBD in forest interiors and lower AGBD at edges (Fig. 16), outperforming Yang in variability representation, while Zarin, Santoro, and ESA CCI also perform well in this aspect. Coarse-resolution datasets, such as the Hu map (Fig. 16f) and Su map (Fig. 16g) (1000 m), show significantly reduced accuracy in fragmented landscapes, making them unsuitable for precise carbon accounting. Furthermore, the quality of forest thematic maps plays a critical role in AGBD mapping; only CFATD, Zarin, and Santoro effectively detected shelterbelts within agricultural areas, whereas other datasets either partially missed or entirely failed to capture these features (The third scenario of Fig. 16).

Figure 16 Comparison of AGBD product details across different datasets. (a) Google high-resolution imagery (© Google Earth); (b)–(d) represent the CFATD products for the years 2000, 2010, and 2020, respectively; (e)–(j) display AGBD maps from Zarin et al. (2016), Hu et al. (2016), Su et al. (2016), Santoro et al. (2021), Yang et al. (2023), and ESA CCI. Each row corresponds to subregions representing different forest types. Note: The numbers in parentheses indicate the corresponding year of each dataset.

4.1.2 Dynamic and trend comparison

Our AGBD dataset reveals a consistent upward trend from 1985 to 2023, highlighting the steady recovery and effective management of China's forests (Fig. 17). This trajectory is broadly in line with national afforestation and conservation programs, and provides a robust depiction of long-term carbon stock dynamics (Zeng et al., 2023). To benchmark our results, we compared them with four independent AGB time-series datasets (see Sect. S2).

Figure 17 Comparative analysis of temporal dynamics across different AGB products.

The comparison highlights substantial differences in both magnitude and temporal pattern. For instance, the Hengeveld dataset (Hengeveld et al., 2015) exhibits pronounced fluctuations, including a decline during the late 1990s and early 2000s, while the ESA CCI product shows a general decrease in the 2010s. Chen et al. (2023) also report consistently lower AGBD values than our estimates. These discrepancies arise from several key factors: (1) spatial resolution – most global products are provided at ≥ 1 km resolution, which smooths local heterogeneity and suppresses extremes captured by our 30 m dataset; (2) forest definition and masking – different tree-cover thresholds (e.g., Chen et al., 2023 use > 15 % canopy cover) alter the mapped forest extent and introduce mixed pixels, leading to lower mean AGBD. Notably, when Chen's product is aggregated with the 1 km forest cover map used by Hu et al. (2016), its average AGBD increases substantially, ranging from 68.94 Mg ha⁻¹ in 2002 to 75.58 Mg ha⁻¹ in 2021, highlighting the sensitivity of biomass estimates to forest definitions. (3) input data and modeling approaches – radar- or lidar-upscaled products differ in their sensitivity to canopy structure compared with optical or inventory-driven models; (4) many existing datasets rely on coarse LULC or forest masks that cannot effectively capture China's large-scale afforestation efforts (Cai et al., 2024), thereby underestimating biomass gains; and (5) temporal coverage and reference data – the calibration periods and reference samples vary among products, creating temporal mismatches in trend estimation (Quegan et al., 2019; Santoro et al., 2021).

Despite these differences, our dataset shows stronger agreement with independent evidence from national forest inventory statistics and regional studies (Fang et al., 2001; Liu et al., 2015; Piao et al., 2005; Zeng et al., 2023), all of which support a long-term upward trajectory of China's forest biomass. This comparison further underscores the value of high-resolution, long-term mapping for capturing heterogeneous forest dynamics in China.

4.2 Limitations and future improvement

Current AGBD maps often overestimate low values and underestimate high values – a persistent pattern well documented in previous studies (Réjou-Méchain et al., 2019; Rodríguez-Veiga et al., 2019). These discrepancies are largely driven by sensor limitations, such as the saturation of optical and radar signals in high-biomass forests, leading to systematic underestimation in dense tropical or subtropical stands. The overestimation of low AGBD is more complex, often arising from regression models fitted to limited or non-representative plot data. Additional systematic biases can result from measurement errors, incorrect allometric models, and environmental factors (e.g., precipitation, snow cover) that influence remote sensing signals (Santoro et al., 2021). Furthermore, mismatches between the spatial resolution of remote sensing data and the size or geolocation accuracy of field plots can distort pixel-level estimates, underscoring the need for improved plot–pixel matching and more representative field measurements. It is important to note that the validation accuracy in this study was achieved using a strictly curated dataset. Including raw field plots without visual screening for homogeneity would likely result in significantly lower correlations due to the scale mismatch between field measurements and satellite pixels.

Against this backdrop, our estimates indicate a substantial increase in mean forest AGBD from 95.74 to 122.69 Mg ha⁻¹ and in total AGC stock from 5.50 to 13.97 PgC over the study period. These upward trends are broadly consistent with the scale and trajectory of China's large-scale afforestation and forest rehabilitation programmes since the 1980s – especially post-1999 initiatives such as the Grain-for-Green Program – which have simultaneously expanded forest area and increased stand biomass through natural growth and management. Independent comparisons with provincial NFI statistics and several external products show similar upward trajectories (see Figs. 5 and 17), lending confidence to the observed direction of change.

Nevertheless, we acknowledge several sources of systematic uncertainty that could affect the magnitude of these increases. First, optical indices tend to saturate in very high-biomass forests, potentially biasing estimates downward in tropical and old-growth stands. Second, our definition of forest (≥ 20 % tree-cover threshold) directly influences the mapped forest extent and thus the total AGC stock; a higher threshold would yield lower total area and carbon estimates (see Sect. S5). Third, the use of a fixed biomass-to-carbon conversion factor (0.5) ignores species- and region-specific variation in carbon content, which may introduce bias in regional and national totals. Finally, the GEDI footprint AGBD dataset, which informs part of our modelling framework, lacks direct training data from China and neighbouring countries such as Russia (Duncanson et al., 2022). Although the Level 4A algorithm is designed for geographic transferability (Kellner et al., 2023), the absence of local calibration data may still introduce uncertainty in parts of China.

While deep learning approaches have been employed here to reduce spectral saturation effects and improve estimation accuracy, further advances are possible. Integrating multi-source remote sensing data (optical, radar, LiDAR) can better capture forest structural complexity and mitigate the limitations of individual sensors (see Sect. S6). Developing region- and forest-type-specific allometric models will further refine biomass estimation in heterogeneous landscapes. Expanding calibration and validation with more extensive and diverse ground data, particularly from underrepresented regions, will improve both accuracy and transferability. Incorporating time-series data into modelling frameworks can also help distinguish between short-term fluctuations and long-term growth trends (Araza et al., 2022; Duncanson et al., 2019).

To guide proper use of these maps, we note several best-use recommendations: the products are most reliable at regional to national scales rather than for individual plots; caution is warranted when interpreting high-biomass tropical or subtropical forests due to potential saturation, and users are encouraged to cross-validate with other datasets when conducting local analyses; and users should consider forest definition thresholds and local biomass-to-carbon conversion factors when estimating carbon stocks. Where possible, independent validation with local field data is recommended before applying these maps for management or policy decisions. Finally, continued development of deep learning models and integration of multi-source data are likely to further enhance both the accuracy and applicability of future AGB products.

5 Data availability

The China Forest AGB Time Series Dataset (CFATD) generated by this study is accessible via the Google Earth Engine script at this link: https://code.earthengine.google.com/4f8ad8d32ddb84e826e941a95f31f9be (last access: 22 February 2025). Please note that access to the script requires a Google Earth Engine login. Users also can download the CFATD from Zenodo: (Part I: 1985–1993) https://doi.org/10.5281/zenodo.12620984 (Cai et al., 2025a), (Part II: 1994–2001) https://doi.org/10.5281/zenodo.12637101 (Cai et al., 2025b), (Part III: 2002–2008) https://doi.org/10.5281/zenodo.12655492 (Cai et al., 2025c), (Part IV: 2009–2015) https://doi.org/10.5281/zenodo.12658255 (Cai et al., 2025d), (Part V: 2016–2021) https://doi.org/10.5281/zenodo.12742210 (Cai et al., 2025e), (Part VI: 2022–2023) https://doi.org/10.5281/zenodo.12747329 (Cai et al., 2025f).

6 Conclusions

In this study, we developed a comprehensive approach to estimate the spatiotemporal distribution of forest AGB across China by integrating multisource datasets – including spectral, climatic, topographic, and tree cover data – within a deep learning framework. Using over 50 000 GEDI sample points, we trained a ResNet model to produce the China Forest AGB Time Series Dataset (CFATD), a 30 m resolution dataset spanning 1985 to 2023, accompanied by uncertainty estimates. Validation with the 2019–2021 GEDI dataset (R² = 0.91, RMSE = 16.49 Mg ha⁻¹, Bias = 0.50 Mg ha⁻¹) and multi-year ground survey data (R² = 0.63, RMSE = 68.26 Mg ha⁻¹, Bias = −19.87 Mg ha⁻¹) confirmed the robustness of the AGB time series maps, which effectively capture the dynamics of forest biomass. The CFATD not only aligns with existing products in spatial consistency but also provides an improved representation of fine-scale spatial heterogeneity and temporal trends. Our nationwide AGB estimates revealed that as of 2023, China's average forest AGBD was 122.69 ± 13.94 Mg ha⁻¹. The regions with the highest AGB carbon densities were concentrated in the Southwest, Northwest, Northeast, and Taiwan. The total forest AGC stock was estimated at 13.97 ± 0.87 PgC, with 41.5 % concentrated in five provinces, highlighting the uneven distribution of forest resources. Between 1985 and 2023, average AGBD increased from 91.58 ± 9.27 to 122.69 ± 13.94 Mg ha⁻¹, and total carbon stock rose from 5.50 ± 0.23 to 13.97 ± 0.87 PgC, with an annual net accumulation rate of 0.22 PgC. We found that 65.1 % of the increased AGC stock was attributed to forest growth, with the remaining 34.9 % attributed to forest expansion. In terms of carbon emissions, tree cover loss contributed 47.8 % of AGC reductions, while forest loss (forest-to-non-forest transitions) accounted for 52.2 %. The CFATD fills a critical gap by providing the first long-term, high-resolution forest AGB dataset for China. This dataset will be crucial for monitoring forest carbon dynamics, assessing climate change impacts, and guiding forest management and conservation policies.