Forest aboveground biomass stock and resilience in a tropical landscape of Thailand

Half of Asian tropical forests were disturbed in the last century resulting in the dominance of secondary forests in Southeast Asia. However, the rate at which biomass accumulates during the recovery process in these forests is poorly understood. We studied a forest landscape located in Khao Yai National Park (Thailand) that experienced strong disturbances in the last century due to clearance by swidden farmers. Combining recent field and airborne laser scanning (ALS) data, we first built a high-resolution aboveground biomass (AGB) map of over 60 km² of forest landscape. We then used the random forest algorithm and Landsat time series (LTS) data to classify landscape patches as non-forested versus forested on an almost annual basis from 1972 to 2017. The resulting chronosequence was then used in combination with the AGB map to estimate forest carbon recovery rates in secondary forest patches during the first 42 years of succession. The ALS-AGB model predicted AGB with an error of 14 % at 0.5 ha resolution (RMSE=45 Mg ha⁻¹) using the mean top-of-canopy height as a single predictor. The mean AGB over the landscape was 291 Mg ha⁻¹, showing a high level of carbon storage despite past disturbance history. We found that AGB recovery varies non-linearly in the first 42 years of the succession, with an increasing rate of accumulation through time. We predicted a mean AGB recovery rate of 6.9 $\mathrm{Mg}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, with a mean AGB gain of 143 and 273 Mg ha⁻¹ after 20 and 40 years, respectively. This rate estimate is about 50 % larger than the rate prescribed for young secondary Asian tropical rainforests in the 2019 refinement of the 2006 IPCC guidelines for national greenhouse gas inventories. Our study hence suggests that the new IPCC rates, which were based on limited data from Asian tropical rainforests, strongly underestimate the carbon potential of forest regrowth in tropical Asia. Our recovery estimates are also within the range of those reported for the well-studied Latin American secondary forests under similar climatic conditions. This study illustrates the potential of ALS data not only for scaling up field AGB measurements but also for predicting AGB recovery dynamics when combined with long-term satellite data. It also illustrates that tropical forest landscapes that were disturbed in the past are of utmost importance for the regional carbon budget and thus for implementing international programs such as REDD+.

1 Introduction

Tropical forest disturbances and subsequent biomass recovery through time significantly affect the global carbon cycle (Harris et al., 2012). Although secondary forests in the tropics could constitute a major global carbon sink, the magnitude of such a sink remains poorly known (Chazdon, 2014; Lugo and Brown, 1992). A previous study estimated that 40 years of carbon storage in regenerating tropical forests from Latin America offset the past 19 years of carbon emissions from fossil fuels and industrial production in this region (Chazdon et al., 2016). Thus, there has been much interest in quantifying the ability of tropical secondary forests to sequester carbon in order to reduce uncertainties in the global carbon balance (e.g., Chai, 1997; Lohbeck et al., 2015; Stas et al., 2017).

Previous studies have used long-term forest plot surveys along chronosequences to quantify forest carbon dynamics in secondary tropical forests (Chazdon et al., 2007; N'Guessan et al., 2019; Norden et al., 2011, 2015; Poorter et al., 2016a; Rozendaal and Chazdon, 2015). Although long-term forest plots are essential for understanding the dynamics of tropical forests (Losos and Leigh, 2004), they are scarce, inherently labor-intensive, expensive and time-consuming to maintain, and not evenly distributed in the tropics. In addition, most studies of carbon dynamics along tropical forest successions are concentrated in Latin America (Chave et al., 2020; Letcher and Chazdon, 2009; Norden et al., 2015; Poorter et al., 2016a; Rozendaal et al., 2017; Rozendaal and Chazdon, 2015, but see N'Guessan et al., 2019, for Africa). They show high among-site variation in forest carbon recovery rates, suggesting a high context dependence (Chazdon et al., 2007; Norden et al., 2011, 2015), partly depending on climate conditions (Poorter et al., 2016a). A few pantropical studies have shown that the carbon potential of Latin American forests is smaller than that of Southeast Asian and African forests (Feldpausch et al., 2012; Sullivan et al., 2017). However, a recent study based on a compilation of published data throughout the pantropics surprisingly found that the forest carbon sequestration potential of Asian tropical secondary rainforest was in fact much lower than in American and African rainforests. This work led to a recent refinement of the 2006 IPCC guidelines for national greenhouse gas inventories (Requena Suarez et al., 2019; IPCC, 2019). Whether these new estimates are representative of Asian tropical rainforests is highly uncertain, due to a critical lack of data for this region. This issue is especially crucial for Asian tropical forests where half of the forests have been disturbed during the last century, resulting in the dominance of secondary forests throughout the region (Achard et al., 2014; Mitchard et al., 2013; Stibig et al., 2014).

