Doubling of annual forest carbon loss over the tropics during the early twenty-first century
Global map of forest cover change and its validationWe use a high-resolution map of global forest cover change15 (annual intervals at 30 m spatial resolution; version 1.7) to quantify forest loss across the tropics. The GFC dataset maps where and when forests were converted (naturally and anthropogenically) from 2001 to 2019. Trees are defined as all vegetation taller than 5 m in height, and forests are defined with a tree canopy threshold of at least 30%. The definition of forests includes plantations and tree crops such as oil palm. Forest loss is the mortality or removal of all tree cover within a pixel.Our study uses the v.1.7 product that spans the period 2001–2019 for the analysis. Different methods are used for detecting forest cover loss in two periods (2001–2010 and 2011–2019). This change in detection method as well as in satellite data (Landsat 7 and Landsat 8) might result in inconsistencies of data during the two periods. Therefore, we perform an independent assessment of the v.1.7 product throughout the study period (2001–2019) using stratified random-sample reference data. We randomly sample 18,000 pixels, a much larger sample population than the assessment of the original product (v.1.0, 628 pixels across the tropics; ref. 15), and visually interpret forest loss using Landsat imagery. Specifically, we randomly select 50 path/row locations (World Reference System II) of Landsat imagery in each tropical continent (150 path/row locations in total in the three tropical continents; Supplementary Fig. 6). For each path/row location, we randomly select 20 loss pixels and 10 non-loss pixels in each period of 2001–2005, 2006–2010, 2011–2014 and 2015–2019, with total sampling pixels of ~18,000 ((20 loss + 10 non-loss) × 50 path/row locations × 4 periods × 3 continents). Our study compares the increase in forest carbon loss from the start (~2001) to the end (~2019) of the study period. Thus, we divided the whole 19 yr study period into four subperiods (2001–2005, 2006–2010, 2011–2014 and 2015–2019), with the first five years considered as the start period and the last five years considered as the end period for the comparison. Some path/row locations do not have 20 loss pixels in a specific period; for example, there is no loss detected in some locations of the Sahara. Therefore, we sample 11,198 loss pixels and 6,000 non-loss pixels (Supplementary Data). Finally, we download time-series Landsat imagery covering 1999–2020 to visually interpret these pixels as reference data.Following the suggestion of Global Forest Watch19 and per best practice guidance of ref. 18, we use a stratified random-sample approach for area estimation, which is independent of the method and satellite changes in the GFC data. The sampling reference data (Supplementary Data) are used to estimate loss area:$${it{p}}_{hi} = w_hfrac{{mathop {sum}limits_{j = 1}^{150} {n_{hij} times A_{hi}} }}{{mathop {sum}limits_{j = 1}^{150} {n_{hj} times A_i} }}$$
(1)
where phi is the stratified random-sample estimated area for GFC map class h that is classified as reference class i; wh is the proportion of the total area GFC map class h; nhij is the number of pixels in GFC map class h that is classified as reference class i in jth Landsat path/row location; nhj is the total number of pixels in GFC map class h in jth Landsat path/row location; Ahj is the total area in GFC map class h in jth Landsat path/row location; and Aj is the total land area of jth Landsat path/row location.We then calculate the error matrix, which includes overall accuracy (OA), user’s accuracy (UA) and producer’s accuracy (PA), as follows:$${mathrm{OA}} = mathop {sum}limits_{h = 1}^H {p_{hh}}$$
(2)
$${mathrm{UA}}_h = frac{{p_{hh}}}{{p_h}}$$
(3)
$${mathrm{PA}}_j = frac{{p_{jj}}}{{p_j}}$$
(4)
where UAh is the UA for stratum h; ph and pj are the total area in stratum h and j, respectively; and PAj is the PA in stratum j.Using the preceding equations, we estimate OA, UA and PA in the four periods (Supplementary Table 2). In general, OAs are >99%, UAs are >88% and PAs are >72% in each period. The stratified random-sample approach for area estimation is considered the most robust method to investigate loss trends in GFC product and can avoid inconsistencies of the dataset due to changes in detection model and satellite sensors19. Our results show that the forest loss from the stratified random-sample approach is similar to GFC mapped loss, both of which show a consistent increase during the four periods of 2001–2019 (Fig. 1a), confirming the increasing forest carbon loss across the tropics during the early twenty-first century.Forest carbon stocksWe estimate forest (aboveground and belowground) biomass carbon losses by co-locating GFC loss data with corresponding biomass data. Forest biomass maps are not universally reliable, owing to uncertainties and some degree of bias. We use four biomass maps to quantify forest carbon stocks, which helps reduce the uncertainties and bias44. The four maps were developed by refs. 9,48,49,20 and are hereafter referred to as ‘Baccini’, ‘Saatchi’, ‘Avitabile’ and ‘Zarin’ maps, respectively. The Baccini map, derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data, presents aboveground live woody biomass (AGB) across the tropics at 500 m spatial resolution. The Saatchi map, also derived from MODIS data, presents total forest carbon stocks at 1 km spatial resolution across the tropics. The Avitabile map, integrated from the Baccini and Saatchi maps, shows AGB at 1 km resolution across the tropics. The Zarin map, derived from Landsat data, presents AGB across the globe at 30 m resolution.Belowground root biomass (BGB) data are sparse because measurements of BGB are time consuming, laborious and technically challenging50. Thus, we calculate BGB (in Mg ha−1 biomass) from AGB maps (Baccini, Avitabile and Zarin) using an empirical model at the pixel level51:$${mathrm{BGB}} = 0.489 times {mathrm{AGB}}^{0.89}$$
(5)
Total forest biomass is calculated as the sum of AGB and BGB. Finally, total forest carbon stocks (MgC ha−1) in live woody forest are estimated as 50% of total biomass20,50. The Saatchi map provides total forest carbon stocks rather than AGB. The total forest carbon stocks are calculated from AGB using the same method mentioned in the preceding50. Thus, we estimate AGB and BGB from total forest carbon stocks in the Saatchi map using the preceding method to separate aboveground and belowground parts of forest carbon stocks.There are inconsistencies in the MODIS-derived biomass maps (Baccini, Saatchi and Avitabile) and Landsat-derived GFC data41, which may underestimate forest carbon loss by the three coarse forest biomass maps (Supplementary Fig. 7). The Zarin map is derived from Landsat data and considers tree cover using GFC data, and tree loss can be co-located with the corresponding biomass20, indicating the consistencies of the Zarin map and GFC. Therefore, to correct the three biomass maps with coarse resolution and reduce the inconsistencies, we resample the Zarin map from 30 m to 500 m (the resolution of the Baccini map) and 1 km (the resolutions of the Avitabile and Baccini maps) and calculate forest carbon loss using the two resampled biomass maps. We then estimate the ratios of forest carbon loss derived from the 30 m biomass map to the forest carbon loss derived from resampled biomass maps in each GFC tile (10° × 10°). The ratios are then used as a scale factor to correct the three biomass maps (Baccini, Saatchi and Avitabile). For the resampling, forest carbon density at 500 m or 1 km is averaged from all 30 m pixels in the corresponding locations.Soil organic carbon stocksDeforestation not only causes forest biomass carbon loss, but also results in loss of SOC10. We calculate SOC loss at 0–30 cm soil depth as:$${mathrm{SOC}}_{{{{mathrm{loss}}}}} = {mathrm{OCS}} times theta$$
(6)
where SOCloss is SOC loss resulting from forest loss measured in MgC ha−1, OCS is SOC stocks at 0–30 cm depth measured in MgC ha−1 and θ is the SOC loss rate.We obtain SOC stocks at 0–30 cm depth from SoilGrids (version 2.0), created by the International Soil Reference and Information Centre52. SOC stocks are calculated using a calibrated quantile random forest model at a spatial resolution of 250 m. We further resample the data from 250 m to 30 m using the nearest-neighbour method to match the scale of forest cover loss data.SOC loss rate is affected predominantly by land-use types following forest loss and tree species53. The loss rate data are compiled from a previous meta-analysis, which summarizes the rate of SOC loss resulting from forest loss over the tropics50. SOC losses resulting from primary and secondary forest loss differ (Supplementary Table 3). We use a map of primary humid tropical forests in 2001 developed by ref. 54 to classify primary and secondary forests. The land covers following forest loss are determined according to the driver of forest loss (see Drivers of forest carbon loss). The spatial resolution of the SOC data is coarse compared with GFC data, which may result in inconsistencies in the calculation. However, as SOC loss resulting from forest loss accounts for a small proportion of total forest carbon loss (8%), we ignore the potential inconsistencies.Drivers of forest carbon lossWe determine drivers of tree-cover loss using the dataset generated by ref. 27. This dataset shows the dominant driver of tree-cover loss at each 10 km grid cell for 2001–2019. There are five categories of drivers of tree-cover loss: commodity-driven deforestation, defined as permanent and/or long-term clearing of trees to other land uses (for example, commodity croplands), shifting agriculture, forestry, wildfire and urbanization. Commodity-driven deforestation, shifting agriculture and forestry dominate tropical forest loss27. Thus, wildfire and urbanization are combined and categorized as ‘others’. In tropical Africa, spatial patterns of commodity-driven deforestation are almost similar to that of shifting agriculture, as pointed out by ref. 27, resulting in large uncertainties in separating commodity agriculture from shifting agriculture. Since commodity-driven deforestation is usually for large-scale agricultural plantations, we treat commodity-driven deforestation as loss for large-scale agriculture. Because shifting agriculture is usually smallholder and/or patchy farming systems, we term shifting agriculture as small-scale agriculture, which may include commodity agriculture with similar spatial patterns to shifting agriculture in some regions such as tropical Africa. The driver data are created using decision-tree models trained by ~5,000 high-resolution Google Earth imagery cells, showing overall accuracy of 89 ± 3% from a separate validation of more than 1,500 randomly selected cells. We resample the data from 10 km resolution to 30 m using the nearest-neighbour method to match the scale of forest cover loss data.Interpretation of post-forest-loss land covers in 2020We collect cloudless and very-high-resolution satellite imagery in 2020 from Planet to determine the fate of the agriculture-driven forest loss during 2001–2019. Planet provides two products, RapidEye (at a spatial resolution of 5 m) and Doves (at a spatial resolution of 3 m; 4-band PlanetScope Scene). We randomly sample 500 pixels that show forest loss owing to agriculture expansion during 2001–2005 and 500 pixels that show forest loss owing to agriculture expansion during 2015–2019, then check cloudless satellite imagery in 2020 to visually interpret the land cover of each sampled pixel in 2020. We classify three types of land cover: agricultural land, forest/shrubland and others (Supplementary Fig. 4). Rubber and oil-palm plantations are classified as agricultural lands.Uncertainty and methods for analysisWe use committed emissions of forest carbon, even though some of this carbon will be lost only in later years or transited to other carbon pools or stored as wood products12. Forest carbon loss is defined as gross carbon loss due to forest removal (as indicated by GFC product), including (aboveground and belowground) forest biomass carbon and SOC losses. We calculate only the gross loss of forest carbon stocks while gain of carbon via reforestation and afforestation is not considered.We first estimate reference sample-based forest loss area using sample data:$${mathrm{AS}}_j = {mathrm{AM}}_jmathop {sum}limits_{h = 1}^H {p_{hj}}$$
(7)
where ASj and AMj are reference sample-based forest loss areas and mapped forest loss areas in stratum j, respectively.To estimate forest carbon loss (aboveground, belowground and soil carbon) using sample-based forest loss area, we apply a ‘stratify and multiply’ approach21,22 by assigning mean forest carbon density for each stratum. Aboveground and belowground forest carbon losses are estimated using four biomass maps (Baccini, Saatchi, Avitabile and Zarin), and we report ensemble mean ± s.d. from the four maps as our best estimate.Mountain forest carbon loss at different elevations is calculated by overlaying mountain polygon from the Global Mountain Biodiversity Assessment inventory55 (version 1.2) and the 30 m ASTER Global Digital Elevation Model56 (version 3).Although our validation shows that the GFC (v.1.7) product could accurately map forest loss, we cannot reduce the omission and commission errors. In addition, the accuracy of the disturbance year is 75.2%, with 96.7% of the disturbance occurring within one year before or after the estimated disturbance year in GFC product15. Therefore, we calculate 3 yr moving averages of annual forest and related carbon losses for time-series analysis, following the suggestion of the Global Forest Watch19 and refs. 28,38. We use a non-parametric Theil–Sen estimator regression method53 to detect trends in time-series results and test the significance of the trend by Mann–Kendall test57.Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
