Variance as a Leading Indicator of Regime Shift in Ecosystem Services

Many environmental conflicts involve pollutants such as greenhouse gas emissions that are dispersed through space and cause losses of ecosystem services. As pollutant emissions rise in one place, a spatial cascade of declining ecosystem services can spread across a larger landscape because of the dispersion of the pollutant. This paper considers the problem of anticipating such spatial regime shifts by monitoring time series of the pollutant or associated ecosystem services. Using such data, it is possible to construct indicators that rise sharply in advance of regime shifts. Specifically, the maximum eigenvalue of the variance-covariance matrix of the multivariate time series of pollutants and ecosystem services rises prior to the regime shift. No specific knowledge of the mechanisms underlying the regime shift is needed to construct the indicator. Such leading indicators of regime shifts could provide useful signals to management agencies or to investors in ecosystem service markets.


INTRODUCTION
Regime shifts are substantial reorganizations of social-ecological systems (Scheffer et al. 2001, Carpenter 2003, Walker and Meyers 2004, Brock 2006).Although examples are known from a wide variety of systems, regime shifts are difficult to study (Scheffer and Carpenter 2003) because they usually involve spatially extensive regions, stochastic processes, events that are infrequent relative to a human lifetime, and processes at two or more scales of space and time (Carpenter 2003).Because of these challenges, we are only in the early stages of understanding environmental regime shifts.Nevertheless, regime shifts have important implications for human livelihoods and well-being.The Millennium Ecosystem Assessment (MA 2006) noted that regime shifts, which it called "nonlinear changes," could substantially alter the flow of ecosystem services, i.e., the benefits that people obtain from ecosystems.
Recent research shows that ecosystem dynamics become more variable prior to some regime shifts (Berglund and Gentz 2002a,b, Kleinen et al. 2003, Oborny et al. 2003, van Nes and Scheffer 2003, Brock et al. 2006, Carpenter and Brock 2006).This finding suggests that changes in variance can be a leading indicator of an impending regime shift.Brock et al. (2006) explain how this can occur for one-dimensional systems.Carpenter and Brock (2006) showed that the variance component related to an impending regime shift could be separated from environmental noise using methods that required no knowledge of the mechanisms underlying the regime shift.This finding suggests that it might be possible to detect impending regime shifts using routine monitoring data even in the absence of a detailed ecological understanding of the underlying nonlinearities.It remains to be seen how variance indicators will perform in the real world of managing ecosystem services.
This paper expands the study of variance indicators of impending regime shifts in two main ways.First, we extend the approach to multivariate, spatially distributed systems.For ease of exposition, we analyze a two-dimensional system, but our methods generalize directly to higher-dimensional systems.In this context, we sketch a framework for choosing empirical variance indicators that are likely to be sensitive to impending regime shifts.Second, our model couples ecosystem dynamics to the flow of ecosystem services.This suggests a forward-http://www.ecologyandsociety.org/vol11/iss2/art9/looking management system based on expected future flows of ecosystem services.Although this paper focuses on variance indicators rather than policies for future flows of ecosystem services, our findings can be nested in a more complete analysis of ecosystem service markets and policies.

Overview of the model
We consider a minimal model of two regions connected by the transport of a pollutant (Fig. 1).The pollutant is released as a side effect of activities that yield economic benefits.However, the pollutant has adverse effects on an ecosystem service such as food production or the regulation of water quality.Thus, the net present value to the society in a region is a balance between polluting activity and the ecosystem service.This situation is analogous to many situations in modern society (MA 2006).Some examples are (1) greenhouse gas pollution from burning fossil fuels, which produces the economic benefit of energy production but also has adverse effects on climate and ecosystems; (2) nitrogen emissions from livestock feeding operations, which are linked to the economic benefit of food production but also produce adverse effects on air and water quality; and (3) pollution with persistent organic toxins, which are associated with economic benefits but also cause unwanted harm to ecosystems and human health.Within a region, the pollutant can be permanently removed (the decay process in Fig. 1) or sequestered (the sink process in Fig. 1).The decay process is linear in the pollutant, and the sink process becomes saturated at high pollutant levels.The linear process is rate limited only by the supply of the pollutant.An example is the destruction of organic pollutants by UV irradiance.The saturating process is rate limited by factors other than pollutant supply.An example is sequestration in soils of phosphorus, an important cause of water pollution, which can be limited by the capacity of the soil to bind phosphate ion.Another example is sequestration of CO 2 in ecosystems, which is limited by the capacity of the plants to assimilate carbon.
The pollutant moves between regions from the region of higher abundance to the region of lower abundance (the mixing process in Fig. 1).Many pollutants are transported across regional boundaries.Atmospheric pollutants such as greenhouse gases are the most familiar example.Water pollutants are transferred among regions that share shorelines on a single body of water.Other pollutants are directly transported by people.For example, manure is a waste product of animal production that causes air and water pollution and is also transported from regions of high manure production to regions in which it is needed as a fertilizer.

