Biomass Increases Go under Cover: Woody Vegetation Dynamics in South African Rangelands.

Woody biomass dynamics are an expression of ecosystem function, yet biomass estimates do not provide information on the spatial distribution of woody vegetation within the vertical vegetation subcanopy. We demonstrate the ability of airborne light detection and ranging (LiDAR) to measure aboveground biomass and subcanopy structure, as an explanatory tool to unravel vegetation dynamics in structurally heterogeneous landscapes. We sampled three communal rangelands in Bushbuckridge, South Africa, utilised by rural communities for fuelwood harvesting. Woody biomass estimates ranged between 9 Mg ha(-1) on gabbro geology sites to 27 Mg ha(-1) on granitic geology sites. Despite predictions of woodland depletion due to unsustainable fuelwood extraction in previous studies, biomass in all the communal rangelands increased between 2008 and 2012. Annual biomass productivity estimates (10-14% p.a.) were higher than previous estimates of 4% and likely a significant contributor to the previous underestimations of modelled biomass supply. We show that biomass increases are attributable to growth of vegetation <5 m in height, and that, in the high wood extraction rangeland, 79% of the changes in the vertical vegetation subcanopy are gains in the 1-3 m height class. The higher the wood extraction pressure on the rangelands, the greater the biomass increases in the low height classes within the subcanopy, likely a strong resprouting response to intensive harvesting. Yet, fuelwood shortages are still occurring, as evidenced by the losses in the tall tree height class in the high extraction rangeland. Loss of large trees and gain in subcanopy shrubs could result in a structurally simple landscape with reduced functional capacity. This research demonstrates that intensive harvesting can, paradoxically, increase biomass and this has implications for the sustainability of ecosystem service provision. The structural implications of biomass increases in communal rangelands could be misinterpreted as woodland recovery in the absence of three-dimensional, subcanopy information.

Introduction

Woody biomass is a fundamental expression of terrestrial ecosystem functioning, (e.g. primary productivity, land-atmosphere gas exchange and nutrient regulation), and can be used for the quantification of ecosystem services, such as fuelwood and carbon sequestration. Biomass distribution reflects the spatial pattern of topo-edaphic and climatic gradients [1-3] and responses to disturbance [4-7]. However, biomass estimation remains challenging, particularly in environments with highly variable species composition and structural complexity [8-10].

Savannas, as complex tree-grass ecosystems, are structurally heterogeneous and are best described by three-dimensional metrics [11]. As such, savannas are ideal for examining the biomass dynamics in structurally complex vegetation. While total precipitation sets the upper boundaries on woody cover in savannas [12], their ‘woody cover potential’ is often unrealised [13-14] as a result of disturbances, such as fire [15-19] and herbivory [20-22]. A major driver in savanna ecosystem structure and function is the influence of people on the landscape [15,23], particularly through natural resource use, such as fuelwood harvesting [24]. Yet, the contributions of anthropogenic changes to savanna biomass dynamics are poorly understood.

Millions of people in Africa rely on woody vegetation for energy, extracted from both communal [25-27] and protected areas [28-29]. Within southern Africa, South Africa has a high per-capita use of fuelwood as a primary energy supply; despite having substantial access to electricity (66% of national population) [30]. Within this context, 93% of current fuelwood demands are no longer met by collection of dead wood [31]. Thus, live wood harvesting occurs around settlements and is a major driving force in woodland degradation in semi-arid ecosystems in southern Africa, particularly in the South African Lowveld (low altitude) savannas [7,24,32]. This is concerning because localised fuelwood scarcity is already being experienced, and the situation is unlikely to improve in the future [33]. Indeed, localised fuelwood shortages have facilitated the development of fuelwood markets [34-35], effectively increasing the harvestable area and thus the impacts of fuelwood extraction may become less of a localised phenomenon. Despite fuelwood markets contributing to rural livelihoods [34-35], they have the unfortunate knock-on effect of artificially maintaining perceptions of fuelwood abundance [36]. Although a depletion of woodland biomass was predicted to occur in Bushbuckridge, South Africa, by 2011 [24] and more recently, by 2024, at current extraction rates [32], the interactions between socioeconomic and environmental factors driving natural resource use are complex, non-linear systems that are difficult to quantify [37]. However, the above predictions do raise the concern that woody vegetation harvesting, driven by increased demand and greater extraction amounts is unsustainable [38] and reduces the ability of ecosystems to provide ecosystem goods and services, fuelling the link between rural poverty and environmental impoverishment [39].

Wood harvesting changes not only biomass, but also vertical stratification of vegetation. Vertical vegetation complexity has relevance to ecosystem function as canopy height is related to biomass and productivity [40], biodiversity [41-43] and contributes to structural heterogeneity [44]. We submit that a method of understanding and, potentially, improving biomass change estimations, is to examine the vertical vegetation structure. We believe that by observing the interplay between woody biomass change and subcanopy structural change, drivers of biomass dynamics may be revealed.

