Drought-induced regime shift and resilience of a Sahelian ecohydrosystem

The Sahel (a semi-arid fringe south of the Sahara) experienced a long and prolonged drought from the 1970s to the mid-1990s, with a few extremely severe episodes that strongly affected ecosystems and societies. Long-term observations showed that surface runoff increased during this period, despite the rainfall deficit. This paradox stems from the soil degradation that was induced by various factors, either directly linked to the drought (impact on vegetation cover), or, in places, to human practices (land clearing and cropping). Surface runoff is still increasing throughout the region, suggesting that Sahelian ecohydrosystems may have shifted to a new hydrological regime. In order to explore this issue, we have developed a simple system dynamics model incorporating vegetation–hydrology interactions and representing in a lumped way the first order processes occurring at the hillslope scale and the annual timestep. Long term observations on a pilot site in northern Mali were used to constrain the model and define an ensemble of plausible simulations. The model successfully reproduced the vegetation collapse and the runoff increase observed over the last 60 years. Our results confirmed that the system presents two alternative states and that during the drought it shifted from a high-vegetation/low-runoff regime to the alternative low-vegetation/high-runoff one, where it has remained trapped until now. We showed that the mean annual rainfall deficit was sufficient to explain the shift. According to the model, vegetation recovery and runoff reduction are possible in this system, but the conditions in which they could occur remain uncertain as the model was only constrained by observations over the collapse trajectory. The study shows that the system is also sensitive to the interannual and decadal variability of rainfall, and that larger variability leads to higher runoff. Both mean rainfall and rainfall variability may increase in central Sahel under climate change, leading to antagonist effects on the system, which makes its resilience uncertain.


