Regime shifts in eutrophied lakes: a mathematical study

Ecological Modelling 179 (2004) 115–130

Regime shifts in eutrophied lakes: a mathematical study P.D.N. Srinivasu∗ Department of Mathematics, Andhra University, Visakhapatnam 530003, India Received 3 June 2003; received in revised form 14 April 2004; accepted 3 May 2004

Abstract In this article, a two dimensional coupled differential system is constructed to represent the dynamics of a lake exposed to eutrophication. The aim is to throw some light on the flipping behaviour associated with the lakes by incorporating vital nonlinear interactions between the phosphorus in water and sediments in a simple model. This model reveals important general patterns observed with the lake ecosystems. It explains the surprising dynamics of eutrophied lakes like regime shifts and the reasons for cyclic behaviour observed during the management of lakes, which are not explained by the existing models. The analysis identifies strategies to break away from the cyclic behaviour during the lake management. A necessity to increase the dimension of the model in order to explain the dynamics in irreversible lakes, in which eutrophication cannot be reversed even by severe reductions in the phosphorus input, has been highlighted. Several modifications are suggested to increase the realism of the model and for further contribution. © 2004 Elsevier B.V. All rights reserved. Subj. Class.: 92A17 Keywords: Lake; Eutrophication; Regime; Shifts; Mathematical; Model

1. Introduction Mother nature is extremely complex. Modelling any component of the nature, for example an ecosystem, requires a deep understanding of various interactions between the vital components that constitute the system. Basically, we model a physical system to gain more insight into its behavior which in turn helps in predicting its future and to manage the system efficiently. Although it is very difficult to incorporate all the nonlinear interactions into the mathematical model, we cannot afford to ignore important interactions in the ecosystem in order to make the model sim∗ Corresponding author. Tel.: +91 891 2710016; fax: +91 891 2755324. E-mail address: [email protected] (P.D.N. Srinivasu).

ple. Especially, when we wish to make a mathematical study to manage an ecosystem, we must make sure that the model captures most of the vital interactions of the system. Otherwise, any management decisions made based on a poor model are likely to fail in achieving the objective of the decision during implementation. Alternative steady states and regime shifts in lake ecosystems have fascinated several scientists and it has become a focus of study for many limnologists and ecosystem managers (Scheffer, 1997, 1999; Scheffer et al., 1993; Scheffer and Jeppensen, 1998; Carpenter, 2003; Kohei and Nakajima, 2002; Lisa et al., 2002). Lakes exposed to eutrophication are known to exist in two alternative states, oligotrophy and eutrophy (Scheffer et al., 1993; Smith, 1998). Oligotrophy lakes are characterized by low nutrient inputs, low to moderate levels of plant production

0304-3800/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2004.05.005

116

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

and relatively clear water. Eutrophic lakes have high nutrients inputs, high plant production and murky water (Carpenter et al., 1999a,b; Peterson et al., 2003). Depending on the nature, lakes can be distinguished in the following three ways (Carpenter et al., 1999a). Reversible lakes: eutrophication can be reversed by phosphorus input controls alone. Hysteric lakes: eutrophication can be reduced by severely reducing phosphorus inputs for an extended period of time, or by combining phosphorus input controls with temporary interventions such as chemical treatment. Irreversible lakes: eutrophication cannot be reversed even by sever reductions in phosphorus input, although it may be reversed in some cases by additional costly interventions. Eutrophication of lakes is usually caused by the excessive input of nutrients, primarily phosphorus (Schindler, 1977). Agricultural runoff is the chief source of phosphorus input into the lakes (Carpenter et al., 1998). As the lakes are enriched, the phosphorus from water gets absorbed into the sediments at the bottom of the lake (Schindler et al., 1987). Hence, there is loss of phosphorus to sediments from the water. This accumulated phosphorus in the sediments get recycled into the water again and this process of recycling increases as the accumulaton in sediments increase (Nurnberg, 1984; Soranno et al., 1997). The phosphorus in the water is also reduced due to out flow and sequestration in biomass of consumers or benthic plants (Carpenter et al., 1999a). The recycling of phosphorus takes place in lakes due to several other mechanisms too. For example, the recycling from sediments is often increased due to resuspension by waves or benthivorous fish (Sas, 1989; NRC, 1992; Nurnberg, 1984). Certain proportion of the phosphorus absorbed by the herbivorous zooplankton and fish also get recycled into the water through excretion or egestion. This process is a potential important source of phosphorus which promotes the growth of phytoplankton in lakes (Carpenter et al., 1992; Elser et al., 1998; Schindler and Eby, 1997). Eutrophic conditions of the lake influences growth and regulates stages of development of some of the species living in the environment (Bayraktaroglu et al., 2003). Dynamics of lake eutrophication has been modelled by many scientists. Several models have been

proposed to represent eutrophication. Dimension and complexity varies from model to model. Some models were meant to study current state and understand characteristics of specific lake ecosystems while some others were meant to derive guidelines to manage the eutrophication in lake ecosystems in general. In Kohei and Nakajima (2002), the lake ecosystem is modelled using reaction diffusion advection equations. They show that there are two stable vertical patterns of phytoplankton over a certain range of parameters of their model. This model exhibits multiple stable equilibria and hysterisis effect and catastrophic transition when one of the parameters is changed continuously. Ichiro (1995) obtains some sensitive indices to discriminate between oligotrophy and eutrophy. These are derived based on some experimental data. Christopher et al. (2000) suggest a two component model to account for continuum of oligotrophic to eutrophic conditions in order to model periphyton biogeochemical relationships in the Everglades National Park (ENP). However this model was not structured to simulate shifts in periphyton composition frequently observed as a consequence of phosphorus enrichment. Santanu et al. (2001) developed a five dimensional model to study the impact of eutrophication on the ecosystem dynamics. Fu-Liu et al. (2001) proposed a GIS-based method of lake eutrophication assessment to study the spatial distribution of eutrophication conditions in the lake environment. Carpenter and co-workers (Carpenter et al., 1999a,b; Lisa et al., 2002; Peterson et al., 2003; Janssen and Carpenter, 1999; Ludwig et al., 2003) have developed several models (both single and multi dimensional, linear and nonlinear, discrete and continuous in nature) to represent the lake dynamics under various conditions and study the management of lake ecosystems from different perspectives. These models are developed to understand the regime shifts in the lake ecosystems and to derive more insight into possible outcomes of ecosystem management and provide with a tool box of approaches to be applied in the real life case studies. The dynamic equation representing the state of the lake (Carpenter et al., 1999a) is given by dP rP2 = l − sP + 2 dt m + P2