Vertical subcanopy structure of vegetation canopies, however, cannot be derived from traditional two-dimensional remote sensing methods and top of canopy cover is a poor predictor of subcanopy cover [45]; three-dimensional (3-D) field-based efforts are impractical at landscape scales. Light detection and ranging (LiDAR) is a valuable tool for repeat estimation and monitoring of biomass, whilst providing subcanopy information, over large geographic areas and with fine-scale detail [46]. Repeat LiDAR campaigns have enabled tracking of woody biomass change as well as variation in the 3-D structure of the vegetation, providing the means to test previous fuelwood supply-demand model predictions [24,32], and to make inferences about the sustainability of wood provision under continued wood extraction pressure. The aim of this research is to utilize the power of airborne LiDAR to assess changes in aboveground biomass and subcanopy structure, as a unique window into unravelling vegetation dynamics in structurally heterogeneous landscapes.

Methods

Study Site

Permission to conduct fieldwork in the Bushbuckridge communal rangelands was granted by the local headmen. This study is part of a broad, long-standing relationship with the local community and the University of the Witwatersrand to conduct ecological research in their communal land. The field studies did not involve endangered or protected species. The study sites were located within the Bushbuckridge Municipality in the Lowveld region, a semi-arid savanna in South Africa. Summer rainfall (October to May) usually falls in convective thunderstorms and ranges between >900 mm per annum in the west and 500 mm per annum in the east with an mean annual precipitation (MAP) coefficient of variation of 25%. Summers are hot and humid with mean daily maxima of 30°C and winters are mild and dry with mean daily maxima of 23°C. Droughts can be prolonged and may be experienced every ten years. Within the timeframe of this study (2008–2012), the 2006–2007 and 2007–2008 summer rainfall was below average and the 2011–2012 was a particularly wet summer. Within seasons, notable rainfall peaks occurred in April 2010 (4.1-fold more rain than the monthly 8-year average) and January 2012 (2.4-fold higher than the monthly 8-year average).

The terrain is shallowly undulating and the geology is dominated by granite with local Timbavati gabbro intrusions. Classic catenal sequences are common in areas with shallow, sandy, dystrophic soils on the uplands and deeper, clayey, eutrophic soils on the bottom slopes [7]. The predominant vegetation type is granite lowveld, but the region also contains gabbro grassy bushveld and legogote sour bushveld [47]. Common plant species on the granite Lowveld uplands include: Terminalia sericea, Combretum zeyheri and C. apiculatum; the bottom slopes are characterised by Acacia nigrescens, Dichrostachys cinerea and Grewia bicolor [47]. Other frequently occurring species are Sclerocarya birrea, Lannea schweinfurthii, Ziziphus mucronata, Dalbergia melanoxylon, Peltophorum africanum and Pterocarpus rotundifolius. The majority of the woody biomass in the region is formed from S. birrea, Pterocarpus angolensis and A. nigrescens [7].

Bushbuckridge is surrounded by conservation land (both state-owned and private) [48] which increases the pressure for grazing and harvesting outside of protected areas. An overgrazing land-use legacy exists from intensively stocked, white-owned cattle farms from 1913 onwards [49]. Apartheid followed in 1948, with the Promotion of Bantu Self-Government Act of 1959, which forced black South Africans to live in ‘homelands’ [49]—centralised settlements on farms of 1000–2000 ha. Bushbuckridge Municipality was formed from the joining of Mhala in Gazankulu and Mpulaneng in Lebowa [2], with settlement boundaries defined by the old cadastral borders of the historical cattle ranches [50]. Although Bushbuckridge falls under state control, there is customary communal land tenure controlled by headmen who zone the land into residential, arable and communal areas for grazing of livestock and collection of timber and non-timber products (e.g. thatch, fruit, medicine) [51]. The settlements range from small, isolated villages to larger, dense settlements along major roads [33]. Human population density sharply increased between 1972 and 1994 to approximately 300 people/km2 [49] but these growth rates have declined over the past ten years [35]. Commensurate with human population growth in the area, the spatial footprint of the residential regions has expanded [37,52]. A foreboding of this decline was an observed reduction in the size-class distribution of the woodland vegetation with increasing distance from certain settlements [53].