Model details
Within each region i, the dynamics of the pollutant P and ecosystem service S are assumed to follow (1) where M is the emission rate, D is the transport rate between region j and region i, c is the decay rate, s is the maximum sequestration rate, f(P) scales the sequestration to the amount of pollutant, σ is the standard deviation of shocks to the system, dW is a Wiener noise process that is N(0,dt) and independent for each i, m is the pollutant level at which sequestration is half the maximum rate, q determines the slope of f(P) near P = m, r is the renewal rate of the ecosystem service, k is the impact of the pollutant on renewal of the ecosystem service, and h is a parameter for the effects of extracting or using the ecosystem service.The model for S is chosen to be a simple function that yields S inversely proportional to P. The dynamics studied in this paper result from the behavior of P. Future work should extend the approach to more realistic ecosystem service models.
A rational social planner would consider the net benefits of pollutant-emitting activities, which generate M, and the ecosystem service, which generates hS and is reduced by increasing P, as well as any other economic consequences of P. The expected net benefits, integrated over infinite time with an appropriate discounting factor, would then be considered (Ludwig et al. 2005).Policies for M http://www.ecologyandsociety.org/vol11/iss2/art9/and h would be chosen to maximize these expected net benefits.We will leave the economic analysis of this model for later work and focus on the ecological dynamics in this paper.Here we will assume that economic considerations lead to some appropriate choice of M and h for each region.These choices may change as events unfold and circumstances change over time.We wish to study the dynamics as M 1 slowly rises past a critical point.
The effects of three levels of M within one region are presented in Fig. 2. The blue curve shows sf(P) vs. P.The red straight lines show the graph of the function of P defined by M + D(P j -P) -cP for three values of M. Intersections denote equilibria.When M is low, there is one stable equilibrium at a relatively low value of P. When M is high, there is one stable equilibrium at a relatively high level of P. At the intermediate M, there are three equilibria.These are unstable, whereas the low-P and high-P equilibria are stable.This paper will focus on the case in which M rises slowly and gradually from a value with two stable P equilibria (like the middle red line) past the critical point at which the low-P http://www.ecologyandsociety.org/vol11/iss2/art9/equilibrium disappears, leaving only the high-P equilibrium (like the upper red line).

Model analysis
To build intuition about the dynamics and variability of P, we plot potential surfaces and stationary distributions for the model using equations presented in Appendix 1.A ball-and-cup diagram for the system is given by a plot of the negative potential vs. P 1 and P 2 .Here we plot −log (potential + constant) for convenience, where the constant is chosen to yield positive arguments to the log function.Note that the maxima and minima of −log(potential + constant) are identical to those of −potential.Over a long period of time, the system approaches a stationary probability distribution (Appendix 1).The spread of the stationary distribution is related to the variability of P near steady state.The stationary probability distribution is the limit distribution that can be shown to be unique for the diffusion systems considered here.It can be found by solving for the steady-state density of the Fokker-Planck equation (Horsthemke andLefever 1984, Berglund andGentz 2002b).It is possible to study regime shifts of stochastic systems via bifurcation analysis of the steady-state solutions of the Fokker-Planck equations.We do not follow that path here, because the system has a potential.However, many ecologically interesting systems of two dimensions and higher do not have potentials (Brock et al. 2006).The framework presented in this paper can be generalized to these cases by studying steady-state solutions of the Fokker-Planck equations.
To study the dynamics of variance, we simulated time series of P and S using the Euler method adapted for the Ito solution of stochastic differential equations (Horsthemke and Lefever 1984).We do this by iteratively computing using Eq. 2 below.Each time step dt is of length 1/n.Over each small increment we compute (2) where Z i,t is an independently drawn random number (normal, mean 0, variance 1) for each small time increment and region.Results presented here use n = 36.For each time step, the n values of P 1 , P 2 , S 1 , and S 2 are used to compute a covariance matrix.The variance indicator for the time step is the largest eigenvalue of that covariance matrix.This eigenvalue is strongly affected by variance in the spatial unit that is approaching a critical point, but only lightly affected by spatial units that are merely noisy.Technical details on choosing variance indicators for impending regime shifts are addressed in Appendix 1.