(1)

where P is the concentration of phosphorus in the water, l the rate of P input, s the rate of P loss and the term rP2 /(m2 + P 2 ) represents the phosphorus input

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

into the water due to recycling from sediments. The analysis developed based on the model (1) establishes existence of either one globally stable equilibrium point or three equilibrium points in the P-space. In the later case, two of them are stable and one is unstable. The unstable equilibrium is always positioned in between the two stable equilibria. Thus, all nontrivial solutions with initial values to the left (right) of the unstable equilibrium approach the stable equilibrium existing to the left (right) of the unstable equilibrium in the P-space. The region to the left (right) of the unstable equilibrium is termed as oligotrophy (eutrophy) regime. The unstable equilibrium point has been interpreted as a separatrix between the two regimes. The regime shifts, often observed in the lakes, are interpreted as the movement of the phosphorus (the state variable) from domain of attraction of one stable equilibrium to another (Carpenter et al., 1999a,b). Here, it is important to note that it always requires a sudden impulse or a perturbation of the state to shift from one regime to another. A continuous transition of the state from one side of the unstable equilibrium point to another side (through the unstable equilibrium point) is not possible owing to the nature of equilibrium points. Also, it is always possible to choose the parameters of the model (1) in such a way that the equilibrium point(s) of the system is (are) sufficiently far from zero or closer to zero. For example, if the input rate is high, then for certain values of s, r and m, all the equilibrium points of the model (1) will be sufficiently large. As a result, for such high concentrations of phosphorus, the lake may be eutrophic at all the three equilibrium states in reality. Thus, the unstable equilibrium point in the P-space need not always be the boundary between the oligotrophic and eutrophic regimes. For, even if the inflow rate is high, the rate of phosphorus loss could also be high enough to reduce the concentration of the phosphorus in the water to a low level. Hence, the lake could be oligotrophic in spite of high rate of input. Similarly, even low rate of input could result in driving the lake to eutrophy state if the phosphorus loss is very less from the water as this could result in increased concentration of phosphorus in the lake due to accumulation. Since it is not possible for the state to move from one side of the separatrix to another (for a given a set of parameter values in model (1)), the so called movement of state from one regime into another is not possible

117

if we accept the definitions of regimes as in Carpenter et al. (1999a,b). But practically these regime shifts are observed in the lakes and these shifts take place very quickly. This surprising behaviour in the lake ecosystems call for a new definition of regimes. It is appropriate to define the regime of the lake in terms of the concentration of phosphorus in various components like in water, in sediments and in other relevant biomass in the lake rather than defining it in terms of the influx of phosphorus and relative to the position of the unstable equilibrium point. The recycling process due to accumulated phosphorus in the sediments plays very important role in the dynamics of the lake. Hence, the phosphorus concentration in the sediment deserves separate attention. Hence, involving the dynamics of phosphorus concentration in the sediments would help to obtain an improved picture of the lake ecosystem. This implies that the lake ecosystem is to be represented using two state variables, viz., the phosphorus concentration in water and phosphorus concentration in the sediments. The oligotrophic and eutrophic regimes could be defined in terms of these two concentrations. Therefore, the total concentration of phosphorus in the lake at any time could be defined as the sum of the concentrations of the phosphorus in the water and in the sediments. In Carpenter et al. (1999a), Janssen and Carpenter (1999), the concentration of phosphorus in the mud received its deserving attention and a model has been proposed to represent the lake dynamics by involving the concentration of phosphorus in the mud in order to enhance the realism of the ecosystem. The equation representing the dynamics of the concentration of phosphorus in the mud has been termed as the mud equation and the corresponding variable is called as slow variable. The proposed difference equation model is as follows: Pt+1 = (1 − s − h)Pt + lt + rMt f(Pt ) Mt+1 = (1 − b)Mt + sPt − rMt f(Pt ) f(P) =

Pq mq + P q

where P and M represent the concentration of the phosphorus in the water and in mud, respectively, l the input rate, s and h the proportions of the P lost due to sedimentation and hydrologic outflow, r the proportions of M recycled to the water, b the proportion

118

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

of M buried permanently, m the P level at which the recycling is half maximal and the exponent q controls the steepness of the recycling curve. In the above model, the rate of loss of P to the mud is taken as a constant, s. It would be more realistic if this term is taken as M dependent. The recycling term is given by rMPq /(mq + P q ). This term is increasing with respect to both the state variables P and M. This implies that more the concentration of P in the water more will be the recycling from the sediments. It is well known that any solution offers resistance to dissolve more solute in it when it is nearing saturation. Following this principle it would be realistic to take a concave function of M in the place of s and take the recycling function in such a way that it is non increasing with respect to P and increasing with respect to M. Management of eutrophication for lakes and economic analysis of the eutrophied lakes have been studied by several authors based on the foresaid models (Carpenter et al., 1999a,b; Janssen and Carpenter, 1999; Kenneth et al., 2004; Maler et al., 2004; Brock et al., 2004; Xepapadeas et al., 2004) and some other linear difference models (Peterson et al., 2003). In fact, it was observed that certain behavior like cycling and the return of the lakes to eutrophy state in spite of controlling the input could not be explained by the existing models and a need for creating a new model was felt in order to explain the foresaid behavior of the lakes (Peterson et al., 2003). The model suggested in this article explains the cyclic behaviour (i.e., flipping between oligotrophic and eutrophic regimes) of the lakes, and also the reasons for the lake to return to eutrophic regime in spite of controlling the inputs. The analysis of the model identifies various control parameters in the system using which the state of the lake can be driven to remain in the required regime eventually. In this article, we model lake eutrofication using a simple two dimensional nonlinear coupled differential system. The aim is to capture the flips observed in the lakes and derive some strategies to control these flips so that the ecosystem can be managed efficiently. Though the dimension and complexity of the model can be increased by involving more state variables and nonlinear interactions, we stick to this minimal model as our aim is to capture the patterns observed in reality. Models are nothing but simple theories about the cause of observed patterns and use of simple hypothesis for explaining these patterns in nature. Non-