Within Bushbuckridge, three communal rangelands were chosen to represent different levels of natural resource utilisation. These rangelands are zoned for use by the following villages: Justicia; Croquetlawn, Ireagh and Kildare; Xanthia and Agincourt (Fig 1). The rangelands were classified according to the relative wood extraction pressure assessed using 2008 data on the number of people and households accessing a given rangeland and relative to this corresponding rangeland area: high (9.2 people ha-1, 1.56 households ha-1; using 2155 ha of rangeland); intermediate (1.8 people ha-1, 0.35 households ha-1; using 1815 ha of rangeland); and low (0.21 people ha-1, 0.04 households ha-1; using 4425 ha of rangeland) (see [53] for detailed demographic data). Although each rangeland is used by its corresponding settlements, use is not exclusive to these villages and foreigners (both local and cross-border immigrants) are known to harvest from these areas [38]. The intermediate-use intensity rangeland (Justicia) is the only example of exclusive access, as it is fenced on two sides by private conservation land and its location makes it more difficult to access from other villages [32].

Fig 1

Study sites in Bushbuckridge municipality, located in the South African Lowveld.

Sites are classified (from west to east) as low, high and intermediate wood extraction pressure based on the number of households and people utilising each rangeland. Settlements that utilise each rangeland are shown, including the names of the major settlements, as well as the location of the gabbro intrusions in the predominantly granitic landscape.

Field-derived biomass estimates

All field data were collected concurrently with the airborne LiDAR campaigns in April 2012. Field-plots (total n = 56; high extraction site n = 16; intermediate extraction site n = 20; low extraction site n = 20) of 25 m x 25 m were established within the extent of the communal rangelands LiDAR coverage, and their locations recorded with a differential Global Positioning System (Trimble GeoXH Handheld GPS). All heights and basal stem diameters on stems thicker than 5 cm on trees taller than 1.5 m in height were recorded. A ‘tree’ may refer to a single-stemmed or multi-stemmed individual derived from the same rootstock, whilst ‘stem’ refers to the all branches derived from a single point on the ground. These height and basal stem diameter field data were used to estimate field biomass using allometric relationships from Colgan et al. [9], an extensive harvesting study with the same woody species composition as Bushbuckridge, in the form: where m is dry aboveground stem mass (kg), D is stem diameter (cm), H is height (m) and ρ is a unitless wood-specific gravity constant. The individual stem masses where then summed within each 25 m x 25 m plot to obtain plot-level field biomass, reported in Mg ha-1.

Light detection and ranging (LiDAR) data

