Carbon sink despite large deforestation in African tropical dry forests (miombo woodlands)

For the change analysis, re-measured plots between ILUA I and ILUA II were identified through spatial overlap in a GIS in order to ensure comparison of the same location. Re-measured plots amounted to 829 out of 839 plots measured in ILUA I and included in 212 1 km by 1 km clusters of four plots (figure 2). To compare the two assessments, we selected the trees of the ILUA I survey based on their position along the central transect, namely the trees located within the first 50 m of the 250 m transect. It was not possible to perform a tree-by-tree comparison which would have required trees to be tagged with a unique identifier. The comparison is done in terms of the sampled area, on a plot-by-plot basis.

When the ILUA I and ILUA II species lists were compared, some spelling errors and species known under different scientific names were noticed. To obtain valid results in terms of woody plant diversity, an unique species list was generated from the two surveys where we removed duplications or spelling mistakes that would have affected the evaluation of change in diversity. The final list was verified and validated by a botanical curator and co-author (N.Z.) and included 238 woody species belonging to 154 genera, showing high diversity and endemism over the country.

We calculated woody plant diversity indices and stand structure variables for each of the surveys, at both the cluster and the plot levels. We calculated common alpha (α) biodiversity indices including tree species richness, Shannon's entropy, Simpson's index, the Shannon and the Simpson diversity numbers, as well as the Pielou, Shannon, and Simpson evenness indices [52]. With the FD package in R [53], we calculated functional diversity indices using important life-history traits, wood density and seed mass, gathered from global data sets [54–56]. Stand structure variables included basal area, stump density, and biomass for all trees and shrubs ≥5 cm diameter at breast height (dbh) as well as stem density for all stems ≥1.3 m in height. The plot biomass was calculated by summing the individual tree biomasses estimated with the allometric models developed by Chave, Rejou-Mechain [57] using the dbh, height and wood density (wd) of each individual. The average of the four plots was used at the cluster level. Then, we evaluated the mean difference and variation in metrics between the two repeated measurements.

The difference in aboveground biomass between the repeated measurements, hereafter called productivity, was considered a key variable to assess TDFs resilience. Because these TDFs are subjected to frequent fires as well as other disturbances, the difference in biomass can be either positive (gain) or negative (loss). Changes associated with belowground components, litter and woody debris were not accounted for because these components were not measured during the first survey.

In addition to the survey data from ILUA I and II, we extracted explanatory variables from land-cover change maps and from global data sets. The land-cover change maps were developed by the ILUA II project for years 2000, 2010, and 2014, using Landsat imagery (30 m pixel). The maps are based on the six Intergovernmental Panel on Climate Change land-cover categories (cropland, forestland, grassland, settlement, wetland and other land). It includes classes that show no change (stable class over time) and classes showing the conversion of forests (appendix S1, figure S1.1).

We calculated statistics from the land-cover change maps, using tabulated areas and zonal statistics in ArcGIS in a 1 km radius of the plot center. We also developed patch and landscape metrics to describe landscape heterogeneity with the 2014 land-cover map, using FRAGSTATS [58] in R. These fragmentation metrics (∼50 metrics) were calculated using a 1 km radius from the cluster center and were developed for the forestland category. One example of fragmentation metric is the 'core area of forestland', which represents the area within a forest patch beyond some specified edge distance or buffer width [59]. We also compiled fire area from the MODIS burned area products (MCD45A1) [60] from 2006 to 2014 in order to evaluate the presence or absence of fire, the fire frequency and the time since fire (appendix S1, figure S1.2). Other variables related to human influence were assessed including the distance to roads, population within 5 km radius and distance to towns.

Bioclimatic variables derived from monthly temperature and rainfall data were extracted at the plot center to represent annual trends, seasonality and limiting environmental variables relevant to tree distributions [61]. We also used soil organic carbon, bulk density, and pH collected at the cluster level through the soil assessment performed in ILUA II. The clay/sand ratio and the total sand were extracted from soil GIS layers provided by the ZFD.

To understand the impacts of different drivers on woody plant diversity and ecosystem functions, we used a structural equation modeling, called path analysis. This statistical approach is a straightforward extension of multiple regression, but in addition to merely explaining variation in a response variable, path analysis attempts to describe and quantify causal relationships [62]. This analysis aimed to provide estimates of the magnitude and significance of hypothesized relationships between sets of variables to understand the impacts of different trajectories of change and capture the main drivers. We first developed a directed acyclic graph (figure S1.3), laying out the different sets of variables and the hypothetical directed links. We used correlation tables to select the most important variables for our set of descriptors. Using the lavaan package in R, we developed a series of models for different trajectories of change based on the field information, including biomass gain over time, biomass loss over time, forestland plots remaining forestland, and forestland plots converted to non-forestland, at the cluster and the plot level.