linear models of ecosystems dynamics that incorporates positive feed backs and multiple, internally reinforced states have considerable explanatory power. The fact that intriguing and complex patterns in nature can often be explained by simple models is seductive to many scientists. However, others feel instinctively repulsed by the minimal model approach. This split notoriously pronounced by biologists. Nonetheless, most people agree that scientific progress benefits from the combination of models and experiments (Scheffer, 1999). The model constructed in this article is designed to illustrate general patterns of system behaviour, rather than to make specific predictions. This model is usable, understandable and can easily be modified to accommodate unforeseen circumstances and new ideas. It is to be noted that even this simple model yields complex dynamics that mimic the observed patterns known from many case studies and discussed over and over in the literature. Although simple, this model provides ideas leading towards more realistic management strategies to control the ecosystem. Section wise split up of this article is as follows. In Section 2, the lake eutrophication has been modelled as a two dimensional coupled differential system. In Section 3, the model is analyzed to investigate the existence and nature of the equilibria. A few sufficiency conditions are derived for the existence of multiple equilibria. Based on the analysis, a few control parameters have been identified to help in restoring the lake. The analysis is interpreted in ecological perspective. The flipping behavior observed in the lakes has been explained. A solution to avoid these cycles is suggested. In Section 4, a necessity to increase the dimension of the model is stressed to explain the irreversibility properties of lakes exposed to eutrophication for longer duration. In Section 5, an illustrative model is considered by taking a specific absorption function and by assigning particular values to the parameters. Numerical simulations are performed and corresponding phase portraits are presented to illustrate the behavior of the system under various conditions. The reasons for cyclic behaviour, observed during lake management, is explained. Conclusions and modifications to increase the realism of the model are suggested in Section 6 for future work. Necessary proofs and mathematical analysis are given in the Appendix A.

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

2. The model We assume that the lake is represented by two state variables given by P and Q where P is the concentration of phosphorus in the water of the lake and Q is the concentration of phosphorus in the mud at the bottom of the lake. We have the following (observations) assumptions to assist us to model the lake dynamics. The concentration of phosphorus in the water increases due to input from the watershed and recycling from the sediments and it decreases due to out flow and absorption by zooplankton, plant and fish living in the lake. Similarly, increase in the concentration of phosphorus in the mud is due to absorption of phosphorus from water and decrease is due to recycling from the sediments. Let us call these assumptions as basic assumptions for future reference. Thus, we have the following coupled dynamic equations representing the lake dynamics. dP rQ2 = l − sP − Pf(Q) + 2 dt m + Q2

(2)

dQ rQ2 = Pf(Q) − 2 dt m + Q2

(3)

where l is the rate of phosphorus input from the watershed. The rate of P loss per unit time is represented by s. This loss process is due to out flow and sequestration in biomass of consumers or benthic plants. The function f represents the Q dependent rate of loss of phosphorus to mud due to sedimentation. Alternatively, f represents the phosphorus absorption function by the mud or mud response function to phosphorus. We assume that the rate of loss of phosphorus to mud due to sedimentation is positive when Q = 0 i.e.,f(0) > 0. This response increases until Q reaches a critical value at Q = a and monotonically decreases for Q > a with a zero at Q = K > 0, where K represents the maximum concentration of phosphorous the sediments in the lake can hold. This behavior of f represents that the absorption of phosphorus into the mud initially increases and reaches its maximum at a and there after the sediments offer resistance to the absorption and hence there is a continuous decrease in the absorption as the concentration of phosphorus in the mud exceeds a and there will be no absorption when Q becomes K i.e., when the mud is fully saturated with phosphorus. Hence, we

119

take the function f to be a concave function with a maximum in (0, K). The term rQ2 /(m2 + Q2 ) (an increasing function of Q) represents recycling of phosphorus from the sediments which clearly depends on the concentration of the phosphorus in the mud and also on the history of the phosphorus loading into the mud. r is the maximum rate of recycling of phosphorus from sedimentation. m is the half saturation value of the recycling function (Carpenter et al., 1999a). It is extremely important to observe that, for a given absorption function f(Q), the rate of change of Q decreases with decrease in P and increase in r. Thus, increase in recycling decreases the concentration of the phosphorus in the sediments. We have the following proposition. Proposition 1. All the solutions of the system (2, 3) which initiate in R2+ are uniformly bounded. The above result implies that, given the input rate and other parameter values in the ecosystem the concentrations P and Q will always remain finite and bounded. Hence, with out loss of generality we can assume that the state space of the system (2, 3) to be the unit square.

3. System analysis In this section, we shall investigate the existence and nature of the equilibria admitted by the system (2, 3) and also study the properties of the paths in the phase space. We give an intuitive representation for oligotrophic and eutrophic regimes in terms of the state variables and interpret the results in the ecosystem perspective. Clearly, an equilibrium of the system represented by (P ∗ , Q∗ ) satisfies the relations: P∗ =

l s

P ∗ f(Q∗ ) =

(4) rQ∗2 m2 + Q∗2

(5)

The number of equilibrium solutions admitted by (2, 3) are same as the number of solutions of (5). This number depends on the nature of f and values of the involved parameters. It is interesting to note that, irrespective of the number of equilibrium points admitted

120

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

by the system (2, 3), the P component is the same for all the equilibrium points and only the Q component differs. From the properties of the absorption function f(Q) and the recycling function rQ2 /(m2 + Q2 ) the following propositions can be easily established and hence we state them without proof. In fact, the equation (5) always admits odd number of solutions (one or three). Let g(Q) = rQ2 /(m2 + Q2 )P ∗ and let a be such that f  (a) = 0 with 0 < a < 1. Proposition 2. If g(Q) < f(Q) on (0, a] then the Eq. (5) has a unique solution belonging to (a, 1) and if g(Q) > f(Q) on (a, 1] then all solutions of (5) belong to (0, a]. Proposition 3. The Eq. (5) admits at most one solution in the interval (a, 1). It can be verified that, if the system admits three equilibrium points given by E1 , E2 and E3 then E1 and E3 are locally stable and E2 is a saddle (refer appendix for equilibrium analysis). From the dynamic nature of the equilibrium points and invariance of the system in the unit square, it can be observed that the stable manifold of E2 is the separatrix of the system which divides the unit square into two disjoint invariant sets each containing one of the locally stable interior equilibrium points (E1 or E3 ) in its interior. The two branches of the unstable manifold of E2 connect E2 to E1 and E3 (The stable and unstable manifolds of E2 are represented in Fig. 1). In fact, the invariant regions to which the stable equilibrium points belong are their respective basins of attraction. Thus, solutions initiating in the region to the left (right) of the separatrix will reach E1 (E3 ) eventually. Since we are representing the state of the lake by two variables we define the oligotrophy and eutrophy regimes only in terms of these two variables. For the reasons mentioned in Section 1, the separatrix of the system need not always represent the boundary between the oligotrophy and eutrophy regimes. Since these two states are determined depending on the concentrations P and Q, with out loss of generality, we can assume that there exists a continuous curve in the (P, Q) space which also divides the unit square into two connected sets each representing one of the foresaid regimes. Therefore, these two regimes