The communal rangelands were surveyed with airborne laser mapping as part of a Carnegie Airborne Observatory (http://cao.ciw.edu/) campaign in April 2008 and April 2012, concurrently with the collected fieldwork data in 2012. Small footprint, discrete-return LiDAR is a remote sensing method which estimates 3-D vegetation structure over large areas. The 2008 LiDAR data were collected from 2 000 m a.s.l. with the CAO-Alpha system with a laser pulse repetition frequency of 50 kHz and laser spot spacing of 1.1 m (see [54]); the 2012 data were collected with CAO-2 AToMS with a laser pulse repetition of 100 kHz and laser spot spacing of 1m (see [55]). The LiDAR system also provides accurate geo-locational information generated by a high performance inertial management unit (IMU) and global positioning system (GPS) [54]. The LiDAR product is a 3-D point cloud from which a canopy height model (CHM) was constructed from the difference between the digital terrain model (DTM, interpolated from the last LiDAR returns) and the digital surface model (DSM, interpolated from the first LiDAR returns). Spatial errors on the more coarse of the two products (2008 data) were <0.20 m vertically and <0.36m horizontally [54]. Although different sensors and processing methods were used for the 2008 and 2012 data, errors between corresponding DTM’S were <15cm.

Volumetric pixels (voxels) are formed by aggregating LiDAR laser returns into 1 m height classes [56]. The position of each voxel is taken from the voxel centroid relative to the ground. LiDAR return frequency, within each voxel, are reported as a percentage relative to the total number of LiDAR points in the complete vertical column, including the ground returns. These data are used to quantify subcanopy (i.e. vegetation beneath the canopy cover) structure.

LiDAR-derived biomass estimates

LiDAR-derived metrics of woody vegetation can be used to estimate allometric relationships and infer biomass [2,8,9,32,57-58]. We derived a biomass regression model according to previously established methods by correlating the plot-level field-allometry and a corresponding LiDAR-derived H x CC (height x canopy cover) predictor metric calculated for each 25 m x 25 m grid cell created to correspond to the 25 m x 25 m field plots; H is plot-averaged (mean pixel height values >1.5 m) and CC is the proportion of canopy cover per plot (proportion of pixels >1.5 m in height). Both values were extracted from the CHM (see [9] for details). The H x CC metric is not only ecologically meaningful as it is an approximation of wood volume, but it also gives the best results over more complex metrics [2]. The height mask (>1.5 m) was used to account for the possibility of ground and tall grass being misclassified as vegetation. The LiDAR-derived predictor metrics were trained against field-derived biomass for each rangeland as they all exhibit different vegetation structural patterns, resulting from variable rainfall, different geologies and wood extraction pressures. Not only were these site-specific models able to explain more variation than one general equation; they were also deemed more ecologically valid. Biomass maps were then created by applying the site-specific biomass models to the LiDAR CHM extent (masked at heights >1.5 m) for each rangeland for both 2008 and 2012. Only grid cells that fit the criteria of an average height of >1.5 m (once pixels of <1.5 m were excluded) were used to estimate biomass as this is the vegetation that the fieldwork included. However, the cells that matched these criteria varied in both number and spatial location between 2008 and 2012. For the purposes of biomass change detection, only those cells that met the average height criteria for both years in the same location were considered. Riparian areas adjacent to streams in the rangelands were excluded from the biomass maps as they require separate calibration [2]. Similarly, cultivated fields and built-up areas were excluded.

LiDAR-derived subcanopy analysis

The voxel data (5 m x 5 m x 1 m) were resampled to 25 m x 25 m x 1 m, making the data comparable to the biomass grid cell sizes, and stacked into the following ecologically relevant, vertical height classes: 1–3 m (shrubs and small trees in the ‘fire trap’ [16]); 3–5 m (trees in the ‘elephant trap’ [22]); 5–10 m (tall trees contribute to structural diversity and thus to ecosystem function [59]); >10 m (very tall trees, ‘keystone structures’ [60], are often culturally important trees conserved in the rangelands [61]). These data were used to detect changes in the distribution of the vegetation size classes within the vertical vegetation column. “LiDAR returns” refers the percentage of laser pulses that were emitted from the sensor, hit an object and returned to the sensor. In the results, “Total % LiDAR returns” refers to the returns for the full vegetation column—excluding the ground returns. “% Subcanopy returns” refers to the LiDAR hits within a particular height category. Higher subcanopy returns implies greater density of vegetation in that height class.

Data extraction and analysis

Features of the settlements (e.g. roads, villages, crop fields) and rivers were manually digitised using a combination of SPOT 5 imagery (panchromatic-multispectral merge (480–890 nm), 2.5 m spatial resolution, www.spotimage.com) and aerial photographs (50 cm resolution, www.ngi.gov.za). Biomass estimates were extracted from the maximum number of randomly distributed points with a minimum enforced distance of 50 m to avoid spatial autocorrelation, based on the results of semivariograms (calculated in ENVI v4.7). All data were analysed in R v3.0 (R Core Team), including descriptive statistics, linear regression models and correlations. Biomass estimates were tested with Shapiro-Wilk Normality tests from the “fBasics” package and all sites in both 2008 and 2012 were found to be non-normally distributed (p < 0.001). Thus, a non-parametric Wilcoxon rank sum test was used to analyse differences between means over time within sites.

Results

Biomass models

A strong relationship existed between the field allometry and LiDAR metrics, although the highly heterogeneous rangeland resulted in high root mean square error (RMSE) values in both high and low use sites on granitic substrates (18.6 and 19.1 Mg ha-1, respectively) (Table 1). The increase in variability with increase in biomass indicated (S1 Fig) less agreement between the field allometry and LiDAR metrics at higher biomass values. This is a common phenomenon, termed ‘heteroskedasticity’, of model performance at higher biomass levels where the error variance is not consistent over all the observations [62]. Most typically, modelling the error structure shows a fanning pattern of increasing variance with increasing biomass [62], and this is true of the residual structure for both the high and low wood extraction sites (S1 Fig).

Table 1

Site-specific biomass models derived from field allometry and LiDAR metric linear regression.

Extraction pressure	Model	R2	n	RMSE (Mg ha-1)
high	y = 2312.3x - 157.14	0.78	16	18.6
intermediate	y = 409.57x + 252.74	0.60	20	4.8
low	y = 913.9x + 127.86	0.68	20	19.1

In the model equations, y refers to the plot-level (25 m x 25 m) biomass estimate (kg/625 m2) and x to the LiDAR-derived H x CC predictor metrics, where H is plot-averaged height (> 1.5 m) and CC is the proportion of canopy cover (> 1.5 m in height) per plot. Root mean square error (RMSE) was reported in Mg ha-1 for ease of interpretation and n is number of 25 m x 25 m plots.

Biomass dynamics

Mean biomass (± SD) in 2008 at the high, intermediate and low extraction sites was: 26.99 ± 16.43 Mg ha-1 (n = 102 cells), 9.42 ± 4.13 Mg ha-1 (n = 291 cells), and 21.18 ± 12.04 Mg ha-1 (n = 1654 cells), respectively. Biomass increased significantly at all sites between 2008 and 2012 by an average 18.38 Mg ha-1 (highest use site: W = 3036, p <0.001), 5.45 Mg ha-1 (intermediate use site: W = 16780, p <0.001), and 11.34 Mg ha-1 (low use site: W = 771641, p <0.001) (Table 2).

Table 2

Mean biomass increase (Mg ha-1) at sites under varying wood extraction pressures.

	Extraction pressure
	High (n = 102)	Intermediate (n = 291)	Low (n = 1654)
2008 (mean ± S.D.)	26.99 ± 16.43	9.42 ± 4.13	21.18 ± 12.04
2012 (mean ± S.D.)	45.37 ± 28.37	14.87 ± 6.76	32.52 ± 17.60
Absolute increase	+18.38	+5.45	+11.34
Relative increase (%)	+68.08	+57.80	+53.57

n is the number of 25 m x 25 m grid cells in each rangeland.

Variability increased with increased biomass, particularly in the high and low extraction pressure sites (Table 2). Represented as a rate of biomass change, the mean annual woody biomass productivity (± 95% spatial confidence interval) translates to 14 ± 1.39% p.a, 12 ± 0.08% p.a. and 11 ± 0.00% p.a for the high, intermediate and low wood extraction sites, respectively. These increases were despite ongoing wood harvesting in these rangelands. Relative to the starting biomass, all mean increases were greater than 50% (Table 2). Extreme biomass increases were related to large changes in relative height (Fig 2) and relative canopy cover (e.g. >50% increase in canopy cover results in biomass increases of >20 Mg ha-1, Fig 3). However, the extreme biomass changes (i,e. >40 Mg ha-1) predominantly occurred in the 1–3 m height class (Fig 2A and Fig 3A). Biomass increases of >40 Mg ha-1 did not occur in height classes >5 m (Fig 2C and Fig 3C). The largest increases in biomass occur in the high wood extraction site when compared with the same increases in relative height (Fig 2A and 2B) and canopy cover (Fig 3A and 3B) in the other rangelands. There are no data for the high extraction site for the 5-10m height class as there are no grid cells with an average height >5 m in this rangeland (Fig 2C and Fig 3C).

Fig 2

Height-specific biomass change as a function of relative height change per grid cell.

Height categories are a) 1–3 m, b) 3–5 m and c) 5–10 m for rangelands of high, intermediate and low wood extraction pressure. There were no data for the 5–10 m height class in the high wood extraction rangeland and the >10 m height class for all rangelands as there were no grid cells with an average height over 10 m. Grid cell size: 25 m x 25 m.

