ARTICLE https://doi.org/10.1038/s41467-019-13798-8 OPEN Forest management in southern China generates short term extensive carbon sequestration Xiaowei Tong1,2,13, Martin Brandt 2,13, Yuemin Yue1,3*, Philippe Ciais 4, Martin Rudbeck Jepsen2, Josep Penuelas 5,6, Jean-Pierre Wigneron 7, Xiangming Xiao 8, Xiao-Peng Song 9, Stephanie Horion2, Kjeld Rasmussen2, Sassan Saatchi10, Lei Fan 7, Kelin Wang1,3*, Bing Zhang11, Zhengchao Chen11, Yuhang Wang 12, Xiaojun Li7 & Rasmus Fensholt 2 1234567890():,; Land use policies have turned southern China into one of the most intensively managed forest regions in the world, with actions maximizing forest cover on soils with marginal agricultural potential while concurrently increasing livelihoods and mitigating climate change. Based on satellite observations, here we show that diverse land use changes in southern China have increased standing aboveground carbon stocks by 0.11 ± 0.05 Pg C y−1 during 2002–2017. Most of this regional carbon sink was contributed by newly established forests (32%), while forests already existing contributed 24%. Forest growth in harvested forest areas contributed 16% and non-forest areas contributed 28% to the carbon sink, while timber harvest was tripled. Soil moisture declined significantly in 8% of the area. We demonstrate that land management in southern China has been removing an amount of carbon equivalent to 33% of regional fossil CO2 emissions during the last 6 years, but forest growth saturation, land competition for food production and soil-water depletion challenge the longevity of this carbon sink service. 1 Key Laboratory for Agro-ecological Processes in Subtropical Region, Institute of Subtropical Agriculture, Chinese Academy of Sciences, Changsha, China. 2 Department of Geosciences and Natural Resource Management, University of Copenhagen, Copenhagen, Denmark. 3 Huanjiang Observation and Research Station for Karst Ecosystem, Chinese Academy of Sciences, Huanjiang, China. 4 Laboratoire des Sciences du Climat et de l’Environnement, CEA-CNRS- UVSQ, CE Orme des Merisiers, Gif sur Yvette, France. 5 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Spain. 6 CREAF, Cerdanyola del Vallès, Spain. 7 ISPA, UMR 1391, INRA Nouvelle-Aquitaine, Bordeaux Villenave d’Ornon, France. 8 Department of Microbiology and Plant Biology, University of Oklahoma, Norman, OK, USA. 9 Department of Geosciences, Texas Tech University, Lubbock, TX 79409, USA. 10 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA. 11 Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing, China. 12 State Key Laboratory of Earth Surface Processes and Resource Ecology, Faculty of Geographical Science, Beijing Normal University, Beijing, China. 13These author contributed equally: Xiaowei Tong, Martin Brandt. *email: ymyue@isa.ac.cn; kelin@isa.ac.cn NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 M itigating climate change while securing livelihoods is a deforested areas, afforestation of croplands for conservation major societal challenge in the twenty-first century. purpose, and development of industrial timber and paper pro- Reducing emissions of CO2 from the combustion of fossil duction19. These diverse actions produce contrasted forest types fuels and deforestation, and developing land-based technologies to with different carbon sequestration and storage potentials, which sequester atmospheric CO2 are crucial to limit global warming1. are evaluated in this study. We apply a range of observational Current research indicates that vegetation growth and photo- data based on optical satellite imagery (including a full area synthetic activity are increasing globally2,3, which translates into an coverage of southern China at 2 m) and low-frequency passive increased terrestrial biomass carbon sink and thereby contributes to microwaves to identify different management strategies and mitigate the growth rate of atmospheric CO24. Cropland and pas- biomass carbon sinks from 2002 to 2017 (see Methods; ture land may accumulate carbon in soils2,5, but the storage refs. 23–25). Low-frequency microwave data are also used to potential is uncertain, whereas many studies indicate that adequate evaluate SM trends during the period 2010–2017 (see Methods; forest management can enhance biomass carbon stocks6,7. ref. 24). Results show that new forests are widespread over Forestation (reforestation and afforestation) sequesters carbon southern China representing an extensive carbon sink and off- in biomass, but the scalability of this land management option for setting an amount of carbon equivalent to one-third of the meeting ambitious warming limitation targets has been ques- regional fossil CO2 emissions, but possible soil water depletion tioned due to the sheer amount of land area required1. Further, and the limited availability of arable land challenge the longevity forestation of arable land results into trade-offs with the local of the carbon sink service. food production7,8, and forests tend to have a high water use, so that new plantations decrease soil moisture (SM) and renewable freshwater resources9,10. Results Moderate forest harvesting practices on forested farmlands can Tree cover and CO2 emissions increased between 2002 and generate carbon stocks and at the same time provide economic 2017. Forest cover has increased considerably at the expense of output from timber products11. As such, expansion of managed bare ground and short vegetation (grassland and cropland) in forests provides alternative income for farmers, and an increase in southern China (we studied eight provinces covering China’s the time-averaged standing carbon stock. Most of our current largest forests17,26) from an average tree cover of 21% in 1982 to knowledge of the impact of forest management on climate 27% in 1999 and 38% in 2016 (+0.48% year−1)27 (Fig. 1a; Sup- change, however, is based on models and site-scale evidence, and plementary Table 1). Most of the increase in tree cover and real-world empirical information about the performance of decrease in bare ground/short vegetation occurred after the diverse forest management strategies in relation to carbon storage implementation of forestation policies18 around the year 2000, and SM at larger observational scale is needed12. making this region the largest tropical area of tree cover increase Southern China is an important case region in the sub-tropics, globally in recent time (Fig. 1b). as it hosts a large number of carbon-emitting economic activities Fossil fuel carbon (C) emissions from the eight provinces of and because it experienced intensive land use changes, being now southern China amounted to 0.21 Pg C in 1997 and increased to among the most dynamic areas of managed forests in the 0.63 Pg C in 2012 where after they stabilized28 (Fig. 1c). The world13–17. Since 2000, government-funded efforts to combat carbon sequestration in aboveground vegetation was estimated land degradation and fight poverty have been put in place18,19. using moderate resolution imaging spectroradiometer (MODIS) The main land use decisions have been to expand the forested and SM and ocean salinity (SMOS) low-frequency passive area and reduce agriculture on marginal sloping lands, and to microwave data (see Methods; refs. 23–25). We used a static intensify crop cultivation on more fertile/less erodible soils. benchmark map of aboveground biomass to train a machine Planted trees have generated income19 and woody vegetation learning algorithm applied on annual MODIS imagery to cover has greatly increased18,20. The results of this land use estimate changes of aboveground carbon density from 2002 to transformation are relevant to other subtropical regions that are 201723, the period in which extensive forest management has experiencing fast tree-cover loss and continued expansion of taken place in the region. Aboveground biomass retrieved from croplands in hill slopes21,22. SMOS-passive microwave observations were calibrated by a set of Management practices for forested land in southern China benchmark maps, which25 was used as an independent evaluation include protection of existing (old-growth) forests, recovery of of changes in carbon density calculated from MODIS optical data a b c 100 0.75 75 0.50 Pg C year –1 Cover (%) 50 0.25 25 0.00 0 –0.25 1990 2000 2010 2000 2005 2010 2015 Short vegetation Tree cover Bare ground Tree cover increase > 0.5% year –1 Tree cover decrease > –0.5% year –1 Net C change CO2 emissions Fig. 1 Dynamics in forests and carbon stocks in southern China. a AVHRR estimated27 dynamics in tree cover, and cover fraction of short vegetation, which is grass and cropland, as well as bare ground from 1990 to 2016 in southern China. b The same data source is used to show global increases and decreases in tree cover (larger 0.5% year−1) from 2000 to 2016. The study region is marked with a black rectangle. c Fossil CO2 emissions (2000–2015)30 from the eight provinces in southern China (error bars are standard deviations) compared to net C changes (2002–2017) in aboveground vegetation calculated from MODIS data for the same area (shaded area is the uncertainty based on the RMSE between MODIS- and SMOS-based estimates; see text). Pixel number for eight provinces = 8,573,649. Source data are provided as a Source Data file. 2 NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 ARTICLE a b 100°E 110°E 120°E 54 r 2 = 0.90 C density n=8 Mg ha–1 year–1 MODIS C density 52 (Mg C ha–1) –5 30°N 50 –1 48 0 1 f 46 60 Mg C ha–1 5 55 44 44 46 48 50 52 54 50 20°N SMOS C density (Mg C ha–1) c 0.25 m3 m–3 100°E 110°E 120°E 10 d Soil 0.2 moisture m m–3 year–1 3 0.15 Change (%) 5 2010 2012 2014 2016 –0.04 0 30°N –0.01 0 –5 0.01 –10 0.04 g Mg C ha–1 Rainfall C density Soil moisture 50 20°N e 0 300 600 km 45 0.4 Density Sign. decrease (p < 0.05) m3 m–3 0.4 0.2 0.3 0 2010 2012 2014 2016 –2 0 2 4 6 Mg C ha–1 year–1 Fig. 2 SMOS-based carbon density and soil moisture changes 2010–2017. a Comparison of independent MODIS- and SMOS-based aboveground C density estimates (annual averaged values of the study area). b SMOS aboveground C density trends during 2010–2017. c Changes (in %) of CHIRPS rainfall, SMOS C density and SMOS soil moisture during 2010–2017. Boxplots show whiskers, quartiles and median across 25 × 25 km2 pixels. d Spatial distribution of SMOS soil moisture trends during 2010–2017 (n = 2279 pixels). Transparent areas have been disregarded for insufficient data quality. Significant negative trends (p < 0.05) are shown in purple. e Density plot showing the distribution of SMOS aboveground C density trends in the purple areas of d with a significant (p < 0.05) decrease in soil moisture (n = 172). f Example of an area with a significant (p < 0.05) decrease in soil moisture and a strong increase in C density. g Example from Guangxi province showing an area with large-scale tree plantation and massive harvest in 2013 and a significant (p < 0.05) decrease in soil moisture. after the year 2010. The SMOS instrument is able to sense the were predominantly positive in these areas (Fig. 2e, f). A hotspot entire vegetation layer without saturation in densely vegetated of decreased SM was found in Guangxi, which is known for areas. The MODIS estimate of the carbon sink was 0.11 (±0.05 extensive Eucalyptus sp. plantations30. Here, the example shown temporal uncertainty quantified as root mean square error in Fig. 2g points towards a large-scale harvest (followed by (RMSE) between MODIS and SMOS) Pg C year−1 (0.61 ± regrowth) in 2013. 0.2 Mg C ha−1 year−1) for 2002 to 2017, which represents an To further study the impact of different management types on amount of carbon equivalent to 20% of the regional annual fossil C sinks, we stratified land and forest management over the whole CO2 emissions from 2002 to 2015, or 33% of the emissions after region into eight different types (Table 1), including dense/ 2012, a year after which the emissions have stagnated (Fig. 1c). A minimally managed forest, persistent non-forest, areas with fast/ single drought year, however, can offset this annual mitigation of slow forest gains, deforested areas and two intensities of forest emissions, as for example in 2011, where the ratio of carbon sink rotation (Table 1; Supplementary Figs. 1–6). The classification is via forest growth to fossil CO2 emissions dropped close to zero. based on mapping of the duration, magnitude and direction of SMOS low-frequency passive microwave data (available since human-induced disturbances (Supplementary Fig. 7) from annual 2010 at a resolution of 25 × 25 km2) showed a high agreement of forest probability time-series data (2002–2017, 500 × 500 m2 aboveground carbon changes with the MODIS estimates, both resolution), trained and compared with a full area coverage of spatially (r = 0.7) and temporally (r = 0.9) (Fig. 2a, b). The spatial 2-m GF-1 satellite imagery. The association of the different land distribution of SM trends derived from SMOS observations (fully use types with carbon sequestration is shown in Fig. 3, and the independent of the biomass estimations, due to the multi-angular spatial distributions are shown in Fig. 4a, b (close-ups are shown and dual polarization of the sensor) is different from trends in the in Supplementary Figs. 5, 6). Due to their coarse spatial resolution aboveground biomass C sink during 2010–2017 (Fig. 2c, d). The (25 × 25 km2), biomass carbon change data from SMOS could not SMOS SM data showed an overall decrease, while rainfall29 and be downscaled into the different land use types, but they were biomass C density increased (Fig. 2c). Statistically significant (p < used as independent check of the spatially detailed MODIS-based 0.05) decreases in SM were found in 8% of the area (n = 172 grids estimates at 500 m resolution. Similarly, SM trends were assessed with 25 × 25 km2); however, the SMOS-based C density trends at 25 km resolution. NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 Table 1 Types of land use and land use change in southern China. Type Management Tree plantation/ Tree cover Forest loss Mean C density Net C sink (Pg C year−1)/ harvest intensitya (%)b (%)b (Mg C ha−1)c % contributionc 1 Dense forest Persistent, often protected old Very low/very low 56 ± 14 2.8 ± 7 105 ± 11 0.005/4% forests 2 Forest Persistent, semi-managed Low/low 44 ± 16 3±8 75 ± 16 0.023/20% 3 Non-forest Persistent, for example, farmland, Low/low 14 ± 12 0.7 ± 3.4 22 ± 16 0.032/28% sugarcane and fruit trees 4 Recovery Non-forest to forest (slow) Medium/none 28 ± 15 0.5 ± 2.5 47 ± 12 0.016/14% 5 Afforestation Non-forest to forest (fast) Very high/very low 27 ± 14 0.8 ± 3.3 44 ± 12 0.021/18% 6 Deforestation Forest to non-forest None/very high 40 ± 16 12 ± 17 64 ± 13 −0.00007/−0.06% 7 Rotation Medium-scale forestry, changes High/high 32 ± 16 4 ± 10 53 ± 11 0.017/15% between forest and non-forest 8 RotationL Large-scale forestry, low recovery Low/very high 37 ± 15 9 ± 15 61 ± 11 0.0008/0.8% Dense forest, forest and non-forest are land use types with minor disturbances. The other types include land that experienced changes between forest and non-forest. The percentage of tree cover in each land use type (±standard deviation) is from a Landsat-based tree-cover map in 201026. Forest loss is also derived from26 and represents the average forest loss per 500-m grid (±standard deviation) from 2000 to 2017. Mean C density and net C sequestration are from our MODIS-based estimates from 2002 to 2017 (±RMSE from calibration biomass data as measure for spatial uncertainty, Supplementary Table 3). See Supplementary Figs. 1–6 for illustrations and Supplementary Table 2 for more details aSee Supplementary Fig. 7 bTaken from ref. 26 cMODIS estimates 2002–2017 (see Methods) a C stocks b C density c CO2 emission sequestration 5 0.03 1: Dense-forest 1.5 4 2: Forest Mg C ha–1 y–1 3 5 0.02 3: Non-forest Percent y–1 1 Pg C y–1 4: Recovery Pg C 2 5: Afforestation 1 0.5 2.5 0.01 6: Deforestation 0 7: Rotation –1 0 8: RotationL 0 0.00 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 Type of land use Type of land use Type of land use Stock 2002 Gross loss Gross gain Eight provinces China Fig. 3 Dynamics of carbon stocks for different types of land use. a C stocks estimated from MODIS data are shown for the first year of study in 2002 in grey, with consecutive gross gains and losses for 2002–2017 in green and orange, respectively. b Net changes in C density (2002–2017) from MODIS data. Error bars are the temporal uncertainty based on the RMSE with SMOS estimates. c Annual increase in C stocks for each land use type in Pg C year−1 on the right-hand axis and as fraction of fossil CO2 emissions for the eight provinces and for China on the left-hand axis (2002–2015). Pixel number for eight provinces = 8,573,649. Source data are provided as a Source Data file. Dense forests and persistent forests. Dense forests are protected Recovering forest and afforestation areas. If an area changed remnants of natural forests or dense secondary unharvested from non-forest to forest, the velocity (the ratio of magnitude to forests with a dense tree cover and no disturbances (line 1 in duration) of the change determined the type of forestation (lines 4 Table 1). Dense forests covered only 8.8% of the region, but and 5 in Table 1). The recovery type is associated with a slow stored 20.5% of the total aboveground C stocks (1.71 Pg C) in the increase in forest cover, likely caused by natural regeneration or beginning of the study period (2002) due to their high C density planted slow-growing species, and covered 5.4% of the area. In (105 Mg C ha−1 from 2002 to 2017). C dynamics in dense forests contrast, the afforestation type was associated with a rapid remained small during 2002–2017, and contributed 4% (0.005 Pg increase in forest area and cover, suggesting plantations of fast- C year−1) to the region’s net C uptake, which removed 0.9% of growing species (7.4% of the area). Recovery and afforestation did the annual fossil CO2 emissions (Fig. 3a–c). not include large-scale harvesting, but often included forestation Persistent forests (20.5% of the area) are less dense forests with of former agricultural land (Fig. 4c). minor disturbances, but no major harvests during the observation The C stocks of these two types of forestation land use were period (line 2 in Table 1). This land use type stored 2.76 Pg C low in 2002 (recovery: 0.37 Pg C; afforestation 0.48 Pg C) (Fig. 3a), (33%) in 2002 (mean C density of 75 Mg C ha−1) and showed a but increased rapidly from 2002 to 2017 (recovery: 0.016 Pg C substantial stock increase (0.023 Pg C year−1) from 2002 to 2017 year−1; afforestation: 0.021 Pg C year−1) (Fig. 3a, b), and (Fig. 3a), contributing 20% to the region’s C sequestration, contributed 14% (recovery) and 18% (afforestation) to the equivalent to 4.2% of the provincial CO2 emissions (Fig. 3c). region’s C sequestration, thus storing on average 7.5% of the annual industrial CO2 emissions (Fig. 3c). Persistent non-forest land. Persistent non-forests are typically farmland or grassland, and include small forest patches and fruit trees (line 3 in Table 1). Non-forests covered 43.8% of the region Forest extraction. We identified three forest land use types and contained 1.61 Pg C (in 2002), despite a mean C density involving timber extraction at large scales (lines 6, 7 and 8 in of 22 Mg C ha−1. From 2002 to 2017, the C uptake was Table 1). The intensity of both the forest removal and replanting 0.032 Pg C year−1, caused by the large area occupied (Fig. 3a, b). defines each type. Conventional (short-)rotation forestry involves Net changes in C contributed 28% to the region’s C sequestration, clear-cutting of mature stands and subsequent regrowth (i.e., storing 5.2% of the annual CO2 emissions (Fig. 3c). dynamic changes between forest and non-forest) and occurred in 4 NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 ARTICLE a 100°E 110°E 120°E b 106°E 106°E 30°N 25°N 0 10 20 km c 0 250 500 km Dense forest Afforestation Forest Deforestation 24.9°N Non-forest Rotation Recovery RotationL 104.7°E 0 10 20 km d e 12.5 50 CO2 sequestration (%) 10.0 7.5 Pg C 25 5.0 0 2.5 –25 0.0 2004 2008 2012 2016 2005 2010 2015 f g 50 Timber harvest [Pg C] 0.012 40 Areal ratio (%) 0.008 30 20 0.004 10 0.000 0 2002 2006 2010 2014 Types of land use Fig. 4 Land use in southern China. a The area classified into the eight types of land use. b Forestry is particularly common in northern Guangxi. This close- up is shown as a black rectangle in a. c Forestation of farmlands on sloping hills is a major reason for the increase in forests. This example image (GF-1, 2-m resolution) is also from northern Guangxi (not the same area as b). d Contributions of the types of land use to store fossil CO2 emissions for the eight provinces in southern China (considering net C changes) from 2002 to 2015. e C stocks for the types of land use from 2002 to 2017. f Annual timber harvest (from forestry statistics at province level; n = 8) converted to Pg C. g Areal ratios of the types of land use. See Supplementary Figs. 5, 6 for close-up illustrations. Source data are provided as a Source Data file. 10.6% of the area. Harvesting in rotation areas occurred at a rotationL and deforestation areas acted as a C source (Fig. 3c) rather small spatial scale. Areas that experienced larger-scale from 2002 to 2015. harvesting without adequate regrowth were termed rotationL Forestry statistics at province level showed that on average (1.9% of the area). Forest removal without any regrowth during 0.008 Pg C were harvested for timber production each year, which 2002–2017 was defined as deforestation (1.6% of the area). The C amounts to 46% of the average annual net changes (0.017 Pg year−1) stocks of rotation, rotationL and deforestation areas were 0.99, from rotation areas (Fig. 4f). Before 2009, the harvest was higher 0.22 and 0.2 Pg C in 2002, respectively, with a positive net balance than the net C sequestration, which however reversed in the later in rotation and rotationL areas (0.017 and 0.0009 Pg C year−1, years. Timber harvest accounted for 0.004 Pg C in 2003 and 0.012 Pg respectively) and a negative net balance in deforestation areas C in 2017. (−0.0001 Pg C year−1) between 2002–2017 (Fig. 3a). It has to be noted that the C sink here is calculated from forest growth without considering the full life cycle of extracted wood. It also Discussion ignores changes in litter, coarse woody debris and soil carbon. Within the past 20 years the subtropical and mountainous These three land use types with forest extraction make a con- landscapes in southern China have transitioned from traditional tribution of 15%, 0.8% and –0.06% to the regional C sequestration agriculture towards managed forests mainly dedicated to the and represent an amount of carbon sequestered equivalent to production of wood products. Our study showed that newly 2.7 % of regional fossil CO2 emissions in rotation areas, whereas planted forests considerably increased C sequestration while more NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 than tripling timber harvest. The forests that were already present considered (life-cycle analysis). If the wood enters a short-life at the beginning of the study period contributed only marginally product such as paper and swiftly progresses to decomposition or to the observed large increase in C sequestration, while the combustion, the C sequestration in the forest is counterbalanced plantation of non-forest areas had the highest contribution. at best or represents a C source when accounting for the C Although unharvested forest plantations showed the highest C emissions during processing of the wood43,44. If the wood is used uptake, harvested forests also had a strongly positive impact on C for construction or other long-term uses, the standing C stock is stocks, as long as harvest was followed by forest regrowth. extended for several decades45,46. Finally, using the biomass for A decline in SM in areas of intensive forestation and increased bioenergy production while ensuring regrowth following harvest C stocks has been shown by model simulations10, and the SMOS could also help mitigate climate change. satellite data indicate a similar development in several areas of Uncertainties in the present study are mainly related to the southern China. A possible reason is the increased water demand spatial resolution of the data sets and the benchmark maps of tree plantations and the associated increased evapotranspira- used to calibrate the carbon density maps. The resolution of 25 × tion. Although products on actual evapotranspiration do not 25 km2 of the SMOS data is too coarse to distinguish different agree on the magnitude of increase for the study area (Supple- types of land use. The resolution of 500 × 500 m2 of the MODIS mentary Fig. 8), the data support that forest plantations increase data also represents mixtures of different processes, since neither evapotranspiration, while lowering the albedo and the land- forest expansion/harvesting nor stand increments or degradation surface temperature31,32. However, while hotspot areas of usually cover the full extent of a MODIS grid cell. Frequent decreased SM were found to overlay extensive commercial tree changes between tree planting and harvesting at small scales pose plantations (e.g., in Guangxi), statistically significant (p < 0.05) challenges for accurate classification, which tends to primarily decreases in SM were limited to small parts of the area. Moreover, reflect major disturbances measurable at larger scales. Finally, C droughts (e.g., in 2011) may have increased the region’s evapo- dynamics estimated with optical data are restricted to the upper transpiration, with implications for the water use efficiency (gross green canopy layer, and cannot directly measure the non-green primary productivity/evapotranspiration), which was found to wood part of the vegetation. Larger uncertainties are however decrease for southern China over the period 2000–201133. The mainly found in tropical rainforests47 with a low spectral varia- continued monitoring of both C stocks and SM dynamics is thus bility both in the spatial and temporal domain. This is less of an essential to understand the impact of land use policies and issue in southern China, where the vegetation is less dense, and droughts on ecosystems. Here, our study demonstrated the rapid canopy changes occur, enabling our machine learning applicability of low-frequency passive microwaves as a tool ful- model that translates MODIS optical data into biomass to be filling these demands34,35. trained by observations with a large dynamic range. Moreover, no While 30% of the region is covered by protected forests or saturation was found in our carbon estimates (Supplementary forests without major active management during the study per- Fig. 10), and the relatively good agreement of our assessment iod, 27% is actively managed as afforestation or from various based on MODIS and high-resolution land cover data with forest extractive uses. Those managed forests contributed 47% to independent data sets on forest dynamics28,29 and biomass the region’s carbon sequestration, while generating income for the (SMOS)27 provide confidence in our estimates (see Methods). local population through the marketing of wood, other timber The numbers on C stock increases presented (0.11 Pg C year−1 products and fruits19,36. About one-third of the area’s fossil CO2 for southern China) are considerably higher than previous esti- emissions were stored by aboveground vegetation biomass mates, also considering that our numbers are only based on increase during the past 6 years, and based on this number, it can aboveground vegetation biomass. For example, Fang et al.48 used be extrapolated that an additional forestation area of ~3 million inventory data to estimate the C sink of forests in the entire China km2 would be required to reach net zero emissions (100%). This at 0.075 Pg C year−1 for the 1980s and 1990s. Piao et al.38 found is impossible, since only ~1 million km2 are currently non- the average terrestrial carbon sink of all China (excluding crop- forested, of which 51% are farmlands37. Areas that were not lands) to be 0.177 Pg C year−1 for the same period, of which 58% classified as forest also played an important role in our assess- can be attributed to vegetation biomass. The period of analyses of ment, as the general C uptake in non-forested areas was high. above-mentioned studies ends, however, around 2003. Here we However, the annual C turnover and losses in drought years were showed that new forests, which were not present in 2003, con- considerable in non-forest areas (with only little impact on for- tributed with 47% to the region’s C sequestration. Moreover, non- ests), suggesting limited C sequestration abilities of these areas if forested areas, commonly not considered as a sustainable C sink, tree cover is not increased38–40. contributed further 28%, whereas only 24% of the C sink could be The spatial extent of dense forests (largely old forests) was attributed to forests, which were already present in 2003, small, yet our study showed that those forests were an important explaining the large differences to previous studies. and stable carbon stock that needs to be preserved. A total of While the increase in C stocks over the period studied is sig- 1.71 Pg C (20.5% of the total C stock) would be lost if all dense nificant, in absolute terms as well as relative to the CO2 emissions forests in southern China were cut, corresponding to ~9 years of of the eight provinces, it is evident that further enhancement of fossil CO2 emissions. However, the limited capacity of dense this sink in the future is unlikely. The sloping lands that offer the forests to store fossil fuel CO2 emissions highlights an inherent best opportunity for new forest plantations replacing low pro- dilemma of forest management in relation to mitigating climate ductivity croplands are now to a great extent already converted change; undisturbed forests maximize the standing C stock and and further increase in C storage will require improved man- provide a permanent reduction of the atmospheric C content, but agement of currently planted areas to reach a higher C density. once these forests reach a climax state, the net aboveground C Whether this may be seen as a viable option for the region is sequestration is low41,42. Managed forests store less C in the determined by a range of biophysical, economic and political standing biomass and have a smaller potential as permanent C factors. Further conversions from low to high C density forest stock, but sequester substantial amounts of C annually18,42. To types are unlikely, unless economic incentives (e.g., in the form of evaluate the potential of managed forest for mitigating climate compensation payments, high market prices, or ‘Payments for change, two aspects are critical: first, the harvested areas should Environmental Services’) are provided. Thus, future sustainable be systematically replanted to sustain long-term carbon uptake, development will notably depend on the extent to which man- and second, the subsequent use of harvested C should be agement strategies and economic factors will be able to combine 6 NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 ARTICLE substantial and increasing C stocks with a favourable economic indicated the likeliness that an area showed characteristics similar to those of the output. This requires a more refined analysis including different forest training points, implying that temporally increasing forest probability represented the growth of forests50. Likewise, a decrease in forest probability was rotation cycles, plantation types (that is species and harvest linked with harvesting. The assumption that forest probability was associated with purpose) and price/demand scenarios. the fraction of a grid covered by forest was evaluated with independent data (see Although our study found a large increase in C sequestration of below). which about half was caused by newly planted and managed forests, our numbers also show that storing fossil CO2 emissions from forest plantations alone is not a viable strategy for miti- Types of land use. We applied Landsat-based Detection of Trends in Disturbance and Recovery (LandTrendr) to the maps of annual forest probability to identify the gating anthropogenic carbon emissions, due to the sheer amount types of land uses51. We used similar settings as suggested by refs. 52,50 because of of land areas needed. It should be kept in mind, however, that a comparable data sources and time periods. We set the significance threshold for the substantial part of the industrial CO2 emissions in southern model fitting to 0.05 and the maximum number of trajectory segments to six. The China are caused by the production of goods that are not con- LandTrendr algorithm is used for partitioning of the annual data into fitted seg- ments. A segment ends if a major disturbance is detected. The outputs of the sumed in this area, but are likely nationally and globally traded. algorithm are the number and duration of the segments, and the magnitude of the Additionally, large-scale forestation programmes can have an disturbances. We found no evidence of natural causes of the disturbances to forest adverse impact on the ecosystem by decreasing SM and thus probability (only 7% of the disturbances were detected in 2009–2011, when the water availability for local organisms, showing that carbon most severe drought of the past century occurred), suggesting that the disturbances were due to human management. We used the LandTrendr output data to calculate sequestration and economic benefits should not be the only two indices: managed forest increase and managed forest decrease (i.e., forest aspects considered when debating the sustainability of land and increases/decreases caused by human management), by dividing the magnitudes of forest management. the strongest positive/negative changes by the durations of the changes (Supple- mentary Fig. 7). The indices ranged from 0 to 1, with high/low values implying rapid and strong/slow and weak changes. A high value in the managed forest Methods increase map was thus interpreted as large-scale plantations of fast-growing trees, Overview. This study estimated dynamics in forested areas and C density at a scale and a high value in the managed forest decrease map was interpreted as a clearcut. of 500 × 500 m2 for 2002–2017 (see flowchart in Supplementary Fig. 9). We pro- The indices were classified into several intervals: 0–0.2 (very low), 0.2–0.4 (low), duced annual forest maps using MODIS, which were compared with independent 0.4–0.6 (intermediate), >0.6 (high). tree-cover data at scales of 30 m28 and 5.6 km27. We derived a typology of forest We subsequently used the maps of annual (2002–2017) forest probability to continuity and dynamics from the annual forest maps. Since forest management in derive nine classes of forest dynamics. The forest classes included (1) dense forests this area is a very dynamic process, the typology requires annual maps at a rea- (forest probability ≥0.8), (2) forests (≥0.5) and (3) non-forest (<0.5). The threshold sonable spatial scale, which is not provided by any publicly available product. As number for dense forests is based on Wang et al. 50, who used field observations to analyses were conducted at a spatial scale of 0.25 km2, the land-surface types identify old forests (above 0.8). Random Forest commonly uses a probability inevitably included cases of multiple processes and the final classification was threshold of 0.5 to distinguish if a pixel belongs to a class or not. We used this determined by the majority of these processes within a grid cell, implying that the threshold to define if an area belongs to the forest (probability above 0.5) or non- data likely over- or underestimated the real areal extents, particularly for forest forest type. The classes dense forest, forest and non-forest were assigned if the dynamics. For example, an area of 500 × 500 m2 reported to belong to the affor- probability scores remained within these intervals during the full period; if not, the estation type means that afforestation was the major process occurring in this grid type of forest dynamic was determined by the threshold between forest and non- cell, but does not imply that trees had been planted on the entire area of 0.25 km². forest (0.5). Class (4) represents forest increases changing from <0.5 in 2002 to >0.5 in 2017 without negative disturbances, and class (5) represents the same but with MODIS satellite data. Daily BRDF (bidirectional reflectance distribution func- disturbances. Class (6) represents forest decreases when an area changed from >0.5 tion) corrected MODIS imagery (MCD43A4 collection 6)49 from 2000 to 2017 was in 2002 to <0.5 in 2017 without positive disturbances (e.g., regrowth), and class (7) the basis for this study. We used all seven bands. Cloud cover is severely impacting represents the same but with disturbances. Class (8) represents an area fluctuating the number of good quality observations from optical remote sensing sensors in between forest and non-forest, but remaining non-forest (2002 compared to 2017), southern China, and we chose the annual median of daily MCD43A4 images to and class (9) represents the same but for remaining forest. retrieve an annual image. Due to the fact that even the annual median images were Finally, a grid cell is characterized by three values determining the final land use not free of noise, we additionally applied a 3 years moving median window to the type: (a) the strongest managed forest increase, (b) the strongest managed forest annual time series, shortening the available period of data to 2002–2017. This decrease and (c) the forest dynamic (nine possibilities, previous paragraph). The aggressive cleaning was necessary because all kind of noise would be mapped as combination of these three values determined the final type of land use of an area. “clearcut” or other forest dynamic, which would alter our results (see later We also determined the number and duration of the segments for each possible sections). combination to include the velocity and duration of tree growth before harvest to potentially gain insight into the type of tree species (slow- and fast-growing species) (Supplementary Table 2). The final eight types (Table 1; Supplementary Table 2; Forest probability. We created annual forest/non-forest probability maps from Supplementary Figs. 2–6) included: (1) dense forests (probability always ≥0.8, no annual MODIS data. Here we used forest/non-forest training points manually major disturbances), (2) persistent forests (probability always ≥0.5, no major selected from ~10,000 GF-1 satellite images (pansharpened true-colour composites disturbances) and (3) persistent non-forests (probability always <0.5, no major at 2 × 2 m2 for 2013–2017) to train a Random Forest classifier (using the standard disturbances), which are continuous types with no major human forest setting of 500 trees) with MODIS MCD43A4 reflectance (seven bands) and a management (managed forest index values <0.4) and long-lasting segments Shuttle Radar Topography Mission (SRTM) digital elevation model for predicting (a segment is a period without disturbance). These three types do not exclude the annual forest probability. “Probability” is the standard output of a Random planting and harvesting in forest/non-forest areas, but these activities were minor. Forest classification and shows how likely a pixel belongs to a certain class. If only (4) The recovery class implies a low but positive development in the managed two classes are available (here forest and non-forest), the probability shows the forest index (i.e., no major disturbance, a slow forest increase, no decrease, long- likeliness if a pixel belongs to the forest (1) or non-forest (0) class. A total of 110 lasting segments) and changes from non-forest to forest. These areas typically had forest and 323 non-forest training points were manually selected from the GF-1 greatly reduced human disturbance throughout the period of analysis that allowed images, supported by Google Earth. The manual selection of the points followed the tree cover to gradually increase. This type does not exclude either tree planting strict rules: First, the areas should be clearly interpretable in the GF-1 imagery. or harvesting in a given grid cell, but the majority of cells had slow and steady Second, the areas should be as spatially uniform as possible, that is, either a dense increases in tree cover without disturbance. (5) Afforestation represents tree forest or clearly no forest. Third, we used the 500 × 500 m2 grid from the MODIS plantations of medium- to fast-growing species and changed rapidly from non- resolution to ensure that a minimum of nine pixels (the point was set in the centre forest to forest. No major decreases were detected (low managed forest index). pixel) are covered by the class. Finally, the forest training points should mainly (6) Deforestation represents areas changed from forest to non-forest, with all represent dense forests that were stable during 2002–2017 (historical Google Earth segments being negative. (7) Rotation represents areas with frequent changes from images were used when possible), so we used the average of all years as input for non-forest to forest (and vice versa), without the requirement of a particular type as the model for each of the seven bands, which further reduced noise and guaranteed starting and ending years. This type mainly represents the rotation of planting and a stable model50. A three-fold cross-validation of the model (excluding random harvesting activities, which could vary in velocity and magnitude. The average subsets of the training data that were then predicted) was satisfactory (r = 0.93). duration of a segment in this type was 4.6 ± 1.3 year−1, indicating the typical However, since the chosen points were rather homogeneous, we applied compar- growing cycle of a forest patch53. (8) RotationL represents fast and rapid decreases isons with independent forest cover maps to assess the quality of the forest in forest probability and changes from forest to non-forest, characterizing areas probability maps (see section Intercomparison of independent data sets). We with large-scale harvesting. The difference between the types of rotationL and applied this model to the annual MODIS bands to generate annual probability deforestation is that segments of positive slopes can occur in rotationL areas; maps of the forest/non-forest classification. The probability ranged from 0 to 1 and however, regrowth does not reach the state of forest (probability <0.5). NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 Carbon density. We used a static unpublished global benchmark map of C density types were very low (Table 1), and only the rotation and deforestation types of woody vegetation for 2015 to study forest carbon changes in China54,55. The showed substantial forest losses (Supplementary Fig. 11). map was constructed using three types of data: (1) more than 8 million waveforms Finally, we compared the MODIS-derived maps of C density with the SMOS from the Geoscience Laser Altimeter System (GLAS) onboard ICESat-1 (Ice, Cloud L-VOD (L-band vegetation optical depth) data set (currently being the only and land Elevation Satellite) that reflect the vertical structure of forest and woody spatially explicit data set on dynamics of aboveground biomass C)24,25 as derived vegetation using a semi-systematic sampling approach, (2) a large number of from low-frequency passive microwaves. This data set is available at a resolution of ground data from national inventory and research biomass plots to develop models 25 × 25 km2 for 2010–2017, and the conversion was done using both Saatchi’s and for converting the GLAS lidar data to C density for different forest types globally Baccini’s biomass benchmark maps following25, and taking the average of both and (3) satellite imagery from Landsat, ALOS Phased Array L-band Synthetic calibrations. At the pixel level, the spatial correlation between the MODIS Aperture Radar data for 2015, and land elevation data from SRTM. Using GLAS- C density and the L-VOD C density was satisfactory (r = 0.7). For temporal predicted biomass as the dependent variable, satellite imagery as independent comparisons, we averaged for each year both the MODIS C density and L-VOD variables and a Bayesian maximum entropy algorithm54,56, C density and its maps over the region and the temporal dynamics for both data sets agreed well uncertainty were predicted and were mapped globally at 1-ha spatial grid cells. The (r = 0.90), providing confidence in the dynamics of C density. The uncertainty in machine learning algorithm uses 8 million GLAS derived C density as training data the temporal dynamics was assessed by the RMSE between MODIS- and SMOS- to predict the probability of each pixel falling into a range of C density. The derived changes in C density. For an uncertainty assessment based on different probability maps are then combined to calculate the mean and variance of C benchmark calibration maps, we refer to ref. 25. Due to the high resolution required density for each pixel54. Details of the prediction algorithm and comparisons with to distinguish different forest types, the temporal uncertainty could only be other machine learning approaches have been discussed in earlier publications56. provided for the study area as a whole, and the uncertainty of the different forest The map is an improvement over earlier pan-tropical maps by including recent and types was limited to the spatial level. Here we calculated the RMSE between the advanced remote sensing data, and new models that include variations of wood benchmark map used for calibration and the MODIS-derived C density maps for density of trees across tropical regions. each forest type. The uncertainty of the map at each grid cell included uncertainty associated The negative trend found in SMOS SM is contrasting other studies showing a with C density models and spatial correlation54,57. Although, the previous maps positive trend in SM this period24. However, whereas the SMOS satellite is able to were criticized for dense tropical forests due to use of optical satellite data and the retrieve SM data in densely vegetated areas59, traditional SM products based on impact of the wood density variations58, dense rainforests are not present in our high-frequency passive and active microwave data do not agree in China60 and study region and no saturation between our data and the benchmark map was often fail in densely vegetated areas59. Model simulations and field measurements observe. While normalized difference vegetation index (NDVI) clearly saturates instead confirm the negative trend in SM10. over forests (Supplementary Fig. 10), satellite imagery from higher-resolution radar sensors provide improved sensitivity to forest structure and biomass across CO2 emissions. The data set from ref. 28 contains fossil CO2 emissions for each regenerating and lower C density forests across the study region. province for 1997–2015. We used the sectoral approach and converted CO2 to C by The C density map was used to train a model of gradient boosted decision trees dividing the value by 3.67. with MODIS (MCD43A4 7 bands; NDII, EVI2, MCD43A3 shortwave albedo) and an SRTM elevation model, following23. Model training, prediction and validation were done only using data covering the study area. We used 50% of the available Forestry data. The timber production data were collected from the Chinese pixels (n = 286,824) to train the model, and the remaining 50% were used for Forestry Statistical Yearbooks, including the timber volumes from 2002 to 2017 for validation. The quality of the model was generally high (r = 0.86, slope = 0.74, each province. We applied the accumulation method to estimate timber carbon RMSE = 11 Mg C ha−1), and the C density of the original map could be based on the volume of timber. The main species are pine, cedar and eucalyptus, reproduced (estimates and spatial uncertainty for each class are presented in and we used the average wood basic density of these three species to calculate the Supplementary Table 3) without saturation over forests (Supplementary Fig. 10). biomass of the timber production, which was then converted to carbon by mul- The model was applied to the data sets for each year to derive maps of annual tiplying the biomass with 0.5. C density. The dynamics were compared with low-frequency passive microwave satellite data (next section), which also includes an estimation of the temporal SMOS microwave data. L-VOD was used as independent data set to assess the uncertainties. dynamics of C density, and SM was used to monitor the patterns of SM dynamics. Both L-VOD and SM products were retrieved simultaneously from the SMOS observations35. The use of multi-angular and dual-polarization SMOS observation Intercomparison of independent data sets. Since no field data on land use types ensures that L-VOD and SM are independent of each other. was available, we conducted a visual comparison between the classification and the The SMOS L-VOD and SM data sets used here were produced using the SMOS-IC full areal coverage of GF-1 satellite imagery at a resolution of 2 m (~10,000 images; algorithm (version 105) at a resolution of 25 × 25 km2 for 2010–201761–63. This see Supplementary Figs. 2–4). However, no statistics could be derived from this algorithm was selected here, as it does not rely on ancillary data such as modelled SM comparison, thereby requiring additional result comparisons to be taken based on and vegetation optical indices (such as LAI or NDVI), making the SMOS-IC product independent data sets. The first step was a comparison with the Landsat-based tree- more robust for applications monitoring both biomass carbon and SM dynamics62. cover map26 for 2010, available at a resolution of 30 m (Supplementary Fig. 1). The SMOS-IC data were processed as follows: first, scene flags were applied and map was aggregated to 500 m to show the average tree cover of a 500 m pixel from both the L-VOD and SM data for ascending (ASC) and descending (DESC) orbits our analysis. The agreement between the aggregated tree-cover map and our forest were filtered according to flags for strong topography, frozen soils (soil temperature probability map was very high (r = 0.88, n = 8,796,179), confirming that our <273.5 K), urban areas and water bodies. We then applied a filter for the effects of probability maps capture well the average tree cover per grid cell. radio frequency interference (RFI) by excluding data with a TB-RMSE index higher We then used the annual long-term data set from ref. 27 for 1982–2016 to than 10 K. We further applied specific filtering to the L-VOD data: a moving provide an overview of the long-term dynamics and to compare with our average smoothing (window size = 30 days) was applied on the filtered ASC and classification of forest dynamics (Supplementary Fig. 1). The data set is based on DESC L-VOD time series. Outliers, corresponding to data lower (higher) than the long-term data record advanced very high-resolution radiometer (AVHRR) 10th (90th) percentile of residues (defined here as the differences between raw imagery and contains the fractions of tree cover, short vegetation and bare ground L-VOD data and smoothed L-VOD data) were excluded. To further filter out RFI at a resolution of 0.05°. The average tree cover (2001–2016) spatially agreed well effects, the filtered ASC and DESC L-VOD data were combined into one L-VOD with our forest probability map (aggregated to 0.05° and averaged over the same data set by keeping only 30 observations with the lowest TB-RMSE values for period) (r = 0.81, n = 111,954), confirming that our probability maps could also each year. indicate the fraction of a given grid cell with forests (the higher the probability, the higher the forested fraction). A change in the probability could thus reflect a fractional increase of forest cover within the grid cell. We then compared the Data availability temporal dynamics of the tree-cover maps with our classification for the MODIS MCD43 is available from Google Earth Engine. Tree-cover data from Landsat overlapping period (2002–2016). The slopes (increase in forest cover per year) for can be downloaded at https://earthenginepartners.appspot.com/science-2013-global- the forest types (dense forest, forest) (+0.60% year−1), non-forest (non-forest) forest/download_v1.2.html. Tree cover from AVHRR (1982–2016) can be downloaded at (+0.30% year−1) and forest increase (+0.75% year−1) (afforestation and recovery) https://search.earthdata.nasa.gov. CO2 emission data can be downloaded from https:// clearly differed (Supplementary Table 1). The forest types (dense forest, forest) had www.nature.com/articles/sdata2017201. SMOS data are available from https://www. high fractions of forest cover in 1999 (44%) before start of the MODIS time series catds.fr/. The C density map is available from S.S. The forest management classification (as used here). The fraction of forest cover for the forest-increase types started at and carbon density maps will be made available upon acceptance of the manuscript. 26% in 1999, but increased to 41% in 2016. The forest cover for the non-forest type was 14% in 1999 and remained at 21% in 2016. Deforestation (and replanting) typically occurs at a fine spatial scale, so we Code availability The LandTrendr code can be downloaded at http://emapr.ceoas.oregonstate.edu/tools. applied the Landsat-based forest-loss map from26 identifying forest losses from 2000 to 2014. We aggregated the 30-m forest-loss map to 500 m and compared the html. The IDL version, which was used for this study, is available at https://github.com/ average loss rate (per 500-m grid) with the land use types defined here. The forest KennedyResearch/LandTrendr-2012. Random Forest was run in a GRASS GIS losses in the forest, non-forest and forest-increase (afforestation and recovery) environment with the r.randomforest package. 8 NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 ARTICLE Received: 18 July 2019; Accepted: 26 November 2019; 34. Rodríguez-Fernández, N. J. et al. An evaluation of SMOS L-band vegetation optical depth (L-VOD) data sets: high sensitivity of L-VOD to above-ground biomass in Africa. Biogeosciences 15, 4627–4645 (2018). 35. Kerr, Y. H. et al. The SMOS soil moisture retrieval algorithm. IEEE Trans. Geosci. Remote Sens. 50, 1384–1403 (2012). 36. Wu, T. et al. Carbon sequestration by fruit trees—Chinese apple orchards References as an example. PLoS ONE 7, https://doi.org/10.1371/journal.pone.0038883 1. Baldocchi, D. & Penuelas, J. The physics and ecology of mining carbon dioxide (2012). from the atmosphere by ecosystems. Glob. Change Biol. 25, 1191–1197 (2019). 37. Hollmann, R. et al. The ESA climate change initiative: satellite data records for 2. Huang, K. et al. Enhanced peak growth of global vegetation and its key essential climate variables. Bull. Am. Meteorol. Soc. 94, 1541–1552 (2013). mechanisms. Nat. Ecol. Evol. 2, 1897 (2018). 38. Piao, S. et al. The carbon balance of terrestrial ecosystems in China. Nature 3. Zhu, Z. et al. Greening of the Earth and its drivers. Nat. Clim. Change 6, 458, 1009–1013 (2009). 791–795 (2016). 39. Hua, F. et al. Tree plantations displacing native forests: the nature and drivers 4. Keenan, T. F. et al. Recent pause in the growth rate of atmospheric CO2 due to of apparent forest recovery on former croplands in Southwestern China from enhanced terrestrial carbon uptake. Nat. Commun. 7, ncomms13428 (2016). 2000 to 2015. Biol. Conserv. 222, 113–124 (2018). 5. Chen, C. et al. China and India lead in greening of the world through land-use 40. Jia, B., Wang, Y. & Xie, Z. Responses of the terrestrial carbon cycle to drought management. Nat. Sustain. 2, 122 (2019). over China: modeling sensitivities of the interactive nitrogen and dynamic 6. Erb, K.-H. et al. Unexpectedly large impact of forest management and grazing vegetation. Ecol. Model. 368, 52–68 (2018). on global vegetation biomass. Nature 553, 73–76 (2018). 41. Sharma, T., Kurz, W. A., Stinson, G., Pellatt, M. G. & Li, Q. A 100-year 7. Luyssaert, S. et al. Trade-offs in using European forests to meet climate conservation experiment: Impacts on forest carbon stocks and fluxes. For. objectives. Nature 562, 259 (2018). Ecol. Manag. 310, 242–255 (2013). 8. Mehrabi, Z., Ellis, E. C. & Ramankutty, N. The challenge of feeding the world 42. Chazdon, R. L. et al. Carbon sequestration potential of second-growth forest while conserving half the planet. Nat. Sustain. 1, 409 (2018). regeneration in the Latin American tropics. Sci. Adv. 2, e1501639 (2016). 9. Jackson, R. B. Trading water for carbon with biological carbon sequestration. 43. Fang, J. et al. Forest biomass carbon sinks in East Asia, with special reference Science 310, 1944–1947 (2005). to the relative contributions of forest expansion and forest growth. Glob. 10. Li, Y. et al. Divergent hydrological response to large-scale afforestation and Change Biol. 20, 2019–2030 (2014). vegetation greening in China. Sci. Adv. 4, eaar4182 (2018). 44. Earles, J., Yeh, S. & Skog, K. E. Timing of carbon emissions from global forest 11. Hogarth, N. J., Belcher, B., Campbell, B. & Stacey, N. The role of forest-related clearance. Nat. Clim. Change 2, 682–685 (2012). income in household economies and rural livelihoods in the border-region of 45. Dymond, C. C. Forest carbon in North America: annual storage and emissions southern China. World Dev. 43, 111–123 (2013). from British Columbia’s harvest, 1965–2065. Carbon Balance Manag. 7, 8 (2012). 12. Bellassen, V. & Luyssaert, S. Carbon sequestration: managing forests in 46. Stockmann, K. D. et al. Estimates of carbon stored in harvested wood products uncertain times. Nat. N. 506, 153 (2014). from the United States Forest Service Northern Region, 1906-2010. Carbon 13. Viña, A., McConnell, W. J., Yang, H., Xu, Z. & Liu, J. Effects of conservation Balance Manag. 7, 1–16 (2012). policy on China’s forest recovery. Sci. Adv. 2, e1500965 (2016). 47. Hansen, M. C., Potapov, P. & Tyukavina, A. Comment on “Tropical forests 14. Piao, S. et al. Detection and attribution of vegetation greening trend in China are a net carbon source based on aboveground measurements of gain and over the last 30 years. Glob. Change Biol. 21, 1601–1609 (2015). loss”. Science 363, eaar3629 (2019). 15. Hua, F. et al. Opportunities for biodiversity gains under the world’s largest 48. Fang, J. Y., Guo, Z. D., Piao, S. L. & Chen, A. P. Terrestrial vegetation carbon reforestation programme. Nat. Commun. 7, 12717 (2016). sinks in China, 1981–2000. Sci. China Ser. D. 50, 1341–1350 (2007). 16. Ahrends et al. China’s fight to halt tree cover loss. Proc. R. Soc. Ser. B 284, 49. Schaaf, C. B. et al. First operational BRDF, albedo nadir reflectance products 20162559 (2017). from MODIS. Remote Sens. Environ. 83, 135–148 (2002). 17. Fang, J. Changes in forest biomass carbon storage in China between 1949 and 50. Wang, Y. et al. Major forest increase on the Loess Plateau, China (2001–2016). 1998. Science 292, 2320–2322 (2001). Land Degrad. Dev. https://doi.org/10.1002/ldr.3174 (2018). 18. Tong, X. et al. Increased vegetation growth and carbon stock in China karst 51. Kennedy, R. E., Yang, Z. & Cohen, W. B. Detecting trends in forest via ecological engineering. Nature. Nat. Sustain. 1, 44–50 (2018). disturbance and recovery using yearly Landsat time series: 1. LandTrendr— 19. Delang, C. O. & Yuan, Z. China’s Grain for Green Program (Springer temporal segmentation algorithms. Remote Sens. Environ. 114, 2897–2910 International Publishing, 2015). https://doi.org/10.1007/978-3-319-11505-4. (2010). 20. Brandt, M. et al. Satellite-observed major greening and biomass increase in 52. Sulla-Menashe, D. et al. Detecting forest disturbance in the Pacific Northwest South China Karst during recent decade. Earth’s Fut. 6, 1017–1028 (2018). from MODIS time series using temporal segmentation. Remote Sens. Environ. 21. Zeng, Z. et al. Highland cropland expansion and forest loss in Southeast Asia 151, 114–123 (2014). in the twenty-first century. Nat. Geosci. 11, 556 (2018). 53. Li, X. et al. Effects of successive rotation regimes on carbon stocks in 22. Davis, K. F., Yu, K., Rulli, M. C., Pichdara, L. & D’Odorico, P. Accelerated eucalyptus plantations in subtropical China measured over a full rotation. deforestation driven by large-scale land acquisitions in Cambodia. Nat. Geosci. PLoS ONE 10, e0132858 (2015). 8, 772–775 (2015). 54. Saatchi, S. S. et al. Benchmark map of forest carbon stocks in tropical regions 23. Baccini, A. et al. Tropical forests are a net carbon source based on across three continents. Proc. Natl Acad. Sci. USA 108, 9899–9904 (2011). aboveground measurements of gain and loss. Science https://doi.org/10.1126/ 55. Carreiras, J. M. et al. Coverage of high biomass forests by the ESA BIOMASS science.aam5962 (2017). mission under defense restrictions. Remote Sens. Environ. 196, 154–162 24. Brandt, M. et al. Satellite passive microwaves reveal recent climate-induced (2017). carbon losses in African drylands. Nat. Ecol. Evolution 2, 827–835 (2018). 56. Xu, L. et al. Spatial distribution of carbon stored in forests of the Democratic 25. Fan, L. et al. Satellite observed pantropical carbon dynamics. Nat. Plants 5, Republic of Congo. Sci. Rep. 7, 15030 (2017). 944–951. 57. McRoberts, R. E. et al. Using a finer resolution biomass map to assess the 26. Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover accuracy of a regional, map-based estimate of forest biomass. Surv. Geophys change. Science 342, 850–853 (2013). 40, 1001–1015 (2019). 27. Song, X.-P. et al. Global land change from 1982 to 2016. Nature 560, 639–643 58. Mitchard, E. T. et al. Uncertainty in the spatial distribution of tropical forest (2018). biomass: a comparison of pan-tropical maps. Carbon Balance Manag. 8, 10 28. Shan, Y. et al. China CO2 emission accounts 1997–2015. Sci. Data 5, 170201 (2013). (2018). 59. Ma, H. et al. Satellite surface soil moisture from SMAP, SMOS, AMSR2 and 29. Funk, C. et al. The climate hazards infrared precipitation with stations—a new ESA CCI: a comprehensive assessment using global ground-based environmental record for monitoring extremes. Sci. Data 2, 150066 (2015). observations. Remote Sens. Environ. 231, 111215 (2019). 30. Zhang, H. et al. Biomass and carbon storage of Eucalyptus and Acacia 60. Yin, Z. et al. Evaluation of ORCHIDEE-MICT-simulated soil moisture over plantations in the Pearl River Delta, South China. For. Ecol. Manag. 277, China and impacts of different atmospheric forcing data. Hydrol. Earth Syst. 90–97 (2012). Sci. 22, 5463–5484 (2018). 31. Yao, Y. et al. The effect of afforestation on soil moisture content in 61. Wigneron, J.-P. et al. Modelling the passive microwave signature from land northeastern China. PLoS One 11, e0160776 (2016). surfaces: a review of recent results and application to the L-band SMOS & 32. Tan, Z.-H. et al. Rubber plantations act as water pumps in tropical China. SMAP soil moisture retrieval algorithms. Remote Sens. Environ. 192, 238–262 Geophys. Res. Lett. 38, https://doi.org/10.1029/2011GL050006 (2011). (2017). 33. Liu, Y. et al. Water use efficiency of China’s terrestrial ecosystems and 62. Fernandez-Moran, R. et al. SMOS-IC: an alternative SMOS soil moisture and responses to drought. Sci. Rep. 5, 13799 (2015). vegetation optical depth product. Remote Sens. 9, 457 (2017). NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13798-8 63. Al-Yaari, A. et al. Assessment and inter-comparison of recently developed/ Additional information reprocessed microwave satellite soil moisture products using ISMN ground- Supplementary information is available for this paper at https://doi.org/10.1038/s41467- based measurements. Remote Sens. Environ. 224, 289–303 (2019). 019-13798-8. Correspondence and requests for materials should be addressed to Y.Y. or K.W. Acknowledgements Peer review information Nature Communications thanks the anonymous reviewer/s for X.T. and Y.Y. were funded by the National Key Research and Development Program of their contribution to the peer review of this work. Peer reviewer reports are available. China (no. 2016YFC0502400, 2018YFD1100103), the National Natural Science Foun- dation of China (41930652), and the Strategic Priority Research Program of Chinese Reprints and permission information is available at http://www.nature.com/reprints Academy of Sciences (XDA23060100). X.T. was funded by a Marie Curie fellowship, grant number 795970. M.B. received funding from the AXA research grant. R.F. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in acknowledges the funding from the Danish Council for Independent Research (DFF) published maps and institutional affiliations. Grant ID: DFF-6111-00258. X.X. was funded by the NASA Land Cover and Land Use Change Program (NNX14AD78G). P.C. and J.P. were supported by the European Research Council Synergy grant ERC-2013-726 SyG-610028 IMBALANCE-P. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give Author contributions appropriate credit to the original author(s) and the source, provide a link to the Creative X.T. and M.B. contributed equally to this manuscript. X.T., M.B., R.F., M.R.J., S.H., Y.Y., Commons license, and indicate if changes were made. The images or other third party P.C., J.P. and K.W. designed the study. X.T., M.B., Y.Y and Y.W. conducted the analyses material in this article are included in the article’s Creative Commons license, unless with support by P.C., J.P., X.-P.S. and L.F. The data were provided by J.P., L.F., S.S., indicated otherwise in a credit line to the material. If material is not included in the X.-P.S., B.Z. and C.Z. M.B. and X.T. drafted the manuscript with contributions by P.C., article’s Creative Commons license and your intended use is not permitted by statutory J.P., R.F., K.R., M.R.J., J.-P.W., X.X., X.-P.S., S.H. and all other authors. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Competing interests The authors declare no competing interests. © The Author(s) 2019 10 NATURE COMMUNICATIONS | (2020)11:129 | https://doi.org/10.1038/s41467-019-13798-8 | www.nature.com/naturecommunications