Identifying Ecosystem Function Shifts in Africa Using Breakpoint Analysis of Long-Term NDVI and RUE Data

: Time-series of vegetation greenness data, derived from Earth-observation imagery, have become a key source of information for studying large-scale environmental change. The ever increasing length of such series allows for a range of indicators to be derived and for increasingly complex analyses to be applied. This study presents an analysis of trends in vegetation productivity—measured using the Global Inventory Monitoring and Modelling System third generation (GIMMS3g) Normalised Difference Vegetation Index (NDVI) data—for African savannahs, over the 1982–2015 period. Two annual metrics were derived from the 34 year dataset: the monthly, smoothed NDVI (the aggregated growth season NDVI) and the associated Rain Use Efﬁciency (growth season NDVI divided by annual rainfall). These indicators were then used in a BFAST-based change-point analysis, allowing the direction of change over time to change and the detection of one major break in the time-series. We also analysed the role of land cover type and climate zone as associations of the observed changes. Both methods agree that vegetation greening was pervasive across African savannahs, although RUE displayed less signiﬁcant changes than NDVI. Monotonically increasing trends were the most common trend type for both indicators. The continental scale of the greening may suggest global processes as key drivers, such as carbon fertilization. That NDVI trends were more dynamic than RUE suggests that a large component of vegetation trends is driven by precipitation variability. Areas of negative trends were conspicuous by their minimalism. However, some patterns were apparent. In the southern Sahel and West Africa, declining NDVI and RUE overlapped with intensive population and agricultural regions. Dynamic trend reversals, in RUE and NDVI, located in Angola, Zambia and Tanzania, coincide with areas where a long-term trend of forest degradation and agricultural expansion has recently given way to increases in woody biomass. Meanwhile in southern Africa, monotonic increases in RUE with varying NDVI trend types may be indicative of shrub encroachment. However, all these processes are small-scale relative to the GIMMS NDVI data, and reconciling these conﬂicting drivers is not a trivial task. Our study highlights the importance of considering multiple options when undertaking trend analyses, as different inputs and methods can reveal divergent patterns.


Introduction
The grasslands and savannahs of Africa are hot-spots of contemporary environmental change. Rainfall across the continent varies on decadal time scales, with prolonged droughts impeding vegetation in water-limited regions [1]. In addition, the population of Africa has grown by roughly

NDVI Data
We used the GIMMS3g NDVI dataset for the period 1982-2015, which is derived from the AVHRR sensors carried on-board the seven National Oceanic and Atmospheric Administration (NOAA) satellites. NDVI is calculated using the standard equation: where, ρN IR is reflectance in the near infrared range, and ρRed is reflectance in the visible red part of the electromagnetic spectrum [26]. The GIMMS3g dataset is generated from AVHRR imagery. The raw daily 1 km pixels are binned into optimum 15-day, 8 km composites. The pre-processing routines include empirical mode decomposition and Bayesian methods to remove non-vegetation related artefacts, such as orbital drift, solar zenith angle effects, and sensor transitions. A radiative transfer model is used to correct for the volcanic aerosol loading caused by El Chichón and Mount Pinatubo eruptions [27].

Productivity Indicators
We processed the raw 15-day GIMMS data into maximum monthly value composites and discarded the months for the incomplete 1981 year, resulting in a 34-year (408 months) time-series. We interpolated missing months and smoothed the series using a Savitzky-Golay filter (parameters: order = 1, length = 3, scaling factor = 30). We then generated an annual aggregation of the NDVI values within the growing season (i.e., the small integral). If a pixel contained two growth seasons, mainly areas in East Africa with bimodal rainfall, both seasons were included.
RUE was calculated as the ratio of growth season NDVI over total annual rainfall, using the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) dataset [28]. The CHIRPS data was resampled to match the resolution of the GIMMS NDVI layers using nearest neighbour resampling.
To spatially refine our area of analysis to savannah regions, we calculated the median NDVI using the full monthly time-series and masked pixels with a median NDVI < 0.15 or > 0.8 from further analysis. This removed areas with very low vegetation cover (e.g., the Sahara and Namib deserts) and dense forests (e.g., the Congo and Guinean forests), and is in line with thresholds used in other studies [20].