Fig 3

Height-specific biomass change as a function of relative change in canopy cover per grid cell.

Height categories are a) 1–3 m, b) 3–5 m and c) 5–10 m for rangelands of high, intermediate and low wood extraction pressure. There were no data for the 5–10 m height class in the high wood extraction rangeland and the >10 m height class for all rangelands as there were no grid cells with an average height over 10 m. Grid cell size: 25 m x 25 m.

Vegetation structural dynamics

Total % canopy returns increased between 2008 and 2012 in all rangelands, but up to 79% of the total change in canopy returns was attributable to the increase in the 1–3 m height category within the subcanopy (Fig 4). Losses in subcanopy returns were only found in the high wood extraction rangeland, and only in the 5–10 m height class (Fig 4A). There was little contribution to total change in % subcanopy returns from the >10 m height class (Fig 4). Although the high and low extraction rangelands had fairly similar overall increases in % total canopy returns, this was not the case with relative change (from 2008), where the highest extraction site was far greater (e.g. relative canopy returns for height class 1–3 m: 425%, 387% and 90% for high, intermediate and low extraction, respectively). Thus, the order of relative change in % canopy returns followed the gradient of wood extraction levels at the different sites.

Fig 4

Height-specific subcanopy returns (%) (mean ± standard deviation) for 2008 and 2012.

Wood extraction levels are: a) high (n = 102 cells), b) intermediate (n = 291 cells), and c) low wood extraction (n = 1654 cells). Contribution of height class change (subcanopy returns) to total change (total vegetation column) (%) is the black bar represented by values on the secondary axis. e.g. In the high wood extraction rangeland, 79% of the change in the total vegetation column was attributable to the 1–3 m height class.

Another indicator of shrub level increase in the rangelands is the change in the number of cells that remained after an average height mask was applied (i.e. that fulfilled the average height criteria threshold to be included in the biomass analysis), expressed as a percentage of each rangeland. The high extraction rangeland changed from 10% of the rangeland that met the average height (>1.5 m) criteria mask in 2008 to 15.9% of the rangeland in 2012 (χ2 1 = 107.6; p <0.001); the intermediate use site doubled in the percentage of rangeland that met the average height criteria from 8.5% to 17.4% (χ2 1 = 780.8; p <0.001); and the low use rangeland increased from 54.2% in 2008 to 63.8% of the rangeland in 2012 (χ2 1 = 220.7; p <0.001).

Association between biomass change and vegetation subcanopy returns

