Catastrophic Regime Shift in Water Reservoirs and São Paulo Water Supply Crisis
Abstract
The relation between rainfall and water accumulated in reservoirs comprises nonlinear feedbacks. Here we show that they may generate alternative equilibrium regimes, one of high watervolume, the other of low watervolume. Reservoirs can be seen as socioenvironmental systems at risk of regime shifts, characteristic of tipping point transitions. We analyze data from stored water, rainfall, and water inflow and outflow in the main reservoir serving the metropolitan area of São Paulo, Brazil, by means of indicators of critical regime shifts, and find a strong signal of a transition. We furthermore build a mathematical model that gives a mechanistic view of the dynamics and demonstrates that alternative stable states are an expected property of water reservoirs. We also build a stochastic version of this model that fits well to the data. These results highlight the broader aspect that reservoir management must account for their intrinsic bistability, and should benefit from dynamical systems theory. Our case study illustrates the catastrophic consequences of failing to do so.
1 Introduction
Complex socioecological and socioenvironmental systems often involve interactions between natural elements and human action in a nonlinear and adaptive way (Folke, 2006). Such systems have received much attention (Gordon et al. , 2008; Österblom et al. , 2013) recently and may display many of the common dynamical features present in natural systems. Here we will be interested in the existence of alternative stable states, representing different possible dynamical regimes, and transitions between them. These tippingpoint transitions have been well studied and are common in natural systems (Scheffer et al. , 2001). Desertification (Klausmeier, 1999; Meron et al. , 2004; Scheffer & Carpenter, 2003) and lake eutrophication (Scheffer & van Nes, 2007; Carpenter, 2005) are only the most notable ones among a plethora of cases. Characterization of a transition from timeseries data can be assessed using indexes related to either returntime to equilibria or variability near the tipping point (Dakos et al. , 2012). These techniques, together with models that encompass the main features of the system, may be extended to systems where the interaction with the human factor creates the possibility of a certain level of control or, at least, of action attempting to bring the system to a situation that is considered desirable. Our attention will be focused on the dynamics of water reservoirs, which has not been studied under this light, and on the characterization of alternative stable states and associated transitions. We show a case in which these dynamic properties had been neglected in the management of one of the largest water reservoirs in the world, the Cantareira system, which serves the Metropolitan Area of São Paulo (MASP).
MASP faces an unprecedented crisis of water supply. This exceptional situation can be readily seen from the plot of the volume of water stored the Cantareira system, the main reservoir serving the area, over the past twelve years (Fig.(1)). The volume of water decreased sharply from mid2013 and the operational capacity of the reservoir was depleted in July 2014. Since then water withdrawal is done by pumping of the socalled “strategic reserve” or “dead volume”. The São Paulo Water Company, SABESP, began to reduce withdrawals in January 2014, and by May 2015 the total outflow was 40% of the average values. Additionally, the last rainy season (October 2014 – March 2015) provided more rainfall compared to the previous two years. Despite that, only of 15% of total volume had been recovered, and the reservoir remains operationally exhausted. This situation seems typical of catastrophic regime shifts, driven by a tippingpoint transition. In this work, we will establish this regime transition on firm grounds and discuss its consequences. For sake of brevity, we will use the term reservoir to mean a system of interconnected reservoirs.
We propose that reservoir watervolume is subject to drastic regime changes due to the underlying bistability of the system. Bistability is the condition in which there are at least two possible stationary states of a system for the same set of parameters. In our case, the parameter is rainfall and the variable that may assume two possible alternative stable states is the volume of water stored in the reservoir. If the reservoir is at high levels the catchment has more stored water and thus more of the rainfall will flow into water the reservoir. On the other hand, for low water levels the most obvious change is that much of the water is absorbed by dry soil and thus a lower proportion of rainfall becomes stored water in the reservoir. Many other factors may be play and the exact relation between rainfall and the amount of water stored in the reservoir is an instance of the rainfallrunoff problem in hydrology (Beven, 2011). What matters for our discussion is that these feedback mechanisms are drivers of critical transitions.
Both states, highlevel and lowlevel regime, are resilient in ranges of rainfall and outflow parameters that may overlap, leading to hysteretic behavior: the paths to and from the lowlevel state are not the same. Resilience of the highlevel water volume regime may be lost, for instance, if the rainfall index falls below a certain critical value, or if the reservoir is overexploited. The system then transits to the lowlevel regime, in a tippingpoint transition entailing a regime shift. The dramatic side of this fact is that the lowlevel regime is also resilient, so either rainfall index has to increase beyond the mentioned critical value, or outflow has to be reduced drastically in order to produce a backwards transition to the highlevel regime. This state of affairs prompts us to characterize transitions by indexes related to the critical slowing down phenomenon (Dakos et al. , 2015). These indexes are usually connected to research of earlywarning signals of critical transitions. Discussions in recent literature have raised doubts about how early the signals are (Boettiger et al. , 2013), but this is not actually relevant in the present case, as we are concerned mainly with characterization and not anticipation of the transition.
In what follows we proceed by the following steps: we first inspect the data to establish the existence of a feedback mechanism; next we look for indicators of critical transitions in two steps, first by a qualitative assessment of the data, and second by demonstrating the existence of a peak in the variance characteristic of near tipping point regions. We then turn to a further step by proposing a mathematical model for the dynamics of the system, in two versions. First, a deterministic version where we can understand the dynamical behavior of reservoirs. Next, we propose an stochastic extension of the model, via a stochastic differential equation which can be fitted to actual data accurately and used for shortterm extrapolations.
2 Model Independent Results
The feedback mechanism
Given a certain amount of rainfall, how much water will effectively flow into the reservoir? This is connected to the rainfallrunoff problem, central in hydrology. The big picture is that the wetter the catchment the more of the rainfall is converted in runoff water (Beven, 2011), that in this case will flow to the reservoir. A complete assessment of this situation involves modeling the rainfallrunoff relation in an explicit spatial setting, considering soil heterogeneities, terrain and many other variables. However, in order to proceed in a simpler way, we postulate that the volume of water stored in the reservoir is a surrogate for the drainage basin condition. To verify that this is a sound assumption, we show in Fig.(2) the ratio of water inflow to the rainfall in the Cantareira reservoir as a function of the volume, yielding a clear trend: the higher the volume, the more efficiently the rainfall becomes water flowing into the reservoir.
Assessment of a critical transition
In Fig.(3) we show the data for the Cantareira reservoir in the rainfall vs. storedwater plane. In normal situations, we would expect to have nearly closed curves in the upper part of the plane. These curves come from the seasonality of rainfall, which induces similar oscillations in the water stored in the reservoir. However, what is seen in the plot is that from 2013 on, the system spins down to a new cycle on the lower part of the plane. This lower cycle represents, again, seasonal variations, but the volume now oscillates around much smaller values. In the inset of Fig.(3) we show what would be the result obtained from an ordinary differential equation displaying bistable behavior, in which the rainfall parameter oscillates in time and has a decreasing trend, inducing a regime change. Similarities are evident. The specific equations that generate Fig.(3) will be discussed in the Mathematical Models section.
The previous discussion gives a qualitative view of a transition, but we would like to take a step forward and provide a more quantitative assessment. In recent years, a toolbox for detecting earlywarning signals of critical transitions has been developed (Dakos et al. , 2012; Dakos, 2015). These are based on the critical slowing down phenomenon, implying that the lagautocorrelation function and variance indicators should increase near the transition. As discussed in Boettiger & Hastings (2012), a series of underpinnings exist concerning the reliability of these indicators, which may fail to anticipate the occurrence of a transition. However, in this work we are interested in characterizing the regime shift rather than anticipating it. In order to do so, we evaluated evidence of a critical transition with a diffusiondriftjump model (ddj). If the time series is long enough and has high frequency of observations, this model gives an accurate description of the underlying process that generates the data (Carpenter & Brock, 2011). Variations along time are described as a combination of instantaneous changes by deterministic trends (the drift), instantaneous variation (diffusion) and occasional uncorrelated changes (jumps). The model also allows to estimate the variance conditioned to these components. In theory, this conditional variance goes to infinity at bifurcation points of the underlying dynamics (Dakos et al. , 2012). We found a very clear peak in the conditional variance near the beginning of 2014 (Fig.(4)). This gives the first characterization of an indicator of critical regime shifts for water reservoirs. Runningwindow statistics like variance and autocorrelation at lag one (Dakos et al. , 2012) followed the same trend.
As discussed in Dakos et al. (2015), critical transitions are better assessed when, together with a statistical indicator, like the conditional variance used in Fig.(4), a mechanism leading to bistability is evident. The feedback discussed above is the basic mechanism we hypothesize. A further step, therefore, is to build a mathematical model that takes into account the feedback and may be fitted to the data, and show that it exhibits bistability. This is what we present in the next section.
3 Mathematical models for reservoir dynamics
Deterministic model
The change in time of storedwater volume in a reservoir is determined by inflow and outflow rates. Outflow rates are a consequence of management decisions. Inflow comes ultimately from rainfall through a set of hydrological processes. This leads to an equation of the form:
(1) 
where is the volume of water stored in the reservoir, is the outflow rate and is the inflow rate, with a measure of the rainfall. We now discuss these terms separately.
Outflow.
We assume that withdrawal of water is managed and saturates to a maximum, , when the stored volume is high, but is reduced as water becomes scarce (which could be done both to prevent complete disaster or due to the difficulty in pumping water from low levels). Such assumption is realistic for the Cantareira system: while stored volume was above 33% water withdrawals oscillated around 36 , according to the maximum quota set by the Brazilian Water Agency, ANA (Fig5). As the volume dropped outflow decreased because ANA pushed the main reservoir operator SABESP to cut down withdrawals. To model this situation we use the following simple form:
(2) 
where is the volume of water in which the withdrawal is reduced to . It is a measure of how cautious the water management is.
Inflow.
As already implied in the notation , inflow depends not only on rainfall, but also on volume. As discussed above, we propose that inflow is higher for higher volumes. A simple assumption is thus:
(3) 
where , and are phenomenological parameters. We also introduced a maximum volume for the reservoir, . Once it is attained, no increase in stored water is possible and excess inflow is released downstream through a spillway.
Putting everything together, we come to the following differential equation:
(4) 
if . If , the second term on the righthand side is zero. Equation (4) is not meant to be a full model of all hydrological processes, but rather a description of the basic dynamical factors affecting the volume of water stored in a reservoir. In spite of its simplicity, a stochastic version of it fits very well to observed data, as will be seen below.
If the rainfall is constant over time, Equation (4) has either one or two stable fixed points, depending on the value of . Fig.(6) shows the general situation, assuming specific values for the parameters. The existence of a bistability region is clear.
The existence of the bistability region depends essentially on the fact that inflow is dependent on the volume. The form of the outflow term is not essential: all that is required is that it goes to zero as goes to and that it saturates at some value. The existence of a maximum volume is, from the mathematical viewpoint, a form of avoiding arbitrarily large volumes.
Obviously, in real situations is not constant. It has clear seasonal variations and may be subject to longterm trends. Solutions, in this case, will not tend to a constant value. Bistability will be reflected in the existence of different regimes, corresponding to oscillations around the fixed points of the timeindependent problem, as long as the seasonal variations and trends are not too large. In Fig.(7) an example is shown, with . We also plot the solution as a curve in the xplane (inset of Fig.(3)), showing the distinctive feature of oscillations in the upper part, spinning down on a short time scale to the lower part. This pattern is present also in the observed relationship between rainfall and stored volume in the reservoir (main panel of Fig.(3).
Stochastic model
In order to compare theory to data we need to take a further step and allow for stochastic fluctuations. This leads us to consider the following stochastic differential equation:
(5) 
where is the outflow of water (), is the mean rainfall in the previous 30 days (mm) and is the standard Wiener process. The term expresses an instantaneous stochastic Gaussian noise in the stored volume, with zero mean and standard deviation proportional to the stored volume. We also assume that the observed values of stored volume follow a Gaussian distribution with mean value equal to the expected value of the stochastic equation and an unknown standard deviation . The reason to take as a 30day mean is empirical: volume does not respond instantaneously to rainfall, and complex hydrological processes tend to spread effects over time. This nonlocality in time is encompassed by taking a mean rainfall over a certain period. The period itself (30 days) is justified a posteriori as the one giving the best fit to data.
Equation (5) keeps the deterministic skeleton from the previous discussion, at the same time allowing us to take , , and inflow from timeseries of recorded data for the Cantareira reservoir. We recall that our model predicts catastrophic shifts because water inflow depends on the stored volume, which creates a feedback mechanism. To evaluate the support provided by data to this hypothesis, we evaluated a competing model in which inflow depends only on rainfall:
(6) 
The dataset are the time series of daily inflow, outflow, rainfall, and stored water volume. The model provided good fits for time spans up to one year. In this paper we focused on the last 365 of the time series (May 31 2014 to May 31 2015). Because simultaneous fitting of all five parameters ( , , , and ) did not converge or provided unreasonable estimates, we callibrated both models in 3 steps:

The observation error was conservatively estimated from trajectory matching (King et al. , 2009).

The exponent was estimated as the slope of a Gaussian linear regression of reservoir area in function of reservoir volume in log scale.
Equation 5, which describes expected water inflow as a function of rainfall and current stored volume, provided a much more plausible fit to the observed time series (loglikelihood ratio: 12.11, Tab. 1). The model that does not take into account the effect of volume on inflow (Eq.(6)) underestimated the stored volume in most of the period, and did not predict the increase and further stabilization of stored volume since February 2015 (Fig.(8)). The better fit of equation 5 supports the hypothesis that the ratio of inflow to rainfall depends on the volume. This feedback is caused by the interaction between rainfall and stored volume, a surrogate of the hydrological state of the catchment. This, in turn, substantiates our statement about the existence of alternative states due to a feedback process.
Model  LL  

Eq.5:  5.998  1.043  0.590  0.00231  6676.0  
Eq.6:  392494.3  0.926  –  0.00948  6688.1 
4 Concluding remarks
We have demonstrated that the dynamics of the volume of water stored in a reservoir is subject to abrupt regime shifts and hysteretic behavior. These results follow from: (i) a qualitative assessment of the dynamics; (ii) a quantitative evaluation of indexes characterizing regime shifts; (iii) a modeling approach that gives a mechanistic view of the dynamics; and (iv) a fit of a corresponding stochastic model to the data.
These results have serious significance for water reservoirs management. Resilience of the normal operating conditions may be lost under drought conditions in favor of a new low water levels regime, which is itself resilient under a range of external conditions. Recovery of the normal condition can, thus, be difficult, and may demand potentially extreme measures. These measures are those that allow for decrease in the outflow of the reservoir, representing a possible burden to populations served by the reservoir. This implies that cautious management of a reservoir should avoid the regime shift. The major question then becomes: how to achieve such management? We have taken a first step here, identifying the need to acknowledge the potential for such catastrophic transitions. Under this perspective, early warning signals of catastrophic regime shifts are a possible avenue to tackle the problem in the absence of detailed hydrological data, but the extent of their predictive power has yet to be ascertained.
Studies in water resources have mainly favored fullblown models to depict the ever increasing understanding of hydrological processes (Beven, 2011). Nevertheless simple mathematical descriptions of the first principles that operate in catchments also proved to be useful both in fitting data and understanding the underlying mechanisms (Beven, 2011; Young et al. , 2004). We followed this synthetic approach to propose a model akin to those used in system dynamics. Our models are a baretothebones description on how rainfall turns into stored water in a reservoir. In doing that we depict many hydrological processes in a few fixed parameters. What could remain from such a phenomenological approach are very general dynamic properties. In this sense, we show how alternative states and regime shifts result from a very simple first principle, namely, nonlinear feedback mechanisms. We show that a model that has these feedbacks – and not much more than that – provides a more plausible description of observed data. In doing so we aim to highlight key properties that water reservoir share with other socioenvironmental systems, which can contribute to fullyfledged models.
In this study we focused on a specific reservoir, the Cantareira system in the metropolitan area of São Paulo. Let us now elaborate on this specific case. Signals of below the average rainfall were present since mid2013, with corresponding decrease in water inflow. The conditional variance increased steadily in the very beginning of 2014, indicating an ongoing regime shift. The average inflow rate in January 2014 was 15.7 , whereas the outflow was kept at 34.3 , a very unusual deficit for the rainy season in the region. Actually, the 30day moving average of inflows and outflows presented deficits continually from May 2013 to February 2015, contrary to the recurrent pattern where the period from DecemberMay always had greater inflows than outflows. The operators began decrementing releases only in March 2014, and at a pace that was not enough to bring back the system to the normal operating conditions. The obvious question is why the management policies could not avoid the collapse of the system. Official documents available to the public (Agência Nacional de Águas & Departamento de Águas e Energia Elétrica, 2015) show that the Brazilian Water Agency (ANA) and the São Paulo State Water Department (DAEE) agreed that operators (SABESP and other water companies) could take up to 36 of water from the Cantareira system. The agreement also included a reservoir rule curve that prescribed maximum allowed withdrawals according to the stored volume of the system. These were not respected since January 2014. Even if the limits were met, they were probably too loose. For instance, the outflow allowed for the rainy months of November/December 2014 would be 27 even if the operational capacity had been completely depleted. Since March 2014, ANA and DAEE abandoned the rule curve and adopted tighter limits. The maximum values have been negotiated periodically, and are currently 17.5 (JuneAugust 2015) and 13.5 (SeptemberNovember 2015). Nevertheless, the National Center for Surveillance of Natural Disasters (http://www.cemaden.gov.br/) forecasts that the Cantareira system will stay below its operational capacity at least until the beginning of 2016.
In summary, our results and our study case show that the management of reservoirs should take alternative regimes into account and avoid a transition to lowvolume regimes. Failing to do so represents a prolongated burden, extending well beyond the period of anomalous rainfall, because outflow has to be kept as low as possible until a backwards transition occurs. Therefore managers should act as another feedback mechanism in the socioenvironmental system that keeps it in the desired regime notwithstanding external forces like climate anomalies. In not doing so, managers of the Cantareira system acted like one more external force that pushed the reservoir to a catastrophic shift.
Data Sets
Data for the time series of stored water volume in the main reservoir system of São Paulo State (Cantareira) and for the corresponding rainfall index were obtained from the São Paulo Water Company (SABESP) database, which can de accessed at http://www2.sabesp.com.br/mananciais/DivulgacaoSiteSabesp.aspx. Data for the inflow and outflow were obtained from http://www2.sabesp.com.br/mananciais/Divulgacaopcj.aspx and processed through Optical Character Recognition (OCR) software. Starting January 15th, 2015, a more convenient tool has been provided by SABESP at http://site.sabesp.com.br/site/interna/Default.aspx?secaoId=553. All data and codes used in the analyses are available at https://github.com/cantareira/plos.
Acknowledgments
The authors thank FAPESP, CNPq and CAPES (Brazil) for partial support.
References
 Agência Nacional de Águas & Departamento de Águas e Energia Elétrica (2015) Agência Nacional de Águas, & Departamento de Águas e Energia Elétrica. 2015 (June). Dados de referência acerca da outorga do sistema Cantareira. http://arquivos.ana.gov.br/institucional/sof/Renovacao_Outorga/DDR_Sistema_Cantareira12Jun15FINAL.pdf.
 Beven (2011) Beven, Keith J. 2011. Rainfallrunoff modelling: the primer. John Wiley & Sons.
 Boettiger & Hastings (2012) Boettiger, Carl, & Hastings, Alan. 2012. Quantifying limits to detection of early warning for critical transitions. Journal of The Royal Society Interface, 9(75), 2527–2539.
 Boettiger et al. (2013) Boettiger, Carl, Ross, Noam, & Hastings, Alan. 2013. Early warning signals: the charted and uncharted territories. Theoretical ecology, 6(3), 255–264.
 Carpenter & Brock (2011) Carpenter, SR, & Brock, WA. 2011. Early warnings of unknown nonlinear shifts: a nonparametric approach. Ecology, 92(12), 2196–2201.
 Carpenter (2005) Carpenter, Stephen R. 2005. Eutrophication of aquatic ecosystems: bistability and soil phosphorus. Proceedings of the National Academy of Sciences of the United States of America, 102(29), 10002–10005.
 Dakos (2015) Dakos, Vasilis. 2015. Early warning signals toolbox. http://www.earlywarningsignals.org. accessed on 29 Jan 2015.
 Dakos et al. (2012) Dakos, Vasilis, Carpenter, Stephen R, Brock, William A, Ellison, Aaron M, Guttal, Vishwesha, Ives, Anthony R, Kefi, Sonia, Livina, Valerie, Seekell, David A, van Nes, Egbert H, et al. . 2012. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PloS one, 7(7), e41010.
 Dakos et al. (2015) Dakos, Vasilis, Carpenter, Stephen R, van Nes, Egbert H, & Scheffer, Marten. 2015. Resilience indicators: prospects and limitations for early warnings of regime shifts. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 370(1659), 20130263.
 Folke (2006) Folke, Carl. 2006. Resilience: The emergence of a perspective for social–ecological systems analyses. Global environmental change, 16(3), 253–267.
 Gordon et al. (2008) Gordon, Line J, Peterson, Garry D, & Bennett, Elena M. 2008. Agricultural modifications of hydrological flows create ecological surprises. Trends in Ecology & Evolution, 23(4), 211–219.
 King et al. (2009) King, Aaron A, Ionides, Edward L, Bretó, CM, Ellner, Steve, Kendall, Bruce, Wearing, Helen, Ferrari, Matthew J, Lavine, Michael, & Reuman, Daniel C. 2009. pomp: Statistical inference for partially observed Markov processes. http://pomp.rforge.rrproject.org.
 Klausmeier (1999) Klausmeier, Christopher A. 1999. Regular and irregular patterns in semiarid vegetation. Science, 284(5421), 1826–1828.
 Liu & West (2001) Liu, Jane, & West, Mike. 2001. Combined parameter and state estimation in simulationbased filtering. Pages 197–223 of: Sequential Monte Carlo methods in practice. Springer.
 Meron et al. (2004) Meron, Ehud, Gilad, Erez, von Hardenberg, Jost, Shachak, Moshe, & Zarmi, Yair. 2004. Vegetation patterns along a rainfall gradient. Chaos, Solitons & Fractals, 19(2), 367–376.
 Österblom et al. (2013) Österblom, Henrik, Merrie, Andrew, Metian, Marc, Boonstra, Wiebren J, Blenckner, Thorsten, Watson, James R, Rykaczewski, Ryan R, Ota, Yoshitaka, Sarmiento, Jorge L, Christensen, Villy, et al. . 2013. Modeling SocialEcological Scenarios in Marine Systems. BioScience, 63(9), 735–744.
 Scheffer & Carpenter (2003) Scheffer, Marten, & Carpenter, Stephen R. 2003. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends in ecology & evolution, 18(12), 648–656.
 Scheffer & van Nes (2007) Scheffer, Marten, & van Nes, Egbert H. 2007. Shallow lakes theory revisited: various alternative regimes driven by climate, nutrients, depth and lake size. Hydrobiologia, 584(1), 455–466.
 Scheffer et al. (2001) Scheffer, Marten, Carpenter, Steve, Foley, Jonathan A, Folke, Carl, & Walker, Brian. 2001. Catastrophic shifts in ecosystems. Nature, 413(6856), 591–596.
 Young et al. (2004) Young, Peter C, Chotai, Arun, & Beven, Keith J. 2004. Databased mechanistic modelling and the simplification of environmental systems. Pages 371–388 of: Wainwright, John, & Mulligan, Mark (eds), Environmental Modelling: Finding Simplicity in Complexity. Chichester, UK: Wiley.