Time-Series Analyses
To detect shifts in the productivity indicators we used the Breaks For Additive Season and Trend (BFAST) process [19]. BFAST is a time-series decomposition technique, modelling a time-series using parameters quantifying the within-season fluctuations, the underlying trend, and residual noise. When applied to annually aggregated series, as in this case, the model contains only the trend and residual components. Potential non-linear changes in the decomposed parameters are represented by incorporating piecewise linear model segments. Where a series is better represented by multiple linear segments, opposed to a single monotonic trend, a breakpoint is established and dated using structural change tests.
We used the BFAST01 implementation of BFAST, which produces a model containing either one or zero breakpoints, for both of our productivity indicators [29]. The rationale for limiting analysis to a single breakpoint is that one (significant) break is likely to represent the most ecologically relevant shift in a time-series [20]. We limit break detection to include segments of seven years or longer. Trend changes were identified by four structural change tests: the Bayesian information criterion (BIC), the Lagrange multiplier (SupLM), moving sum of residuals (OLS-MOSUM), and the F-test (supF).
A change point was set if any single test identified a significant break (p < 0.05). Breakpoints detected by this process are similar but not synonymous with the concept of ecological tipping points which are generally considered irrecoverable [30]. Here, we use the term 'ecosystem function shift' to indicate a significant development in large-scale ecological function, but make no inference on the reversibility of these processes.
Breakpoints identified by the BFAST01 process were classified into change types, based on direction and significance of the constituent linear models and any identified breakpoint. Our classification scheme was largely based on de Jong [29]. As we analyse RUE trends, we do not refer to increasing or decreasing trends as greening or browning, as patterns may be due to rainfall changes, not vegetation condition. We identified the following six change types: • Monotonic: increase (a significant increase with no significant break detected); • Monotonic: decrease (a significant decrease with no significant break detected); • Interruption: increase with negative break (a significant increase with a significant break followed by a significant increase); • Interruption: decrease with positive break (a significant decrease with a significant break followed by a significant decrease); • Reversal: increase to decrease (a significant increase with a significant break followed by a significant decrease); • Reversal: decrease to increase (a significant decrease with a significant break followed by a significant increase).
We applied BFAST01 models on both growth season NDVI and RUE time-series. All analyses were conducted in the R statistical environment [31]. The raster and rgdal packages provided spatial data handling functionality [32,33]. The GIMMS3g data were downloaded and processed using the gimms package [34].Trends were calculated in the bfast package [19]. Finally, generic data handling and plotting functions were provided by the tidyverse and ggplot2 packages [35].

Interpretation of Classified Trends
To analyse the classified trends, we incorporated land cover and climatic zone data to test for broad-scale patterns in the trend types. A land cover map was obtained from the European Space Agency's Climate Change Initiative (ESA CCI) product. This is series of multi-sensor derived 300 m spatial resolution land cover maps for 1992-2015, designed to be consistent with the United Nations' Land Cover Classification System [36]. We used a map for the year 2000, as this is broadly central in our time-period. The raw data were reclassified into a simplified typology of five types: grassland, shrubland, croplands, woodlands, and mosaic (heterogeneous mixture of different types). Smaller classes covering urban and wetland areas were discarded. The 300 m layer was resampled to match the NDVI-derived layers using nearest neighbour resampling ( Figure 1a).
Climate zone classifications were extracted from the United Nations' Food and Agriculture Organization (FAO) Global Agro-Ecological Zone (AEZ) database. The AEZ is a broad-scale climatic classification based on the Climate Research Unit (CRU) precipitation and temperature datasets [37]. The derived zones are designed to be indicative of the climatic factors that limit vegetation growth. The raw data were resampled to match the NDVI-derived layers using nearest neighbour resampling (Figure 1b).

