Carbon storage potential in degraded forests of Kalimantan, Indonesia

The lidar-derived MCH provides strong power-law relationships with AGB for dryland (R² = 0.81) and wetland forests (R² = 0.79). The power-law models are significantly different, showing a higher biomass in wetland forests for a given lidar MCH because of a relatively high tree number density and basal area compared with drylands forests (figure 3). We assumed the two models developed in this study capture the largest structural and allometric differences in the region between wetlands and drylands. Without having adequate plots in all forests types present in our study domain, we were not able to verify if there are other vegetation types such as primary swamp pole forests and mangroves with distinct lidar models. Similar assumptions have also been made about the AGB estimates at the ground plots. We assumed that the pan-tropic allometric model provides reasonable estimates of AGB for all forest types in the region (Chave et al 2014).

Zoom In Zoom Out Reset image size

The wall-to-wall AGB map of Kalimantan derived from RF algorithm captures the imprints of land use activities and fragmentations on the spatial variations of AGB (table 1 and figure 4). Dense canopy forests with tall trees and high AGB are located inland away from coastal regions and with biomass density increasing at higher elevation where more intact forests are located (figures 4(a) and (b)). The average AGB density of the entire study region from the RF map is 178.4 ± 7.5 Mgha⁻¹, ranging from 34.7 ± 10.1 Mgha⁻¹ for young tree plantations to 312.7 ± 20.6 Mgha⁻¹ for intact lowland forests. The four target land cover classes (denoted with * in table 1) account for 90% of total AGB with an average of 274.1 ± 14.0 Mgha⁻¹.

Zoom In Zoom Out Reset image size

The estimates of average AGB and uncertainty from the original lidar samples and the RF map are provided in table 2. The AGB of intact lowland (312.7 ± 20.6 Mgha⁻¹) and montane forests (311.0 ± 51.2 Mgha⁻¹) are similar and together they occupy an area nearly as large as the secondary and degraded forests, which stores in average 51 Mgha⁻¹ less than intact forest. The AGB density of peat swamp forests is even lower (211.8 ± 12.7 Mgha⁻¹) but covers a large spectrum of AGB density (the 5th and 95th percentile equal 30.9 and 403.3 Mgha⁻¹). Peatlands are heterogeneous and occupy a highly fragmented landscape ranging from highly dense primary peat swamp pole forest to the stunted peat swamp padang forest with significantly lower AGB. Our result suggests that some of the high biomass forests in Kalimantan are located in the peat swamps surrounding Lake Sentarum and in the Sebangau National Park (dark red areas in the northwest and south of Kalimantan visible in figure 4(a)) that correspond to primary peat swamp tall pole forests (section 3.5). At the same time, the peat swamp forests also show an overall lower AGB due to large-scale degradation and forest fires, particularly in areas outside and surrounding the Sebangau National Park that include large-scale rice fields (Rose et al 2011).

Among the remaining vegetation types in Kalimantan, the scrublands cover an area of about 6.88 Mha with average AGB of 54.8 ± 5.5 Mgha⁻¹, and swamp scrublands 3.96 Mha with an average AGB of 40.7 ± 4.9 Mgha⁻¹ (figure 4(b)). Tree plantations cover 7.84 Mha and have an average 34.7 ± 10.1 Mgha⁻¹ indicating that the IMF map only considers the most recent plantations in this class and may have confused older and higher biomass plantations with some remaining natural trees as degraded forests. Note that we assign a higher confidence to the results regarding the target land cover classes because the lidar-AGB models have been developed using field inventory collected over those forest types.