Variance near a regime shift: simple case
Here we introduce a simplified heuristic to build intuition about changes in variance near a regime shift (Brock et al. 2006, Carpenter andBrock 2006).
The key is to recognize three time scales: (1) very slow change in an exogenous driver or slowly changing system component, (2) medium-speed change in the state variable subject to regime shift, and (3) fast changes in random shocks to the rate equation for the state variable.Assume that these three time scales are quite different, in the sense that a faster variable can approach equilibrium over a time period of only slight change in a slower variable.For this paper, the slowest speed is the very gradual change in pollutant input M, the focal state variable (medium speed) is P, and the fast shock process is denoted by dW.We introduce the variable a to represent the discrete shocks drawn from a random-number generator to simulate the solution of the P dynamics over very short intervals of time.For a simplified one-dimensional case, define dP/ dt = V(M,P,a), let var(da) denote the variance of the random disturbances (a constant in this case), and let var(dP) represent the variability we observe in P over a time interval during which many small random shocks occur.Suppose that M is moving so slowly that P is near an equilibrium, that is, 0 = V (M,P + dP,a + da).For each shock a + da, we assume the fast dynamics of a are so fast that the system has relaxed to a new steady-state P + dP on the time scale relevant for P. Then it can be shown, by expanding V(M,P + dP,a + da) in a Taylor series about (P,a) that Fig. 2. Rate of input + diffusion -decay (red lines) for three values of pollutant discharge in region 1, M 1 (0.9, 1.5, and 2.1), and the rate of loss from the saturating sink (blue line).For a given minimum pollution rate in region 1, M 1 , equilibria occur at the points at which the red line crosses the blue line.The pollutant level in the other region was set at 1.6.Other parameters: M 2 (pollution discharge rate in region 2) = 0.75; D (transport rate between regions) = 0.1; c (decay rate) = 0.05; s (maximum sequestration rate) = 0.5; m (pollutant level at which sequestration is half the maximum rate) = 5; q, which determines the slope of f(P) near P = m, = 6.P stands for amount of pollutant, and f(P) scales the sequestration to the amount of pollutant.
To understand the implications of Eq. 3 for variance near a regime shift, consider Fig. 3, which plots V vs. P for a case with lower M (blue line) and higher M (red line).Suppose the pollutant level is near the lower stable equilibrium.The red line (higher M) is closer to regime shift.Clearly, dV/dP is lower for the red line, and therefore var(dP) must be larger for the case closer to regime shift.Intuitively, the variance of da gets magnified more and more as dV/ dP, evaluated at steady-state P, gets closer and closer http://www.ecologyandsociety.org/vol11/iss2/art9/ to zero.Earlier papers discuss this point in more detail (Brock et al. 2006, Carpenter andBrock 2006).
Although this heuristic shows how the small shocks are magnified near a regime shift, it is a onedimensional case that makes simplifying assumptions about time scales, proximity to equilibrium, and near-linearity.We now turn to the more complex spatial model by first considering the potential surfaces and stationary distributions and then presenting a simulation example.