General Patterns
In this study we apply a trend break time-series analysis on two complimentary indicators of ecosystem productivity-the growth season integrated NDVI and the corresponding RUE for African savannah and grassland regions. By comparing the spatio-temporal distribution of trends, and changes in the trend direction, shifts in ecosystem function can be inferred [20,29]. The spatial distribution of the trend types, their significance, and the year in which any breakpoint occurred for NDVI and RUE, are shown in Figure 2. Figures 3 and 4 show these layers stratified according to the land cover types and agro-ecological zones of Figure 1.
Our most notable observation is that both NDVI and RUE are highly dynamic, with 38% and 52% of pixels having a statistically significant trend, respectively (Figure 2a,d). This approximately corresponds to an area of 7 million km 2 for NDVI and 5.5 million km 2 for RUE. Monotonic increases form a majority of significant trends for both indicators, constituting 60% of NDVI and 49% of RUE pixels with a trend (37% and 23% of all pixels, for NDVI and RUE, Figure 2d-f). The spatial distribution of RUE trends were more homogenous than for NDVI, with more contiguous areas of distinct trend types occurring at specific years ( Figure 2d).
The number of trend breaks identified in the growth season NDVI increased in the late 1990s (Figures 3c and 4c). This increase was initially concentrated in semi-arid and sub-humid agro-ecological zones, with arid zones featuring a later increase, around 2006 (Figure 4c). This spike in change points was not reflected in the RUE time-series, which remained relatively stable throughout the study period, with the exception of a spike in 2004. Relative stability in RUE and a prominence in water-limited semi-arid zones indicates that the increase in trends breaks around the millennium is precipitation driven. By the late 1990s the Sahel had recovered from the earlier drought, and certain areas were either stabilising or transitioning into a decreasing trend (Figure 2a). A large cluster of NDVI trend breaks (decrease with positive break type) occurred in western Tanzania around 1998, coinciding with the onset of a prolonged drought following a wet El Niño year [38,39]. Again, this cluster does not significantly overlap with RUE trends. Controversially, Zhao and Running [40] postulated a general increase in negative vegetation trajectories since the millennium, driven by increasing drought and temperature in the southern Hemisphere [18,41]. Although we observe an increase in trend change points, these do not appear to relate to ecosystem state shifts, as the corresponding RUE is relativity stable-implying precipitation driven changes. Furthermore, there was no continental (or hemispheric) shift towards negative trends in either RUE or NDVI and the overall trend continues to be of vegetation greening.   Many of the detected change points were only identified by one or two structural change tests, out of a possible four, increasing the likelihood that detected change points include false positives ( Figure 5). There was little overlap between areas of multiple positive tests for NDVI and RUE. For NDVI clusters of multiple break identification occurred in central and East Africa, participially Zimbabwe and Tanzania; for RUE western Angola and Sudan were notable hotspots.
On a continental scale, trend breaks in NDVI increased around the millennium, in line with earlier studies [29]. However, there was no noticeable up-tick during the major El Niño/La Niña events in 1997 and 2002, respectively. These events produce well documented regional anomalies in vegetation, but do not appear to initiate a functional shift [38]. Furthermore, the eruption of Mount Pinatubo in 1991 did not produce a spike of change points. The increased aerosol load from this eruption caused a temporary decline in tropical forest photosynthesis [42], which would not be expected to occur in our study region. Nevertheless, the absence of an apparent effect indicates the competence of the GIMMS pre-processing routines, as uncorrected aerosols would be expected to reduce NDVI values. The dominant trend dynamic for both NDVI and RUE was monotonic increases (Figure 2a,d). This is in agreement with other studies using both NDVI-derived indicators and ecosystem models [6,43,44]. As increases occur almost pan-continentally, large-scale drivers should be considered as contributing factors for globally observed trends [44]. Over the past 40 years, biogeochemical changes have occurred, such as increases in reactive nitrogen deposition and elevated atmospheric CO 2 , which increase the maximum potential vegetation production [44][45][46]. Furthermore, the climatic constraints on plant growth have weakened, as rainfall has generally increased, whilst air temperature changes have been relatively minor [47,48]. Combined, these factors would lead to an assumption of increasing vegetation turnover and production, especially for semi-arid regions, where constrains on plant growth are higher [49]. That RUE is increasing in synchrony with NDVI supports non-climatic causation, as RUE trends are adjusted for changes in rainfall limitations.

