Dryland Vegetation Functional Response to Altered Rainfall Amounts and Variability Derived from Satellite Time Series Data

Vegetation productivity is an essential variable in ecosystem functioning. Vegetation dynamics of dryland ecosystems are most strongly determined by water availability and consequently by rainfall and there is a need to better understand how water limited ecosystems respond to altered rainfall amounts and variability. This response is partly determined by the vegetation functional response to rainfall (β) approximated by the unit change in annual vegetation productivity per unit change in annual rainfall. Here, we show how this functional response from 1983 to 2011 is affected by below and above average rainfall in two arid to semi-arid subtropical regions in West Africa (WA) and South West Africa (SWA) differing in interannual variability of annual rainfall (higher in SWA, lower in WA). We used a novel approach, shifting linear regression models (SLRs), to estimate gridded time series of β. The SLRs ingest annual satellite based rainfall as the explanatory variable and annual satellite-derived vegetation productivity proxies (NDVI) as the response variable. Gridded β values form unimodal curves along gradients of mean annual precipitation in both regions. β is higher in SWA during periods of below average rainfall (compared to above average) for mean annual precipitation <600 mm. In WA, β is hardly affected by above or below average rainfall conditions. Results suggest that this higher β variability in SWA is related to the higher rainfall variability in this region. Vegetation type-specific β follows observed responses for each region along rainfall gradients leading to region-specific responses for each vegetation type. We conclude that higher interannual rainfall variability might favour a more dynamic vegetation response to rainfall. This in turn may enhance the capability of vegetation productivity of arid and semi-arid regions to better cope with periods of below average rainfall conditions.