Potential surfaces and stationary distributions
To understand changes as pollutant emission slowly increases in one region, we plotted negative potential surfaces for the pollutants at four levels of emission in region 1 (Fig. 4).Valleys in the surface represent attractors, as in a ball-and-cup stability diagram.At the lowest emission rate (M 1 = 1.25, upper left), there is a deep attractor at low pollutant levels and a shallow attractor at moderately high levels of the pollutant in both regions.When the emission rate is increased to 1.5, a third attractor appears, with moderately high pollutant levels in region 1 and low pollutant levels in region 2. The attractor with low pollutant levels becomes more shallow, and the attractor with relatively high levels of both pollutants deepens.When the emission rate increases to 1.75, the attractor with relatively high pollutant levels becomes even deeper, and the attractor with low pollutant levels becomes quite shallow.At an emission rate of 2.0, the highpollutant valley appears to be the only attractor.
Stationary probability distributions are roughly the same as the negative potential surfaces turned upside down (Fig. 5).At the lowest emission rate (M 1 = 1.25, upper left), the probability mass is concentrated at low pollutant levels, and thus the overall variance of P i is relatively low.When the emission rate is increased to 1.5, there are three distinct peaks, with the greatest mass of probability at relatively high pollutant levels, suggesting fairly high variance of P i .The dominance of probability peak at high P becomes successively greater when emission is 1.75 and then 2.0, suggesting lower variance of P i .

Dynamics of variance
Although the potential surfaces and probability distributions are useful for building intuition about the attractors for the system and the dispersion of P over long periods of time, a person who is observing the system and making decisions about it may not have this information.It is more likely that the decision maker will observe time series of pollutant levels or the supply or price of the ecosystem service and draw inferences from these.Observable features of these time series can provide signals of impending regime shift in advance of the actual regime shift.To illustrate this point, we simulated dynamics over long periods of time while slowly raising M 1 .
The time series from such a simulation show two regime shifts (Fig. 6).In the first regime shift, region 1 moves to a moderately high pollutant level, followed by the movement of region 2 to a moderately high pollutant level.In the second regime shift, region 2 moves to a high pollutant level, followed by the movement of region 1 to a high pollutant level.Corresponding shifts in ecosystem services are opposite in direction to the shifts in pollutant level.
The variance index spikes before each regime shift (Fig. 7).If we focus on a narrower window of time around the regime shift, it appears that the increase in variance occurs roughly 5-20 time steps prior to the regime shift (Fig. 8).Note that variance is smoothed such that the index at time t = average of index at times t, t -1, t -2.This backward-looking smoothing might degrade the lead time of the indicator.Despite this, there is still a clear rise in the indicator prior to each regime shift.