Introduction
Due to interactions between changing external forcing factors (e.g. rainfall) and internal dynamics, an environmental system (e.g. a forest) may shift into new states where it stays even if the forcing factors return to their initial value. In such a system, two alternative stable states may exist towards which the system is attracted (e.g. forest versus savannah, Staver et al 2011). For certain conditions, both equilibria can coexist and the system can shift abruptly from one to the other when critical thresholds (or tipping points) are passed (review in Scheffer et al 2001). This kind of behaviour has been widely studied in ecological systems (Scheffer and Carpenter 2003, Rietkerk et al 2004, Hirota et al 2011, van Nes et al 2014, Yin et al 2014, among others. Attention has been also paid to regime shifts in human societies driven by environmental factors (Di Baldassarre et al 2013, Broderstad and Eythorsson 2014, Sivapalan 2015, Downey et al 2016, Kuil et al 2016. Global change might increase the potential for the tipping of some key elements of the Earth's system (Lenton et al 2008, Barnosky et al 2012, Steffen et al 2018. Whilst the existence of tipping points at the plenatery scale (Brook et al 2013) and their propagation across scales  is debated, examples exist at the regional scale, such as the the collapse of the 'green Sahara' and the shrinking of Lake Chad around 6000 BP as a response to a gradual change in the Earth's orbital parameters (Claussen et al 1999).
Marked and persistent ecohydrological changes have been observed in the Sahel (a semi-arid belt south of the Sahara) over the past 60 years. This region experienced a severe and prolonged drought from the 1970s onwards (Lebel and Ali 2009), with extremely severe episodes that have strongly and durably affected ecosystems and societies. Over most of the Sahel, greening has been observed by satellite imagery since the 1980s (Fensholt et al 2012. Concurrently, long-term observations have shown that the outflow of Sahelian watersheds has increased from the 1950s onwards (Mahé and Paturel 2009, Gardelle et al 2010, Descroix et al 2012, Gal et al 2016, including during the drought (the so-called 'Sahelian paradox'; Favreau et al 2009). Field studies showed that these apparently conflicting changes resulted from smaller scale, soildependent, processes. Vegetation (mainly herbaceous) recovered on deep sandy soils, driven by post-drought rainfall increase, whereas it decayed on shallow soils, despite the annual rainfall trend, causing erosion and runoff , Trichon et al 2018, Gal et al 2017. In some highly populated areas (e.g. South West Niger), this so-called regreening was not observed due to the impact of human activities (land clearing, crops which are less green than rangelands, fertility losses; Favreau et al (2009), .
As in other semi-arid regions, the Sahelian hydrology is strongly dependent on surface conditions (e.g. Casenave andValentin 1992, Peugeot et al 1997). At the elementary scale, heavy rains favour soil crusting all the more since the vegetation cover is sparse; a loss (gain) of vegetation cover favours (prevents) erosion and fertility losses, which in turn prevents (favours) vegetation expansion. These small-scale mechanisms combine into a positive feedback loop which has been shown to be involved in larger scale desertification dynamics (D'Odorico et al 2013, Wilcox et al 2017, Saco et al 2018. If the feedback is strong enough, catastrophic transitions between the alternative low and high runoff states (and the corresponding high and low vegetation ones) can occur (Mayor et al 2013, Kefi et al 2016. The observed ecohydrological changes on the scale of the whole Sahel show some similarities with these small-scale mechanisms. They suggest that the 1970-1994 drought could have triggered a shift to a high runoff state at the elementary scale, which would have resulted in the new hydrological regime observed at the regional scale. The outflow of Sahelian watersheds has continued to increase in recent years , which is associated with an intensification of rainfall (Taylor et al 2017 and floods (Wilcox et al 2018, Tamagnone et al 2019, and these trends are expected to persist with climate change (Giannini et al 2013, Monerie et al 2017, Martin 2018. Land conversion to croplands is also expected to continue over the coming decades as a result of one of the highest population growths worldwide (UN 2017). How these changes will impact future water resources and living conditions is still largely unknown. In particular, the possibilities for societies to adapt to adverse changes will be radically different, whether the ecohydrological response to changing forces is fast or slow  and whether it is driven by gradual versus critical transitions between states (Scheffer et al 2001).
Exploring whether interactions between the hydrological cycle and the vegetation dynamics might explain the recent changes in runoff conditions in the Sahel and involve some potential for tipping is therefore a scientific challenge with a high societal impact. While the existence of bistability and alternative states in semi-arid ecosystems has been investigated either from a conceptual (e.g. Vetter 2005, Turnbull et al 2008, observational (e.g. Hirota et al 2011, King et al 2012, Holmgren et al 2013, and modelling (e.g. Kefi et al 2010, van Nes et al 2014, Cueto-Felgueroso et al 2015 perspective, very few studies combined dynamic modelling and observations (e.g. Yin et al 2014). We have thus developed a model designed to capture the first-order dynamical interactions (including feedback loops) between the hydrological cycle and the vegetation dynamics in the Sahel, inspired by the work of Scheffer and Carpenter (2003) and van Nes et al (2014). The model was tested against data from an experimental site in Mali, with the aim of: (i) demonstrating that it is able to reproduce the observations; (ii) assessing whether the observed evolution might correspond to regime shifts between alternative stable states, triggered by the drought, and (iii) exploring the possible future evolution of the system under climate change and its resilience to rainfall fluctuations.

Study site
The site used in this study is located in northern Mali (Sahel, West Africa), near a place named Ortonde (15.15°N; 1.56°W), located 20 km from Hombori village. The climate of the region is semi-arid. The yearly rainfall averages 375 mm over 1936-2015, with rains only occurring between June and September.
The study site is a typical Sahelian banded vegetation pattern also called tiger bush, made of elongated thickets of vegetation perpendicular to the hillslope alternating with runoff-prone bare soil areas (Hiernaux and Gérard 1999). Tiger bush is a natural system encountered in many semi-arid areas over the world (d 'Herbés et al 2001, Deblauwe et al 2008. In this type of ecosystem, the vegetation bands do not produce runoff while intercepting the runoff generated on the upslope bare soil areas. Capture and infiltration of this extra amount of water are essential for the persistence of the vegetation (Ludwig and Tongway 1995, Valentin et al 1999, Saco et al 2007. Due to low population density, wood cutting and grazing by livestock can be considered negligible on this site, which is also not affected by fires. More details can be found in Trichon et al (2018).
Rainfall was measured at the Hombori meteorological station by the Malian Meteorological Service and the AMMA-CATCH Observatory (Galle et al 2018). The drought period  resulted in a 100mm deficit (-24%) as compared to the predrought  one (Le Barbe et al 2002, Lebel and Ali 2009, Trichon et al 2018; the annual rainfall slightly increased from 1995 but remained lower than before the drought (table 1). The inter-annual variability, assessed by the standard deviation of the annual rainfall, remained fairly constant over the different periods (about 95-110 mm y -1 ). The lag-1 autocorrelation over the full period, equal to 0.45, was used as a measure of the temporal structure of the annual rainfall series.
The evolution of the land cover on the site was assessed by Trichon et al (2018), from field surveys over a transect (since 1985), and aerial and satellite images . This dataset describes the evolution of the perennial woody cover, the areas where seasonal herbaceous vegetation can develop, and the runoff-prone bare soil areas (supplementary section S1 available online at stacks.iop.org/ERL/14/105005/mmedia). Those variables evolve over years, with a negligible intra-seasonal cycle. The well-developed woody vegetation bands have progressively shrunk since 1955, especially since the drought. In 1955, the structure of vegetated/bare areas and the near-absence of rills suggested that the thickets completely infiltrated runoff from bare soil patches, and that hillslope-scale runoff was very low, as expected for a healthy tiger bush. A network of rills and gullies progressively developed, leading to a drainage of about 11% and 26% of the site in 1985 and 2015, respectively (Trichon et al 2018). These changes show an increase in the interconnection of bare soil patches and of hillslope-scale runoff, although quantitative measurements are lacking. This site is typical of the areas which did not recover after the drought (Gardelle et al 2010).

Model development
The model that was developed describes the main interactions between the land cover and the hydrological cycle for this type of system. It works at a yearly time step and hillslope scale. It is a lumped model, with no dimension in space. Only the main vegetation/ hydrology interactions are represented and detailed processes (nutrients/carbon cycle, small-scale erosion, crusting, sediment fluxes...) and their spatial and intra-seasonal variability are implicitly embedded in the lumped model parameters. The two state variables of the model are: (i) the fraction of the surface covered by perennial woody vegetation W (−) and (ii) the runoff-prone bare soil fraction B (−). The remaining fraction H corresponds to areas where annual herbaceous can develop, hence : Thus W and H represent a fraction of the total surface (and not a leaf area); they vary annually and do not have a seasonal cycle. The total W, B and H fractions are simulated without specifying their spatial repartition within the hillslope. The dynamics of each state variable is represented by an ordinary differential equation in time. The annual precipitation P (mm) is the only external force driving the system. In the following, upper case symbols refer to time-varying variables, and lower case symbols to constant parameters.

Water balance
The annual precipitation P (mm) is partitioned into runoff R (mm), which is the net runoff exported out of the system, and I=P−R (mm), which represents the total water available over the year. It is not an actual water storage but an indicator of the soil water amount available for the vegetation during the year; it aggregates evapotranspiration and infiltration.
The hillslope scale runoff coefficient Ke represents the fraction of rainfall P converted to runoff R, which is written: In the field, the runoff produced locally on bare soil areas can reinfiltrate in herbaceous and woody  (2008), the spatial structure of bare soil patches was parametrized by the flowlength metrics (FL), defined as the mean length of uninterrupted runoff paths along the slope. Rodríguez et al (2018) proposed a formulation of FL adapted to the representation of the system on a grid with a known pixel size. Based on their work, we developed a modified formulation of FL which is better adapted to our lumped approach, which writes: where l(−) is a constant parameter representing the effect of the size and spatial organization of patches. FL * ranges from 0 to 1, and represents the hydrological connectivity along the hillslope (Okin et al 2015). The hillslope-scale runoff coefficient Ke (equation (2)) was assumed to be proportional to the flowlength: where ke B (−) is the small scale runoff coefficient of bare soil areas.

Woody vegetation dynamics
The meaning of the parameters of this equation are described in table 2, and details on each term are given in supplementary section S2.

Bare soil dynamics
The bare soil fraction B increases at the expense of herbaceous areas H due to crusting and erosion, and decreases due to vegetation expansion. This dynamics is represented by three main processes: (i) increase due to crusting by rainfall (e.g. small-scale splash effect); (ii) increase due to runoff (crusting, erosion); and (iii) decrease due to recolonization by perennial vegetation. For simplicity, we assumed that rainfall and runoff effects on B were proportional to the yearly amounts P and R, and we did not parametrize the effects of rainfall intensities at the event timescale. The dynamics of B writes: where α p and α r are parameters; the rightmost term is identical to equation (S2) (recolonization) and appears in equation (6). Finally, the model is composed of the differential equations (6) and (7), together with equations (1) to (5).

Model implementation
The model includes 10 parameters (table 2). The range of l was constrained by fitting equation (4) (table 2).
From these ranges, we built an ensemble of 10 6 sets of parameters using a Latin hypercube sampling method with uniform parameter distributions, with the constraint r r r g . For each parameter set, the model was initialized with W and B observed in 1955 (table S1), and run from 1955 to 2015 using the series of annual rainfall from Hombori. The goodness-of-fit between the simulated and observed time series of W Table 2. Variable definitions and ranges of possible values. ( * ) Assuming recolonization of bare areas is slower than growth in herbaceous areas, we imposed r r <r g . ( ** ) From Peugeot et al (1997) and Galle et al (1999 (2018), the calibration was performed without accounting for the uncertainties on W and B, and the impact of this assumption was evaluated (supplementary section S8). The model was implemented under the R environment (R/3.4.3), using the 'deSolve' package (Soetaert et al 2010) for solving differential equations and 'SAFER' for latin hypercube sampling (Gollini et al 2015. The code was parallelized using the 'parallel' package (R Core Team 2018) and run on a 28-core node on the MESO@LR-Montpellier University computing facility. The simulation of 10 6 sets of parameters over 60 years took ≈2 h on this platform.

Exploration of the system dynamics
We used the 1000-member simulation ensemble to analyze the properties of the modeled system. For each member, we first determined the equilibrium states for W, B and Ke from the final state reached by the system when forced by long series of constant rainfall (1000 years) with P ranging from 150 to 600 mm (supplementary section S5). As the state towards which the system is attracted changes with changing rainfall forcing, the equilibrium states are virtually never reached in the real world. Then, in order to evaluate how rainfall variability determines the system final state  (2009), we used a red-noise model in which the autocorrelation decreases exponentially as a function of the lag. For a given initial condition and parameter set, the model was forced over 2000 years with the series of synthetic rainfall, exploring a range of mean P, SD and Corr. The response of the system was assessed by the statistical distribution of W, B, and FL * over the last 1000 years, to get rid of the influence of the initial condition.

Model calibration
The 1000-member simulations ensemble reproduced very well the observed evolution of W and B (figure 1) with a spread lower than the uncertainty on the observations (RMSE<0.002 4). The flowlength FL * (proportional to the runoff, equation (5)) was derived from B and depends on l (equation (4)). The lack of quantitative observations to constrain FL * (hence Ke), the nonlinear dependence between B and FL * (equation (4)) and the additional degree of freedom brought by l explains the large spread of FL * and its increase over time. Despite the spread, the increase of FL * is consistent with the development of rills and gullies observed on the site (Trichon et al 2018).

Equilibrium states
For all members, the system equilibrates to a degraded state (low W, high B and FL * ) for P<350 mm, and to a vegetated state for P>580 mm (figure 2). The fraction of woody cover W at equilibrium increases with the rainfall amount, and B and FL * decrease accordingly. The herbaceous cover fraction H at equilibrium (not shown) is about 0 for the degraded state and ranges from 0.03 to 0.1 for the vegetated one.
For 37% of the members, a single equilibrium state exists (monostability) for each rainfall amount, and the equilibrium value gradually changes with precipitation (figures 2(a)-(c)). For the remaining members (63%), the system is bistable over a range of P values where two equilibria coexist (figures 2(d)-(f)). The bounds of the bistability range correspond to the critical thresholds (tipping points) of the system, which delimit the range of P where the equilibrium state differs depending on whether the system is declining or regreening. The equilibrium curves (figure 2) and the precipitation critical thresholds for mono-and bi-stable members (table 3) are similar along the decline branches, but very different along the regreening ones. Since observations were only available along the decline trajectory of the system, the regreening branches were poorly constrained. Consequently, the value of the critical regreening thresholds and the proportion of bistable members remain uncertain.  Table 3. Annual precipitation thresholds P (mm) for quantiles 10, 50, and 90% of the ensemble members with and without bistability. For monostable cases, the rainfall threshold was defined as the value of P above which a runoff exists (R>0.1). For all members, the degraded state was the only one possible for drought (P = 325 mm) and postdrought (P = 357 mm) mean precipitation (figure 2). For the pre-drought mean precipitation (P = 421 mm), an equilibrium vegetated state with low runoff exists for 73% of the members, of which an alternative degraded equilibrium exists for about one half of them. For the remaining 27%, the attracting equilibrium is the degraded one. All these results were obtained from an ensemble of parameter sets selected without taking into account the uncertainties on the observations, which we could not quantify (section 2.3). We have shown with a sensitivity analysis (supplementary section S8) that accounting for these uncertainties moderately changes the values of the critical thresholds and the proportion of mono-and bi-stable members. Nevertheless, the uncertainties do not affect the qualitative behavior described above and do not call into question the existence of the two alternative stable states, nor the other conclusions of the study.
For most of the members, the system equilibrium is the vegetated/low runoff state before the drought, and the alternative low vegetation/high runoff state since then. However, equilibrium states do not provide information on actual trajectories because the way in which the system converges to these states depends on its evolution rate and on the rainfall variability.

Actual trajectories
When the model was initialized with the 1955 state and forced with the observed rainfall, all the simulations crossed several times the equilibrium curves during the pre-drought period, oscillating around a vegetated state with W≈0.15, B≈0.80 and FL * ≈0.40 ( figure 3). Due to a few very dry years (around 1984) in the drought period, it was pulled far away from this equilibrium (to the left of the equilibrium curve), while W (resp. B, FL * ) was decreasing (resp. increasing). After the drought, the system remained in a degraded state, despite the occurrence of some wet years (P>450 mm) corresponding to a vegetated equilibrium. Although the range of rainfall variability in every subperiod encompassed P values associated with vegetated and degraded equilibrium, the system did not adjust to them (figure 3), which means that its adjustment timescale is larger than the timescale of rainfall variability.
These results suggest that the drop in the mean annual rainfall between 1970 and 1994 was a major driver of the observed decline of the Ortonde tiger bush. However, the inter-annual variability of rainfall may also have contributed to the destabilization of the system.

Sensitivity to rainfall variability
Except if the system is extremely slow (insensitive to forcing variability) or extremely reactive (it adjusts immediately to the equilibrium state corresponding to the forcing), the variability of the forcing and its temporal structure can drive oscillations of the system between states (Bathiany et al 2018) and trigger critical transitions (van der Bolt et al 2018). The sensitivity of the system to rainfall variability was assessed using synthetic annual rainfall series with prescribed mean, standard deviation and lag-1 autocorrelation (section 2.4). This sensitivity is illustrated in figure 4 for the two members corresponding to the median member (supplementary section S6) of the monostable and the bistable simulation subsets, initialized with the 1955 condition.
In the monostable case (figures 4(a1)-(a3)), the distribution of W widens (blue colors) and the mode (most frequent value, crosses) decreases when SD (figure 4(a1)) and Corr (figure 4(a2)) increase: the system gradually shifts from a vegetated to a degraded state when the rainfall variability increases and/or when its temporal structure strengthens ( figure 4(a3)). The increased dispersion of W for high values of SD and Corr means that the system can reach a high vegetation states at some time, but it is attracted back to a low W state.
In the bistable case (figures 4(b1)-(b3)), the same type of behaviour is obtained except that beyond a SD threshold (SD≈50 mm for Corr=0.45 in figure 4(b1)), the distribution of W becomes bimodal (W≈0.15 and W≈0.01). For SD>80 mm, the distribution abruptly shrinks and the mode shifts towards W≈0 ( figure 4(b1)). The shift to the degraded state occurs for a lower variability (SD) when Corr increases ( figure 4(b3)). It means that due to longer dry and wet periods (high Corr), less deviation from the mean rainfall is required to attract the system towards the degraded state (van der Bolt et al 2018). In all cases, the system forced with variable precipitation oscillates around a W (resp. B and FL * ) value which is lower (resp. higher) than the equilibrium reached with a constant rainfall (SD=0 on figures 4(a1), (b1)). Increased autocorrelation, which implies longer dry and wet periods, moves the system further from equilibrium during precipitation anomaly periods, which amplifies the effect of the inter-annual variability (figures 4(a3) and 4(b3)). This result can be interpreted by asymmetric declining versus regreening processes: the system degrades more during dry years than it recovers during wet years, which pushes it to a more degraded state when rainfall variability increases.
These results confirm that no vegetated state was possible in both drought and post-drought periods, and that the system remained trapped in the basin of attraction of the degraded state until now. They also show that the reduction of the mean annual rainfall during the drought period is sufficient to explain the observed decline of this tiger bush.

System resilience
With climate change, the inter-annual variability as well as the length of dry and wet periods may increase in central Sahel in the next decades (Monerie et al 2017, Martin 2018. In this region, there is no consensus on the sign of annual precipitation change between global circulation models (Christensen et al 2014). Recent studies over the region showed that the annual rainfall may increase in central Sahel (Giannini et al 2013, Park et al 2016, Monerie et al 2017, although this scenario is still uncertain. This trend could bring favourable conditions for the recovery of the Ortonde tiger bush. The resilience of the system, defined as its potential to recover a vegetated state with low runoff, was assessed with the same kind of analysis than in the previous section. The model was initialized with the current state (2014), and forced with synthetic rainfall series with increasing mean rainfall: P=425 mm (≈pre-drought), P=450 mm (+6%), and P=475 mm (+12%), and varying SD and Corr.
For the median monostable member, the system response to variable rainfall with pre-drought properties is similar to that obtained with the 1955 initial conditions (figures 5(a) and 4(a)), and the current state is stable (W≈0.025). For higher annual precipitation (figures 5(c), (e)), the system recovers to W values that are larger if P is large, except if both SD and Corr are high (figures 5(c3) and (e3)).
For the median bistable member, the system remains on the degraded state for pre-drought mean precipitation, with little fluctuations ( figure 5(b)), whereas a vegetated state exists with the 1955 initial condition ( figure 4(b)). This difference underlines the effect of the bi-stability: with the 2014 initial condition, the system lies in the basin of attraction of the degraded state, whereas it is attracted by the alternative vegetated state with the 1955 condition. For P=450 mm, the system is bi-stable for 50 mm<SD<130 mm (Corr=0.45). Bi-stability almost disappears for P=475 mm, and the system recovers except when SD110 mm and Corr0.5 (figure 5(f)). Along with vegetation recovery, the model simulates a reduction of the runoff capacity (supplementary section S7), associated with the reduction of bare soil areas. However, the rills and gullies that developed on the site will hardly disappear. These durable physical changes imply that the system cannot recover following the same pathway back. As such changes in the hydrological structure are not represented in our model, our results could overestimate the resilience of the system. For all members the shift between states occurs for higher mean precipitation when the rainfall variability (SD and/or Corr) increases, e.g. the iso-value W=0.05 moves to domains with higher SD and Corr when P increases (figures 5(a3)-(c3)-(e3), and (d3)-(f3)). It implies that under variable precipitation forcing, the regime shift will occur for higher precipitation values than assessed from constant rainfall (table 3 and figure 2). According to the model, a regreening of the system from its current state is possible. It requires the mean annual precipitation to be at least equal to its predrought value (P≈425 mm) which is sufficient if the rainfall variability is low and if the system is not bistable. If the system is bistable, a higher mean precipitation is required, at least equal to the regreening critical threshold for constant rainfall ( figure 2 and table 3). This value is not sufficient considering annual rainfall variability, as illustrated in figure 5(d3) for the median, bistable member: with the current rainfall variability (SD=105 mm, Corr=0.45), P=450 mm is no longer sufficient for regreening, whereas it is the regreening threshold in table 3.
Considering the future rainfall increase as a working assumption, we have estimated from Monerie et al (2017) that the mean annual precipitation may reach about 475 mm by the end of the 21st century in Northern Mali. Similarly, the rainfall variability may increase, with an estimated SD increment of about 50 mm. Thus, according to our simulations, the simultaneous increase in annual totals and interannual rainfall variability would lead to antagonistic effects on the resilience of the system and, therefore, an uncertain future. The red noise model used to assess the effect of rainfall variability poorly represents the strong decadal structure of Sahelian rainfall (Dieppois et al 2013), and the effect of the temporal structure of rainfall may be stronger than estimated. Furthermore, a larger intraseasonal variability is also anticipated in a warmer climate (Martin 2018), which could also affect the potential for resilience.
In addition to the intraseasonal rainfall variability and the hydrological effects of rills, other factors such as air temperature have been omitted in our simple model. We used rainfall as unique forcing because vegetation and runoff changes were more driven by rainfall than temperature changes over the past decades , 2009, Leauthaud et al 2017. Temperature increase is a robust trend globally (Christensen et al 2014), it and could play a major ecohydrological role in the future. Increased temperature-controlled evaporation may increase the water stress (Young et al 2017). However the effects of higher temperature on the physiology of tropical plants are largely unknown and deserve specific studies (Jones 1992, Cavaleri et al 2015. Although direct and indirect temperature effects are uncertain, it is probable that higher temperature will reduce rather than improve the future resilience of Sahelian ecosystems.

Conclusions and prospects
We developed a system dynamics model to represent the ecohydrological evolution of a Sahelian tiger bush over decades, constrained by observation data. The model represents the first order interactions (including retroactions) between hydrology and vegetation. Although based on a simple representation of the studied tiger bush, the model was able to reproduce the observed vegetation decline and the associated increase in surface runoff over the 1955-2015 period. Based on field work, these trends are robust and are iconic of the declining systems in the region. According to our simulations, the system was oscillating around a vegetated state before the 1970-94 drought, as a response to rainfall variability. Then it started to shift to a degraded (low vegetation/high runoff) state around which it currently remains, despite a slight rainfall recovery. We have confirmed that this decline corresponds to a shift between two stable states: a vegetated/low runoff state and a low vegetation/high runoff state. The drop in the mean annual precipitation associated with the drought was found to be sufficient to explain this regime shift. We could not conclude whether the shift occurred as a gradual response to the changing forcing or if it implied critical transitions (hence bistability), both modalities of change being plausible in our simulation ensemble.
We showed that increased variability in annual rainfall (amplitude of variation and temporal structure) pushes the system towards more degraded states. Rainfall increase in central Sahel is a possible trend anticipated in some climate change scenarios, but its effect on vegetation recovery may be offset by the concurrent increase expected in rainfall variability. According to our simulations, the Sahelian tiger bush is resilient. However, as the model was only constrained along the decline trajectory observed on a particular site, the critical rainfall conditions for resilience could not be evaluated precisely. Moreover, the simple model structure and some unaccounted processes (impacts of higher temperature, effect of gullies) might well lead to over-estimate the resilience of the system. Our study has yielded robust qualitative conclusions on the existence of alternative stable states in such systems and on their resilience potential, but due to the remaining uncertainties we could not make quantitative predictions.
The main conclusion of this study is that a drought-induced regime shift can explain the runoff increase observed at the hillslope-scale in some Sahelien landscapes. In a future work, the model will be applied to various sites representative of the diversity of Sahelian ecohydrosystems, including regreening ones . This will quantify the role of hillslope-scale regime shifts in the concurrent regreening and runoff increase trends observed at the regional scale.