The FDI derived degradation map shows a clear gradient of decreasing forest degradation towards inland and mountainous areas (figure 4(c)). The impact of degradation on AGB storage can be estimated using the average density for different FDI classes that range from intact old growth (357.2 ± 12.3 Mgha⁻¹) to severe degraded forests (134.2 ± 6.1 Mgha⁻¹, table 1). The AGB of intact forests increased by about 50 Mg/ha on average after improving the classification of intact and forests using FDI. The AGB distribution shown in violin plots (figure 4(d)) demonstrates a significant overlap among the FDI classes indicating the presence of a large variability of AGB within each class. For instance, intact forest can store low values of AGB on higher elevations and rough topography areas. The intact old growth forests cover the largest area (7.19 Mha), followed by moderate (4.96 Mha) and light (4.59 Mha) degraded forests and with high (2.05 Mha) and severe (3.52 Mha) degraded forests being less abundant.

According to the GIS data layers of logging concessions in Kalimantan (figure SI4), only 20% of the surface covered by secondary and degraded forest class are outside the official logging concession areas. Thus, logging companies are currently managing 22% out of the 27% of Kalimantan's forest biomass stored in secondary and degraded forests. Not surprisingly, the AGB density of forests within the logging concessions (258.0 ± 16.5 Mgha⁻¹) is nearly identical to AGB density in the secondary and degraded forest class (tables 1 and SI3). Other agroforestry classes such as palm oil (9.82 Mha, 59.7 ± 11.3 Mgha⁻¹) and wood fiber (5.5 Mha, 95.4 ± 13.2 Mgha⁻¹) cover large areas and store moderate amounts of AGB. We found the AGB values of forest types in different Kalimantan administrative areas had no significant differences, suggesting that the land use activities and practices are relatively similar across jurisdictions (see SI5.3).

The difference between lidar- and RF-estimated averaged AGB is relatively low for the target land cover classes, with 1.9% for intact lowland forest, −0.9% for intact montane forest, 3.6% for secondary and degraded forests, and 1. 7% for peat swamp forests (the percentage is relative to average AGB and negative values means that RF underestimates lidar estimates). While that difference is also relatively low regarding the FDI classes (>8%), it reaches 49.5%, 24.3%, and −60.5% for swamp scrublands, scrublands and plantations, respectively (table 1). Nevertheless, these land cover classes account only for the 11% of the total area of Kalimantan and their effects on the average difference between lidar and RF derived AGB for the entire region remain bounded at about 6%.

We analyzed the AGB density for the main wetlands present in Kalimantan using AGB estimated directly from the lidar samples (SI4 and table 2). Primary peat swamp tall pole forests (also called tall interior forests, Wösten et al 2008) show a mean AGB density (457.0 ± 13.8 Mgha⁻¹) significantly higher than secondary tall pole forests (207.4 ± 0.7 Mgha⁻¹). These results must be taken with caution because the lidar-AGB model for wetlands does not include plots with AGB exceeding 350 Mgha⁻¹, but the model has been used to predict significantly larger biomass forests. Nevertheless, the AGB map indicates that tall interior forests are within the densest areas of Kalimantan, which agrees with field measurements reported in Verwer and van der Meer 2010 (645.34 Mgha⁻¹). Tall interior forests are highly heterogeneous in terms of forest structure ranging from swamp areas with large tree gaps and low AGB to areas with densely packed and tall trees (5th and 95th percentile in table 2).

Of the remaining wetland vegetation types, peat swamp pandang forest has an AGB density of 95.5 ± 1.2 Mgha⁻¹, because of their short stature and less homogenous cover. Other wetland classes all have low AGB density. The burned peatlands with AGB of 18.5 ± 1.3 Mgha⁻¹ show very little recovery since being classified in the IMF map, and riverine forests, composed largely of woody vegetation, have higher AGB density (54.0 ± 1.8 Mgha⁻¹) than Alang-Alang, which is composed of herbaceous plants.

