Modeling Secondary Forest Regrowth and Carbon

To model the regrowth of secondary forests we applied a space-for-time substitution method. Instead of tracking the associated Aboveground Carbon (AGC) regrowth over time, the regrowth was estimated by considering the available ages of the standing secondary forest area in 2017 and the associated AGC at the same time. Here we explain the methods used to determine secondary forest AGC using the ESA-CCI Aboveground Biomass (AGB) product (100-m) for the year 2017²³ (see Supplementary Notes 1 and 2). All analysis was carried out in the original product units (AGB) but expressed as AGC by assuming a 2:1 ratio of biomass to carbon²⁴ .
The ESA-CCI AGB product was only released in late 2019 and was in its early phases of development at the time of use. However, given that its spatial resolution was high enough to separate areas of only secondary forest and its recent acquisition warranted its use for this research. Only areas of secondary forest greater than 9000 m² were considered for further analysis, an area approximately equal to 1 pixel of the ESA-CCI product. Despite limiting the study to these larger secondary forest polygons, we were still left with just under 2.5 million polygons of secondary forest to analyse. The secondary forest map was laid over the AGC data and the modal AGC was extracted for each secondary forest polygon using the “zonal_stats” function available in the “rasterstats” module for the programming language “Python” (v3.6). We then aggregated the AGC values by the age of secondary forest and used the median AGC value for each age in further analysis. We applied a bias correction to the median AGC values, subtracting the smallest median value from all values to shift the data to begin at or near 0 Mg C ha⁻¹ AGC for a 1-year-old secondary forest.
Following this, we used six remote sensing products of driving variables widely accepted to influence regrowth of forests. The data products included four environmental drivers (1–4) and two anthropogenic disturbance drivers (5–6): (1) Mean annual downward shortwave radiation (for the period 1985–2017)^{26 (link)}, (2) Mean annual precipitation (for the period 1985–2017)²⁷ , (3) the mean Maximum Cumulative Water Deficit (MCWD) (for the period 1985–2017)^{65 (link),66} , (4) Soil Cation Concentration^{30 (link)}, (5) Annual burned areas (between 2001 and 2017)³¹ and (6) Number of times a secondary forest area was deforested between 1987 and 2017 (repeated deforestations) (this study). These products all have different spatial resolutions (Supplementary Table 1) and so had to be resampled to the size of secondary forest pixels (30-m spatial resolution) using the “resample” package in the Geographic Information System programme, ArcMap10.6. We calculated the key zonal statistics of these variables such as the mean value of the driver affecting a specific area of secondary forest, again using the “zonal_stats” function in Python.
The drivers were then grouped according to numerical limits, such as the 25, 50 and 75th percentiles. We then modelled the AGC for the age of secondary forest under these groupings using the commonly used Chapman-Richard model for regrowth^{67 (link)}: $Y_t=A{\left(1-e^{-kt}\right)}^c\pm \epsilon ;A,k\mathrm{and}c>0$
where Y_t refers to the AGC at age t; A is the AGC asymptote or the AGC of the old-growth forest; k is a growth rate coefficient of Y as a function of age; c is a coefficient that determines the shape of the growth curve; and ε is an error term. We applied the “nls” function available in the “nlstools” package for the statistical software R (v4.0.2)^{68 ,69} . We assumed that after a given amount of time, the AGC could return to levels equivalent to old-growth forests, and reach a pre-calculated asymptote. As such, we extracted the median, bias-corrected AGC value of old-growth forests under each variable condition from the ESA-CCI AGC product to represent the value of the asymptote (Supplementary Fig. 6 and Supplementary Tables 8 and 9). From this, we could also determine if and when the modelled AGC of secondary forest regrowth would reach those equivalent to old-growth forest levels. Forcing the models to “fit” to an expected value for the asymptote value naturally increases the error of our model, partly due to heterogeneity in old-growth forest values within each variable condition.