Regional Patterns
Southern Africa, and especially South Africa, featured predominantly monotonic increases in RUE, interspersed mainly with some interrupted increase break type or a decrease to increase reversal. These increases coincide with a variety of NDVI trend types, including: monotonic increases, increase to decrease as well as decrease to increase reversals. That RUE is unaffected by the underlying trend in NDVI suggests an ecological transition: for southern Africa, most likely, shrub encroachment [50,51]. Increase in shrub cover, at the expense of grasses, increases RUE, as woody plants are more effective at exploiting residual soil moisture [52]. Pixels with changing trend directions in NDVI show a strong spatial alignment with areas of increasing maximum foliage coverage that [46] mapped and attributed to CO 2 fertilisation. Shrub encroachment may produce a variety of NDVI trends, as shrublands have a less pronounced seasonal peak, due to maintaining canopy cover longer than the more ephemeral grasslands. Increases in shrub cover is considered a land degradation process in southern Africa, highlighting the need for caution in interpreting indicators of vegetation productivity without consideration of the local ecological context. Central Africa, covering Angola, Zambia and northern Mozambique, was a highly dynamic region, featuring extensive decrease to increase trend reversals, for both NDVI and RUE (Figure 2a,d). Ecologically, this zone marks a transition from southern African savannah to semi-tropical woodland. This region is notable for the range and dynamism of land use/land cover changes. The end of the 1975 to 2002 Angolan civil war, facilitated an increase in agriculture, which materialised through intensification of cropping frequency within shifting cultivation [2,53]. Further agricultural expansion has occurred in Zambia, as part of an up-serge in soybean plantations [54]. These changes have driven declines in NPP and forest cover, as identified by References [55][56][57]. Yet, in recent years, woody biomass has increased, due to rapid regrowth and shrub expansion [3,58]. Reconciling these dynamic processes within 8 km-resolution pixels, is a challenging endeavour. Over the long-term (1980-present), it is plausible that forest degradation, through agricultural expansion and charcoal harvesting, has driven negative trends, with recent regrowth reversing this process. RUE reversals in this region are clustered around the late 2000s, corresponding to a period of increasing biomass mapped by McNicol et al. [3].
The Sahel and West Africa displayed a gradient in NDVI trends, with the northern Sahel a contiguous band of monotonic increase, the southern Sahel predominantly increasing to decreasing reversals, and the Guinean coast monotonically decreasing (Figure 2a). The trend reversals occurred around the early 2000s (Figure 2c), coinciding with early reports of a slowdown in vegetation gains around the millennium [29]. 'Increasing' pixels in the North are in the more arid regions, which have seen the largest increases in rainfall in the post-drought years, leading to increases in woody plant cover [59]. 'Decreasing' pixels show an alignment with regions of high and rapidly increasing human population. Population growth was identified by Brandt et al. [4] as a localised constraint on continental-scale woody cover increases between 1992-2011. The mechanism for this decrease was postulated as increasing demand for fuelwood and agricultural land, driving the degradation of woodlands and clearing of savannahs. As these decreases in NDVI align with decreasing RUE, changes in vegetation-not rainfall-are supported. Pixels featuring trend reversals align with intensive agricultural zones, where greening processes are modulated by underlying soil types [60].
Multiple trend reversals and interruptions were identified in both indicators over East Africa, in particular, Tanzania, but also Kenya and Ethiopia. Analysis in this area may be complicated by the bimodal rainfall regimes. The wider eastern African region, as noted by Vrieling et al. [61], is experiencing significant decreases in the length of the growing season(s), with vegetation senescence occurring up to one month earlier by 2010 than in the early 1980s. The most pronounced declines in growth season length in western Tanzania, show a strong geographic alignment with monotonically decreasing NDVI pixels [61]. Trend reversals may also be driven by a rapid shift from the wet El Niño year of 1997 into a multi-year drought [38,39].

Limitations
This study, as with all GIMMS NDVI studies, has a number of limitations. Firstly, reconciling coarse resolution pixels with on-the-ground studies and observations is non-trivial, as very few, if any, field studies cover the area of even a single 64 km 2 pixel. Secondly, the use of temporal trends in NDVI is based on the assumption that changes in vegetation will be captured by the trend statistic(s). This is not always the case: simulations have shown that in some contexts, major reductions in NDVI are required to produce negative significant trends [17]. Finally, many detected breaks were not identified by multiple structural change statistics, increasing the likelihood for false positives in the trend change layers.

Conclusions
We analysed long-term NDVI data, processed into two complementary indicators of ecosystem productivity, for the 1982-2015 period-the aggregated growth season NDVI and associated RUE. To detect ecosystem functional shifts, we employed a trend breaks approach. Our main observations are as follows: (1) for both NDVI and RUE, monotonic increase was the most common trend type observed; (2) RUE was less dynamic than NDVI, indicating that a large component of NDVI variability is driven by rainfall variability; and (3) areas where trends reverse or are interrupted align with known ecological shifts, including a host of human ecosystem modification processes. These results highlight the benefit, and necessity, of considering multiple approaches to trend analysis of vegetation index data. A reliance on a any single indicator or analysis, may miss important processes, and comparisons of different approaches can produce meaningful information. Funding: TH was supported by an MMU post-graduate fellowship. ES was supported by an EU Marie Currie Career Integration Grant (PCIG12-GA-2012-3374327). We are grateful to NASA and the USGS for the NDVI data.

Conflicts of Interest:
The authors declare no conflict of interest.