Rising variance and regime shifts
Our findings corroborate earlier work showing that ecosystem behavior becomes more variable prior to a regime shift (Berglund and Gentz 2002a,b, Kleinen et al. 2003, Oborny et al. 2003, van Nes and Scheffer 2003, Brock et al. 2006, Carpenter and Brock 2006).This paper shows that rising variance serves as a leading indicator in a situation in which ecosystems are coupled across space.The detection of rising variance requires no special knowledge of http://www.ecologyandsociety.org/vol11/iss2/art9/Fig. 3. Net rate of change for a single region vs. pollutant level for two rates of input to region 1.Blue: the pollution rate for region 1 (M 1 ) = 1.6; red: M 1 = 1.9).Circles denote stable points.The pollutant level in the other region was fixed at 1.6.Other parameters: M 2 (the pollution rate in region 2) = 0.75; D (the transport rate between regions) = 0.1; c (the decay rate) = 0.05; s (the maximum sequestration rate) = 0.5; m (the pollutant level at which sequestration is half the maximum rate) = 5; q (which determines the slope of f(P) near P) = 6.P stands for amount of pollutant, and f(P) scales the sequestration to the amount of pollutant.ecosystem dynamics.A simple indicator based on the dominant eigenvalue of the covariance matrix for observed time series can detect the rising variance.This is so because the dominant eigenvalue is strongly influenced by the spatial unit that is closest to a regime shift.Therefore, the dominant eigenvalue can provide early clues of regime shift even in cases in which the investigator is unsure about which ecosystem variables are most important.Also, the dominant eigenvalue may be relatively insensitive to junk variables that are highly noisy but carry no signal related to an impending regime shift.In Appendix 1, we address technical issues related to constructing variance indicators for regime shifts.Potential surface (negative log potential) vs. pollutant in region 1 and region 2 for four different input rates to region 1.Parameters: M 1 (pollution discharge rate in region 1) = 1.25, 1.5, 1.75, and 2; M 2 (pollution discharge rate in region 2) = 0.75; c 1 (decay rate in region 1) = 0.05; c 2 (decay rate in region 2) = 0.1; D (transport rate between regions) = 0.1; s (maximum sequestration rate) = 0.5; m (the pollutant level at which sequestration is half the maximum rate) = 5; q (which determines the slope of f (P) near P) = 6.P stands for amount of pollutant, and f(P) scales the sequestration to the amount of pollutant.
In practice, the indicator will be sensitive to the time interval over which the covariance matrix is calculated.The time interval should be short enough that the slow driver, M in this case, does not change much, and long enough for the shocks to approach the stable distribution.Carpenter and Brock (2006) comment on some issues related to filtering data for variance signals.In general, however, much work remains to be done in the area of designing filters for variance signals of regime shifts in data.
It is important to note that we studied a "hard" regime shift in which there was a discontinuous change because of the disappearance of an attractor.Fig. 5. Stationary distribution vs. pollutant in region 1 and region 2 for four different input rates to region 1.Parameters: M 1 (pollution discharge rate in region 1) = 1.25, 1.5, 1.75, and 2; M 2 (pollution discharge rate in region 2) = 0.75; c 1 (decay rate in region 1) = 0.05; c 2 (decay rate in region 2) = 0.1; D (transport rate between regions) = 0.1; s (maximum sequestration rate) = 0.5; m (the pollutant level at which sequestration is half the maximum rate) = 5; q (which determines the slope of f(P) near P) = 6.P stands for amount of pollutant, and f(P) scales the sequestration to the amount of pollutant.
For other kinds of regime shifts, variance signals may be weak or absent.For example, the pitchfork bifurcation, in which one regime splits continuously into two regimes, may not be preceded by an increase in variance, depending on the standard deviation of the shocks and the relaxation time of the fast variables in the system (Berglund and Gentz 2002b).A simple ecological example is the paradox of enrichment (Rosenzweig and MacArthur 1963), in which gradual nutrient enrichment destabilizes a resource-consumer system.In the original case, there is a rise in variance before the regime shift (Nakajima and DeAngelis 1989) comparable to the phenomenon described in this paper.However, if the paradox of enrichment is placed in an ecosystem model with nutrient recycling by the consumer, variance decreases over a period of time as the critical point is approached, then increases sharply Fig. 6.Time series of ecosystem services (upper panel) and pollutant (lower panel) in each region vs. time as M 1 (pollution discharge rate in region 1) increases linearly from 2.0 to 2.1.Other parameters: r (renewal rate of the ecosystem service) = 1; h (effects of extracting or using the ecosystem services) = 0.5; k (impact of the pollutant on the renewal rate of the ecosystem service) = 0.5; M 2 (pollution discharge rate in region 2) = 0.75; c 1 (decay rate in region 1) = 0.05; c 2 (decay rate in region 2) = 0.1; D (transport rate between regions) = 0.1; s (maximum sequestration rate) = 0.5; m (the pollutant level at which sequestration is half the maximum rate) =5; q (which determines the slope of f(P) near P) = 6; σ = 0.05.P stands for amount of pollutant, and f(P) scales the sequestration to the amount of pollutant.just before the critical point is reached (Nakijima and DeAngelis 1989).Thus, real systems could feasibly have complexities that cause variance to decrease for a while before arriving at a regime shift.Ecosystems seem to exhibit a considerable diversity of regime shifts, a range of speed differences between fast and slow variables, and diverse