Remote sensing technology has emerged as a promising tool for extrapolating local field carbon estimates over landscapes, regions, or at the global scale (Gibbs et al., 2007; Goetz et al., 2009). However, current long-term (>20-year) satellite data such as Landsat are weakly sensitive to forest carbon, especially in high-biomass forests (Ferraz et al., 2018; Lu, 2006; Meyer et al., 2019; Zheng et al., 2004). Yet, these data may be used to produce reliable land-cover classifications (e.g., forest versus non-forest areas; FAO, 2010). They allow for assessing the dynamics of deforestation and reforestation worldwide (Hansen et al., 2013) and can thus monitor disturbance history, particularly the time since abandonment of agriculture (Cohen et al., 1996; Masek and Collatz, 2006). However, the forest carbon dynamics associated with such deforestation and reforestation events remains highly uncertain due to the large uncertainties of global carbon maps (Mitchard et al., 2013, 2014; Réjou-Méchain et al., 2019).

On other hand, airborne laser scanning (ALS) provides accurate landscape-scale estimates of forest structural parameters (Maltamo et al., 2005; Næsset, 2002; Wulder et al., 2012). When calibrated with field-based estimates of aboveground biomass (AGB), ALS metrics can be used to produce high-resolution forest carbon maps, even for highly carbon-dense tropical forests (Asner et al., 2010; Cao et al., 2016; Ferraz et al., 2018; Kronseder et al., 2012; Labriere et al., 2018; Zhao et al., 2009; Zolkos et al., 2013). Multi-temporal ALS acquisitions may thus provide direct estimates of the carbon balance of tropical forest landscapes (Dubayah et al., 2010; Meyer et al., 2013; Réjou-Méchain et al., 2015). However, due to its relatively recent emergence, ALS technology cannot yet be used to investigate long-term dynamics directly (>10 years).

Combining long-term (>40-year) land cover change assessment from satellite data archives (e.g., Landsat) and contemporary lidar AGB maps may be a promising avenue for understanding long-term forest carbon dynamics. Such an approach has been successfully implemented in temperate and boreal forests (Bolton et al., 2015; Pflugmacher et al., 2012, 2014; White et al., 2018; Zald et al., 2014). However, to our knowledge, it has not been yet used to assess the forest carbon resilience of tropical forests (but see Helmer et al., 2009, who used satellite-based lidar).

In this study, we combined extensive field, ALS, and Landsat time series (LTS) data to assess the spatial variation of AGB and forest AGB dynamics of secondary forests in a Thai landscape. More specifically, we first calibrated a robust ALS-AGB model to produce a fine-scale AGB map at the landscape scale. We then used a random forest machine-learning algorithm to classify historical Landsat images from 1972 to 2017 into forest and non-forest classes. Using this information over time, we generated a cumulative forest gain map over a period of 42 years of succession. We finally combined this chronosequence with our ALS-AGB map to estimate the forest carbon resilience of secondary forests during the 42 first years after land abandonment.

Figure 1 Location of the study area in Thailand (a) and in the Khao Yai reserve (b). The central map (c) illustrates the lidar top of canopy height (TCH) in the study area at 1 m resolution and the location of the 70 studied plots (in black). Examples of the different stand development stages are illustrated (d–f; SIS: stand initiation stage; SES: stem exclusion stage; and OGS: old-growth stage).

2 Materials and methods

2.1 Study area