We report the uncertainty in estimating average AGB density for all land cover types and also at jurisdictional scales (table 2, SI3 and SI4) using error propagation throughout the AGB estimation and considering the spatial autocorrelation of the errors (SI6.4). The uncertainty associated with the lidar-AGB models for drylands and wetlands differ significantly (RMSE of 62.2 Mgha⁻¹ and 19.28 Mgha⁻¹, figure 3) due to more heterogeneous nature of dryland structure and degradation compared to wetland forests. Both models showed small bias through cross-validation with plot data with 2.03 Mgha⁻¹ for drylands −1.19 Mgha⁻¹ for wetlands. Similarly, the RF predictions also had small systematic error from cross-validation (bias = 0.49 ± 10.9 Mgha⁻¹) but a relatively large random error (RMSE = 101.9 ± 11.3 Mgha⁻¹⁾ (table SI5).

The uncertainty of the map also remains bounded (<7%) for every land cover class of interest ensuring that the wall-to-wall map provides precise average and total AGB estimates at the scale of land cover classes or jurisdictions. The only exception is the intact montane forests with average AGB of 311 ± 51.2 Mgha⁻¹ and uncertainty that reaches to about 16% (table SI6). For degradation classes derived from FDI thresholds, the uncertainty remains lower (<4%) than similar classes from the IMF map, suggesting that the average AGB of these classes of degradation can be estimated with higher confidence because they are structurally more homogeneous than the IMF map (table 1).

The RF predicted AGB at 1 ha pixels shows a relatively good agreement with an independent set of ground-estimated AGB (R² = 0.73, figure 5). The field plots used for comparison have been covered by airborne lidar but were not used for the development of lidar-AGB models (SI6.5). However, these results should be taken with caution because the plots are located within the lidar flights and the RF predictions are prone to overfitting the lidar data used in training the map. Field plots located in areas with forest structure significantly different from those covered by the lidar may have larger uncertainty.

Zoom In Zoom Out Reset image size

The derived AGB distribution in Kalimantan forests is based on a probabilistic airborne lidar sampling that ensures unbiased estimates of mean and total forest AGB subject to the choice of the lidar-AGB model (Ståhl et al 2016), similar to ground inventory sample measurements that can produce unbiased estimates given the right biomass allometry. However, due to requirements for planning airborne flights in the region, the lidar sampling design provided only 29 lidar scenes for a total of 29 000 ha. These samples, although randomly located across the region, are not widespread and do not adequately capture the variations of all types of forests. Selecting smaller sampling areas (∼500 ha) and a larger number of samples could have improved the RF predictions and the uncertainty of estimates for all forest types present in the region. Given additional resources, an improved sampling strategy based on land cover stratification can be used for future sampling design.

We collected a larger number of airborne lidar data over existing plot networks to allow development of reliable and unbiased lidar-AGB models. The lack of any systematic inventory plots in the region and the diversity in existing plot size and quality are considered sources of errors that cannot be readily quantified or mitigated. However, our methodology allows for improvement of lidar-AGB models by acquiring additional ground inventory data within the lidar coverage.

Unlike the previous AGB maps (Saatchi et al 2011a, Baccini et al 2012), we developed forest-type specific lidar-AGB models (for drylands and wetlands) that adjust to differences in the 3D forest structure and tree density in order to reduce the uncertainty associated with lidar estimates of AGB. The average AGB values for the target land cover classes derived from the RF prediction map are in agreement with the lidar-derived estimates, suggesting the map as a reliable tool to estimate the average AGB density for different land cover classes or at landscape scales (100–10 000 ha). One of the limitations of the RF prediction algorithm is the problem of overfitting the distribution that may lead to underestimation of the range of biomass variability, a so-called dilution bias. Although we included a bias correction approach that improves the uncertainty of the RF prediction bias (Xu et al 2016), the map still shows differences at the 5th and 95th percentiles with the lidar-AGB distribution (table 1).

Furthermore, the map may slightly underestimate high biomass density due to the lack of sensitivity of satellite imagery used in RF models. Upcoming satellite missions (NASA-ISRO NISAR mission, NASA GEDI mission and ESA BIOMASS mission) are expected to improve the sensitivity of measurements to forest structure and AGB across the entire biomass range, particularly in the dense tropical forests and across the complex terrains (LeToan et al 2011, Yu and Saatchi 2016). However, the underestimation is in general much smaller than the predicted sensitivity of the satellite data used as input layers to develop the machine learning algorithms. Both the existing radar and optical sensors are known to be saturated over dense tropical forests when used directly in parameteric models to estimate forest biomass (Saatchi et al 2011b). The non-parametric machine-learning approach has a significantly more efficient way of using the spatial variations of the data layers to extrapolate or estimate forest biomass from the training data. The approach includes the use of spatial distribution of training data from lidar-derived biomass to increase the probability of predicting biomass values across landscapes. For example, most remaining high biomass forests are located across upland and higher elevation landscapes and the RF model uses the SRTM (elevation and surface ruggedness) to predict the high biomass forests in Kalimantan (S13). In addition, other layers such as short wavelength and near infrared bands of multi-spectral Landsat imagery potential help to separate high biomass forests due to shadows from dense layered canopy structures and older aged leaves from open and younger aged leaves of degraded and secondary forests (Steininger 1996).

The average AGB estimates for drylands intact forests reported in the literature in recent years are relatively higher (436 Mgha⁻¹ in Slik et al 2010, 426 Mgha⁻¹ in Qie et al 2017) than the estimates from the RF map (357.2 Mgha⁻¹) or lidar sampling (386.7 Mgha⁻¹). Similarly, the RF AGB prediction of the peat swamp forests of 218 ± 12.7 Mgha⁻¹ is relatively lower than those published by Verwer and van der Meer (2010) of 278.85 Mgha⁻¹ and, Murdiyarso et al (2010) of 269.55 Mgha⁻¹. We expect our results from either the RF map or the lidar sampling to be more realistic because of the probabilistic sampling approach used in our study. Most research plots in these studies do not follow a systematic sampling and are often located in areas of high biomass forests that lead to an overestimation of the AGB in Kalimantan, the so-called majestic forest effect (Sheil 1995). The high-resolution lidar observations used in our study capture the natural variability of structure and biomass in intact forests caused by topography, soil variations, wind and other natural disturbance. By providing the 5th and 95th biomass values for different forest types, the range of biomass in intact forests can be readily compared with published results. Furthermore, by adjusting the values of mean canopy height and percent of large trees, the spatial products from this study can be used to delineate areas of high biomass density characteristic of the majestic forests.

Following the CDM of the Kyoto Protocol guidelines (UNFCCC 2002), we define forest as an area with vegetation higher than 5 m and with cover that exceeds 30% of the 1 ha grid cells. We developed the FDI to map different degrees of forest degradation and improve the IMF based on Landsat imagery. The average biomass loss from degradation based on our study is much larger (up to 50%) than what is reported in the literature for selectively logged forests in Kalimantan (3%–15%) (Sasaki and Putz 2009, Kronseder et al 2012, Pearson et al 2014). Using the FDI index, we were able to delineate the lightly and moderately degraded forests within the Landsat land cover map that closely capture the selectively logged forests. However, other types of degradations due to different logging intensity (Pearson et al 2014), fragmentations, loss of biomass from edge effects, and fire, reported from our study, must be considered as degradation in accounting for carbon loss and emissions in the region.

Based on the RF map, the intact and light degraded forests cover (11.78 Mha) approximately the same size area of moderate to severe degraded forests combined (10.53 Mha), suggesting a large area of dryland forests with significant potential for carbon sequestration. By focusing on dryland forests, the total carbon stored in about 22.3 Mha of forests is approximately 3.1 PgC (using average AGB of 284.5 Mgha⁻¹ from table 1, and the carbon fraction of 0.48). The emission factors associated with degradation can be calculated by the difference between the average carbon storage of intact and degraded forests that multiplied by the area occupied by the degraded forests provide the carbon sequestration or storage potential of degraded forests. We find the storage potential of degraded forests to be about 0.8–1.1 PgC depending on emission factors calculated from lidar or RF map. Kalimantan degraded forests have significantly larger storage potential per ha (70.2 MgCha⁻¹) than the entire Latin America second-growth forests of age 1 to 60 years (35.33 MgCha⁻¹) and over a shorter period of time (Chazdon et al 2016).