To compare plot woody plant diversity and stand structure metrics between ILUA 1 and ILUA 2, we used paired t-tests when the variables could be normalized to approximately obey a Gaussian distribution (e.g. with log-transformation); otherwise we used the Wilcoxon paired-sample test (or Wilcoxon Signed-Rank Test).

The fact that the individual trees were not tagged is a limitation of the data set, and we adopted a conservative approach to account for this. This limitation made it more difficult for us to identify potential measurement errors in the data set (e.g. inclusion or exclusion of a tree at edge of a plot). Because we have quantified a biomass increment over the time period, we have purposefully removed sites showing a higher gain that could potentially be due to the inclusion of trees in the second survey that were not measured during the first survey. To ensure robust results for statistical analysis, we have established a threshold to remove the plots with biomass gain >50 Mg ha⁻¹ over the period, which is equivalent to an increment of about 3.5 Mg C ha⁻¹ yr⁻¹. We evaluated that the plots with >50 Mg ha⁻¹ biomass gain have a higher probability to be affected by measurement error. We report summary statistics for the whole data set. For statistical analysis, including path analysis, t-tests, and Wilcoxon paired-sample tests, we have excluded the sites with higher biomass gain for a final sample size of 754 plots. We used growth rates with and without applying the threshold to test for sensitivity of this parameter on carbon absorptions.

To identify which species contributed the most to biomass changes, we ranked the 50 most dominant tree species based on the total biomass for ILUA 1 and for ILUA 2 and compare the difference in biomass for each species for all the plots (n = 754) and also by separating between plots showing a biomass increment (n = 400) and those with no change or a loss (n = 354) during the period.

We compared emissions from deforestation and absorptions by forestlands, contrasting different sensitivity scenarios for biomass increment in forestlands. We used error propagation to account for the different sources of error in land areas, carbon stocks and growth rates and quantify the overall uncertainty around final estimates. For carbon emissions, we used estimates for the period 2000–2010 and 2010–2014 published by the Republic of Zambia [49], in table 6, that we converted from CO₂ into C. To account for carbon sequestration from forestlands, we used the stable forestland area from 2000–2014 published in ILUA II [63] that we multiplied by the biomass increment or growth rate values calculated from the difference between ILUA I and ILUA II. We used sensitivity analysis to test the effect of different growth rates estimates. We included the mean and the median biomass increment for plots of 'forestland remaining forestland' between the two surveys and the mean and median by removing all sites with a productivity higher than 3.5 Mg C ha⁻¹ yr⁻¹. These growth rates include forestland plots with degradation or disturbances, but excluded plots that have undergone conversion from forestland to other land use. Then, we calculated the difference between emissions and absorptions for the 2000–2014 period.

The comparison over time showed statistically significant differences between the ILUA I and II, with a gain in richness, basal area, biomass, and stump density (figure 3). The annual difference between those metrics was positive, with a mean difference of 0.15 ± 0.02 and 0.16 ± 0.03 species per year for richness, 0.12 ± 0.02 and 0.13 ± 0.03 m² ha⁻¹ yr⁻¹ for basal area, 1.00 ± 0.17 and 1.17 ± 0.27 Mg ha⁻¹ yr⁻¹ mean annual biomass gain and 0.26–0.31 Mg ha⁻¹ yr⁻¹ median annual biomass gain, for the plot and the cluster level respectively (table 1). We also noticed a statistically significant decline in evenness (figure 3) and found no significant difference in stem density $(V=109\,830,p \mbox{-} {\rm{v}}{\rm{a}}{\rm{l}}{\rm{u}}{\rm{e}}=0.8229).$

Zoom In Zoom Out Reset image size

Trajectories of change were statistically explained by different human drivers and natural factors. Models results were consistent at the cluster and plot levels. We found persistent patterns amongst trajectories. One commonality was the positive relationship between basal area in the second survey and species richness or productivity. Richness was negatively affected by the temperature in the warmest quarter of the year. Productivity was negatively affected by species evenness (or positively with dominance) and by functional richness. Climatic and edaphic variables had no significant impact on basal area.

When all plots were included, the core area of forestland (a fragmentation metric) and the distance to road variables had a positive effect on basal area and richness, which means that plots located in continuous, non-fragmented forestland and away from roads showed higher richness and basal area, also having an indirect positive effect on productivity (figure 4(a)).

Zoom In Zoom Out Reset image size

In plots exhibiting biomass gain between the two surveys, we identified more significant factors than in other trajectories, including environmental, diversity and anthropogenic (or human-proxy) variables (figure 4(b)). Fire occurrence decreased basal area of plots gaining biomass but was not identified to be a key variable under other trajectories. Core area of forestland was positively linked to richness. The distance to roads was related to higher basal area and productivity. We saw a positive effect of functional richness on basal area. Evenness was most strongly negatively associated with basal area and productivity than in any other trajectory of change. Soil organic carbon was also linked to higher productivity. The explained variation in species richness was lower for this trajectory, indicating that other factors not accounted for are also important.

