Carbon storage through China’s planted forest expansion

China’s extensive planted forests play a crucial role in carbon storage, vital for climate change mitigation. However, the complex spatiotemporal dynamics of China’s planted forest area and its carbon storage remain uncaptured. Here we reveal such changes in China’s planted forests from 1990 to 2020 using satellite and field data. Results show a doubling of planted forest area, a trend that intensified post-2000. These changes lead to China’s planted forest carbon storage increasing from 675.6 ± 12.5 Tg C in 1990 to 1,873.1 ± 16.2 Tg C in 2020, with an average rate of ~ 40 Tg C yr−1. The area expansion of planted forests contributed ~ 53% (637.2 ± 5.4 Tg C) of the total above increased carbon storage in planted forests compared with planted forest growth. This proactive policy-driven expansion of planted forests has catalyzed a swift increase in carbon storage, aligning with China’s Carbon Neutrality Target for 2060.


Supplementary Note 2 | Definitions of vegetation, temporal and textural indices.
The vegetation indices used in this work included the bare soil index (BSI), enhanced vegetation index (EVI), normalized difference vegetation index (NDVI), modified soil-adjusted vegetation index (MSAVI), and soil-adjusted vegetation index (SAVI).
They were calculated using the equations below: are the coefficients of the trigonometric components with frequencies   , and  0 is the coefficient when the frequency is zero.
Textural features were calculated using the gray-level co-occurrence matrix approach, which is widely used in textural information extraction in remote sensing classification [4][5][6] .This method can be easily implemented in google earth engine (GEE) platform and enables the extraction of 18 textural indices for the input band 7 .
The textural information between natural forests and planted forests showed considerable differences in the remote sensing images 8 .All 18 textural indices of each spectral band and vegetation indices were retained as the feature library (Supplementary Table 5).Each textural index is explained in detail on the website (https://developers.google.com/earth-engine/apidocs/ee-image-glcmtexture).
As terrain also influences the distribution of planted forests 9,10 , elevation, aspect, and slope features were derived according to the digital elevation model (DEM) data of the Shuttle Radar Topography Mission in GEE.In total, 220 features were constructed for mapping planted forests (Supplementary Table 5).
Given the computation resource limitation of GEE, we divided China into six geographic zones and implemented mapping at regional scales.The random forest (RF) classifier was used to map planted forests owing to its ability to handle highdimensional data, tolerance to sample errors, and robustness to missing data 2,11 .After tests were conducted with different numbers of trees, from 10 to 500, the number of trees was set to 100, according to the exported overall classification accuracies.
To reduce the number of features used in the classification process and economize computational resources, we used the recursive feature elimination with crossvalidation (RFE-CV) feature selection approach on a local computer to construct the optimal feature set for each region and period 12 .This approach involves searching for a subset of features, starting with all features in the training dataset and eliminating features that negatively influence accuracy through cross-validation.RFE-CV was implemented using Python 3.8 on a local computer, with the k-fold parameter in CV set to 5. The optimal number of features and selected features for each region were determined according to the output standard deviation of accuracy (Supplementary Fig. 12 and Table 4).Using the constructed feature set for each region and period, we implemented planted forest mapping using RF and the training data on the GEE cloud-processing platform.

Supplementary Note 4 | Accuracy assessment of planted forest maps
The accuracy of the planted forest maps was assessed using four data sources: (1) the validation data from field surveys (35,074 samples), (2) the validation data from stratified random sampling, (3) the statistical data of China's planted forest resource inventory between 2000 and 2020 for each province, and (4) intercomparison with forest inventory map.The validation data from field surveys and stratified random sampling were used to construct a confusion matrix, in which the overall accuracy (OA), F1 score, producer's accuracy (PA), and user's accuracy (UA) were calculated using the following equations to assess the accuracies.
where   is the number of testing samples mapped correctly, Sum is the total number of testing samples, PA is the producer's accuracy, and UA is the user's accuracy.The macro-average of the F1 score is the average F1 value for all classes.
During the calculation of the average F1 value, all classes are assigned equal weights, and the actual sample frequency of occurrence is ignored.

Validation data from field surveys
A total of 35,074 field samples were used to validate our planted forest maps (Supplementary Fig. 12b).The evaluation results (Supplementary Fig. 14) indicated that the OA of the resultant maps ranged from 78.4% to 82.7% (Supplementary Fig. 4a).The UA (a measure of commission error) for planted forests was between 78.0% and 83.0%(Supplementary Fig. 4a).The PA (a measure of omission error) ranged from 78.0% to 88.0%.The F1 score (a balance of UA and PA) ranged from 0.88 to 0.93 (Supplementary Fig. 4a).Accuracy increased with time, with the planted forest map for 2020 yielding the highest accuracy.

Validation data from systematic random sampling
The systematic random sampling approach, which is regarded as one of the most robust methods for assessing the accuracy of LULC maps in remote sensing, was also applied to validate the resultant maps 13,14 .We first generated a grid with a resolution of 0.1°×0.1° in China's forest regions to collect samples for accuracy assessment.For each grid, one point was randomly generated, resulting in the generation of 15,578 random points (Supplementary Fig. 13).Each point was interpreted using the planted forest inventory map and google earth high resolution images by vegetation ecological experts.The systematic random validation samples for each map were generated and used to calculate the OA, F1, PA, and UA.The validation results revealed that the OA and F1 scores of the planted forest map for 2020 were 81.8% and 0.9, respectively, while the UA and PA of planted forests were 82.0% and 90.0% (Supplementary Fig. 4b), demonstrating the high reliability of the planted forest dataset for revealing the spatiotemporal dynamics of planted forests.

Supplementary Table 2 | Changes in C storage caused by conversion of other LULC to
planted forests at the regional scale.
1 +   ) − (  +   ) ( 1 +   ) + (  +   ) is the reconstructed NDVI time series, m is the number of harmonics and set to 1 in this study.  is the time when the original value of  ̃ was observed, where j=1, 2, …, N with N as the maximum number of observations in a time series.  , 3here   ,   ,   , and   are the surface reflectance values of Blue, Red, near-infrared, and shortwave infrared of the Landsat images.Temporal features were extracted through analysis of the 1990-2020 NDVI and EVI time-series data derived from Landsat 5 TM Collection 1 Tier 1 8-Day NDVI/EVI composite, Landsat 7 Collection 1 Tier 1 8-Day NDVI/EVI composite, and Landsat 8Collection 1 Tier 1 8-Day NDVI/EVI composite imagery.Harmonic analysis, a popular reconstruction method, was applied using the equation below to reduce the random noise of the NDVI/EVI time series and ensure the reliability of the extracted temporal features3.Then, temporal features, including amplitude, phase, magnitudes of the fitted time series, and the root mean squared error (RMSE) between the original and fitted values, were all extracted as the temporal feature set.