Fig. 1. This figure represents the stable and unstable manifolds of E2 . The broken arrows represent paths of the system. Note that paths with initial values on either side of the stable manifold approach the corresponding stable equilibrium lying in that side.

can be identified with reference to a smooth curve in the (P, Q) space which can be considered as the boundary between these two regimes. Let this curve be given by M(P, Q) = 0 which is continuous, satisfying M(0, 0) < 0 and M(1, 1) > 0. It is also realistic to assume that M(P, Q) is increasing with respect to both its arguments. Thus, M(P, Q) = 0 divides the unit square into two disjoint sets. The oligotrophy regime is characterized by M(P, Q) < 0 and eutrophy regime is characterized by M(P, Q) > 0. Refer Fig. 2 for characterization the regimes. Note that, for the considered system (2, 3), the P ∗ value of an equilibrium point is always a constant given by l/s. Thus, if M(l/s, 0) > 0 then all the equilibrium points of the system lie in the eutrophy state. Hence, irrespective of the initial concentrations and the number of admitted equilibrium points, the lake will eventually become eutrophic and it will remain in that state as long as M(l/s, 0) remains positive. This is due to the fact that the  - limit set of the system (2, 3) is completely contained in the eutrophic regime. If the (P ∗ , Q∗ ) is the unique equilibrium of the system satisfying M(P ∗ , Q∗ ) < 0(M(P ∗ , Q∗ ) > 0) then the eventual state of lake will be oligotrophy (eutrophy) (refer Figs. 3 and 4).

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

Fig. 2. This figure represents the oligotrophic and eutrophic regimes in the PQ-space.

From the above observations, we see that reduction in inputs (l) and increasing the out flux (s) would decrease the P component (i.e., l/s) in an equilibrium state of (2, 3). The decrease in the P component

121

Fig. 4. This figure represents the existence of a unique eutrophic globally stable equilibrium point and flow direction of solutions towards this equilibrium point, a situation where the lake will turn eutrophic even though it is oligotrophic initially.

also decreases the Q component (due to the dynamic relation between them). Further decrease in Q can be achieved by either reducing the absorbing capacity of the mud or by increasing recycling. Hence, some or all the above four controls can be used to reduce the total concentration of phosphorus in the lake. If the reduction in these two components yields an unique oligotrophic equilibrium point to the system then the lake will turn oligotrophic gradually. If the system (2, 3) admits multiple (three) equilibria (i.e (P ∗ , Q1 ), (P ∗ , Q2 ), (P ∗ , Q3 )) then we have the following cases. A : M(P ∗ , Q1 ) < 0, M(P ∗ , Q2 ) < 0, M(P ∗ , Q3 ) < 0 B : M(P ∗ , Q1 ) < 0, M(P ∗ , Q2 ) < 0, M(P ∗ , Q3 ) > 0 C : M(P ∗ , Q1 ) < 0, M(P ∗ , Q2 ) > 0, M(P ∗ , Q3 ) > 0 D : M(P ∗ , Q1 ) > 0, M(P ∗ , Q2 ) > 0, M(P ∗ , Q3 ) > 0

Fig. 3. This figure represents the existence of a unique oligotrophic globally stable equilibrium point and flow direction of solutions towards this equilibrium point, a situation where the lake will turn oligotrophic even though it is eutrophic initially.

In the cases A and D, all solutions of the system will eventually enter the oligotrophy regime and eutrophy regime, respectively. In the cases B and C, we may have situation where the state is oligotrophic (eutrophic) initially, becomes eutrophic (oligotrophic) for some time and later again returns to oligotrophy (eutrophy) regime and remain there for all future times. This illustrates the regime shifts where the lake flips

122

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

4. Reversible and irreversible lakes Now, following definitions of the reversible and irreversible lakes, let us assume that the input of phosphorus from watershed is completely shut off (i.e., l = 0 according to our basic assumptions). Hence, the dynamics of the lake get modified to rQ2 dP = −sP − Pf(Q) + 2 dt m + Q2

(6)

dQ rQ2 = Pf(Q) − 2 dt m + Q2

(7)

We have the following proposition. Proposition 4. All the solutions of the system (6, 7) approach (0, 0) asymptotically. Fig. 5. This figure represents the case B in the text where the system (2, 3) admits three equilibrium points given by E1 , E2 and E3 . The broken line arrows represent paths of the system with different initial values. Observe that the solutions lying in the region of attraction of E1 initiate in the eutrophic regime, change the regime in a finite time and reach the oligotrophic stable equilibrium point E1 eventually. Those initiating in the region of attraction of E3 undergo flips from eutrophic to oligotrophic and then eutrophic again and finally reach E3 , the eutrophic stable equilibrium point.

from one regime to another with out any interventions. Case B is represented in Fig. 5. In all cases, the eventual value (E1 or E3 ) of the solutions will depend on its initial position relative to the separatrix of the system. The only way to change a limit of a solution to E3 (E1 ) with initial value in the domain of attraction of E1 (E3 ) is to impulsively push the state of the system across the separatrix of E2 at some time t ≥ 0 into the region of attraction of the required limiting value. Ecologically, this means that the management actions like decreasing inflow (l), increasing flushing (s) (both these actions decrease the P component) and manipulations like dredging, etc., which decrease the Q component (NRC, 1992; Cooke et al., 1993) should be undertaken in such a way that the values of P and Q get reduced and shift the state, i.e., (P, Q) to the left of the separatrix. Once the state is pushed to the left of the separatrix, then automatically the lake will return to oligotrophic state provided there are no externalities which affect the parameters of the system (refer Fig. 5).