There was a positive correlation between change in biomass and change in % subcanopy returns (Fig 5); particularly in the 1–3 m height class in the high extraction sites (high extraction: r = 0.22, p <0.0001; intermediate extraction: r = 0.58, p<0.0001) and the 3–5 m height class (high extraction: r = 0.62, p <0.0001; intermediate extraction: r = 0.64, p <0.0001; low extraction: r = 0.56, p <0.0001). Although this relationship was also present in the 5–10 m height class at all extraction levels (r >0.31), it degraded at heights >10 m (r < 0.10) (Fig 5). It is interesting to note that the strength of the relationship between change in biomass and change in % subcanopy returns across all height categories was strongest at the intermediate wood extraction site (Fig 5).

Fig 5

Height-specific correlation (p < 0.001) between change in biomass (%) and subcanopy returns (%).

Wood extraction levels for each rangeland are listed per column as high, intermediate and low.

Changes in biomass and height-specific subcanopy returns were spatially associated (S2 Fig). However, these changes were more apparent at <5 m (S2 Fig). Almost no change in % subcanopy return for vegetation >10 m is evident (S2 Fig). The same biomass values for a given grid cell can manifest as different structural profiles. As such, structural profiles could change in different ways whilst maintaining the same overall biomass value outcome. For example, if the site was dominated by grasses with several trees >5 m, that site could, theoretically, show no change in biomass value by 2012, but the structural profile may have changed to predominant shrub cover and fewer tall trees.

Discussion

Large increases in biomass at all sites (Table 2) are in contradiction to previous fuelwood supply-demand models which predicted biomass depletion [24,32,63]. Biomass increases in Bushbuckridge rangelands were attributable (>80%) to vegetation in the 1–3 m height class within the subcanopy (Fig 4), with extreme biomass gains (>20 Mg ha-1) associated with vegetation that gained >25% in height (Fig 2A) or >50% in canopy cover (Fig 3A). This agrees with an observed increase in the number of thinner, taller stems within Bushbuckridge rangelands [35] and more grid cells meeting the average height criteria in each of the rangelands between 2008 and 2012. These low height class increases probably reflect local-scale dynamics of harvesting—more harvesting drives coppicing (resprouting from the stem or roots) in the intermediate and high extraction sites (Fig 2A and Fig 3A)—but the relationship appears more pronounced in the intermediate site as less of the coppice is harvested. It is likely that wood harvesting is acting as a ‘bush thinning’ mechanism, changing the size specific growth rates, particularly in resprouting from stumps with fully-developed root systems [64]. Indeed, thick stands of small-stemmed trees can yield more woody biomass than a few, large trees as a result of divergent, size-specific growth rates [65]. However, low height class increases in biomass could also be a result of newly established bush encroachers which characteristically invade overgrazed and degraded rangelands [66-68]. Biomass estimations for different height classes in a savanna woodland reveal, collectively, greater biomass quantities are located below 4.5 m in height than above; a disparity more prominent immediately after a disturbance [69]. Harvesting has been found to increase the density of smaller stems without changing the height structure of the woodland [70]. Unfortunately, there is a dearth of data on the preferred height of harvested species, only preferred diameter size which ranges, location dependent, between 2–6.5 cm [26,36]. There are records of stems >1 cm being taken, with preference for those >4 cm and almost no stems harvested >20 cm [71]. Extrapolating 1 cm and 5 cm diameter size into available coppice diameter-height allometry relationships [72] suggests pre-harvested heights of 0.74 m and 2.92 m in Dichrostachys cinerea, 0.63 m and 2.07 m in Acacia harveyi, and 0.77 m and 2.44 m for Combretum collinum, respectively. Although the relationship between harvested stem diameters and regrowth shoot length is variable, we can infer that stems harvested for fuelwood are generally <3 m. Therefore, preferred ‘harvesting heights’ coincide with height class with the most subcanopy gains (Fig 4).

Subcanopy biomass increases at low heights in a rangeland context are likely a combination of woody regrowth-response (harvesting effects) [71-74] and bush encroachment (overgrazing effects) [15,75-76], here collectively referred to as ‘bush thickening’. However, these are not mutually exclusive events and can occur together. Low height-class increases occur in Bushbuckridge both as standalone shrubs as well as occurring underneath the canopies of tall trees [45]. Resprouting rates and the subsequent influence on communal rangeland dynamics have been underestimated in the earlier research in this region [77]. Although the Wessels et al. [32] supply-demand model did include resprouting estimates of 89 kg ha-1 yr-1 which is significantly higher than the 20 kg ha-1 yr-1 that the Banks et al. [24] model used; these rates are only from one species, T. sericea, and thus, may underestimate the growth rates for the other predominant coppicing species, e.g. D. cinerea. Previous data suggest that even during a poor rainfall period, in just five months there was coppice of 989 kg ha-1 (6.6% of the total post-harvest biomass) and harvested trees recovered two thirds of their preharvest biomass, with no harvest-induced mortality [71]. T. sericea coppice shoots from established stumps gained between 1–2 m in height over 3 years [78], whilst coppice stands in Malawi and Kenya gained 3m [79] and 2m [80], respectively, over 4 years. This is evident in the annual productivity suggested by the LiDAR-derived estimates of well over 10% p.a. (especially when we consider that this is over and above the biomass removed for wood energy) which exceeds the previous woodland productivity value of 4% [24,32,81]. The disparity in the growth rates is likely a result of higher productivity in the low height classes [69] and a significant contributor to the Wessels et al. [32] underestimation of biomass production rate. Growth rates could also have been affected by the drier than normal conditions in 2007 and, likewise, the high rainfall in 2010 and early 2012. As data collection was subsequent to these events, it is likely that biomass estimates were affected.