Introduction
The importance of rainfall and water availability for dryland vegetation functioning has long been recognized and shown in numerous studies [1][2][3][4]. That is, strong interannual fluctuations in rainfall amounts drive an analogous variability in vegetation dynamics from year to year [5].
In particular, vegetation net productivity is a key variable in terrestrial ecosystem functioning [6], as it represents the underlying driving force of all trophic ecosystem levels [7] and carbon budgets [8] and thus connects plant-scale metabolic processes with ecosystem structure [9].
Drylands cover nearly 40% of the global terrestrial surface and make up a similar number of the global terrestrial share of net primary production [10]. Thus, climate model predictions of globally altered water availability in the upcoming century [11] make it necessary to enhance the understanding of dryland vegetation response to rainfall and future changes thereof [12,13].
Vegetation response to rainfall (β) may be defined as any change in vegetation activity and, consequently, production or cover over time with a corresponding change in rainfall over the same time period for a given location [14,15]. The response is regularly estimated from linear models as slope or linear coefficient ingesting rainfall as explanatory and vegetation productivity (or a proxy thereof) as the response variable [3,6,[16][17][18][19][20]. β thus translates water availability during a period of net growth into any ecosystem or plant-specific functional response as a function of the integral of other constraints potentially influencing plant growth [21]. β can be interpreted as the Water Use Efficiency (WUE) (i.e., the ratio of net photosynthesis to transpiration) under preconditions that the ineffective precipitation (runoff and soil evaporation) was accounted for leading to an intercept = 0. In this case, β becomes a measure of WUE [5,22]. Deriving from WUE, the concept of Rain-Use Efficiency (RUE), defined as the ratio of ANPP (aboveground net primary productivity) to annual precipitation [1] and temporal changes therein was suggested as an integral measure for evaluating land degradation and desertification. It has been widely used to assess non-precipitation related land degradation-or the reverse-in global drylands [23][24][25]. As emphasized in [14,17,23], this may only yield reliable results in cases where the intercept of the relationship between vegetation and water metrics used equals zero (otherwise the β and efficiency will be different entities) which is seldom obtained by using satellite based proxies of vegetation productivity such as annually summed NDVI (normalised difference vegetation index) and annually summed precipitation.
Moreover, the basic assumption involved in the use of RUE is that vegetation productivity is linearly related to precipitation in the absence of human-induced land degradation. If this assumption of linearity does not hold, the normalization for precipitation variability in the spatial or temporal domain does not apply and, consequently, the interpretation of changes in RUE as non-precipitation related land degradation is not justified. Several studies have questioned this linearity or constant slope to apply across drylands characterized by different plant functional types or hydrological conditions. Along a gradient of mean annual precipitation (MAP) stretching across arid to semi-arid regions, linear coefficients usually form a unimodal response curve [13,[26][27][28] following a shift of the relative importance of the factors limiting plant growth and, consequently, any response of vegetation to rainfall [6]. In the arid-most parts (approximately 0-300 mm MAP) physiological constraints such as adaptions of stomatal conductance to high vapour pressure deficit and low specific leaf area affect the relative growth rate [29] and, consequently, limit the plant-specific β. Approaching the semi-arid regions around 300-600 mm MAP, β increases towards a region-specific local maximum [1,6,27] marking the transition zone from physiological to rather edaphic constraints on plant growth [3]. A rapid decrease in β is often observed after the peak ascribed to the increasingly limiting power of nutrients and a decrease of the governing power of precipitation [6,30]. Moreover, local fire regimes and land use practices increasingly shape vegetation dynamics [15,31,32].
There have been several attempts to model the effect of below and above average rainfall on β using a space for time approach, e.g., [13,18]. Other studies have split time series of vegetation productivity and rainfall into two separate periods a priori assuming differences in the prevailing rainfall conditions, e.g., [33,34]. However, no attempt to describe any ecosystem-specific reaction to above or below average rainfall conditions with respect to ecohydrological traits such as β has been made on a regional scale based on gridded data. Additionally, there is evidence that greater environmental heterogeneity (such as rainfall variability) affects the way ecosystems react with changing resource availability [35,36]. Lázaro-Nogal et al. [37] found greater phenotypic plasticity being associated with higher rainfall variability in dryland plant communities which may in turn affect the adaptability of net production to altered annual rainfall. Thus, it may be hypothesized that higher rainfall variability promotes a more dynamic vegetation functional response to rainfall. Consequently, integrating the unimodality framework with proxies for environmental heterogeneity over wet and dry periods will help to reveal potential regional-scale vulnerabilities to changes in rainfall.
In this study, we present results from a novel approach to analyze the spatiotemporal relationship between vegetation and rainfall in drylands providing a basis for predictions of future dryland vegetation sensitivities to short-and long-term changes in rainfall. The approach involves shifting linear regression models (SLR) producing a time series of linear slope coefficients based on gridded time series of vegetation productivity proxies and rainfall estimates covering nearly three decades . We explore the spatial response functions of β along gradients of MAP during relatively dry and wet conditions for two dryland regions in subtropical Africa (West Africa, WA, and South West Africa, SWA) characterized by distinctly different patterns of interannual variability of annual rainfall.
Specifically, (1) we test the importance of interannual rainfall variability on the response to changes in below or above average rainfall conditions (dry or wet) assessed by the shape and amplitude of the β curve over a MAP gradient for the two study regions; (2) We analyze β for the three major vegetation types of African drylands (savannas, grass and crop type vegetation, and shrubs) to reveal vegetation-specific differences.

Study Regions
The importance of rainfall amount and variability for vegetation response was assessed in two major African dryland regions, namely South West Africa (SWA) and West Africa (WA) (Figure 1a,b). The MAP gradients in these regions range from approximately 20 mm to >900 mm. However, only regions with MAP <900 mm were considered as above this threshold rainfall increasingly loses its predetermining influence on vegetation [15,38]. The SWA study region (Figure 1a) is characterized by a southern-centric summer rainy season which lasts approximately from November to April with the south-western parts experiencing significantly shorter seasons compared to the north-east. Vast areas of this region are covered by shrub and savanna vegetation and only minor parts are characterized by grassland ( Figure 1c). The WA study region ( Figure 1b) is characterized by a Northern hemisphere summer precipitation regime receiving nearly all annual rainfall during the rainy period between July and October. Large parts of this region are characterized by grassland or savanna vegetation ( Figure 1d). A more detailed description of the study regions can be found in the Appendix A.

Vegetation Productivity Proxies
Global Inventory Modelling and Mapping Studies normalized difference vegetation index (GIMMS NDVI3g, 8 km spatial resolution, biweekly temporal resolution) [39] data were used to derive annual proxies of vegetation productivity from 1983 to 2011. The NDVI is closely coupled to the fraction of photosynthetically active radiation absorbed by vegetation (fAPAR) in dryland areas [40] and has been shown to be closely linked to ecosystem scale carbon fixation [41] and, consequently, net biomass production [42,43]. Particularly, the cyclic part of the greenness signal measured throughout the growing season is most strongly coupled to vegetation productivity [44,45].
We used a phenological parameterization model [46] to derive the cyclic part of the annual vegetation signal (Appendix B). We subsequently used this annual "cyclic fraction" as a proxy for vegetation productivity. Due to the unitless character of the NDVI, we normalized annual productivity proxies to z-scores as z i = (x i − µ)/σ where z i represents the z-score of the vegetation proxy value x i of time step i of the time series and µ and σ are the time series mean and standard deviation, respectively, for any given grid cell. Insets are enlarged maps of the respective gradients of mean annual precipitation (MAP) derived from African Rainfall Climatology 2 (ARC2) data. White areas on the rainfall maps indicate areas, which were not considered as they were either masked as water bodies or had MAP values > 900 mm. (c-d) Vegetation types found in the South West African (SWA, (c)) and West African (WA, (d)) study regions derived from MODIS MCD12C1 Type 3 land cover data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). White areas depict pixels excluded. Pixels have been nearest-neighbour-resampled to match the GIMMS spatial resolution.

Rainfall Data
Data of the African Rainfall Climatology version 2, (ARC2, 0.1° spatial resolution, available from 1983 to present) [47] were obtained and bi-linearly resampled to match the GIMMS NDVI spatial resolution. The daily rainfall estimates were then summed annually, where one year in the WA region lasts from January to December (Northern hemisphere summer centric) and from August to July in the SWA region (southern hemisphere summer centric). In months where no ARC2 data were available, a monthly estimate from the Global Precipitation Climatology Project (GPCP) v2p2 product [48] was resampled to GIMMS NDVI resolution and scaled to monthly rainfall sums assuming a month of length 365/12 days. To maintain comparability, annual rainfall sums were transformed to z-scores following the procedure described for vegetation productivity proxies.
To analyze β response functions along MAP gradients, we calculated pixelwise mean annual precipitation (MAP (mm)) and mean rainy season length (MSL (months)). A month was considered as belonging to the rainy season when it received more than 20 mm of rainfall. Although this is below the value often assumed for "effective rainfall", e.g., [5,49], it is well in the range commonly used to differentiate between rainy and dry season [50]. To account for the variability of annual rainfall and season length, we computed the pixel-wise coefficient of variation from MAP (CVP) and MSL (CVS) as CV = σ/µ. Insets are enlarged maps of the respective gradients of mean annual precipitation (MAP) derived from African Rainfall Climatology 2 (ARC2) data. White areas on the rainfall maps indicate areas, which were not considered as they were either masked as water bodies or had MAP values > 900 mm. (c-d) Vegetation types found in the South West African (SWA, (c)) and West African (WA, (d)) study regions derived from MODIS MCD12C1 Type 3 land cover data (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). White areas depict pixels excluded. Pixels have been nearest-neighbour-resampled to match the GIMMS spatial resolution.

Rainfall Data
Data of the African Rainfall Climatology version 2, (ARC2, 0.1 • spatial resolution, available from 1983 to present) [47] were obtained and bi-linearly resampled to match the GIMMS NDVI spatial resolution. The daily rainfall estimates were then summed annually, where one year in the WA region lasts from January to December (Northern hemisphere summer centric) and from August to July in the SWA region (southern hemisphere summer centric). In months where no ARC2 data were available, a monthly estimate from the Global Precipitation Climatology Project (GPCP) v2p2 product [48] was resampled to GIMMS NDVI resolution and scaled to monthly rainfall sums assuming a month of length 365/12 days. To maintain comparability, annual rainfall sums were transformed to z-scores following the procedure described for vegetation productivity proxies.
To analyze β response functions along MAP gradients, we calculated pixelwise mean annual precipitation (MAP (mm)) and mean rainy season length (MSL (months)). A month was considered as belonging to the rainy season when it received more than 20 mm of rainfall. Although this is below the value often assumed for "effective rainfall", e.g., [5,49], it is well in the range commonly used to differentiate between rainy and dry season [50]. To account for the variability of annual rainfall and season length, we computed the pixel-wise coefficient of variation from MAP (CVP) and MSL (CVS) as CV = σ/µ.

Land Cover Data
To determine the vegetation-type-specific response to rainfall we used Moderate-Resolution Imaging Spectroradiometer (MODIS) MCD12C1 (0.05 • spatial resolution, annual temporal resolution, 2001-2011) Type 3 (LAI/fAPAR) land cover data (nearest-neighbour resampled to match GIMMS NDVI spatial resolution). This product provides land cover information and sub-pixel frequencies of each class per grid cell. Information about long-term stable vegetation type was obtained as follows [51]: The sub-pixel land cover frequency was obtained per pixel, land cover class and year. For each year, each pixel was categorized into the land cover class with the highest sub-pixel frequency if the frequency was above 80%. If it was below this threshold, the pixel was disregarded from further analysis. We only considered land cover classes "Grasses/Cereal crops", "Shrubs" and "Savanna" as those classes comprise the vegetation types most frequently found in the study regions.

Shifting Linear Regression Models
To obtain estimates of vegetation response to rainfall we computed shifting linear regression models (SLRs) and derived the linear slope coefficients β. The coefficients were obtained using ordinary least square methods (OLS). The SLR is computed for a given set of years of the time series with temporal window length L, where rainfall is the explanatory variable and the vegetation productivity proxy is the response. In the next step, the window is shifted by one year and the model is recalculated to derive a new linear coefficient (β). For a more detailed description of the SLR procedure see Appendix C. For each grid cell, thus a time series of linear slope coefficients of length n + 1 − L was obtained where n is length of the input data time series and L is the temporal window length over which the SLRs are calculated. Typically, this length should at least span 5-10 years [13,52]. We used L of length 11 years. It was consequently chosen to be long enough to capture the typical interannual vegetation response to rainfall yet short enough to reflect the effects of dry and wet periods on β.
All pixels reporting negative slope values were omitted from further analysis as assumed to be a product of poor data quality [27] or containing vegetation not controlled by rainfall. This amounted to a removal of 16.8 per cent and 17.2 per cent of the pixels on average per model time step for WA and SWA, respectively.

Spatial Analysis
All β coefficients were consequently assigned one of two classes; "dry" or "wet" representing the above or below average rainfall during the respective model time step i. Thus, if MAP i > MAP for any given grid cell that pixel for time step i was assigned the class "wet", otherwise "dry". All coefficients were then averaged over 1 mm MAP steps and the classes dry and wet.
We used Generalized Additive Models (GAMs) [53,54] implemented in the R computing environment [55] to derive smooth response curves with MAP as the explanatory variable and the linear coefficients (averaged over 1 mm MAP steps) as the response. The models were parameterized assuming normal error distributions and thus identity link functions.
In a next step, we analyzed differences in peak β with respect to absolute β values (β max ), MAP position. For this, the upper 90th percentile of the fitted response curve of β (β max ) was extracted. We further analyzed the specific β of the three main vegetation types found in the two study regions (savanna, shrub and grass/cereal crop type vegetation) with respect to above and below average rainfall. Figure 2 provides an overview of all processing steps involved.

Figure 2.
Flowchart of the data processing and analysis. MSL refers to mean season length, CVS is the coefficient of variation of season length, MAP represents mean annual precipitation, CFR refers to cyclic fraction, SLR abbreviates shifting linear regression models, β represents vegetation functional response to rainfall, dry/wet refer to below and above average rainfall conditions, respectively, GAM represents generalized additive models. The colouring scheme refers to the figure in which the respective analysis step is presented.

Results
An important difference in the rainfall variability regime exists between the two regions of South West Africa (SWA) and West Africa (WA) ( Figure 3). Interannual variability of annual rainfall strongly increases with decreasing MAP (Figure 3a) for both regions. SWA generally experiences higher interannual variability of annual rainfall for any given MAP as compared to WA. However, the inverse pattern is observed for season length where WA generally is characterized by a higher interannual variability of the wet season length as compared to SWA (Figure 3b).  Flowchart of the data processing and analysis. MSL refers to mean season length, CVS is the coefficient of variation of season length, MAP represents mean annual precipitation, CFR refers to cyclic fraction, SLR abbreviates shifting linear regression models, β represents vegetation functional response to rainfall, dry/wet refer to below and above average rainfall conditions, respectively, GAM represents generalized additive models. The colouring scheme refers to the figure in which the respective analysis step is presented.

Results
An important difference in the rainfall variability regime exists between the two regions of South West Africa (SWA) and West Africa (WA) ( Figure 3). Interannual variability of annual rainfall strongly increases with decreasing MAP (Figure 3a) for both regions. SWA generally experiences higher interannual variability of annual rainfall for any given MAP as compared to WA. However, the inverse pattern is observed for season length where WA generally is characterized by a higher interannual variability of the wet season length as compared to SWA (Figure 3b). Flowchart of the data processing and analysis. MSL refers to mean season length, CVS is the coefficient of variation of season length, MAP represents mean annual precipitation, CFR refers to cyclic fraction, SLR abbreviates shifting linear regression models, β represents vegetation functional response to rainfall, dry/wet refer to below and above average rainfall conditions, respectively, GAM represents generalized additive models. The colouring scheme refers to the figure in which the respective analysis step is presented.

Results
An important difference in the rainfall variability regime exists between the two regions of South West Africa (SWA) and West Africa (WA) ( Figure 3). Interannual variability of annual rainfall strongly increases with decreasing MAP (Figure 3a) for both regions. SWA generally experiences higher interannual variability of annual rainfall for any given MAP as compared to WA. However, the inverse pattern is observed for season length where WA generally is characterized by a higher interannual variability of the wet season length as compared to SWA (Figure 3b).  The response functions of β along MAP gradients based on generalized additive models (GAMs, Figure 4) exhibit unimodal shapes in both regions during dry and wet periods, although particular inter-regional differences exist. All fitted models are significant (p < 0.001, Table 1) and the overall fit as indicated by the R 2 value is higher in the SWA region as compared to WA. Both regions have higher model R 2 during dry periods compared to wet periods (Table 1).
the shape of the curve as compared to SWA. Generally, β values for a given MAP are equal or lower in WA than in SWA and there is no consistent difference in the response curves of dry and wet periods in WA.
Peak β values (βmax) in the SWA region are on average 34% higher compared to the WA region (Table 1, Figure 4c). Differences between the regions account for the largest part of the total data variability. Yet, when each region is considered separately, above or below average rainfall strongly affects βmax in the SWA region compared to WA with peak β values decreasing 21% on average from dry to wet. On the contrary, above or below average rainfall conditions have a small impact on βmax in WA.
Inter-region differences in MAPβ-max values are even more pronounced with MAP values at β peaks in WA being on average 48% higher compared to SWA (Table 1, Figure 4c). In SWA, MAPβ-max values slightly decrease from dry to wet conditions. In WA, on the contrary, MAPβ-max decreases on average 13% from dry to wet conditions.  Peak values β max and the corresponding mean annual precipitation (MAP β-max ). Colouring represents above and below average rainfall conditions over which data have been aggregated. All data have been averaged over 1 mm MAP steps. Table 1. Model summaries of the fitted generalized additive models (GAMs) with n referring to the number of data points used to fit the models; p indicates the significance level derived from the F-statistics of the fitted smoothing terms; EDF is the estimated degrees of freedom of the respective GAM; R 2 is the variability in β explained by MAP with respect to the fitted GAM. β max corresponds to mean peak β; MAP β-max is the mean annual precipitation average of the peak position of β.

Region
Clim. The curves in SWA (Figure 4a) show clearer unimodality for both, wet and dry periods, compared to WA. Moreover, there is a stronger difference between the response curves of dry and wet conditions in South West Africa. While β tends to be lower during wet periods for MAP < 600 mm this reverses for higher MAP as indicated by the intersecting curves. The unimodal shape is less apparent in the WA region (Figure 4b) resulting from a lower dynamic range and a broader peak in the shape of the curve as compared to SWA. Generally, β values for a given MAP are equal or lower in WA than in SWA and there is no consistent difference in the response curves of dry and wet periods in WA. Peak β values (β max ) in the SWA region are on average 34% higher compared to the WA region (Table 1, Figure 4c). Differences between the regions account for the largest part of the total data variability. Yet, when each region is considered separately, above or below average rainfall strongly affects β max in the SWA region compared to WA with peak β values decreasing 21% on average from dry to wet. On the contrary, above or below average rainfall conditions have a small impact on β max in WA.
Inter-region differences in MAP β-max values are even more pronounced with MAP values at β peaks in WA being on average 48% higher compared to SWA (Table 1, Figure 4c). In SWA, MAP β-max values slightly decrease from dry to wet conditions. In WA, on the contrary, MAP β-max decreases on average 13% from dry to wet conditions.

Vegetation Type-Specific Response to Rainfall
The vegetation-specific response to rainfall ( Figure 5) indicates that grass type vegetation shows no significant difference between dry and wet periods irrespective of study region. Moreover, grass/crop β is consistently lower in WA compared to SWA although this difference is small. In West Africa, grass/crop type vegetation has the overall highest β of all vegetation types. Table 1. Model summaries of the fitted generalized additive models (GAMs) with n referring to the number of data points used to fit the models; p indicates the significance level derived from the Fstatistics of the fitted smoothing terms; EDF is the estimated degrees of freedom of the respective GAM; R 2 is the variability in β explained by MAP with respect to the fitted GAM. βmax corresponds to mean peak β; MAPβ-max is the mean annual precipitation average of the peak position of β.

Vegetation Type-Specific Response to Rainfall
The vegetation-specific response to rainfall ( Figure 5) indicates that grass type vegetation shows no significant difference between dry and wet periods irrespective of study region. Moreover, grass/crop β is consistently lower in WA compared to SWA although this difference is small. In West Africa, grass/crop type vegetation has the overall highest β of all vegetation types.
For savanna type vegetation β values in WA are hardly different between wet and dry periods. In SWA, β values are slightly higher during wet periods. The difference of savanna β to grass/crop vegetation in this region is small compared to WA. In both regions, savanna vegetation has the overall lowest β of all vegetation types.
In both regions, shrub type vegetation has a significantly higher β during dry periods compared to wet periods. While shrubs have the overall highest β in SWA they range between grass/crop and savanna type vegetation in WA. In accordance with the overall higher vegetation response to rainfall in SWA, all vegetation types have higher β to in SWA compared to WA. Major vegetation types of the study regions were analyzed as a function of mean annual precipitation (MAP) ( Figure 6). Vegetation types clearly differ in their ranges of MAP. Shrubs dominate for lower rainfall amounts around 250 mm/year. Grasses are mainly found around 400-500 mm MAP and savanna vegetation dominates for MAP > 550 mm. The distribution of vegetation types is relatively similar between the two regions. However, grasses/crops cover a larger gradient in WA and Savanna is mainly found in higher MAP regions in WA as compared to SWA. For savanna type vegetation β values in WA are hardly different between wet and dry periods. In SWA, β values are slightly higher during wet periods. The difference of savanna β to grass/crop vegetation in this region is small compared to WA. In both regions, savanna vegetation has the overall lowest β of all vegetation types.
In both regions, shrub type vegetation has a significantly higher β during dry periods compared to wet periods. While shrubs have the overall highest β in SWA they range between grass/crop and savanna type vegetation in WA. In accordance with the overall higher vegetation response to rainfall in SWA, all vegetation types have higher β to in SWA compared to WA.
Major vegetation types of the study regions were analyzed as a function of mean annual precipitation (MAP) (Figure 6). Vegetation types clearly differ in their ranges of MAP. Shrubs dominate for lower rainfall amounts around 250 mm/year. Grasses are mainly found around 400-500 mm MAP and savanna vegetation dominates for MAP > 550 mm. The distribution of vegetation types is relatively similar between the two regions. However, grasses/crops cover a larger gradient in WA and Savanna is mainly found in higher MAP regions in WA as compared to SWA.

Discussion
This study shows how the differential regional-scale vegetation response to rainfall (β) varies between two African dryland regions based on gridded satellite-derived vegetation productivity proxy and rainfall data. We have shown that a shifting linear regression model can successfully be applied to determine the local vegetation response to rainfall in dry and wet periods, respectively. Moreover, we have demonstrated how the response can be regionalized as a function of mean annual precipitation (MAP), changing water availability, interannual rainfall variability and vegetation type.

Vegetation Response to Rainfall Along MAP Gradients
Our results confirm that vegetation response to rainfall (β) follows the expected unimodality along MAP [6] in both study regions. β is generally higher or at similar levels for any given MAP in SWA compared to WA (Figure 4). Moreover, SWA shows a stronger difference between dry and wet periods, has clearer unimodal shapes along the MAP gradient and possesses higher variability along the MAP gradient. This agrees with the hypothesized higher β dynamics in regions with higher interannual variability of annual rainfall [36,37], such as SWA (Figure 3a). In addition, this suggests that a higher rain use efficiency, which has been observed during dry periods [18], might be particularly favoured by greater rainfall variability.
GAM R 2 values are higher during dry periods compared to wet periods (Table 1). This indicates that during relatively dry conditions MAP is the main determinant of vegetation response to rainfall, while during wetter periods other limitations such as nutrients [30] or land use and fire [38] might become relatively strong factors. Moreover, higher GAM R 2 values in SWA indicate that β in this region is more strongly determined by MAP. A possible fact which might lead to this pattern is the higher population density in WA compared to SWA [56]. More densely populated areas are likely to be characterized by higher land use intensity which in turn might make MAP a slightly less important determinant of β. Additionally, in regions with MAP < 500 mm, SWA is characterized by farm-based livestock keeping whereas WA is characterized by (semi-) nomadic livestock keeping by pastoralists. The latter may lead to spatially diverse effects on the vegetation making MAP a less strong determinant of β.
There are small but apparent deviations of the actual data from the fitted curves such as in SWA around 100 mm MAP during wet periods. Those deviations can be caused by local terrain effects that can lead to higher (e.g., local depressions or valleys through water accumulation) or lower (e.g., ridges through increased runoff) β values [57,58]. Moreover, differences in land use and land use intensity at the local scale may cause deviations [59]. Consequently, β over MAP gradients should preferably be considered over larger scales as in Figure 4. Those small deviations, however, are

Discussion
This study shows how the differential regional-scale vegetation response to rainfall (β) varies between two African dryland regions based on gridded satellite-derived vegetation productivity proxy and rainfall data. We have shown that a shifting linear regression model can successfully be applied to determine the local vegetation response to rainfall in dry and wet periods, respectively. Moreover, we have demonstrated how the response can be regionalized as a function of mean annual precipitation (MAP), changing water availability, interannual rainfall variability and vegetation type.

Vegetation Response to Rainfall Along MAP Gradients
Our results confirm that vegetation response to rainfall (β) follows the expected unimodality along MAP [6] in both study regions. β is generally higher or at similar levels for any given MAP in SWA compared to WA (Figure 4). Moreover, SWA shows a stronger difference between dry and wet periods, has clearer unimodal shapes along the MAP gradient and possesses higher variability along the MAP gradient. This agrees with the hypothesized higher β dynamics in regions with higher interannual variability of annual rainfall [36,37], such as SWA (Figure 3a). In addition, this suggests that a higher rain use efficiency, which has been observed during dry periods [18], might be particularly favoured by greater rainfall variability.
GAM R 2 values are higher during dry periods compared to wet periods (Table 1). This indicates that during relatively dry conditions MAP is the main determinant of vegetation response to rainfall, while during wetter periods other limitations such as nutrients [30] or land use and fire [38] might become relatively strong factors. Moreover, higher GAM R 2 values in SWA indicate that β in this region is more strongly determined by MAP. A possible fact which might lead to this pattern is the higher population density in WA compared to SWA [56]. More densely populated areas are likely to be characterized by higher land use intensity which in turn might make MAP a slightly less important determinant of β. Additionally, in regions with MAP < 500 mm, SWA is characterized by farm-based livestock keeping whereas WA is characterized by (semi-) nomadic livestock keeping by pastoralists. The latter may lead to spatially diverse effects on the vegetation making MAP a less strong determinant of β.
There are small but apparent deviations of the actual data from the fitted curves such as in SWA around 100 mm MAP during wet periods. Those deviations can be caused by local terrain effects that can lead to higher (e.g., local depressions or valleys through water accumulation) or lower (e.g., ridges through increased runoff) β values [57,58]. Moreover, differences in land use and land use intensity at the local scale may cause deviations [59]. Consequently, β over MAP gradients should preferably be considered over larger scales as in Figure 4. Those small deviations, however, are expected to have minor effects on the overall pattern observed (Figure 4) as MAP is the main determinant of vegetation dynamics along the MAP gradients considered in this study [32].
The spatially differential response of vegetation to above or below average rainfall questions the idea of a constant RUE in space and time. We consequently raise some concern about RUE as an integral measure for evaluating land degradation and desertification in various semiarid and arid zones across the globe [1,25,60,61]. Similarly, the more precautious suggestion of a unifying maximum RUE during dry periods [18,19] is challenged by the current results. That is, some MAP regions in fact have notably higher β during drier periods (SWA, MAP < 600 mm). However, vast areas of WA show a limited difference between dry and wet periods or even higher β under wet conditions. This pattern of higher β values during wet periods may point to a positive effect of increasing atmospheric CO 2 on semi-arid vegetation productivity being most evident during periods of above average precipitation [41,62]. This interpretation, however, implies as well that the degree to which an ecosystem responds to rising CO 2 levels may depend on the rainfall variability regime. The higher β curves during wet periods for MAP >600 mm in SWA may point to the existence of a vegetation type-specific transition zone. Below this value, SWA is characterized by grasses and shrubs whereas above savannas dominate the landscape ( Figure 6). Thus, given the savanna-specific β under wet and dry conditions ( Figure 5), higher β under wet conditions in SWA for MAP > 600 mm can be expected.
The phenology of plants in dry environments is closely related to the occurrence and timing of wet seasons. Woody deciduous species, for example, regularly unfold their leaves already before the first rains of a wet season [57] whereas herbaceous plants are triggered by and follow rainfall events very closely with a short delay [63]. For either strategy, however, a consistently recurring rainy season is important. Higher season length variability (Figure 3b, WA) may thus impede an optimal use of water which may ultimately lead to overall lower β values as seen in WA.
Generally, variability of annual rainfall is suggested to have a strong effect on region-scale differences in the vegetation response to rainfall during wet and dry periods. The within-region expression is, however, suggested to be shaped by several other variables such as effects of topography, land use, soil type, variability of wet season length (CVS) and rising atmospheric CO 2 levels.
The ecosystems of the Sahel region covering WA have been reported to have recovered from the severe 1970s and 1980s drought periods [34,42,44,64,65] (see description Appendix A). However, considering the relatively low β during dry periods, it might be suggested that the adverse effects of future dry periods on the vegetation in this region should not be underestimated. Contrastingly, the relatively high β below 600 mm MAP in SWA during dry conditions suggests that vegetation in this region is less vulnerable during dry periods. Above 600 mm, however, the savannas of SWA might be expected to be relatively prone to drought periods ( Figure 4).
We note that rainfall from previous years may well contribute to explaining the patterns of β we have observed despite the fact that soil moisture of drylands returns to the same low level after a long dry season regardless of the amount of rainfall during the previous rainy season. To which degree rainfall from preceding years contributes to the response to rainfall of any given year, however, depends on the annual/perennial composition and age structure of the vegetation in any given dryland ecosystem. This is due to the fact that previous years' growth indeed affects the growth of perennials in any given year, which will in turn depend on the rainfall which fell during those previous years [66]. Further on, there are accumulating effects of increasingly favourable microclimatic, edaphic growing conditions and amount of grass seeds in the topsoil during periods with higher rainfall. We thus accept this influence of previous years' rainfall as leading into any β observed and thus producing an ecosystem-specific functional response to rainfall. An additional analysis (Appendix D), moreover, confirmed that the SLR model only considering rainfall of the same year supersedes models accounting for rainfall of previous years in most instances, spatially as well as temporally.

Dynamics of Peak Vegetation Response to Rainfall
The 90th percentile of fitted β values (β max ) were extracted for SWA and WA as a region-wide characteristic feature to quantify differences in the peaks of vegetation response to rainfall and how the peak correspond to mean annual precipitation. β max is significantly higher in SWA than in WA (Figure 4c). It moreover significantly differs between wet and dry episodes in SWA which does not apply for WA. Both can be expected from the overall higher β in SWA (Figure 4a). It thus further supports the hypothesis that greater rainfall variability (CVP, Figure 3a) may lead to a temporally more dynamic vegetation response to rainfall [37,67]. MAP β-max values (the positions of β max on the rainfall gradients) are significantly lower in SWA compared to WA (Figure 4c). This indicates that the transition zone from rather physiological limitations to more edaphic constraints on β [3] is highly variable between regions.
Although soils certainly shape the vegetation response β to rainfall at the local scale [30], it is clear from Figure 4 that they may have rather an "envelope function" (individual points in Figure 4a,b) in shaping the boundaries of regional scale response to altered rainfall. If soils had a stronger effect, we would not expect a distinct peak β along MAP as found in this study. Although the previously mentioned explanations for the existence of a peak in β (physiological adaptions below the peak, limiting nutrients above the peak [6]) certainly hold true, the MAP range of peak β should not be considered stable within a region.

Vegetation Type-Specific Response to Rainfall
Only small differences in β were observed for grass/crop type vegetation and this vegetation type showed minimal response to below or above average rainfall conditions irrespective of study region. The reason for this is likely that grasses usually regrow their entire above ground biomass annually with only limited possibility to adjust β to dry or wet periods.
Savanna vegetation in both regions shows little difference between wet and dry episodes with a slight tendency towards higher β under wet conditions, particularly in SWA.
The higher β for shrubs in both regions during dry periods might be explained by the fact that vegetation particularly in dry regions is able to increase water use efficiency during dry periods by reducing nitrogen use efficiency and thus maintaining overall photosynthesis rates at similar levels [68]. However, the absence of higher β during dry periods for grass/crop type and savanna vegetation leads to two possible conclusions. Firstly, whether or not a particular vegetation type has higher β and consequently ecosystem scale water use efficiency during dry periods depends on the relative limiting degree of rainfall and thus on the position along the MAP gradient. Savanna and grass type vegetation are both found at considerably higher MAP values than shrubs in both regions ( Figure 6). This might lead to the assumption that the vegetation adaptability to dry periods is most strongly pronounced in arid regions [37]. Secondly, how strongly β increases from wet to dry periods in arid regions depends on the respective variability of annual rainfall with greater interannual variability promoting relatively higher β during dry periods.
All vegetation types show overall higher β in SWA compared to WA. This further supports the hypothesis that higher CVP might promote a more dynamic functional response to rainfall.
Arguably, response curves (Figure 4a,b) are modified by the local vegetation composition with its distinct sensitivity to rainfall. However, MAP is a first order determinant of vegetation composition in arid and semi-arid regions [57] and thus implicitly accounts for the vegetation-specific response to rainfall. Hence, the differences of vegetation-specific β between regions must be a result of physiological and structural adaptions to differing environmental heterogeneity rather than being caused by different positions on the MAP gradient. This conclusion is further supported by the fact that all vegetation types are found in approximately the same regions of MAP in both regions ( Figure 6). The distribution of the major vegetation types, however, indicates differences in the relative abundance between regions ( Figure 6). While WA is rather dominated by grass/crop type vegetation, SWA is strongly dominated by savanna and shrub type vegetation. Thus, rainfall variability may have an indirect effect on vegetation composition with higher CVP promoting shrub type vegetation whereas lower CVP favours grass type vegetation [69].

Conclusions
A novel approach of shifting linear regression models (SLRs) was used to derive satellite-based time series of vegetation functional response to rainfall in two African subtropical dryland regions.
The results presented show a distinct difference between vegetation functional response to rainfall (β) during dry and wet periods for the two semi-arid regions (South Western Africa (SWA) and Western Africa (WA)). The results suggest that interannual variability of annual rainfall (higher in SWA) has a strong influence on β (also higher in SWA). Thus, future studies addressing questions of water availability in ecosystems should not only focus on overall water availability (mm/year) but should as well account for changes in rainfall variability. As shown here, variability in both rainfall total amounts and season length may be expected to affect regional-scale dryland ecosystem dynamics and thus the impact of dry periods and droughts. The ecosystems of the Sahel region (WA) with their overall lower β during dry periods are expected to be vulnerable to future droughts. The vegetation productivity of SWA (MAP < 600 mm) is likely to be less susceptible to changes in water availability given its widespread, relatively high β values during dry periods. However, above this MAP, SWA may be considered to be prone to periods of below average rainfall. The region-specific implications of ecosystem-scale β found here are remarkable also since the results are consistent irrespective of vegetation type. Author Contributions: Gregor Ratzmann designed the concept of the study, analyzed the data and compiled the manuscript with support of all co-authors. Ute Gangkofner developed the phenology model.

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

Additional Information on the Study Regions
The South West African (SWA) study region is characterized by Southern hemisphere summer precipitation. The long-term mean rainy season in the region depending on the position lasts approximately from November to April with the southern parts experiencing significantly shorter seasons compared to the north [70]. The study region in the southern quarter is part of the semi-desert landscape Nama-Karoo with annual low precipitation (<150 mm MAP) [70]. Its vegetation consists mainly of perennial grasses which, however, have been widely replaced by poorly palatable annuals and dwarf shrubs [71]. Virtually the entire east of the study area comprises parts of the western branches of the Kalahari semi-desert. Most of the Kalahari vegetation consists of a variety of tree-grass savanna landscapes with the dominating tree species being of the Acacia genus. West of the Kalahari, former savanna landscapes have been widely replaced by shrubland with most of the shrub species being of the Acacia and Dichrostachys genus [72]. The westernmost part of the study region comprises parts of the highly arid and only sparsely vegetated Namib Desert. Vegetation in this region is mostly characterized of enduring grasses [73].
The West African (WA) study region is characterized by Northern hemisphere summer precipitation regime. The long-term mean duration of the rainy season is from approximately July to October [33]. The WA region recurrently is faced with periods of droughts and famines and is regularly identified as area being highly affected by land degradation [74]. It has suffered severe droughts during the 1970s and 1980s following dramatically below long-term average precipitation [75]. The consequence was a wide spread loss in available forage which triggered a decline of livestock populations of up to 40% [30] and severe famines [76]. In the northern part the WA study region can be characterized as Saharan-Sahelian transition zone with the regional plant community being dominated by perennial grasses [33]. Further south, the Sahelian zone comprises a savanna type vegetation consisting of scattered trees of the Acacia genus and an annual grass understory [77]. The Sudano-Sahelian zone connects the Sahelian zone in the north with more humid regions in the south. It mainly consists of a savanna which is dominated by species of the Combretacae family in the tree layer and by annual grasses in the herbaceous story [33].
In both study regions, recurring wild fires are an important contributor to vegetation dynamics [78] and can ultimately affect the functional composition particularly of savanna type ecosystems [31].

Phenological Parameterization Model
The algorithm of the phenological model [46] used in this study is designed to handle vegetation index data in a bi-weekly temporal resolution. In a first step, it ingests a set of 3 years of data with the year of interest being the centermost. This window is subsequently shifted along the temporal axis as soon as per time step all processing steps have been accomplished. Therefore, the following descriptions are generic and applicable to any given time step. In order to avoid false alarms for the calculation of the start of season, negative outliers are linearly interpolated; considered a negative outlier is a value with a difference of 5% to the preceding and following bi-weekly value relative to the value calculated as max−min of the respective two-year period. The start of the season (SOS) is found by means of an ancillary variable derived as moving sum of the differences between each value of the time series and its predecessor. The moving sum is built over each six consecutive differences, and can be understood as the average change over time during each 3-month window. In cases where no SOS can be found the average SOS of the other years is taken. SOS in general refers to the start of the so called "vegetation year". A vegetation year is composed of commonly one (or two and more in regions with several rainy seasons) vegetation peak(s) (cyclic fraction, CFR) followed by the dry season(s) (see Figure B1). To determine the onset and the end of the CFR of any given year, a baseline is derived, which separates the dry season values from the CFR. Values above this baseline are part of the CFR, values below belong to the dry season. The baseline is derived from the average amplitude between the dry season values and the preceding and respectively subsequent vegetation peak. The CFR is consequently determined as sum of the distances of all time-series values above the baseline to the latter. Thus, the CFR constitutes the integral of the vegetation season values above the (local) "dry season greenness" level expressed by the baseline.
Remote Sens. 2016, 8, 1026 13 of 19 [75]. The consequence was a wide spread loss in available forage which triggered a decline of livestock populations of up to 40% [30] and severe famines [76]. In the northern part the WA study region can be characterized as Saharan-Sahelian transition zone with the regional plant community being dominated by perennial grasses [33]. Further south, the Sahelian zone comprises a savanna type vegetation consisting of scattered trees of the Acacia genus and an annual grass understory [77]. The Sudano-Sahelian zone connects the Sahelian zone in the north with more humid regions in the south. It mainly consists of a savanna which is dominated by species of the Combretacae family in the tree layer and by annual grasses in the herbaceous story [33].
In both study regions, recurring wild fires are an important contributor to vegetation dynamics [78] and can ultimately affect the functional composition particularly of savanna type ecosystems [31].

Phenological Parameterization Model
The algorithm of the phenological model [46] used in this study is designed to handle vegetation index data in a bi-weekly temporal resolution. In a first step, it ingests a set of 3 years of data with the year of interest being the centermost. This window is subsequently shifted along the temporal axis as soon as per time step all processing steps have been accomplished. Therefore, the following descriptions are generic and applicable to any given time step. In order to avoid false alarms for the calculation of the start of season, negative outliers are linearly interpolated; considered a negative outlier is a value with a difference of 5% to the preceding and following bi-weekly value relative to the value calculated as max−min of the respective two-year period. The start of the season (SOS) is found by means of an ancillary variable derived as moving sum of the differences between each value of the time series and its predecessor. The moving sum is built over each six consecutive differences, and can be understood as the average change over time during each 3-month window. In cases where no SOS can be found the average SOS of the other years is taken. SOS in general refers to the start of the so called "vegetation year". A vegetation year is composed of commonly one (or two and more in regions with several rainy seasons) vegetation peak(s) (cyclic fraction, CFR) followed by the dry season(s) (see Figure B1). To determine the onset and the end of the CFR of any given year, a baseline is derived, which separates the dry season values from the CFR. Values above this baseline are part of the CFR, values below belong to the dry season. The baseline is derived from the average amplitude between the dry season values and the preceding and respectively subsequent vegetation peak. The CFR is consequently determined as sum of the distances of all time-series values above the baseline to the latter. Thus, the CFR constitutes the integral of the vegetation season values above the (local) "dry season greenness" level expressed by the baseline. Figure B1. Overview of the principle and terminology used in the phenological parameterization model to derive vegetation productivity proxies (cyclic fraction). One time step represents 14 days. Figure B1. Overview of the principle and terminology used in the phenological parameterization model to derive vegetation productivity proxies (cyclic fraction). One time step represents 14 days.

Shifting Linear Regression Model (SLR)
The SLR is computed on a per-pixel basis. Thereby, the linear regression model is calculated over a given set of years denoted by the window length, L, starting with the first set, and subsequently shifted along a temporal axis in one year steps until the end of the time series. Hence, for a given window length, n + 1 − L models are computed, where n denotes the number of years in the time series. For example, assuming a length of the time series of n = 28 and a window of L = 11, sets of 28 + 1 − 11 = 18 models and linear slope coefficients are derived.

Shifting Linear Regression Model (SLR)
The SLR is computed on a per-pixel basis. Thereby, the linear regression model is calculated over a given set of years denoted by the window length, L, starting with the first set, and subsequently shifted along a temporal axis in one year steps until the end of the time series. Hence, for a given window length, n + 1 − L models are computed, where n denotes the number of years in the time series. For example, assuming a length of the time series of n = 28 and a window of L = 11, sets of 28 + 1 − 11 = 18 models and linear slope coefficients are derived. Figure C1. Schematic drawing of the SLR approach. L denotes the window length.
The principle of the SLR is depicted by Figure C1. The first model in the time series in this schematic is Model 1 ingesting data from year 1 to year L. The model is subsequently shifted by one year and becomes Model 2 etc. Per model, the linear slope coefficient (β) is extracted. Slope coefficients are found based on ordinary least squares (OLS) techniques. The procedure is similar to moving correlation functions [79]. However, the SLR derives quantitative slopes rather than correlation coefficients (r-values).

Complementary SLR Analysis
To determine a potential effect of rainfall of previous years, we conducted the entire analysis using the SLR framework but with two additional model setups. (1) We extended the original model by rainfall of the previous year resulting in a multiple linear regression model with two explanatory variables (termed "Multiple"). (2) We extended the latter model by incorporating an interaction term between rainfall of the concurrent year and rainfall of the previous year (termed "Interactive"). To determine which model supersedes the other models over time and in space we computed time series of Akaike's Information Criteria (AIC) for each of the three model types per pixel. We then determined the model with the lowest AIC per pixel and time step. In a final step we derived the mode of the model with the lowest AICs. Maps of this mode are depicted by Figures D1 and D2. As is apparent from both figures as well as from the relative contribution of the model mode to the total pixel count ( Figure D3) it is by far the model with only one variable (rainfall of concurrent year, termed "Single" in Figure D1-D3) which supersedes both other models spatially and temporally. The principle of the SLR is depicted by Figure C1. The first model in the time series in this schematic is Model 1 ingesting data from year 1 to year L. The model is subsequently shifted by one year and becomes Model 2 etc. Per model, the linear slope coefficient (β) is extracted. Slope coefficients are found based on ordinary least squares (OLS) techniques. The procedure is similar to moving correlation functions [79]. However, the SLR derives quantitative slopes rather than correlation coefficients (r-values).

Complementary SLR Analysis
To determine a potential effect of rainfall of previous years, we conducted the entire analysis using the SLR framework but with two additional model setups. (1) We extended the original model by rainfall of the previous year resulting in a multiple linear regression model with two explanatory variables (termed "Multiple"). (2) We extended the latter model by incorporating an interaction term between rainfall of the concurrent year and rainfall of the previous year (termed "Interactive"). To determine which model supersedes the other models over time and in space we computed time series of Akaike's Information Criteria (AIC) for each of the three model types per pixel. We then determined the model with the lowest AIC per pixel and time step. In a final step we derived the mode of the model with the lowest AICs. Maps of this mode are depicted by Figures D1 and D2. As is apparent from both figures as well as from the relative contribution of the model mode to the total pixel count ( Figure D3) it is by far the model with only one variable (rainfall of concurrent year, termed "Single" in Figures D1-D3) which supersedes both other models spatially and temporally.