Rates of woody encroachment in African savannas reflect water constraints and fire disturbance

The aims of this study were to (1) estimate current rates of woody encroachment across African savannas; (2) identify relationships between change in woody cover and potential drivers, including water constraints, fire frequency and livestock density. The found relationships led us to pursue a third goal: (3) use temporal dynamics in woody cover to estimate potential woody cover.


| INTRODUCTION
Change in woody cover in sub-Saharan Africa is characterized by decreases in populated areas and a slow but steady increase in remote areas, generally referred to as woody (or shrub) encroachment (Brandt, Rasmussen, et al., 2017). Woody encroachment in grasslands and savannas is a global phenomenon and its causes and consequences have been frequently debated over the last decades (D'Odorico, Okin, & Bestelmeyer, 2012). Consequences of encroachment have long been described as malign, being detrimental to livestock productivity by reducing grass cover and incurring expenses for shrub removal (van Auken, 2009). On the other hand, increasing carbon stocks in drylands can offset carbon losses of deforestation, and enhance the provision of fuelwood and timber resources (Birhane, Treydte, Eshete, Solomon, & Hailemariam, 2017). Higher woody canopy cover also affects biodiversity and local hydrological processes, including evapotranspiration, runoff generation and rainfall interception (Archer et al., 2017;Honda & Durigan, 2016). Since it impacts critical ecological processes and peoples' livelihoods, we need a thorough understanding of the extent and rate at which encroachment occurs, and improved understanding of why it occurs, to inform rational management and adaptation strategies.
Several potential causal factors have been linked to encroachment, including climate change (Brandt, Rasmussen, et al., 2017), changes in precipitation intensity (Kulmatiski & Beard, 2013), livestock grazing (Roques, O'Connor, & Watkinson, 2001), fire suppression (Dalle, Maass, & Isselstein, 2006) and increasing atmospheric CO₂ concentrations (Bond & Midgley, 2012;Buitenwerf, Bond, Stevens, & Trollope, 2012). In Africa, the disappearance of free-roaming elephants may also have contributed to the overall trend (Daskin, Stalmans, & Pringle, 2016;Stevens, Erasmus, Archibald, & Bond, 2016). Given the ubiquity of encroachment observations in different parts of the world, the two potential main drivers are increasing CO₂ concentrations and the expansion of livestock grazing. Elevated CO₂ levels can increase the efficiency of photosynthesis and plant wateruse, leading to enhanced growth among many woody species (Ainsworth & Long, 2005;Kgope, Bond, & Midgley, 2010;Morgan, Milchunas, LeCain, West, & Mosier, 2007). The expansion of livestock grazing has been associated with reductions to fire frequency and intensity, and often the removal of wild browsers (Archer et al., 2017). Browsing and fire both suppress woody seedlings and saplings, and are necessary components of higher rainfall savannas where they can maintain a mixed tree-grass state under rainfall amounts well above levels sufficient to support a closed canopy (Bond, 2008;Sankaran, Augustine, & Ratnam, 2013).
Irrespective of the dominant driver, it is clear that encroachment interacts with rainfall which defines the resource-based potential woody cover of tropical savannas. Sankaran et al. (2005) found a linear relationship between potential woody cover and mean annual precipitation (MAP) up to c. 650 mm/year in African savannas, beyond which canopy closure was possible (albeit not always realized due to natural and anthropogenic disturbances). An extension of this relationship can be found in tropical and temperate forests where the upper bound of above-ground biomass is constrained by water deficits (Stegen et al., 2011). It is, however, not clear how rainfall influences tree dynamics other than constraining the upper limit, or what role rainfall may play in regulating rates of woody encroachment.
Most analyses of encroachment focus on specific areas and may be biased in favour of sites where encroachment has been observed.
The most comprehensive analyses of encroachment in African savannas to date include Brandt, Rasmussen, et al. (2017) which used continent-wide passive microwave remote sensing, and Stevens, Lehmann, Murphy, and Durigan (2017) which compiled a metaanalysis of past woody cover change studies. Both these studies point to widespread encroachment in African drylands, with the latter study finding a mean increase of 0.25% woody cover per year.
Our study uses very high resolution (VHR) satellite imagery, sampled at 596 sites in 38 areas across drought-seasonal Africa, to (1) provide unbiased estimates of current trends in woody cover change in the savannas, and (2) explore how rates of change relate to environmental (climate, soil, herbivory and fire) variables and estimates of potential woody cover based on local rainfall amounts.
2 | DATA AND ME TH ODOLOGY 2.1 | Acquiring and sampling satellite imagery We used the same dataset, image preprocessing steps, and classification framework as in a recent study on patterns in African woody vegetation structure (Axelsson & Hanan, 2017) and refer the reader to this paper for a more thorough description of the methodology.
Our sampling frame was sub-Saharan African savannas with a minimum of anthropogenic disturbances. This study used a subset of sites from the previous study where we had overlapping older (2002-2006) and newer (2011-2016) imagery. In all cases, we only acquired imagery from seasons when trees were in full leaf in the particular area. Most of the newer imagery were WorldView 2 while the older images were from the Quickbird 2 satellite. Due to limited image coverage, it is currently not possible to do a continental change analysis with VHR data from a single satellite. In total, we acquired overlapping imagery in 38 areas with a mean time difference of 9.9 years between image pairs. Within the overlapping imagery, we sampled 240 9 240 m sites guided by a 0.04°longitude/ latitude grid which served as a base for site locations. The location of sites was, however, often adjusted to avoid locations with clear anthropogenic imprints such as settlements, roads, and agriculture, geomorphological features such as gullies and rock outcrops, or because of clouds and cloud shadows in the imagery. Where we needed to adjust site location, we selected the closest suitable location to the original point to reduce the potential for observer-bias in selecting the new site. Later, during the classification phase, some sites were eliminated due to difficulties in identifying woody plants because of a lack of green leaves or similar spectral signatures of trees and grasses. We ended up with 596 sites, and due to varying image quality, anthropogenic disturbances and area of image overlap, there was variation in the number of sites per area ( Figure 1).

| Preprocessing and classification of satellite data
Preprocessing steps included Gram-Schmidt pan-sharpening of the blue, green, red and near-infrared bands, and orthorectification using embedded RPC-information and an SRTM 2 DEM (Farr et al., 2007).
The orthorectified imagery were resampled to 0.6 m ground resolution using the nearest neighbour method. Centred at the site coordinates, we then extracted a 400 9 400 pixel (240 9 240 m) area to represent each site. For classification of the woody pixels, we applied an unsupervised (ISODATA) classification approach dividing each image into 18 spectrally different classes. The results were smoothed using a kernel size of 3 pixels. Up to this point, processing was performed in ENVI 5.2 but the final classification step was done using a custom-built software. This step involved manually assigning the unsupervised classes into woody and nonwoody cover types.
The custom-built software helped us streamline the final classification step, and perform it in an efficient and consistent manner.
Examples of two sites with estimates of woody cover for the older and newer imagery are shown in Figure 2.

| Data used in the regression analysis
After image classification, we had estimates of woody cover from the older imagery (initial cover), the newer imagery (final cover), as well as acquisition dates for the imagery. Annual woody cover change was calculated as (final coverÀinitial cover)/time in years.
We investigated relationships between woody cover change and several datasets that could have a causative impact on tree dynamics: 2. MAP change identifies change in rainfall relative to the years prior to the first image acquisition. It was calculated as the difference in MAP between the observation period (described above) and the period from January 1998 to December of the year prior to the observation period.
3. (PotentialÀInitial) cover is the difference between rainfall-based potential woody cover and the initial cover, and measures how much woody cover can expand before it is limited by water deficits. Potential cover was calculated as 0.14 9 MAPÀ14.2 as described in Sankaran et al. (2005), but 30km F I G U R E 1 Location of the 38 study areas, containing 596 study sites, spread out over African rangelands (defined using the Anthropogenic biomes product; Ellis & Ramankutty, 2008). The map to the right shows a study area on the border between Somalia and Ethiopia to illustrate the sampling strategy for study sites (white rings). The placement of sites was guided by a 0.04°longitude/latitude grid (green lines) in areas with overlapping older and newer satellite imagery (blue lines) with local adjustments to avoid anthropogenic influence, geomorphological features and clouds To avoid registering fires identified in adjacent months as separate fires, we counted fire events in consecutive months as a single fire.

| Regression analysis and considerations for spatial autocorrelation
We started by using ordinary least squares (OLS) regression (lm function in R) to identify significant correlative relationships with woody cover change. Due to the clustered spatial pattern of our sites, there was a clear risk of spatial autocorrelation that should be investigated and addressed. Spatial autocorrelation was evaluated by first constructing a neighbourhood matrix where neighbour relations were set for sites closer than 90 km apart, then adding row standardized spatial weights, and finally calculating Moran's I (Moran, 1950) based on the spatial weights matrix and residuals from the OLS regression model. The distance 90 km was defined where spatial autocorrelation for the y-variable (woody cover change) reached its maximum. Calculations were made using the spdep package in R (Bivand et al., 2017). We found significant positive autocorrelation among residuals which leads to a higher risk of falsely identified relationships (Quesada et al., 2012

| RESULTS
The mean annual change in woody cover for all 596 sites was 0.27% per year, with 394 positive, 198 negative, and four with no change.
The mean woody cover change was significantly greater than zero (p < 2.2e-16). When we first averaged cover change for each area and then took the mean of the 38 areas, we got 0.25% per year, with 30 positive and eight negative. The latter approach is likely more reliable because it removes much of the spatial autocorrelation among the individual sites. A density plot of annual woody cover change ( Figure 3) shows that the distribution peaks close to zero.
There was a higher concentration of woody cover losses in East Africa (Figure 4), likely due to the 2011 drought (Lyon & DeWitt, 2012) which impacted the whole region, and in particular southern Somalia, southern Ethiopia and northern Kenya. Our estimates showed woody cover declines in the areas most heavily affected by the East African drought and a lack of large (>0.5% per year) increases in the overall East African region.

| Drivers of change in woody cover
Standard OLS regression between woody cover change and potential drivers resulted in an R 2 of .15 (Table 1), signalling that most of the variation in woody cover change remained unexplained. Initial woody cover and MAP were not significant individually, but the difference between potential cover and initial cover was highly significant (p < .001). As potential woody cover is defined based on mean annual rainfall (per Sankaran et al., 2005), this implies that rates of increase in woody cover are positively correlated with available water resources. The other significant variables were change in MAP (p < .01), which had a positive influence, and fire frequency (p < .01) and cattle density (p < .05), both of which had a negative effect on woody cover change. The positive sensitivity of woody cover to rainfall is consistent with the importance of MAP in determining potential cover. Similarly, the negative sensitivity to fire frequency is expected given the role of fire in suppressing sapling survival (Hanan, Sea, Dangelmayr, & Govender, 2008). However, the negative effect of cattle density is inconsistent with expectations that cattle both reduce competition with herbaceous vegetation and reduce frequency of ground fires (Asner, Elmore, Olander, Martin, & Harris, 2004). The livestock densities were derived via disaggregation of census summary data for administrative areas (Wint & Robinson, 2007), perhaps limiting their utility in this analysis.
The residuals of the OLS regression were spatially autocorrelated (Moran's I: 0.23, p < 2.2e-16) which could compromise model reliability. We therefore complemented the OLS regression with results from a SAR model, which explicitly accounts for spatial autocorrelation by including a spatial structure among the predictors (Table 2).
In this case, (potentialÀinitial) woody cover and fire frequency remained significant, while change in MAP and cattle lose their significance. It is difficult to determine which of the two models is more correct because inclusion of a spatial structure in the regression model can decrease the role of spatially autocorrelated variables whether these are causative or not (Quesada et al., 2012).
We proceeded to further explore the relationship between woody cover change, rainfall, initial and potential cover ( Figure 5).
Here, a clear trend emerged with higher rates of encroachment where the initial cover is low and higher rates in high rainfall systems. As woody cover approaches its climatic maximum, the mean encroachment rate slows and eventually turns negative. Sites that have over-shot the potential cover, i.e. negative (potentialÀinitial) cover, are more likely to suffer mortality. 3.2 | Estimating potential woody cover based on woody cover change trajectories The relationship between woody cover change and initial cover ( Figure 5a) opens up the possibility to use it to derive estimates for potential woody cover, similar to the relationship presented by Sankaran et al. (2005). In Figure 5a, woody cover change for the middle rainfall category (400-700 mm/year) decreases with increasing initial cover and turns negative at c. 40% cover. We infer that 40% cover is a point of attraction for woody cover change that is equivalent to the potential woody cover for this rainfall interval.
By identifying the point of attraction for MAP intervals of 100 mm/year up to 800 mm/year, we derived a linear relationship for potential woody cover ( Figure 6). The remaining sites with MAP above 800 mm/year were placed in a single category. The derived linear relationship, Woody cover = 0.1415 9 MAPÀ23.43, has a similar slope to that presented by Sankaran et al. (2005; Woody cover = 0.14 9 MAPÀ14.2). The category with rainfall above 800 mm/year had a point of attraction above 100% cover, indicating that the average growth of woody cover is dependent on water availability also in mesic savannas, and only slows due to space restriction as it nears canopy closure. We noted that a num- Significance codes: ***p < .001, **p < .01. Nagelkerke pseudo-R 2 : .262, p-value: <2.22e-16.
and variation in environmental factors (Brandt, Tappan, et al., 2017). Potential cover − Initial cover Annual woody cover change (%) F I G U R E 5 Relationship between woody cover change, initial woody cover and rainfall. (a) Boxplot with woody cover change across intervals of initial woody cover and mean annual precipitation. Each sample represented the average woody cover change among sites in unique combinations of area, rainfall category and initial cover category. The red points denote the mean change for each box. Initial woody cover (%) Annual woody cover change (%)

(b)
F I G U R E 6 Inference of resource-limited potential woody cover from woody cover growth rates estimated across gradients of initial cover and rainfall. (a) The derived potential woody cover (dashed red line; Woody cover = 0.1415 9 MAPÀ23.43) calculated using linear regression of the points of attraction in each rainfall interval (red points). MAP for the points of attraction was calculated as the average MAP of the sites in each interval. We highlighted sites from Somalia and Ethiopia (yellow and blue dots) as these are overrepresented among sites with actual cover>potential cover. (b) Illustration of how the point of attraction was derived from the intersection between the regression line (blue) and the zero change line (black). The graph is for the rainfall interval 500-600 mm/year where the intersection lands at a woody cover of 51% identical to the estimate in the meta-analysis by Stevens et al. (2017) even though the data used in that analysis covered other time periods and were more concentrated to the southern part of the continent. This reinforces the conviction of a consistent trend across African savannas. While our estimates were lower in some parts of East Africa, this was likely a result of the 2011 drought occurring within our relatively brief window of observation. In both the OLS and SAR regression analyses, water availability (potential coverÀinitial cover) and fire frequency were significant variables (Tables 1 and 2). Among the investigated variables, these are the ones that can be most clearly linked to current encroachment rates.
As woody cover approaches its climatic maximum, water constraints slow growth and woody communities are more likely to suffer mortality through density-dependent regulation (Browning, Archer, Asner, McClaran, & Wessman, 2008). Because of high competition from neighbouring plants, individuals are more sensitive to droughts and interannual variation in rainfall. The role of fire in suppressing woody encroachment is well established. Frequent fires kill or topkill smaller woody plants, thereby reducing the number of plants that reach a mature fire-tolerant size (Bond, 2008). The effectiveness by which fires can suppress woody cover also depends on woody growth rates and thus water constraints. Higher growth rates lead to a higher likelihood for juvenile trees to reach a fire-tolerant size and escape the fire-trap (Wakeling, Staver, & Bond, 2011). As our analysis avoided areas with high anthropogenic impacts, we found no relationship with population density (Tables 1 and 2). However, according to Brandt, Rasmussen, et al. (2017), there are substantial woody cover decreases in areas with high population growth due to deforestation and land clearing.
Change in precipitation and cattle density proved significant in the OLS regression, but these relationships could have been caused by spatial autocorrelation ( Table 2). The low accuracy of livestock data at continental scales also makes it difficult to accurately estimate their influence in this type of analysis. While this analysis shows what factors control encroachment rates among individual sites, an explanation for the general trend of encroachment is still elusive. Woody encroachment is widespread across African savannas and neither strongly dependent on climatic zone nor soil texture (Tables 1 and 2). It is very possible that the hidden factor causing an overall increase in woody cover is the higher atmospheric CO₂ level, which has been proposed in several previous studies (Donohue, Roderick, McVicar, & Farquhar, 2013;Higgins & Scheiter, 2012). If CO₂ acts as a global fertilizer, this could also increase potential woody cover (Stevens et al., 2016).

| Derivation of potential woody cover based on change trajectories
The strong relationship between initial woody cover, MAP, and rates of woody cover change ( Figure 5) shows that rainfall universally influences tree dynamics and growth, and not only by constraining the potential woody cover. Mean rates of woody cover change are positively correlated to the difference between current and potential woody cover, suggesting that competition among trees is a key variable determining growth rates. A large gap between actual and potential cover leads to a surplus of water that boosts growth rates.
While mesic savannas (MAP > c. 750 mm/year) are not considered water-limited in the sense that canopy closure is possible, their woody growth rates are nevertheless dependent on the difference between potential and initial woody cover. Conversely, a site with woody cover above its potential (possibly due to recent above-average rainfall) is more likely to experience water shortage and tree mortality. As this trend was consistent across intervals of MAP and initial woody cover, we were able to use it to derive a linear relationship for potential woody cover ( Figure 6). While Sankaran et al. (2005) used a 99th quantile piecewise linear regression to derive potential cover, this method uses the woody cover change trajectories. The relationship we found had a similar slope to that in Sankaran et al. (2005) but was shifted to the right, possibly due to the influence of fire and browsing which can suppress woody cover growth rates. This approach also recognizes that episodic recruitment and rainfall variability can cause woody populations to temporarily over-shoot the potential cover based on long-term MAP, as opposed to the more rigid statistical method (the 99th quantile) used by Sankaran et al. (2005). Many sites in the dataset fell outside of both our and Sankaran's line, the majority of which came from areas of Somalia and Ethiopia (Figure 6a). The divergent pattern for these sites could be related to the bimodal rainfall distribution in the Horn of Africa. Plants are not only affected by the total rainfall amount but also its intra-annual distribution (Axelsson & Hanan, 2017;Good & Caylor, 2011), which can determine levels of infiltration and runoff as well as windows of opportunity for seedling growth.
As opposed to the potential cover estimated in Sankaran et al. (2005), we extended the relationship up to 100% woody cover. The cap at c. 80% cover in Sankaran et al. (2005) emerged as a result of their methodology and the exclusive focus on savanna sites. The growth rate among sites with rainfall above 800 mm/year pointed to a point of attraction above 100% (Figure 6a). Growth in woody cover is likely only slowed by space restriction as it nears canopy closure. Beyond canopy closure, water availability continues to impact forest communities by constraining above-ground woody biomass (Stegen et al., 2011).

| CONCLUSION S
We found that woody encroachment is ongoing and widespread across African savannas with a mean rate of 0.25% per year.
Although our results do not explain the overall trend, they show that change in woody cover is significantly related to a proxy for water availability (potential coverÀinitial cover) and fire frequency. The trend of higher encroachment rates with higher water availability (defined by the difference between actual and potential cover) was consistent across a range of MAP from arid to mesic systems. Competition for water thus universally constrains tree growth, also in higher rainfall (MAP >800 mm/year) savannas. Our approach to derive estimates of potential woody cover based on rates of canopy growth, and initial cover, demonstrates that temporal dynamics in woody populations can be used to infer resource limitations. Competition for water constrains tree population growth rates as canopy cover approaches potential cover, and leads to higher mortality among populations that have over-shot their potential cover based on long-term average rainfall. The results of this study highlight the central role of water availability and competition among trees in determining woody cover growth rates and potential canopy cover.

DATA ACCESSIBILI TY
The data used in the analysis will be provided as a supplement of this paper.