The study area of ca. 6400 ha is part of Khao Yai National Park in central Thailand (14^∘25^′20.4^′′ N, 101^∘22^′36.9^′′ E; Fig. 1). Khao Yai is the first national park established in Thailand, in 1962. It is home to numerous endangered plant and animal species (Kitamura et al., 2004). The area receives approximately 2200 mm of precipitation annually, with a dry season of 5 to 6 months (precipitation below 100 mm month⁻¹) from November to April (Brockelman et al., 2011; Chanthorn et al., 2016). The annual mean temperature is about 22–23 ^∘C (Jenks et al., 2011), and the altitude of the study area varies from 650 to 870 m. Before establishment of the park, some areas were used for low-intensity agriculture activities that likely started at the end of the 19th century (Brockelman et al., 2011, 2017) and then naturally reforested at different times depending on when burning ceased (Chanthorn et al., 2016). As a consequence, the landscape constitutes a mosaic of secondary forests of different ages amidst old-growth forests (Chanthorn et al., 2016).

2.2 Field data

We used three sets of forest inventory plots with a total sample area of 35 ha (Fig. 1). First, a large 30 ha contiguous (500 m×600 m) forest dynamics plot, named Mo Singto, was established in old-growth forest after 1998 and completely censused in 2004–2005, 2010–2011 and 2016–2017. The census method follows the protocol of the Center for Tropical Forest Science (CTFS) network to which the plot belongs since 2009 (Brockelman et al., 2011). The second set of plots included eight separate 0.48 ha plots (60 m×80 m) that were established from March to May 2013 and re-censused from November 2017 to January 2018 (Chanthorn et al., 2017). These plots are set along a successional gradient varying from near stand initiation to old-growth forest. Lastly, a 1 ha plot (100 m×100 m) located near the north border of the 30 ha Mo Singto plot was established in a secondary forest in 2005 and then re-censused in 2010 and 2017. In all plots, trees ≥1 cm in diameter at breast height (dbh) were tagged, identified to species, mapped and measured for their diameter, except in the 0.48 ha plots where the minimum dbh was 4 cm. A total of 184 239 individual trees were measured across all the plots, from which 517 trees were measured for height using a pole for short trees (<5 m), a laser range finder (Nikon Forestry 550) for medium height trees (5–7 m) and a Vertex III hypsometer for tall (>7 m) trees (Chanthorn et al., 2017). In this paper, we used the 2017 census data, concomitant with the ALS campaign, to estimate AGB and multiple censuses to estimate the AGB dynamics of secondary plots. For the sake of homogeneity in tree measurements, we used only trees ≥5 cm in dbh in the whole data set.

In order to homogenize plot size, we subdivided all plots ≥1 ha into 0.5 ha subplots. This resulted in 70 plots of either 50 m×100 m (n=62) or 60 m×80 m (n=8) that we classified in three successional stages from young- to old-growth forests following the classification from Chanthorn et al. (2017): stand initiation (early) stage (SIS; n=3); stem exclusion (intermediate) stage (SES; n=5); and old-growth stage (OGS; n=62). Based on interviews of senior park rangers and using Landsat remote sensing data, Chanthorn et al. (2017) estimated that the ages were approximately 15–20 years for SIS forests, 35–40 years for SES forests and unknown but probably older than 200 years for OGS forests. This classification into successional stages followed the framework of Oliver and Larson (1996) who studied successional gradients in temperate forests. Although the original framework considered four successional stages, we did not find any area corresponding to the understory re-initiation stage in the study landscape, i.e., the stage following SES and preceding OGS. Most second-growth forests have regenerated since the Park was established about 50–60 years ago so that old second-growth forests, where understory re-initiation occurs, are very rare in this area. Note also that our study period (1972–2017; see below) cannot account for forests from the SES stage older than 40 years, e.g., that directly started regenerating at the establishment of the park in 1962, as suggested by some hand-drawn historical maps (Cumberlege and Cumberlege, 1963; Smitinand, 1968).

2.3 ALS data

The airborne laser scanning (ALS) campaign was conducted on 10 April 2017 over ca. 64 km² (Fig. 1). The Asian Aerospace Services Limited company (Bangkok) acquired the ALS data with a RIEGL LMS Q680i installed into a Diamond Aircraft “Airborne Sensors” DA-42 fixed-wing airplane. The flying altitude was about 500–600 m above ground level with a 60^∘ field of view, and a pulse repetition frequency of 400 kHz, for which the aircraft maintained an average ground speed of 185 km h⁻¹ capturing the area of interest in 50 overlapping laser strips. We discretized the full waveform data for subsequent analyses resulting in an average point density of ca. 22 points m⁻².