Under the trajectory 'Forestland remaining forestland', we found causal links for environmental factors with richness and productivity as well as the effect of evenness and functional richness on productivity, but no human-related drivers were significant (figure 4(c)). We found that productivity was enhanced on sites with higher clay to sand ratio.

The opposite was identified for 'Forestland converted to non-forestland' where only a human-proxy variable was significant and where the core area of woody vegetation was negatively related to biomass loss (figure 4(d)). We also developed models for plots with biomass loss but the models did not succeed in converging, possibly due to lower sample size (n = 90).

For plots characterized by a gain in biomass, we looked more closely into what species were responsible for this biomass increment. We found that dominant species, those having the largest biomass stock, also show the largest absolute biomass increment. The species showing the highest biomass gain are Julbernardia paniculata, Brachystegia spiciformis, Brachystegia longifolia, Erythrophleum africanum and, Colophospermum mopane (figure 5(a)). Since we found no significant change in stem density in the re-measured plots, it can be considered that this biomass increment is dominantly linked to the growth of individual stems rather than recruitment.

Zoom In Zoom Out Reset image size

On the other hand, we also looked at plots showing no change or biomass loss to determine which species are the most affected by human activities or other disturbances (figure 5(b)). The greatest losses, in order of importance, are for B. spiciformis, C. mopane, Isoberlinia angolensis, J. paniculata, and B. boehmii. Then, overall plots, C. mopane, I. angolensis, and B. spiciformis reveal a biomass loss despite of the gains, while B. longifolia and J. paniculata show an overall gain.

We found that for the 2000–2014 period and accounting for different sources of uncertainty, the carbon absorptions from stable forestlands largely surpass the carbon emissions from deforestation, thus creating a carbon sink for all the sensitivity scenarios tested (figure 6). This carbon sink from forestlands in Zambia goes from −136 618 to −607 172 GgC (table S1.1), showing that these estimates are very sensitive to the growth rate used to quantify carbon absorptions from stable forestlands. We can also note a high uncertainty around emission estimates, which is driven by the uncertainty around the deforestation area estimates. Emission and absorption estimates are presented in table S1.1 in appendix S1 as well as the distribution of the carbon stock change for the subset of plots of stable forestlands (forestlands remaining forestlands) used in sensitivity analysis figure S1.4.

Zoom In Zoom Out Reset image size

The biomass accumulation rates in 'forests remaining forests' quantified in this study (1.19 ± 0.13 Mg C ha⁻¹ yr⁻¹ mean and s.e.m. or 0.39 ± 0.10 Mg C ha⁻¹ yr⁻¹), is higher than the one quantified by Chidumayo [17] of 0.10 Mg C ha⁻¹ yr⁻¹ for old growth permanent plots of Zambia. These are expected to have lower productivity, however the author indicates that the rates could have been up to about 2.0 Mg C ha⁻¹ yr⁻¹ if degradation due fire impacts and harvesting was avoided. Other studies have quantified similar or higher rates but are not strictly comparable in their methodology and in the type of vegetation they measured. For example, Chidumayo [64] found an increment of 0.97 and 0.82 Mg C ha⁻¹ yr⁻¹ (range of 0.20–1.42 Mg C ha⁻¹ yr⁻¹) for coppiced plots of age 3–29 and 48–49 years old in Zambia, respectively. McNicol, Ryan [14] found an average of 0.83 Mg C ha⁻¹ yr⁻¹ for the miombo woodland following shifting cultivation abandonment in Tanzania. In Zambia, Kalaba, Quinn [65] obtained 0.98 and 1.42 Mg C ha⁻¹ yr⁻¹ in agriculture and charcoal fallows respectively.

Woody plant diversity, and more specifically species dominance, appears to play a key role in this carbon sink. With our path analysis results, we found that productivity is negatively associated to species evenness, which has significantly decline between ILUA I and II. Chidumayo [17] indicates that an increase in evenness in miombo woodlands is a good indicator of forest degradation, so our findings may be indicative of overall resilience. Dominant species of the miombo and mopane woodlands played a particularly influential role in this biomass gain, and thus carbon absorptions. Our results point to biomass increment by individual growth in these species and/or species replacement rather than recruitment, since stem density has not significantly changed for the period.

Our study provides evidence for a significant gain in biomass and a carbon sink across Zambia, but we cannot determine if there is an enhancement in carbon uptake, as we have only two points in time. The most likely interpretation for these results is that Zambia TDFs cover a large area compare to the deforested area, explaining why the absorptions surpass the emissions. While many studies have quantified deforestation in Zambia, we do not know of any study that has looked at the absorptions from TDFs at national-scale. Other countries, including Mexico [66], also show high rates of deforestation but because of their large forests extent, absorptions surpass emissions.