Ecosystem services
Our analysis considered changes in an ecosystem service as well as a biogeochemical variable such as a pollutant.Clearly, a more complete economic analysis of this system could be carried out by specifying the form of the expected net present value over future time, although this is not the purpose of the present paper.However, we note that the variance signals discussed here could readily be transferred to economic signals.For example, a jump in pollutant P at some random time in the future will cause ecosystem service S to drop.Furthermore, individuals who have a better knowledge of the system may have more accurate expectations of future ecosystem services than those observers who are not as well informed about the details and structure of the system.Hence, if there were a stock market for trading of claims to the future flow of ecosystem services, then the equity value of such an entity, assuming that equity shares are the only outstanding claims against it, would be given by the expression for discounted net present value over future time (Duffie 1988).Grossman (1989) explains how markets aggregate individual specialized information about assets into market values for those assets.6 with the time frame narrowed to focus on the regime shifts.In this case, the variance index is smoothed by averaging the previous three observations, i.e., variance index at time t = mean of the variance index at times t, t-1, and t-2.
We believe that a market price equation may lead to the development of more sensitive indicators of impending regime change than the variance indicators we develop in this paper.Furthermore, in the case of corporations that have claims on ecosystem services for which there are options markets available on the stocks of such corporations, an even more powerful indicator of impending regime change is available (Bates 1991).Bates compares the prices of "put" options, i.e., options to sell a certain quantity of stock at a specified target price up to a specified date, vs. "call" options, i.e., options to buy a certain quantity of stock at a specified target price up to a specified date.His result applies to "out-of-the-money" options, that is, a call option with a target price higher than the market price of the stock or a put option with a target price lower than the market price of the stock.If out-of-the-money put options are selling for a lot more than out-of-the-money call options for a given stock, then this is strong evidence for an impending sharp drop in the value of the stock.In our context, this would be strong evidence for an impending regime change, i.e., a sharp increase in http://www.ecologyandsociety.org/vol11/iss2/art9/P that harms S. We believe that this approach can be generalized to settings in which the pollutant discharge M or harvest function h(S) is chosen to optimize the expectation of future discounted net present value at each point in time (Carpenter et al. 1999, Carpenter andBrock 2004).
However, markets are available for only a small fraction of the world's ecosystem services, and few of these have options markets (MA 2006).In most cases, markets are absent or severely distorted by government intervention.Thus, in the majority of cases, indicators based on direct monitoring of ecological variables or ecosystem services may be the only means of addressing regime shifts.We believe that research will reveal useful leading indicators of regime shifts in ecosystem services.Such research should use a wider range of models than the aggregated ones presented in this paper, as well as diverse time-series data sets.For example, individual-based models that simulate diverse behaviors of individual human stakeholders and ecosystem components, and their responses to pollutants and services, may yield new insights about social-ecological dynamics prior to regime change (DeAngelis and Mooij 2003).

CONCLUSIONS
Although variance provides a signal of impending regime shift, it remains to be seen whether the message can evoke adaptive social responses.It is plausible that variance signals could trigger appropriate responses by some individual investors in ecosystem service markets.However, by the time the variance signal is clear, the regime shift may be "an accident waiting to happen."Momentum in slowly changing ecosystem components may make it difficult to reverse course (Carpenter 2003).Many factors in societies impede adaptive responses to impending environmental breakdowns (Scheffer et al. 2003).Even if ongoing research verifies the reliability of variance indicators of regime shifts, substantial institutional reforms may be needed to use such indicators in the forward-looking management of ecosystem services.

Fig. 1 .
Fig. 1.Box-arrow diagram of the model showing flows of the pollutant and economic effects in the two regions.

Fig
Fig. 4.Potential surface (negative log potential) vs. pollutant in region 1 and region 2 for four different input rates to region 1.Parameters: M 1 (pollution discharge rate in region 1) = 1.25, 1.5, 1.75, and 2; M 2 (pollution discharge rate in region 2) = 0.75; c 1 (decay rate in region 1) = 0.05; c 2 (decay rate in region 2) = 0.1; D (transport rate between regions) = 0.1; s (maximum sequestration rate) = 0.5; m (the pollutant level at which sequestration is half the maximum rate) = 5; q (which determines the slope of f (P) near P) = 6.P stands for amount of pollutant, and f(P) scales the sequestration to the amount of pollutant.

Fig. 7 .
Fig. 7. Pollutant in region 1 and variance index (upper panel) and pollutant in region 2 and variance index (lower panel) for the simulation shown in Fig. 6.

Fig. 8 .
Fig. 8. Results from Fig.6with the time frame narrowed to focus on the regime shifts.In this case, the variance index is smoothed by averaging the previous three observations, i.e., variance index at time t = mean of the variance index at times t, t-1, and t-2.