The above result implies that, if the input of phosphorus from external activities is completely stopped then the concentrations of phosphorus in water as well as in mud reduce to zero gradually. From the above analysis, we infer the following. If the state dynamics is governed by only the basic assumptions, then we find that it is always possible to get rid of the phosphorus in the lake by arresting the inflow of phosphorus from the watershed into the lake, in which case, all the phosphorus in the lake, i.e., the phosphorus in the water as well as in the mud will get washed away and the lake will be free from phosphorus eventually. But there were observations where the state of the system did not reach a phosphorus free state even after the inflow of phosphorus is completely stopped. This was observed in lakes which were exposed to eutrophication for very long time (Carpenter, 2003). The model (6, 7) does not reflect this behavior. This indicates that there are a few more vital interactions in the lake which our model does not represent. This is contribution of phosphorus to the water due to other recycling processes like excretion or egestion by some consumers like fish in the lake. This adds an additional term to the P dynamic equation which represents the phosphorus recycling due to the consumers in the water. Prolonged exposure to eutrophication force the consumers in the water to absorb phosphorus into their bodies and get accumulated in the organisms like fish (Carpenter et al., 1999a). This has considerable

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

impact on the recycling process in the lake and also on the life cycle of the consumers (Bayraktaroglu et al., 2003). Recycling of phosphorous is critical for maintaining plant production of lakes. In many lakes, primary producers obtain more of their phosphorous from recycling than from inputs. Phosphorous is recycled by several mechanisms in the lakes. Recycling by animals (excretion or egestion) can supply a significant fraction of the phosphorous demand for phytoplankton growth during summer. Both the herbivorous zooplankton and fishes are potentially important sources of phosphorous to phytoplankton (Carpenter, 2003). In order to incorporate the phosphorus recycling due to consumers in the lakes, we assume that the total biomass of consumers in the lake is constant and the recycling rate of phosphorus due to this biomass is proportional to the biomass. Definitely, this is a crude assumption as the consumers biomass is dynamic. But, we assume this only to show that reduction in the inflow need not result in elimination of phosphorus in the lakes if the recycling due to the consumers is incorporated in our model. To get a better understanding of the lakes exposed to eutrophication for longer duration, i.e., irreversible lakes, it is essential to incorporate the dynamics of the consumers also in the model by taking the biomass of the consumers as an additional state variable. Under the above assumption (i.e., the biomass of the consumers is a constant), the term l in the system (2, 3) can be considered as sum of two terms l1 and l2 where l1 is the rate of inflow of phosphorus from watershed and l2 is the phosphorus recycling rate due to the biomass of consumers in the lake. Now, it is easy to observe that even if the inflow of the phosphorus from watershed is completely restricted, i.e., l1 = 0 the state of the system need not return to a phosphorus free state as l2 is positive and the system still admits nontrivial equilibrium point(s) with P component as l1 /s. The change in the dynamics of the system would be the consequence of reduction of input from watershed. This may result in a change in the regime of the lake depending on the magnitude of reduction in the input and the position of equilibrium point(s) prior to reduction. If the original system admits a unique eutrophic equilibrium, then after reduction, the system may still admit one eutrophic equilibrium state (if the reduction is less), or may admit three equilibrium

123

points (if the reduction is moderate) or a unique oligotrophic equilibrium (if the reduction is large enough). Hence, the eventual state of the lake after reduction of the input will depend on the position of the equilibrium(s) after the reduction and current concentrations of phosphorus in the lake. Thus, if the reduction in the inflow (i.e., making l1 zero) results in an unique oligotrophic equilibrium point for the lake then such lakes are reversible lakes. If the lake ecosystem admits three equilibrium points with one stable equilibrium lying in the oligotrophic regime and the other stable equilibrium in the eutrophic regime, then such lakes are hysteric lakes as they require extreme reduction of phosphorus inputs and some temporary interventions to drive the state of the lake to the region of attraction of the stable oligotrophic equilibrium. If the reduction in the inputs fails to create any oligotrophic equilibrium point for the system then such lakes are irreversible lakes.

5. An illustrative example In this section, we consider the lake model with a simple absorbing function and study the behavior of the paths using numerical simulation. We take the absorbing function f(Q) to be η(Q + v)(1 − Q), −1 < v < 1. Clearly, f has a unique maximum at a = (1−v)/2 in (0, 1) and it satisfies all the necessary conditions like f(0) > 0, f  > 0(< 0) in [0, a)((a, 1]). We assume the following values for some of the parameters. l = 0.2,

r = 0.35,

η = 4.5,

v = 0.08.

m = 0.1,

We take three values for s given by 0.8, 0.76, 0.55 (representing high, moderate and relatively low flushing rates, respectively) and study the system for these three cases. Let M(P, Q) = P 2 +Q2 −0.25 = 0 represent the boundary between the oligotrophy and eutrophy regimes in the phase space. This boundary, which connects the points (0, 0.5) and (0.5, 0), is shown in all the Figs. 6 to 10. Practically, this boundary curve can be guessed only by experiments and observations. Here, we are taking the above form for the boundary for the sake of illustration only.

124

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

Fig. 6. Figure showing existence of a unique globally stable oligotrophic equilibrium point (0.1083, 0.25) for the chosen parameters. Note that the high flushing rate (s = 0.8) reduces the concentration of phosphorus in the water as well as in sediments and turn the lake to oligotrophic state even though it is in the eutrophic state initially.

The considered system with s = 0.8 or s = 0.55 admits a unique equilibrium point (0.1083, 0.25) (oligotrophic) or (0.3636, 0.7455) (eutrophic), respectively. The phase portrait of these two cases

are represented in Figs. 6 and 7 respectively. If s = 0.76 the system admits three equilibrium points. These equilibrium points are given by (0.1232, 0.2632), (0.2999, 0.2632) and (0.5371, 0.2632). Let

Fig. 7. Figure showing existence of a unique globally stable eutrophic equilibrium point (0.3636, 0.7455) for the chosen parameters. Note that the relatively low flushing rate (s = 0.55) is unable to reduce the concentration of phosphorus in the water and in sediments sufficiently which turns the lake to eutrophic state even though it is in the oligotrophic state initially.

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

