Rapid tree carbon stock recovery in managed Amazonian forests
by Rutishauser et al. (Current Biology 25:R787–R788)
Supplemental Information
Forword
This page provide the R code and supplemental information regarding the publication Rapid tree carbon recovery in Amazon managed forests published in Current Biology on 21 September 2015 editor’s link
The code below enables to replicate the analysis, produce all graphs and summary tables. A supplementary analysis is also carried out and provided at the end.
Summary data are freely available here Raw data are available upon request: contact@tmfo.org
1 Site selection and biometric data collection
Ten (10) sites spread across the Amazon basin and the Guiana Shield were selected based on the following criteria: (i) located in tropical forests with a total area inventoried ≥ 1 ha; (ii) mean annual rainfall ≥ 1000 mm (Fig. S1); (iii) consistent and detailed information about logging treatments (e.g. number of stems harvested and correspondent biomass removal) and logging impacts (e.g. logging damages assessment); (iv) at least one pre-logging and (v) at least two post-logging censuses. As sites were generally established by different organizations, there is no standardized protocol for data collection among sites, but all sites comply with generally agreed standards [1]. A general description of the sites can be found in [2]. In all plots, trees ≥ 20 cm DBH (diameter at breast height) had their girth measured at 130 cm or above buttresses/deformations, and were tagged and identified to the lowest taxonomical level.
Figure S1: Representativity of TmFO’s Amazonian network (coloured circles) at regional level (grey-scale circles: GS: Guyana shield, CA: Central Amazon, EA: East Amazon, WA: West Amazon, SA: South Amazon): a) Rainfall gradient, b) soil fertility (CEC) and d) soil texture.
2 Data quality checking and biomass computation
To avoid bias due to discrepancies in data quality (e.g. difference in botanical identification or tree species wood density information), a standardized protocol was applied to each site. At first, botanical identification was checked to match the Global Wood Density Database (GWDD) classification [3]. Tree species present in GWDD were assigned correspondent dry wood density (WD, gr.cm-3). When only the genus was present, genus-average WD was assigned and for unidentified species and species not present in the GWDD, plot-average WD was attributed. In the absence of tree height measurements, tree above-ground biomass (AGB) was estimated using the generic allometric model developed by Chave et al. [4] and including WD, DBH and a synthetic climatic index (E).
Above-ground carbon density (ACS) was obtained by multiplying tree biomass by 0.47 [5]. ACS stock of each plot was further computed as the sum of ACS of live trees DBH ≥ 20 cm divided by the plot surface and expressed in Mg C ha-1.
3 Definition of logging intensity and biomass recovery
The same definition of logging intensity was applied at all sites. Due to varying interval length (1 to 4 years) between pre- and post-logging censuses and application of silvicultural treatments (i.e. poisoning, girdling, understorey clearing) at three sites (Paracou, Tapajos and la Chonta), we estimated the minimum carbon stock (ACSmin, Fig. S2) attained at last within 4 years after logging and computed the difference with initial carbon stock (ACS0). This initial ACS drop off, referred to as ACSloss, is due to both timber harvest and mortality of damaged trees (that can affect up to 46% of remaining trees[6]). As residual mortality peaks within the first years preceding logging [7, 8], this approach allows most of logging-induced mortality to be accounted.
Figure S2: Theoretical post-logging biomass recovery. t0 = year of pre-logging inventory, ACS0 = carbon stock at t0, ACSmin = minimum post-logging carbon stock attained within 4 years after logging, tmin = year at ACSmin, tfinal = year of last census and trec = estimated time of recovery.
We found no evidence of strong deviation from linearity of ACS recovery (Figure S3); therefore, we estimated the annualized ACS recovery rates (Mg C ha-1 yr-1) per plot using linear models among all post-logging censuses spreading between tmin and tfinal. Recovery time (trec in years) refers to the estimated time needed to recover initial ACS stock, and is given by the product of initial ACS loss divided by average recovery rate.
4 Explanatory variables
Several explanatory variables were calculated at each site: (1) average pre-logging ACS stock (ACS0 in Mg C ha-1); (2) Basal Area-weighted wood density (or community wood density, WDBA in g.cm-3); (3) stem density (ha-1); (4) average annual rainfall (mm yr-1) that arose from local weather stations; (5) rainfall seasonality (annual standard deviation) were extracted at each site using WorldClim data (Hijmans et al., 2005) using highest resolution (30 arc-seconds or ~1 km). Due to lack of information at all sites, soil properties were extracted from the Harmonized World Soil raster at a resolution of 30 arc-seconds (9). Information on top soil (0-30 cm) quality was extracted at each site: texture, drainage, available water content (range), clay, silt and sand content (%), cation-exchange capacity (CEC, cmol/kg) and bulk density (kg/dm3).
To test for possible circularity between the synthetic climatic index (E) used to compute ACS and the climatic explanatory variables, all analysis were recomputed with another generic allometric model [9], based on local WD and DBH only. All pattern and variables significance remained unchanged (data not shown).
5 Plot selection and weighings
To ensure that observed biomass recovery was mainly related to logging and avoid bias due to stochastic natural mortality (e.g. the 2005 drought and fires), we selected only plots (79 out of 90) with positive recovery rates (e.g. that gain biomass/carbon over the monitored period), as a detailed checking revealed that those 11 plots suffered from wildfires and droughts. As our sample plots and sites vary in both total area and length of time monitored for, the contribution of each site was weighted by the monitoring effort (number of censuses x plot size), as recommended by [14]. Hence, sites with longer and larger monitoring (more prone to capture and depict forest recovery) are given more weight. To avoid artificial inflation of the variance of random effects, the sum of weights was set to 1. Table S1 provides information on initial and final ACS, ACS loss, recovery rate and recovery time for each plot (N=90).
6 Variable selection
Our main point was to understand generic drivers that led recovery time and recovery rate among all sites. We developed a linear mixed model (LMM, package lme4 [15]) in which recovery time and rate were tested over the different biometric response variables defined above. To account for the site effect, we introduced a random site effect. Indeed, most sites are constituted of several contiguous plots in which silviculture treatments (e.g. logging, girdling or understorey clearing) of varying intensities were applied. Such experimental design ensures a relative homogeneity in environmental conditions and forest structure, but might also induce pseudo-replication. Pseudo-replication occurs when multiple samples from a single treatment unit are analyzed, as if they were independent replicates and embed to distinguished the effect due to treatment from other sources of variation [16]. To avoid this bias, a “site-effect” was introduced in the LMM and pre-logging forest structures were accounted for as explanatory variables in the analysis.
The best models are found through conducting an exhaustive screening and ranking using Bayesian Information Criterion (BIC) (package lmerTest [17]). Instead of picking a single “best” model, we averaged the fits of a number of “good” models (model averaging) based on Bayesian Information Criterion (BIC) weights, thereby stressing prediction over precision [18]. Very good fits were effectively found at each site (Figure S3). To reduce residual heteroscedasticity, recovery time along with two explanatory variables (ACS logged and number of trees harvested) was log-transformed.
All analyses were carried out with R language and environment [19].