Home | Archives | About | Login | Submissions | Notify | Contact | Search

 E&S Home > Vol. 11, No. 2 > Art. 9

Copyright © 2006 by the author(s). Published here under license by The Resilience Alliance.
Go to the pdf version of this article

The following is the established format for referencing this article:
Brock, W. A., and S. R. Carpenter. 2006. Variance as a leading indicator of regime shift in ecosystem services. Ecology and Society 11(2): 9. [online] URL: http://www.ecologyandsociety.org/vol11/iss2/art9/


Research

Variance as a Leading Indicator of Regime Shift in Ecosystem Services

William A. Brock1 and Stephen R. Carpenter1


1University of Wisconsin-Madison



ABSTRACT


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.


Key words: ecological economics; ecosystem service; indicators; markets; pollution; regime shifts; thresholds; variance



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-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.



MODEL AND METHODS

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 CO2 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
Equation 1
(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 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 M1 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(Pj - 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 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. P1 and P2. 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 and Lefever 1984, Berglund and Gentz 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
Equation 2
(2)

where Zi,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 P1, P2, S1, and S2 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.



RESULTS

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 and Brock 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

Equation 3
(3)

(Brock et al. 2006, Carpenter and Brock 2006).

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 to zero. Earlier papers discuss this point in more detail (Brock et al. 2006, Carpenter and Brock 2006).

Although this heuristic shows how the small shocks are magnified near a regime shift, it is a one-dimensional 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 (M1 = 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 high-pollutant 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 (M1 = 1.25, upper left), the probability mass is concentrated at low pollutant levels, and thus the overall variance of Pi 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 Pi. The dominance of probability peak at high P becomes successively greater when emission is 1.75 and then 2.0, suggesting lower variance of Pi.

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 M1.

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.



DISCUSSION

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 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.

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. 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 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 distributions of shocks (Scheffer et al. 2001, Carpenter 2003, Walker and Meyers 2005). More research is needed to ascertain which kinds of ecosystem regime shifts are likely to be signaled in advance by increases in variance.

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.

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 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 and Brock 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.



RESPONSES TO THIS ARTICLE


Responses to this article are invited. If accepted for publication, your response will be hyperlinked to the article. To submit a response, follow this link. To read responses already accepted, follow this link




ACKNOWLEDGMENTS

Both authors made equal contributions to this paper. We thank the referees and editors for helpful suggestions. We are grateful for the support of the Environmental Protection Agency and the National Science Foundation.



LITERATURE CITED

Bates, D. S. 1991. The crash of '87: Was it expected? The evidence from options markets. Journal of Finance 46:1009-1044.

Berglund, N., and B. Gentz. 2002a. Metastability in simple climate models: pathwise analysis of slowly driven Langevin equations. Stochastic Dynamics 2:327-356.

Berglund, N., and B. Gentz. 2002b. Beyond the Fokker-Planck equation: pathwise control of noisy bistable systems. Journal of Physics A 35:2057-2091.

Brock, W. A. 2006. Tipping points, abrupt opinion changes, and punctuated policy change. Available online at: http://64.233.161.104/search?q=cache:7fCVq1OqO7UJ:www.ssc.wisc.edu/econ/archive/wp2003-28.pdf+Brock %22Tipping+points%22&hl=en&gl=ca&ct=clnk&cd=1.

Brock, W. A., S. R. Carpenter, and M. Scheffer. 2006. Regime shifts, environmental signals, uncertainty and policy choice. In J. Norberg and G. Cumming, editors. A theoretical framework for analyzing social-ecological systems. Columbia University Press, New York, New York, USA, in press.

Carpenter, S. R. 2003. Regime shifts in lake ecosystems: pattern and variation. Excellence in Ecology Series, Volume 15. Ecology Institute, Oldendorf/Luhe, Germany.

Carpenter, S. R., and W. A. Brock. 2004. Spatial complexity, resilience, and policy diversity: fishing on lake-rich landscapes. Ecology and Society 9(1): 8. [online] URL: http://www.ecologyandsociety.org/vol9/iss1/art8/.

Carpenter, S. R., and W. A. Brock. 2006. Rising variance: a leading indicator of ecological transition. Ecology Letters 9:311-318.

Carpenter, S. R., D. Ludwig, and W. A. Brock. 1999. Management of eutrophication for lakes subject to potentially irreversible change. Ecological Applications 9:751-771.

DeAngelis, D., and W. Mooij. 2003. In praise of mechanistically-rich models. Pages 63-82 in C. D. Canham, J. J. Cole, and W. K. Lauenroth. Models in ecosystem science. Princeton University Press, Princeton, New Jersey, USA.

Duffie, D. 1988. Security markets: stochastic models. Academic Press, New York, New York, USA.

Grossman, S. 1989. The informational role of prices. MIT Press, Cambridge, Massachusetts, USA.

Horsthemke, W., and R. Lefever. 1984. Noise-induced transitions. Springer-Verlag, New York, New York, USA

Kleinen, T., H. Held, and G. Petschel-Held. 2003. The potential role of spectral properties in detecting thresholds in the earth system: application to the thermohaline circulation. Ocean Dynamics 53:53-63.

Ludwig, D., W. A. Brock and S. R. Carpenter. 2005. Uncertainty in discount models and environmental accounting. Ecology and Society 10 (2): 13. [online] URL: http://www.ecologyandsociety.org/vol10/iss2/art13/.

Millennium Ecosystem Assessment (MA). 2006. Ecosystems and human well-being. Our human planet: summary for decision makers. Island Press, Washington, D.C., USA.

Nakajima, H., and D. L. DeAngelis. 1989. Resilience and local stability in a nutrient-limited consumer system. Bulletin of Mathematical Biology 51: 501-510.

Oborny, B., G. Meszena, and G. Szabo. 2005. Dynamics of populations on the verge of extinction. Oikos 109:291-296.

Rosenzweig, M. L., and R. H. MacArthur. 1963. Graphical representation and stability conditions for predator-prey interactions. American Naturalist 97:209-223.

Scheffer, M., and S. R. Carpenter. 2003. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends in Ecological Evoluation 12:648-656.

Scheffer, M., S. R. Carpenter, J. Foley, C. Folke, and B. Walker. 2001. Catastrophic shifts in ecosystems. Nature 413:591-596.

Scheffer, M., F. Westley, and W. Brock. 2003. Slow response of societies to new problems: causes and costs. Ecosystems 6:493-502.

van Nes, E. H., and M. Scheffer. 2003. Alternative attractors may boost uncertainty and sensitivity in ecological models. Ecological Modeling 159:117-124.

Walker, B., and J. A. Meyers. 2004. Thresholds in ecological and social-ecological systems: a developing database. Ecology and Society 9(2): 3, [online] URL: http://www.ecologyandsociety.org/vol9/iss2/art3/.


Address of Correspondent:
Stephen R. Carpenter
Center for Limnology
University of Wisconsin-Madison
Madison, Wisconsin 53706 USA
srcarpen@wisc.edu

Home | Archives | About | Login | Submissions | Notify | Contact | Search