125

Fig. 8. Figure showing existence of three equilibrium points E1 , E2 and E3 for the intermediate flushing rate s = 0.76. The broken line arrows indicate the direction in which the solutions approach the stable equilibrium points E1 and E3 . It clearly illustrates that the limiting value of a path depends on its initial position relative to the stable manifold of E2 .

us denote these points as E1 , E2 and E3 respectively. The phase portrait of this system is represented in Fig. 8. Clearly, E1 and E2 lie in the oligotrophic region and E3 lies in eutrophic region. Note that two solutions with initial values (0.115, 0.6), (0.12, 0.6) move towards E2 for some time and eventually approach the equilibrium points E1 and E3 respectively (refer Fig. 8), under the influence of the unstable manifold of E2 , which connects the equilibrium point E2 to E1 and E3 . This clearly shows that one branch of the stable manifold of E2 lies between the foresaid solutions and the other branch lies between the solutions initiating at the points (0.46, 0.0), (0.465, 0.0) as these two solutions also move towards E2 initially and eventually approach E1 and E3 respectively. This stable manifold of E2 is the boundary between the regions of attraction of E1 and E3 . From Fig. 8, we clearly observe that all solutions initiating in the region to the left (right) of this manifold end up in E1 (E3 ). It has been observed that the lake reverts from eutrophic state to oligotrophic state when the managers reduce the P loading, which alters the state of the lake to oligotrophic abruptly and gradually it again returns to eutrophic state and this cycle continues (Peterson et al., 2003). The model (1) does not explain this flipping behaviour for the reasons mentioned in Section 1.

But this behaviour can be explained by the considered model (2, 3). These cycles could take place whenever the system (2, 3) admits unique eutrophic equilibrium point or three equilibrium points. Consider the above system with s = 0.76 and refer to the phase diagram of this system presented in Fig. 9. Observe that the path with initial value (0.55, 0.0) in the eutrophic regime enters the oligotrophic regime after some time, remains oligotrophic for some time as it approaches E3 it becomes eutrophic again. This situation explains the cycling that is taking place in the lake management. For the reduction in P and Q, due the management action, could be driving the state to the oligotrophic region in the region of attraction of E3 (which lies to the right of the separatrix of the system). In which case, the state alters to oligotrophic abruptly and will turn eutrophic gradually as the state moves out of the oligotrophic region in a finite time. Similar behavior could be observed from paths with initial values (0.2, 0.6) and (0.12, 0.6). Thus, to break away from these cycles we must ensure that the management action drives the state of the system from one side of the separatrix to the other so that the solution will move into the region of attraction of the oligotrophic equilibrium point which will automatically drive the lake to oligotrophic state eventually.

126

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

Fig. 9. This figure illustrates the regime shifts observed in lakes. Note that the paths with initial values (0.55, 0), (0.12, 0.6) and (0.2, 0.6) are initially in the eutrophic regime, which gradually move into oligotrophic regime and after some time they move into eutrophic regime to reach E3 . It is easy to note that this flipping behaviour can be observed for some paths initiating in the region of attraction of E3 . This also indicates that the lake management action should drive the state of the lake to the region to the left of the stable manifold of E2 to turn the lake to oligotrophy failing which the lake may turn eutrophic after some time.

Similarly, if we consider the Fig. 7, representing the phase portrait of the system with s = 0.55, it is easy to observe that as long as the intervention does not alter the regime of the existence of the equilib-

rium point, all the efforts in bringing the lake to oligotrophic state eventually will fail as all the paths will cross the boundary at some time and enter the eutrophic regime and will remain there. Thus, a man-

Fig. 10. This figure represents the eventual elimination of phosphorous in the lake if l = l1 + l2 = 0 i.e all the inputs except recycling from sediments are stopped.

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

agement will succeed in reverting an eutrophic lake to oligotrophic through any action which sufficiently decreases the concentration of phosphorus in the water or reduces the phosphorus concentration in the sediments or both so that the altered system due to management action admits a unique equilibrium point in the oligotrophic regime. In case the altered system admits multiple equilibria, then the management action must ensure that one of its stable equilibrium points lie in the oligotrophic regime and the state is impulsively pushed into the region of attraction of the oligotrophic equilibrium. Finally, the paths of the system with l = l1 +l2 = 0 are presented in Fig. 10, where we clearly see the eventual elimination of the phosphorus in the lake.

6. Conclusions In this article, the dynamics of lake exposed to eutrophication is modelled using a two dimensional coupled differential system involving the concentration of phosphorus in the water and in the sediments as two state variables. At any time, the lake is represented by these two variables. The model explains some key observations made with respect to the eutrophied lakes which are not explained by the existing models, like causes for the regime shifts and the cyclic behaviour observed during the lake management. The model indicates that the lakes that have been eutrophied for only a short time, where phosphorus contribution due to recycling from consumers is less, may be restored by increasing flushing, decreasing inflow and increasing recycling. It is important to note that it is not possible to restore a eutrophied lake by either decreasing recycling or accelerating sedimentation. These interventions may drive the state to oligotrophy immediately but soon the state will turn eutrophic as these actions tend to increase the concentration of the phosphorus in the sediments. Therefore, the treatment should be to increase recycling, decrease the sedimentation and increase flushing so that the phosphorus recycled into the water gets washed away due to increased out flow. Incase of irreversible lakes, phosphorus recycling due to excretion or egestion of consumers in the lake and other processes contributing phosphorus to water need to be given attention to explain the irreversible behavior. In this case, the model indicates the neces-

127