Other hypothetical factors that could explain this biomass gain Zambia's TDFs include changes in local drivers such as land use and management, regional changes in rainfall, in disturbance regimes (fire, grazing), or CO₂ fertilization effect. Based on recent estimates, we know that land conversion has increased between 2000–2010 and 2010–2014. The rainfall amounts in Zambia were variable over the period, showing increases in some areas and decreases in others [67], and therefore are unlikely to be the cause of observed biomass gain. Future study should explore the potential effect of inter-annual atmospheric variability (e.g. ENSO) and intra-seasonal rainfall variability on woody productivity.

Recent human-driven decline in global burned area, including in the African woodland [68], may potentially be linked to the biomass increment that we measured. However, a recent study in Tanzania, a country with comparable vegetation and drivers, shows a consistent fire regime over time for the past ∼40 years [69]. More evidence would be needed to link our results to a potential change in fire patterns.

CO₂ fertilization could be promoting the dominant species growth leading to biomass gain. Species-specific responses to elevated CO₂ concentration have been demonstrated [70], and may depend on the functional traits of woody plants [71]. Some savanna woody species have shown a strong yet different response to elevated CO₂ in experimental studies [41]. Benefits of CO₂ for tree growth and recruitment promote growth of root systems, providing stored energy for re-sprouting after fire and herbivory [72] and favoring clonal reproduction through coppicing and root suckering [73]. Most miombo dominants have the ability to coppice [74, 75]. This capacity is supported by extensive root systems, which spread vertically, with tap roots reaching depths of at least 5 m, and producing root suckers horizontally once the aboveground parts have been removed [74]. CO₂ fertilization generally enhances water use efficiency, possibly lengthening the growing season [70]. This increased water use efficiency has been recently shown for African trees in drier areas [76].

The possibility exists that dominant trees of the miombo woodlands are responding positively to rising atmospheric CO₂ level. However, Stevens, Lehmann [34] suggest that plant nitrogen-fixing ability is critical in predicting woody encroachment in African savannas which is likely caused by rising CO₂ levels. From the encroaching species list in African dry savanna provided by Stevens, Lehmann [34], only Colophospermum mopane is common to the dominant trees gaining biomass in Zambia. The top weight gainers in Zambia are dominated by non-nodulating species. Julbernardia-Brachystegia-Isoberlinia species are characterized by ectomycorrhizal symbiosis [77–79], providing better capacity to respond to adverse conditions through enhanced nutrients and water uptake [79] and promoting the growth of host species [80–82]. Future research should use long tree-ring chronologies to directly assess the potential effect of atmospheric CO₂ on growth of dominant tree species.

The new national-scale repeated measurements conducted at more the 800 sites across Zambia highlights the unique species diversity and functional traits of the miombo woodlands relative to drier African savanna. Most prior research in the African savanna biome has taken place in South Africa, and researchers have commonly considered this biome to be homogenous [30]. Findings from the drier savanna biome are not generalizable to the mesic part where key differences in ecological processes and drivers exist.

Impact from human activities on TDFs were noticeable. We found higher productivity at sites located further away from roads and in relatively un-fragmented forests. Ninety plots of forestland were cleared between the two surveys. Road expansion stimulates deforestation and forest degradation [83–87] by providing access to agricultural expansion and markets.

The effect of fragmentation in this region showing heterogenous landscape patterns has rarely been addressed. The effect of fragmentation in closed humid tropical forest is well understood, with the edge effect increasing tree mortality caused by microclimatic changes, wind damage and invasive species near the forest edges, with associated impacts on biomass [88–90]. In the miombo woodlands, the exact mechanisms that can affect ecological processes other than the direct impact of land-use change are unknown.

Fire had only a significant negative impact on basal area in plots showing biomass gain, indicating that other factors are more important for degrading stands. Under controlled experiments, fire impact on aboveground biomass has been demonstrated by several studies in the miombo woodlands [16, 17, 91]. Fire would act as a demographic bottleneck by disproportionately affecting smaller trees (<2 m in height) and favoring grass layer development. Trees that have survived and reached a certain height are less susceptible to fire-induced mortality [72]. For both ILUA I and II survey, there is an overwhelming presence of surface fires, with very little presence of more damaging crown fire [45]. Future work should build a more fine-resolution long-term fire history (e.g. Landsat 30 m resolution, see Tarimo, Dick [69]) in order to assess the effects of fire and predict its long-term influence on diversity and ecosystem functions.

Heavy wood extraction for charcoal production and construction materials appear to show a negative balance for some of the most productive and dominant species in Zambia. This decline in biomass for B. spiciformis, C. mopane and I. angolensis indicates that despite their high productivity, these species are heavily extracted, a practice that could be unsustainable over time.