Although lower height classes within the subcanopy showed increases across all wood extraction sites (Fig 4), this was not true for subcanopy returns in the 5–10 m class in the high wood extraction site (Fig 4A). Large, fruiting trees are normally conserved by villagers as they are used for a variety of non-timber uses [82-83]. Despite cultural practices against live-wood harvesting of large fruiting trees, villagers acknowledge that they do cut trees, like marula (Sclerocarya birrea), as they feel they have no alternatives in the face of high electricity prices and localised shortages of fuelwood [83]. We observed several felled and pollarded marula trees in the highest wood extraction site and can assume, together with the lack of data for grid cells of average height >5m (Fig 2C and Fig 3C), that the loss of vegetation returns in the 5–10 m height class reflects a localised lack of fuelwood of sufficient quality and quantity in this rangeland. The reduced number of tall trees and abundance of short subcanopy vegetation in the high use rangeland results in a more homogeneous stand structure (Fig 4A), a possible explanation for the stronger relationship between field and LiDAR data in this site (Table 1). Most fuelwood supply-demand models that predicted loss of biomass are not spatially explicit and did not capture the fine scale variation at village level [84-85] or the mismatch between the spatial variability in fuelwood supply relative to centres of demand [35], especially considering vehicles are increasingly being used to transport larger amounts of wood from more distant locations [39,86]. Yet, the Wessels et al. [32] fuelwood model focused on one “best-case scenario” communal rangeland, exclusively utilised by one village and still predicted losses. However, fuelwood demand is not a linear system and people’s responses to changes in their socio-economic and natural resource environment are complex and difficult to quantify [37]; consequently, the community’s adaptive responses are not incorporated in these models. Global and national studies highlight the lack of adaptive capacity of people in the developing world [37,87-88]; however, the strategies people adopt on local and regional scales often reveal surprising resourcefulness in response to change [89-91]. Within the fuelwood context in Bushbuckridge and elsewhere in Africa, responses to localised fuelwood shortages have included: changes in the preferred size class of fuelwood [29,35,86]; switching preferred fuelwood species [25,33,91]; more frequent trips or more time spent per trip to collect fuelwood [31,92]; travelling further from home [37]; use of wheelbarrows and vehicles to collect more wood per trip [33,38,86,93]; development of fuelwood markets [33,36]; and collecting from neighbouring private land [35]. Socio-economic factors also play a role in fuelwood demand dynamics. High dependence on government social grants and migrant worker remittances is characteristic of rural areas [33,94-95]; changes in these economic flows will affect household cash flow and, thus, alter household-level demand for natural resources. These adaptive strategies and socio-economic factors are difficult to capture in a supply-demand models and are a contributing cause to the disparity between predicted and measured biomass in communal rangelands.

Biomass values range between 9 Mg ha-1 (on gabbro) to 27 Mg ha-1 (on granite) which is comparable to the range for field-based allometry studies in the greater Bushbuckridge area (18.9–23.1 Mg ha-1) [7], and the LiDAR-estimates for the conserved Lowveld region (11.9–92.3 Mg ha-1) [2]. The intermediate wood extraction site has had previous estimates of LiDAR-derived biomass for 2008 of 12 Mg ha-1 [32], but this used allometry from Nickless et al. [96] and field–LiDAR biomass regression relationships derived from the regional landscape. Most studies on allometry have focused on temperate zone and deciduous forests (e.g. [58,97-98]) or tropical forest monitoring (e.g. [8,99-101]). Very few have focused on savanna systems (e.g. [2,4,96]). Both Chave et al. [8] and Colgan et al. [9] stress the importance of allometric equation choice on error as even field-allometry had 16% RSE (Residual Standard Error); these errors often compound with averaging. Although Colgan et al.’s [9] plot-averaged LiDAR-derived biomass estimates had 9% more relative error (difference between predicted and measured biomass) than field-harvested biomass, the bias (mean error) was only -3% (compare to Nickless et al. [96] allometry with 15% more relative error and 50% bias) [9]. Our study also excluded all cells that were below 1.5 m in average height in both 2008 and 2012, cutting out a large proportion of the area relative to the portion used in Wessels et al.’s [32] study. Although our biomass model has fewer field-calibration sites than the Wessels et al. [32] study, our calibration sites were specific only to the area the biomass models were applied to.