sity to increase the dimension of the model to three or more by including the dynamics of the consumers like fish, zooplankton, plants etc., in order to explain the observed behavior like increased fish kills (Smith, 1998), algal blooms etc. Optimal management policies can be worked out to find an optimal strategy which is cost effective and also drives the system to the required state based on the considered model (2, 3). These management studies are likely to yield results which will be an improvement over the existing ones. There are several modifications, apart from increasing the dimension of the model, that can be made to the model (2, 3) to make it more realistic. Firstly, the steepness of the recycling curve can be made q (to be fixed based on experimental evidence) with q ≥ 2. The maximum recycling value r can be taken as a function of P, say r(P), a positive decreasing function on R+ , which represents the resistance offered by water to dissolve more phosphorus coming from the sediments. In Peterson et al. (2003), it was mentioned that the recycling starts only after the concentration of the phosphorus in the sediments exceeds certain critical value. This can also be incorporated in the dynamical system. In this case, the positive quadrant (the phase space) has to be divided into two regions, as in Srinivasu and Gayatri (2004), with dynamic equations being different in each region. In Carpenter et al. (1999b), Janssen and Carpenter (1999), it was assumed that certain proportion of phosphorus in the mud gets buried. This process was represented by a term −bMt in the M-dynamic equation (Section 1) which represents decrease of phosphorus concentration in mud. Here, phosphorus is not lost from the ecosystem. Therefore, this process can be more realisticalyy represented in the model by assuming that certain portion of the phosphorus in the mud plays a passive role and does not contribute to the recycling even though it contributes to the total concentration in the mud. Thus, the lake dynamics has immense potential for further research which will help restoration of lakes and also its bio-economics.

Acknowledgements The lecture given by Prof. S.R. Carpenter at the First School on Ecological Economics, held at Abdus

128

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

Salam International Centre for Theoretical Physics, Trieste, Italy during 27 January–28 February 2003 has motivated the author to write this article. The author thanks Prof. Carpenter for the illuminating lecture. The author also thanks Prof. Sir Partha Dasgupta, Prof. Karl-Goran Maler and Matteo Marsili for giving him the opportunity to participate in the foresaid school and for their encouragement. The author is extremely grateful to the referees and Editor-in-Chief for their constructive criticism and helpful comments.

Appendix A

Proof of Proposition 1. We define a function W = P + Q. The time derivative of W along the solutions of the system (2, 3) is given by dW = l − sP. dt For each η > 0, we have dW + ηW = l − sP + η(P + Q). dt ≤ l + ηK + (η − s)P. Now, if we choose η to be smaller than s, we have

Applying theory of differential inequality (Birkhoff and Rota, 1982), we obtain 0 < W(P, Q) 
l + K (1 − e−ηt ) + W(P(0), Q(0))e−ηt < η

0 < W(P, Q) = P + Q <

B = (P, Q) ∈ R2+ : P + Q
  l = + K + ζ for any ζ > 0 , η for all t ≥ T , where T depends on the initial values (P(0), Q(0)). Thus, the set B is invariant set which contains the Ω-limit set of all the paths of the system (2, 3) that initiate in the positive quadrant. 䊐 Equilibrium analysis: The following matrix defines Jacobian of the system (2, 3) evaluated at an equilibrium point (P ∗ , Q∗ ).   2rm2 Q∗ ∗ ∗  ∗  −s − f(Q ) −P f (Q ) + (m2 + Q∗2 )2      2 Q∗   2rm ∗ ∗  ∗ f(Q ) P f (Q ) − 2 ∗2 2 (m + Q ) The associated characteristic equation is given by  2rm2 Q∗ ∗ λ2 − λ P ∗ f  (Q∗ ) − 2 − s − f(Q ) (m + Q∗2 )2 
2rm2 Q∗ ∗  ∗ =0 − s P f (Q ) − 2 (m + Q∗2 )2 From the above equation, it can be observed that if f  (Q∗ ) ≤ 0 then Q∗ ∈ [a, 1) and the equilibrium point is locally stable. On the other hand, if f  (Q∗ ) > 0 then Q∗ ∈ (0, a) and the equilibrium is locally stable if f  (Q∗ ) ≤

dW + ηW < l + ηK. dt

Thus, at t → ∞, we have






l +K η



Hence, all the solutions (P(t), Q(t)) of (2, 3) that initiate in the positive quadrant are confined in the region

2P ∗ rm2 Q∗ (m2 + Q∗2 )2

else it is a saddle. Thus, we observe that whenever the system (2, 3) admits three equilibrium points given by E1 = (P ∗ , Q1 ), E2 = (P ∗ , Q2 ) and E3 = (P ∗ , Q3 ) with Q1 < Q2 < Q3 then we have either {Q1 , Q2 , Q3 } ∈ (0, a] or {Q1 , Q2 } ∈ (0, a] and Q3 ∈ (a, 1). Also it can be established using Bendixon du-Lac criterian (Paul, 1994) that the system does not admit any periodic solutions if P ∗ f  (Q) −

2rm2 Q − s − f(Q) < 0 (m2 + Q2 )2

(A.1)

is satisfied for all Q ∈ [0, a]. Thus, if the the system admits a unique interior equilibrium, then it is globally asymptotically stable if (8) is satisfied. Henceforth, we

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

shall assume that (8) is satisfied throughout the article. Simple calculus can be used to show that if the system (2, 3) admits a unique equilibrium then it is locally asymptotically stable and if it admits multiple equilibrium points then locally stable and saddle natured equilibrium points occur alternatively with the first one being locally stable. Since all the equilibrium points have the same P component (l/s), we order the equilibrium points as per the Q component. Proof of Proposition 4. Observe that the system (6, 7) admits only the trivial equilibrium, i.e., (0, 0). It can be verified that the eigen values associated with this equilibrium point are −(s + f(0)) and 0. Hence, it is locally stable. In fact, this equilibrium point is globally stable. For, consider the following Liapunov function defined on the unit square V(P, Q) = P + Q. The time derivative of V along the system (6, 7) is given by dV dP dQ = + = −sP < 0. dt dt dt This implies that all nontrivial solutions of the system 䊐 (6, 7) approach (0, 0) eventually.

References Bayraktaroglu, E., Legovic, T., Velasquez, Z.R., Cruzado, A., 2003. Diatom Thalassiosira weissflogii in oligotrophic versus eutrophic culture: models and ultrastructure. Ecol. Modell. 170, 237–243. Birkhoff, G., Rota, G.C., 1982. Ordinary Differential Equations. Ginn, Boston. Brock, W.A., Maler, K.G., Perrings, C.A. Resilience and Sustainability: The Economic Analysis of Dynamic Systems. Beijer Discussion Paper 133. Carpenter, S.R., 2003. Regime Shifts in Lake Ecosystems: Pattern and Variation. Volume 15, Excellence in Ecology Series. Ecology Institute, Oldendorf/Luhe, Germany. Carpenter, S.R., Caraco, N.F., Correll, D.L., Howarth, R.W., Sharpley, A.N., Smith, V.H., 1998. Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecol. Appl. 8, 559–568. Carpenter, S.R., Kraft, C.E., Wright, R., He, X., Soranno, P., Hodgson, J.R., 1992. Resilience and resistance of a lake phosphorus cycle before and after food web manipulation. Am. Nat. 140, 781–798.