Post-processing of ALS data and point cloud classification into ground, vegetation or noise was done using TerraScan of Terrasolid Version 14. Points classified as ground were used to build a digital terrain model (DTM) at 1 m resolution using a k nearest-neighbor kriging approach implemented in the LidR R package (Roussel and Auty, 2017). A 1 m resolution canopy height model (CHM) was then computed from the height of the normalized vegetation points, discarding outliers classified as air or noise. Finally, we used the CHM and the normalized vegetation point cloud to derive different forest height metrics at the plot level (Table S1 in the Supplement).

2.4 Landsat data

We retrieved Landsat images (MSS, TM, OLI and TIRS products) for the study area from the Landsat archive (http://glovis.usgs.gov, last access: 17 April 2019) in the 1972–2017 period (WRS-1 138/50 and WRS-2 path/row: 129/50). To minimize the impact of clouds and potentially varying phenology within years, we mostly selected images acquired during the dry season, from November to March. We thus collected Landsat 1–3 MSS data (1972–1983), Landsat 4–5 TM (1984–2011) and Landsat 8 OLI & TIRS (2013–2017) data. We did not consider Landsat 7 ETM+ images due to the failure of the Scan Line Corrector, leading to data gaps. All Landsat images were already orthorectified and displayed an accurate co-registration with ALS data. Before 1984, Landsat MSS collected data at 60 m×60 m spatial resolution in most bands. Thus, to have consistent time series data, we aligned all the post-1983 Landsat data using a reference image from 1972 and aggregated each image to 60×60 m. Over the 44 years, we selected a total of 34 high-quality images, each representing 1 year. For the 11 missing years, we could not find cloud-free images and no image was available in 2012 since we discarded Landsat 7 ETM+ data.

2.5 Field aboveground biomass calculation

We estimated tree aboveground biomass (AGB) using a pantropical allometric model (Eq. 4 from Chave et al., 2014). This model uses the diameter (D), total tree height (H) and wood density (WD) as the predictors and was shown to hold across tropical vegetation types and regions. Wood density was estimated using species- (47 % of stems), genus- (50 %) or stand-averaged (3 %) values from the global wood density database (Chave et al., 2009; Zanne et al., 2009). Tree height was estimated through locally adjusted height–diameter (H-D) models of the form given in Eq. (1):

$$\begin{array}{}\text{(1)}& \ln \left(H\right)=a+b\times \ln \left(D\right)+c\times \ln {\left(D\right)}^{\mathrm{2}}+\mathit{\epsilon },\end{array}$$

where a and b are model parameters to be adjusted and ε is a normally distributed error with mean 0 and standard error σlogH. Tree height was subsequently estimated using the back-transformation formula including a known bias correction (Baskerville, 1972) using following Eq. (2):

$$\begin{array}{}\text{(2)}& H=\exp \left(\mathrm{0.5}\times {\mathit{\sigma }}_{\log H}^{\mathrm{2}}+a+b\times \ln \left(D\right)+c\times \ln {\left(D\right)}^{\mathrm{2}}+\mathit{\epsilon }\right).\end{array}$$

Because the H-D relationship varies along the successional gradient (Chanthorn et al., 2017), we fitted three independent H-D models for the three different successional growth forest stages using 177 measured trees for SIS plots, 159 for SES plots and 181 for OGS plots.

AGB at the plot level was then estimated in Mg ha⁻¹ by summing individual tree AGB for all trees with dbh ≥5 cm belonging to the plot. We did all these analyses using the R BIOMASS package (Réjou-Méchain et al., 2017).

2.6 Lidar AGB model

We relied on a log–log model form given in Eq. (3) to model AGB from ALS data (Asner et al., 2012; Réjou-Méchain et al., 2015):

$$\begin{array}{}\text{(3)}& \ln \left(\mathrm{AGB}\right)=a+b\times \ln \left(\mathrm{L}\mathrm{1}\right)+c\times \ln \left(\mathrm{L}\mathrm{2}\right)+\mathrm{\dots }+\mathit{\epsilon },\end{array}$$

where L1, L2, … are the lidar metrics to be selected (see Table S1) and ε the error term assumed to be normally distributed with zero mean and residual standard error σlogL. Fitting the model with log-transformed variables allows us to model a multiplicative error and thus to account for higher model prediction error with larger AGB values (Zolkos et al., 2013). Using this model, we selected the most predictive lidar metrics from our full set of lidar metrics using a leave-one-out cross-validation (LOOCV) scheme nested within a forward selection procedure. The LOOCV consists of fitting models with all observations except one, and then using the model to predict the value of the observation held out of model calibration. The process is repeated for all observations so that model prediction accuracy, here the root mean squared error (RMSE), can be independently assessed from all observations. This LOOCV approach was repeated for different models following a forward procedure that begins by selecting the most discriminant variable according to the LOOCV-RMSE criterion. The procedure then continues by selecting the second most discriminant variable and so on. To minimize the problem of model overfitting, we only kept explanatory variables that contribute to a decrease in relative RMSE (RMSE divided by the mean observed AGB) by more than 1 %. The selected lidar-AGB model was then used to predict AGB values over the landscape at 60 m resolution, to match the resolution of Landsat images.

2.7 AGB recovery analysis

2.7.1 Forest and non-forest classification

To classify areas as forest or non-forest, we applied the random forest (RF) algorithm independently on each Landsat image to minimize inter-image classification error that may otherwise arise from instrumental (e.g., differences in sensors spectral characteristics) and phenological effects. We used all Landsat bands and their ratios as predictors in our RF classification models, i.e., the four raw bands for Landsat 1–3 MSS data (1972–1983), the seven raw bands for Landsat 4–5 TM (1984–2011) and the nine raw bands for Landsat 8 OLI & TIRS (2013–2017). The normalized difference vegetation index (NDVI) was additionally used as a predictor for all the sensors while the normalized burn ratio (NBR) was only used for Landsat 4–5 and Landsat 8 due to non-availability of short-wave infrared (SWIR) bands in MSS sensors. Thus, we used 18 predictors for MSS, 51 predictors for TM and 83 for OLI & TIRS as an input for the RF algorithm. RF model for each year of the study period was then trained on the same set of pixels that likely remained either forested or non-forested from 1972 up until 2017. This training data set was built using the 2017 ALS data. We first aggregated the 1 m lidar-derived CHM at the same resolution as the Landsat images (60 m resolution) and defined non-forest pixels as pixels having a mean top of canopy height of <5 m (FAO, 2012; Sasaki and Putz, 2009). Because 60 m-scale deforestation is unlikely to have occurred in the area since the establishment of the national park in 1962, areas that were classified as non-forest with the 2017 lidar data very likely corresponded to non-forested areas during the whole study period. By contrast, we defined as forested areas all 60 m pixels that had a lidar mean top of canopy height of >30 m because these tall forests very likely corresponded to forested areas during the whole study period. We thus used a reference set of 400 60 m pixels classified as non-forest and 110 as forest. This data set was then randomly divided into a training data set (60 %) to calibrate the RF models and a validation data set (40 %) to assess the accuracy of the forest and non-forest classification. We only considered classified pixels that had a post-probability of assignment >90 % in the RF outputs (Pickell et al., 2016; White et al., 2018) and calculated the classification accuracy as the proportion of pixels that were correctly classified as forest or non-forest in the validation data set. This statistical analysis was done using the randomForest R package (Liaw and Wiener, 2002).

2.7.2 Forest AGB recovery analysis

We combined time-series-classified Landsat images with the 60 m resolution lidar AGB map to quantify AGB recovery as a function of time. We used classified time series data to assign to each pixel the last date at which a shift from a non-forest to forest status occurred during the study period. Thus, all pixels that did not experience any shift, i.e., that remained non-forested or forested during the whole study period, were discarded from this analysis. To minimize false detection of land cover change due, for example, to atmospheric pollution, we only considered shifts that entailed land cover change for at least two consecutive images. Thus, we did not consider any shift before 1975 because, to be considered, the non-forest to forest shift of a pixel should occur after being classified as non-forest in the two previous images (in our case in 1972 and 1973). Finally, we also discarded pixels that underwent more than four different shifts during the whole study period because numerous shifts are likely to indicate areas prone to forest degradation, e.g., close to human occupancy areas such as roads, introducing a bias in our inferences on the natural successional dynamics. We thus assigned to each pixel the year of the last non-forest to forest shift, if any, and considered this year as the forest recovery starting time. The AGB predicted by the lidar AGB map in 2017 was then used to estimate how much AGB was stored between the forest recovery starting time and 2017 through a non-linear power model.

Figure 2 Lidar-AGB model showing the relationship between field-derived plot AGB and the lidar top-of-canopy height (TCH) at 0.5 ha resolution. The power model is illustrated by the red line, and the points represent the field plot AGB estimates at different successional stages: stand initiation (early) stage (SIS; n=3); stem exclusion (intermediate) stage (SES; n=5); old-growth stage (OGS; n=62) according to the classification by Chanthorn et al. (2017).

3 Results

3.1 Forest biomass stocks

Field plots AGB ranged from 80 to 577 Mg ha⁻¹ (mean of 315 Mg ha⁻¹), with a mean AGB of SIS, SES and OGS plots of 87, 291 and 328 Mg ha⁻¹, respectively. Among all the lidar metrics, the mean of top-of-canopy height (TCH, defined as the maximum height of 1 m resolution pixels) was the best predictor of field AGB estimates with a relative RMSE of 14 % (RMSE=45 Mg ha⁻¹; R²=0.85) at 0.5 ha scale (Fig. 2). Adding a second predictor did not reduce the relative RMSE by more than 1 % (Table S2). We thus kept TCH as a single predictor for our analyses resulting in the Eq. (4) for the lidar-AGB model:

$$\begin{array}{}\text{(4)}& {\text{AGB}}_L=\mathrm{4.30}\times {\text{TCH}}^{\mathrm{1.39}}.\end{array}$$

Using this lidar-AGB model, we predicted AGB over the whole landscape (Fig. 3a). The distribution of AGB values over the landscape was not normally distributed due to an over-representation of pixels with low AGB values. At the landscape scale, predicted AGB ranged from 0 to 681 Mg ha⁻¹ with a mean of 291 Mg ha⁻¹ (Fig. 3b), close to the mean AGB of the field plots.

Figure 3 Lidar-AGB map and the distribution of AGB values over the landscape at 60 m resolution. (a) Spatial distribution of AGB predicted from the lidar-AGB inversion model over the study area; (b) density distribution of predicted AGB over the landscape.

Figure 4 Landsat time-series-derived map showing non-forest-to-forest change over the study area. (a) Map showing spatialized-selected pixel shifts from non-forests to forests over the years. The shade gradient represents pixels that did not experience any shift (permanently forested or deforested) and pixels that experienced a shift but that did not pass our quality procedure during the study period (not selected). (b) Density distribution of selected pixel shifts over the landscape during the study period.

3.2 AGB recovery analysis

Our forest and non-forest classification through time was highly accurate, with 90 % to 99 % of well-classified validation pixels in individual classified images (Table S3, Figs. S1–S2 in the Supplement). Figure 4a illustrates the resulting spatialized time series of non-forest-to-forest shifts over the study area and showed that most (83 %) of the landscape did not experience such a shift at 60 m resolution. About 78 % and 5 % of the study area remained permanently forested and non-forested over the 42-year study period, respectively. Most of the stable non-forested areas correspond to National Park building areas, including tourist shops and guesthouses or to continuously cleared areas such as camping locations. Over the 17 % remaining pixels that experienced a shift, we concentrated our analyses on the 4 % pixels (n=550; ca. 198 ha) that passed our selection procedure and that were well distributed over the landscape (Fig. 4b).

Figure 5 Relationship between forest biomass estimated from a lidar-AGB model and forest recovery time estimated from a time series of classified Landsat images (grey dots). The fitted power model is represented by the red line. Blue lines and dots represent the AGB directly estimated from eight field plots (same plots are joined by a line) in 2013 and in 2017/18 and for which the forest recovery time was inferred from Landsat-derived forest age (Fig. S5).

Considering the selected pixels that experienced a shift from non-forest to forest, we found that AGB accumulated non-linearly through time during the 42 first years of the succession (Fig. 5). A simple power model led to a pseudo-R² of 0.66 and a power exponent of greater than 1, indicating an increase in the rate of AGB accumulation with recovery time. This model predicts an AGB gain of 143 Mg ha⁻¹ after 20 years of recovery and of 273 Mg ha⁻¹ after 40 years (spatialized gain in AGB is shown in Fig. S3). Using field AGB estimates at two different census dates from eight secondary forest plots that started regenerating during the study period (see Fig. S5), we showed that the observed rate of AGB accumulation was similar to the one predicted by our model and also tended to increase with forest age (in blue dots in Fig. 5). Finally, focusing on the 17 % pixels that experienced at least one shift from non-forest to forest since 1972, we estimated that the study area has stored a minimum AGB of 455 Gg, as observed in the 2017 lidar AGB map, equivalent to 214 Gg C during the study period.

4 Discussion

In this study, we showed that the integration of field inventory, Landsat archives and lidar data provides a powerful approach for characterizing the spatio-temporal dynamics of aboveground biomass in tropical forests. While the carbon stocks and recovery potential of Southeast Asian tropical forests are globally poorly known, our approach contributes to a better understanding of the role of these forests in global carbon dynamics. We specifically showed that our study site stores a large amount of carbon, despite its disturbance history, and acts as a strong carbon sink, through secondary succession pathways.

4.1 Spatial variation in AGB

Using extensive field data, we have shown that forest AGB can be predicted with an error of 14 % at a 0.5 ha resolution using a single lidar metric, the mean top-of-canopy-height (TCH), a metric previously identified as a robust predictor of AGB (Asner and Mascaro, 2014). This error typically falls within the range of expected errors at this resolution (Zolkos et al., 2013). Using a robust metric selection approach, we showed that adding any other lidar metrics did not bring any additional information and that our single predictor did not show any saturation for large AGB values. Many studies have used a combination of several lidar metrics selected through less robust approaches, i.e., not through independent validation approaches such as our LOOCV procedure, potentially generating overfitting problems (Junttila et al., 2015). We here confirm, similarly to Asner et al. (2012) and Réjou-Méchain et al. (2015), that simple parsimonious models should be preferred, at least within a given tropical forest landscape. Due to a limited number of field plots in low-biomass areas, we were, however, unable to test whether model prediction error varied with forest stand AGB.

Using this lidar model, we predicted a mean AGB over the landscape of 291 Mg ha⁻¹, corresponding to a carbon density of 137 Mg C ha⁻¹ (using a ratio of biomass to carbon conversion of 0.47; Thomas and Martin, 2012). Using a large network of field plots, a recent pantropical study suggested that Southeast Asian and African forests store significantly more carbon than Amazonian forests (Sullivan et al., 2017). However, in this latter study, Southeast Asian forests were only represented by field data from Indonesia and Malaysia where trees are known to be particularly tall (Coomes et al., 2017; Feldpausch et al., 2011; Jucker et al., 2018). Here, we found that our study forests stored significantly less carbon than forests in Indonesia and Malaysia, where the mean carbon density reached ca. 200 Mg C ha⁻¹ (Sullivan et al., 2017), but as much as in Amazonian forests (mean of 140 Mg C ha⁻¹; Sullivan et al., 2017), even when considering only old-growth forest plots. Whether the relatively low carbon density of our study site, compared to other Southeast Asian forests, is specific to our study area or representative of other Southeast Asian forests should be further investigated.

We found a very high spatial heterogeneity of AGB at the landscape scale with an apparent over-representation of low AGB values. This is most probably the consequence of past human activities in this area up to the establishment of the park that led to the present mosaic of secondary and mature forests. This result indicates that this area is currently likely to be a net carbon sink.

4.2 AGB recovery through time

Combining classified images obtained from LTS and lidar data, we quantified the recovery rate of forests after land abandonment. As expected, we showed a significant increase of AGB with recovery time. After 20 years of recovery, our model predicts an AGB accumulation of 143 Mg ha⁻¹, an estimate slightly higher than the one predicted by Poorter et al. (2016a) in Neotropical secondary forests (122 Mg ha⁻¹). However, this difference can partly be explained by the inclusion of trees between 5 and 10 cm dbh in our study, contrary to the study of Poorter et al. (2016b). AGB accumulation in our study corresponds to a net carbon uptake of 3.4 $\mathrm{Mg}\mathrm{C}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ for the first 20 years. This rate of carbon accumulation is close to the pantropical estimate from Silver et al. (2000) and is similar to the default continent recovery rates given by the previous 2006 IPCC guidelines for national greenhouse gas inventories (IPCC, 2006). However, the 2019 refinement of these guidelines halved the recovery rate estimate for young Asian secondary rainforest (≤20 years; Requena Suarez et al., 2019; IPCC, 2019), suggesting that young secondary forests in Asia store carbon at a much lower rate than in Latin America or in Africa. This new estimate is derived from a very limited data set (seven chronosequences) that may not be representative of Asian tropical rainforests. Besides, these data included very small field plots (≤0.01 ha in size; Hiratsuka et al., 2006; Ewel et al., 1983), potentially leading to important sampling errors (Réjou-Méchain et al., 2014). Given the serious implications of these updated IPCC default rates for Asian countries, we here call for further testing of these new IPCC rates across tropical Asia.

Our model showed that a non-linear power model with an exponent of >1 best fit our data, suggesting an increase in the rate of carbon accumulation during the first 42 years of succession. Contrary to the results found by Feldpausch et al. (2007), the rates of AGB accumulation inferred with our approach provided estimates similar to those obtained from long-term field plot surveys (Fig. 5), validating the chronosequence approach in our study area. Assuming that the carbon recovery rate rapidly decreases after 50–60 years (Brown and Lugo, 1990; Silver et al., 2000), our result suggests a sigmoid relationship of AGB accumulation with time in our study area. Previous studies have shown different models of AGB accumulation with forest age. Saldarriaga et al. (1988) showed that the AGB of Neotropical forests from the upper Rio Negro increased linearly with stand age during the 40 years, while Jepsen (2006) reported a sigmoidal increase in AGB accumulation in Sarawak, Malaysia, as is likely the case in our study area. Finally, working on 41 Neotropical sites, Poorter et al. (2016a) assumed a logarithmic trend in the AGB accumulation over time, hence a decrease of the rate of carbon accumulation through time, probably because they investigated a longer time period. Selecting the sites of Poorter et al. (2016a, b) that had at least 10 observations over the first 44 years (n=21 out of 28 sites, i.e., excluding 7 sites for which model fitting was not possible), site-specific power models revealed that two-thirds of the sites displayed a power exponent of <1 and one-third showed an exponent of >1 (Fig. S4). Thus, the accumulation of AGB with age follows different trends across sites, as already highlighted in previous studies (Kennard et al., 2002; Poorter et al., 2016a; Ray and Brown, 2006; Ruiz et al., 2005; Silver et al., 2000; Toledo and Salick, 2006). Understanding how these trends vary according to abiotic factors (e.g., soil type, rainfall), species assemblage and diversity, and prioritizing effects such as types of land use and land management existing before forest recolonization, constitute an important avenue of research (Chazdon, 2014; McMahon et al., 2019).

Our analysis was based on a forest/non-forest classification through time and our independent validation suggested high overall accuracy (90 % to 99 %), similar to that reported by other studies using Landsat data classification in boreal systems (Bolton et al., 2015; White et al., 2018). Furthermore, our estimate of forest age using this approach was highly consistent with our expectations. Indeed, using our forest plots, we found that the SES and SIS forest stages lasted on average 40 years (range 38–42) and 13 years (range 8–20), respectively, hence very close to suggestions of Chanthorn et al. (2017; Fig. S5). However, our overall approach cannot be replicated easily in human-occupied areas. Indeed, human disturbances lead to forest degradation that, in contrast to deforestation, is not captured by the Landsat signal, so that, when combined with a reference AGB map, natural carbon recovery potential could be seriously underestimated. Because our study area was protected from human disturbances during the study period, we were in very favorable conditions to estimate forest carbon recovery rates and strongly encourage researchers benefiting from similar conditions to replicate our analyses in other study sites.

5 Conclusions

Our study demonstrates that combining field, lidar and long-term satellite data provides an efficient way to assess forest carbon recovery rates during secondary successions. We showed that it produces similar estimates as those inferred from long-term field plots, but at a much lower cost and within a much shorter time frame. Replicating this approach in other protected tropical landscapes, notably in the Asian subcontinent, would thus considerably increase the representativeness of forest carbon recovery rates. This would improve our understanding of the environmental and historical drivers of these varying rates between ecological zones and continents. This is especially important in Southeast Asian forests that constitute a hotspot of biodiversity and carbon, and that are under threat due to the fast changing of both the environment and socioeconomics in this region. Quantifying the rates at which different forest types accumulate carbon should thus stay at the forefront of the research agenda and would greatly benefit the Earth system model community and international policy initiatives such as REDD+.