While we are confident that our biomass estimates are reflecting a true increase, the shortcomings of using this method have the potential to exaggerate increases, particularly error in canopy cover measurements over time. This is of concern when considering leaf area index (LAI) in LiDAR change detection metrics, as both the voxel and the CHM data may be influenced, affecting the biomass estimates as well as the subcanopy LiDAR returns. Although this was controlled for as much as possible by collecting the LiDAR data in the same month each campaign, LAI varies with phenology and with local climatic changes, such as differential rainfall between years, or heavy winds [102]. The relatively high predictive uncertainty (RMSE range: 4.8–19.1 Mg ha-1) in the biomass models occur in the high and low wood extraction rangelands, both of which are situated on granitic geology (Fig 1) which are more heterogeneous in both topographic relief and stand structure, as well as in the resultant biomass (Table 1). In landscape-scale approximations of biomass, errors are introduced and often propagated. The assumption is that individual plant measurement errors will average out over the plot level, provided the plots are large enough and the measurement process is unbiased. There is also an effect of plot size on error; increasing plot size increases the predictive power of the model [10]. However, there is a trade-off between the cost and logistic realities of sampling large plots and the need to sample a large number of plots, as plot number also affects landscape-scale error [9]. Although relative uncertainty in the biomass models was high and may have been reduced by object based image analysis (OBIA) methods applied to single tree crowns to counter vertical structural heterogeneity errors, plot-level averaging methods have a positive trade-off in their simplicity and their ability to average out within-plot variation, particularly the horizontal canopy cover heterogeneity characteristic of savannas.

Conclusions

Savanna-based biomass studies have considerable scope to rectifying the underestimation of carbon sinks and sources, elucidating the woody encroachment problem in savannas and untangling the interactions between bush encroachment/thickening and wood extraction by rural communities. Without high resolution, 3-D vegetation data covering a large area, the landscape-scale increases in biomass over the Bushbuckridge rangelands could erroneously be interpreted as woodlands recovering to an “unaltered” state. Users of two-dimensional, remotely-sensed biomass estimates should remain aware of structural implications in the landscape to make informed conclusions on vegetation dynamics, particularly in the context of increasing savanna bush encroachment in a CO2 rich environment [103-104]. Indeed, it is the low height class vegetation within the subcanopy which determines future woodland structure. Moreover, most carbon cycle studies in Africa neglect domestic emissions from wood harvesting [105] despite knowing the contribution of deforestation and land degradation to carbon dynamics [106]; a recent carbon model has demonstrated the importance of vegetation increases in the southern hemisphere’s semi-arid regions to terrestrial carbon sinks [107]. The repercussions of bush thickening in communal rangelands will have implications for the direct-use values of ecosystem goods and will affect household vulnerability to shocks [39]. Our research suggests that wood harvesting can, paradoxically, exacerbate bush thickening as many of the harvested savanna species have strong regenerative responses [71-72,79-80,108-109]. Not only is coppice an important survival strategy for regenerating woodlands, the resprouted stems may provide a valuable source of future harvestable biomass [74,78,110-112]. There is, however, little information on regrowth rates and response to continued harvesting as well as whether the coppice is of appropriate quality for fuelwood.

Biomass model data.

Data include 2012 LiDAR-derived average height and canopy cover extraction metrics, as well as field-work based allometry. Each line item is per 25 m x 25 m grid cell. Metadata are included in the dataset.

(XLSX)

Click here for additional data file.

Biomass and subcanopy data.

Data include 2008 and 2012 biomass estimates derived from biomass models as well as % subcanopy returns for voxel data for the height class categories: 1-3m, 3-5m, 5-10m and >10m. Each line item is per 25 m x 25 m grid cell. Data are organized per land extraction category into separate worksheets. Metadata are included in the dataset.

(XLSX)

Click here for additional data file.

Biomass changes (Mg ha-1) in relation to relative height and canopy cover change.

Data include biomass change estimates (2008–2012), percentage height and canopy cover changes for each 25 m x 25 m grid cell. Each height class (relative to height in 2008) are shown on separate worksheets. Metadata are included in the dataset.

(XLSX)

Click here for additional data file.

Site-specific biomass model residuals.

The residual spread demonstrates heteroskedasticity with increasing biomass fitted values for rangelands with a) high, b) intermediate and c) low extraction pressure.

(TIFF)

Click here for additional data file.

Biomass changes (%) relative to height-specific change in subcanopy returns (%).

Height categories are: 1–3 m, 3–5 m, 5–10 m and >10 m.

(TIF)

Click here for additional data file.