129

Carpenter, S.R., Ludwig, D., Brock, W.A., 1999a. Management of eutrofication for lakes subject to potentially irreversible change. Ecol. Appl. 9 (3), 751–771. Carpenter, S., Brock, W., Hanson, P., 1999b. Ecoligocal and social dynamics in simple models of ecosystem management. Conserv. Ecol. 3 (2), 4 [online] URL:http://www.consecol.org/ vol3/iss2/art4. Christopher, P.B., Childers, D.L., Dong, Q., Jones, R.D., 2000. Simulation of periphyton phosphorus dynamics in Everglades National Park. Ecol. Modell. 134, 103–115. Cooke, G.D., Welch, E.B., Peterson, S.A., Newroth, P.R., 1993. Restoration and Management of Lakes and Reservoirs. Lewis Publications, Boca Raton, FL, USA. Elser, J.J., Chrzanowski, T.H., Sterner, R.W., Mills, K.H., 1998. Stoichiometric constraints on food web dynamics: a whole lake experiment on the Canadian Shield. Ecosystems 1, 120–136. Fu-Liu, X., Tao, S., Dawson, R.W., Li, B.-G., 2001. A GIS-based method of lake eutrophication assessment. Ecol. Modell. 144, 231–244. Ichiro, Aoki., 1995. Flow-indices characterizing eutrophication in lake-ecosystems. Ecol. Modell. 82, 225–232. Janssen, M.A., Carpenter, S.R., 1999. Managing the resilience of lakes: a multi-agent modeling approach. Consrev. Ecol. 3 (2), 15, [online] URL: http://www.consecol.org/vol3/iss2/art15. Kenneth A., Daily, G., Dasgupta, P., Levin, S., Mäler, K., Maskin, E., Starrett, D., Sterner T., Tietenberg, T., Managing Ecosystem Resources, Beijer Discussion Paper 122. Kohei, Y., Nakajima, H., 2002. Catastrophic transition in vertical distributions of phytoplankton: alternative equilibria in a water column. J. Theor. Biol. 216, 397–408. Lisa, D., Graeme, C., Cumming, S., Crapeneter, S.R., 2002. Multiple states river and lake ecosystems. Philos. Trans. R. Soc. Lond. B 357, 635–645. Ludwig, D., Carpenter, S., Brock, W., 2003. Optimal phosphorous loading for a potentially eutrophic lake. Ecol. Appl. 13 (4), 1135–1152. Maler, K.G., Xepapadeas A., de Zeeuw, A., The Economics of Shallow Lakes. Beijer Discussion Paper 131. NRC (National Research Council), 1992. Restoration of Aquatic Ecosystems: Science, Technology and Public Policy. National Academy Press, Washington, DC, USA. Nurnberg, G.K., 1984. Prediction of internal phosphorus load in lakes with anoxic hypolimnia. Liminol. Oceanogr. 29, 135–145. Paul G., 1994. Stability, Instability and Chaos: An Introduction to the Theory of Nonlinear Differential Equations. Cambridge University Press. Peterson, G.D., Carpenter, S.R., Brock, W.A., 2003. Uncertainty and the management of multi-state ecosystems: an apparently rational route to collapse. Ecology 84 (6), 1403–1411. Santanu, R., Berec, L., Straskraba, M., Jorgensen, S.E., 2001. Optimization of exergy and implications of body sizes of phytoplankton and zooplankton in an aquatic ecosystem model. Ecol. Modell. 140, 219–234. Sas, H., 1989. Lake Restoration by Reduction of Nutrient Loading: Expectations, Experiences, Extrapolations. Academia Verlag, Richarz, St. Augustin, Germany. Scheffer, M., 1997. Ecology of Shallow Lakes. Kluwer Academic Publishers, Dordrecht, The Netherlands.

130

P.D.N. Srinivasu / Ecological Modelling 179 (2004) 115–130

Scheffer, M., 1999. Searching explanations of nature in the mirror world of Math. Conserv. Ecol. 3 (2), 11, [online] URL: http://www.consecol.org/vol13/iss2/art11. Scheffer, M., Hosper, H., Meijer, M.L., Moss, B., Jeppesen, E., 1993. Alternative equilibria in shallow lakes. Trends Ecol. Evol. 8, 275–279. Scheffer, M., Jeppensen, E., 1998. Alternative stable states. In: Jeppensen, E., Sondergaard, M., Kristoffersen, K. (Eds.), Strcturing Role of Submerged Macrophates in Lakes, vol. 131. Springer-Verlag, New York, USA, pp. 397– 406. Schindler, D.E., Eby, L.A., 1997. Stoichiometry of fishes and their prey: implications for nutrient recycling. Ecology 78, 1816– 1831. Schindler, D.W., 1977. Evolution of phosphorus limitation in lakes. Science 195, 260–262.

Schindler, D.W., Hesslein, R.H., Tuner, M.A., 1987. Exchange of nutrients between sediments and water after 15 years of experimental eutrophication. Can. J. Fish. Aquat. Sci. 44, 26– 33. Smith, V.H., 1998. Cultural eutrophication of inland, eustuarine and coastal waters. In: Pace, M.L., Groffmann, P.M. (Eds.), Successes, Limitations and Frontiers of Ecosystem Science. Springer-Verlag, New York, USA, pp. 7–49. Soranno, P.A., Carpenter, S.R., Lathrop, R.C., 1997. Internal phosphorus loading in lake Mendota: response to external loads and weather. Can. J. Fish. Aquat. Sci. 54, 1883–1983. Srinivasu, P.D.N., Gayatri, I.L., Influence of prey reserve capacity on Predator-Prey dynamics, submitted for publication. Xepapadeas, A., de Zeeuw, A., Maler, K.G., Feed Back Equilibria for a Class of Nonlinear Differential Games in Resource Economics. Beijer Discussion Paper 152.