Lake Restoration in Terms of Ecological Resilience: a Numerical Study of Biomanipulations under Bistable Conditions

An abstract version of the comprehensive aquatic simulation model (CASM) is found to exhibit bistability under intermediate loading of nutrient input, supporting the alternative-stable-states theory and field observations for shallow lakes. Our simulations of biomanipulations under the bistable conditions reveal that a reduction in the abundance of zooplanktivorous fish cannot switch the system from a turbid to a clear state. Rather, a direct reduction of phytoplankton and detritus was found to be most effective to make this switch in the present model. These results imply that multiple manipulations may be effective for practical restorations of lakes. We discuss the present results of biomanipulations in terms of ecological resilience in multivariable systems or natural systems.


INTRODUCTION
Several mathematical models of ecosystems have been proposed and used for understanding complex aspects of ecosystems (Walker et al. 1981, DeAngelis et al. 1989, Hulot et al. 2000, Scheffer 2001).Among these models, the comprehensive aquatic simulation model (CASM) (DeAngelis et al. 1989) has been found to be useful in applications for several aquatic ecosystems, and has been used for ecological risk assessment (Bartel et al. 1999, Naito et al. 2002).A simplified, more abstract version of CASM introduced by DeAngelis et al. (1989), abbreviated as abs.-CASM hereafter, is more suitable for a theoretical analysis; it was previously employed for studying stability of ecosystems (DeAngelis et al. 1989).Abs.-CASM retains important characteristics of aquatic ecosystems: material circulation, species interactions among three trophic levels, and the effects of human activities on aquatic ecosystems.DeAngelis et al. (1989) have analyzed abs.-CASM in terms of resilience, sensu Pimm (1984), i.e., the return time following perturbations ("engineering resilience," Holling 1996).For their analyses, monostable states were examined as a function of nutrient loading.
Many ecosystems, including shallow lakes, are known to exhibit bistable states or alternative stable states (Scheffer et al. 2001).Scheffer (2001) and Beisner et al. (2003) have shown bistable behavior in two-or three-variable aquatic models.We show that abs.-CASM also exhibits bistable behavior under a realistic range of parameter values, and that it is useful to study the properties of alternative stable states.
Ecological resilience, below just referred to as resilience, is introduced to describe the properties of alternative stable states in ecosystems (Holling 1973).Resilience concerns large and drastic changes in ecosystems.It characterizes "the amount of change the system can undergo and still retain the same controls on function and structure" (Holling 1973(Holling , 1996)), and is used as a novel concept for ecosystem management (Scheffer et al. 2001, Bellwood et al. 2004, Carpenter and Brock 2004).Resilience is usually illustrated by using a sigmoidal bifurcation curve for a simple one-variable system, as shown in Fig. 1.The horizontal axis is an environmental factor, i.e., control parameter, such Ecology and Society 10(2): 3 http://www.ecologyandsociety.org/vol10/iss2/art3/as nutrient loading in lakes, and the vertical axis is a state variable such as the biomass of algae.Resilience is defined by the distance of the unstable state (dotted line) from the stable states (solid lines) in the direction of the vertical axis (Scheffer et al. 2001).Two kinds of resilience are shown in Fig. 1: resilience of a clear state (bottom arrow) and of a turbid state (top arrow), if we consider the case of eutrophication of lakes.The magnitude of resilience changes according to changes in the environmental factor.A graphical representation such as Fig. 1 makes it easy to understand the concept, however, it should be noted that resilience cannot always be represented by a single number (distance) such as in the bifurcation curve shown in Fig. 1.For multivariable systems or natural systems, several complications arise.Because resilience is embedded in a multivariable-parameter space for such systems, the separatrix, which denotes the points of no return in the state space, may be a highly complex object.We will come back to this point when discussing the problem of choosing a graphical representation of resilience in relation to biomanipulation.
Biomanipulation was introduced as a lakerestoration technique that focused on trophiccascade top-down effects (Shapiro et al. 1975).Much work has been done to confirm the effects of biomanipulation (results are summarized in Jeppesen et al. 1990b, Demelo et al. 1992, Meijer et al. 1994, Hansson et al. 1998, Drenner and Hambright 1999).It is now thought that the interaction of pelagic food chain trophic levels (fish-zooplankton-algae) is not the only mechanism involved in biomanipulation, and that the reduction of benthivorous fish that suspend sediments is also a key process (Hansson et al. 1998, Drenner andHambright 1999).For a long-term effect of biomanipulation, the external nutrient (phosphorus) loading should be reduced below 0.5-2.0g m -2 year -1 (Jeppesen et al. 1990b), suggesting the existence of alternative stable states around that level of nutrient loading in real lakes.Our simulations of biomanipulation by manipulating one (or two) state variable(s) agree with these findings, and give insights into how the concept of resilience can be used in lake restoration.

Mathematical Model
We used the five-variable model of a three-level food chain proposed by DeAngelis et al. (1989) given by (1) where N is the amount of dissolved nutrient in the system, X is the biomass of phytoplankton, Y is the biomass of zooplankton, Z is the biomass of zooplanktivorous fish, and D is the detrital biomass.The parameter I N represents the input rate of the limiting nutrient from the environment.Here we assume that the limiting nutrient is phosphorus, however, it may be nitrogen depending on the N:P ratio (e.g., Takamura et al. 2003).The other parameters are: r N loss rate of the nutrient, γ ratio of nutrient mass to biomass, r 1 maximum growth rate of phytoplankton, k i half-saturation constant of nutrient (i = 1) or biomass (i = 2, 3), d i death rate of organisms (i = 1, 2, 3) or decomposition rate of detritus (i = 4), e i removal rate of organisms (i = 1, 2, 3) or detritus (i = 4), f i feeding rates, and η assimilation efficiency.Here we mention two characteristic points of the model: (i) the uptake kinetics of the nutrient (N) are described by a Michaelis-Menten-Monod function; (ii) both consumption of phytoplankton (X) by zooplankton (Y) and consumption of zooplankton (Y) by fish (Z) are modeled by a Holling type-III function.

Stability Analysis
A numerical linear stability analysis for the fivevariable model was carried out with the input rate (I N ) of nutrient being used as a control parameter.http://www.ecologyandsociety.org/vol10/iss2/art3/Fig. 1.Schematic representation of resilience using a simple bifurcation curve for a one-variable system exhibiting alternative stable states.Solid curves and dotted curve represent stable states and unstable state, respectively.The distance between the stable states and the unstable state represents resilience.Larger changes in the state variable than the resilience cause a state-shift in the ecosystem in the bistable region.
The other parameters were set to the fixed values, as listed in Table 1, with reference to the values reported by DeAngelis et al. (1989).

Numerical Simulation
Simulations of biomanipulation were carried out in a range of I N where there were two stable steady states and one unstable steady state.First, the differential equations were solved numerically under conditions where the initial values were set close to the stable steady states.Second, when the solutions reached the stable steady state, one variable was (or two variables were) manipulated to take values away from the steady state while the other variables remained unchanged, and the differential equations were solved again up to the point where the system reaches a steady state.Finally, the model output was assessed to determine  (DeAngelis et al. 1989), c) a long-term effect of biomanipulation can be expected below this level in real lakes (Jeppesen et al. 1990b), d) the half-saturation constant for phosphorus, e) forage fish, f) omnivorous fish, g) phosphorus composition of algal cells, h) zooplanktons.

Parameters
Values Units Meaning of the parameters Ref.

Stability Analysis
The local stability of the model system is shown in Fig. 2. Bistability was found for the variables N, X, Z, and D in the range of the input rate of nutrient I N from 0.5-1.34mg m -2 day -1 under the conditions listed in Table 1.However, the value of the variable Y was the same for the stable and the unstable points.
The steady-state solutions for the second-top trophic level can be shown to be generally constant as a function of nutrient loading in linear food chain models (DeAngelis et al. 1989, Hulot et al. 2000, DeAngelis 2003).
Similar bistable behavior could be found for different values of the other parameters, for instance, by changing one of the parameter values in the following range k 2 = 5.7-7.2,r 1 = 0.14-0.16,and f 1 = 6.2-6.5, while leaving the other parameters at the original values (DeAngelis et al. 1989).With increasing nutrient input rate (I N ), the amount of nutrient (N) in the system increased in a linear way, the biomass of zooplanktivorous fish (Z) increased gradually, and the biomass of phytoplankton (X) and detritus (D) increased slightly.When I N crossed the bifurcation point, N fell to almost zero, Z increased by nearly double, and X and D abruptly increased more than fourfold.We identify the stable states before and after crossing the bifurcation point with increasing I N (Fig. 2), the nutrient input, with the clear and turbid states of lakes, respectively.

Simulation of Food-Web Manipulation: Relaxation Dynamics
We manipulated each one of the five variables over the whole range of the bistable region, and found that only manipulations of X and D succeeded, whereas the others failed to switch the states.Figure 3a, b shows examples of a successful and a failed case, respectively.In the case of the manipulation of X (phytoplankton, Fig. 3a), all the variables changed the values from the original stable steady state to the alternative stable steady state under the simulated condition.This corresponds to a transition from a turbid to a clear state.Manipulations of D gave similar results.On the other hand, such a state shift did not occur when manipulating Z (zooplanktivorous fish); all the variables returned to the original stable state after a spike-like response at the beginning of the manipulation (Fig. 3b).Manipulations of N and Y gave similar results as those for Z.

Simulation of Food-Web Manipulation: Resilience Diagrams
Similar simulations were also carried out by manipulating two of the variables; each of the ten possible combinations ((X, Y), (Y, Z), etc.) of the variables N, X, Y, Z, D was selected, and the two variables were manipulated synchronously.The system could be switched from turbid to clear or vice versa only when the variable X or D, or both, were involved in the manipulations.Typical results are shown in Fig. 4. The diagram spanned by X and Z (Fig. 4a) shows a "quasi"-basin of attraction when X and Z were manipulated from the turbid state and all the other variables were fixed at the steady state.One can see that the quasi-basin of attraction of the turbid state is much larger than that of the clear state, and that manipulation of only variable Z cannot switch the system from turbid to clear.The location of the transition from the stable to the unstable states obtained by the linear stability analysis differs from that obtained by the present simulations as indicated by arrows in the figure.This discrepancy arises from the fact that bifurcation curves are obtained by the local stability of steady states whereas the biomanipulations are conducted away from the steady states.
The diagrams in the case of X and D manipulation are shown in Fig. 4b (manipulation from the turbid to the clear state), and Fig. 4c (manipulation from the clear to the turbid state).The quasi-basins of attraction were different in the two cases; that of the clear state was larger when the manipulation was started from the turbid state than when it was started from the clear state.This simulated result may differ from our expectation that the quasi-basin of attraction of the clear state should be larger when the manipulation is started from the clear state than when it is started from the turbid state.In terms of the lake restoration, i.e., when the manipulation was started from the turbid state, the quasi-basin of attraction of the clear state was the largest in the case of the combination (X, D).All the above results reflect that the full vector of the manipulated variables in the five-dimensional space determines the quasi-basins of attraction.The present simulations show that manipulation of only one variable (or two variables) will not always switch the states.In the case of manipulating only one variable (or less than five variables) in the present model, the success of switching the states depends on the structure of the basin of attraction in the five-variable space.In multivariable systems, resilience should ideally be represented by the basin of attraction in the full multi-dimensional space.Of course, such a representation would be beyond the limits of our imagination.As a way out of this dilemma, we propose representations of resilience by quasi-basins of attractions, such as those in Fig. 4.These represent resilience in terms of only those variables that are candidates for direct external manipulation or perturbation in the present system.This kind of representation of resilience should be considered for real ecosystem management as well (Scheffer et al. 2001, Bellwood et al. 2004, Carpenter and Brock 2004).
The resilience diagrams as shown in Fig. 4 suggest that the turbid zone of attraction is much larger than the clear zone.Thus, it is quite important for managers and landowners to adopt measures to try to maintain the clear state once it occurs, rather than attempting to switch the lake from turbid to clear.
Restoring a clear state may require lots of effort and expense.
In real lakes, if a long-term effect of a biomanipulation is to be expected, a reduction in the external loading of nutrients is recommended (Hansson et al. 1998).This recommendation implies that bistability or alternative stable states may be achieved before biomanipulation.The range 0.5-1.34mg m -2 day -1 of the nutrient input for the bistable condition in the present analysis is of the same order of magnitude as the range 0.5-2.0g m -2 year -1 , i.e., 1.37-5.48mg m -2 day -1 , recommended for the reduction of external loading of P before biomanipulation (Jeppesen et al. 1990b).
Real-lake ecosystems can be considered to be multivariable systems, and the interaction between organisms and materials is much more complex than that in the present model.The role of submerged macrophytes (Takamura et al. 2003) that affect the long-term success of biomanipulations (Jeppesen et al. 1990b, Hansel-Welch et al. 2003) is not taken into account in the present model.Submerged macrophytes are said to contribute to the stability of the clear-water state, for instance, by absorbing nutrients, stabilizing the bottom sediments, and improving the chances of zooplankton to escape predation and exert increased grazing pressure on phytoplankton (Jeppesen et al. 1990b, Hansel-Welch et al. 2003).For natural biomanipulated systems, there are reports that lake restorations succeeded by reducing the density of zooplanktivorous fish (Jeppesen et al. 1990a, 1990b, Riemann et al. http://www.ecologyandsociety.org/vol10/iss2/art3/1990, Søndergaard et al. 1990).The success is thought to be closely related to the ability of the lakes to establish large distributions of submerged macrophytes (Jeppesen et al. 1990a).The lack of an explicit description of macrophytes in the present model may be a main reason why the shift did not occur when manipulating Z (zooplanktivorous fish).
In addition, benthivorous fish that suspend sediment and contribute to the turbid state (Drenner and Hambright 1999) are not included in the model.The importance of reducing benthivorous fish, e.g., cyprinid fish, was documented by biomanipulation experiments (Hansson et al. 1998).There are also reports that the success in biomanipulations where only cyprinid fish were removed, whereas piscivorous fish caught by trawling were returned to the lakes, is nearly proportional to the amount of the reduced cyprinid fish (Hansson et al. 1998).Not all the cyprinid fish are benthivorous, however, this report suggests that there is a high probability that biomanipulations will fail when benthivorous fish are present.
An implication of the present results is that multivariable manipulations may be effective for practical restoration of lakes that lack benthic fish and display alternative states.There may be cases where significant removal of zooplanktivorous fish is not feasible, however, even in such cases, the present results suggest that manipulations of phytoplankton, detritus, or other species can help restore lakes.In particular, a direct reduction in the abundance of phytoplankton (X) and detritus (D) was shown to be most effective for switching the present system from a turbid to a clear state.This result supports, for example, the regional plan for conservation of lake water quality in Japan (cf.Chiba Prefecture 2001): direct removal of algae and dredging of benthic mud have been and will continue to be carried out to restore eutrophic lakes.

CONCLUSIONS
The alternative-stable-state theory for shallow lakes and the concept of ecological resilience are important for the restoration and management of ecosystems.A graphic representation of resilience using a folded bifurcation curve is useful for understanding the concept.However, a representation of resilience by only one variable is insufficient for multivariable systems.Resilience is embedded in the multivariable phase space.The present study shows that manipulation of only one (or two) of the variables may not switch the system from a stable state to an alternative stable state even though the value of that variable is manipulated beyond the unstable state.This is a characteristic feature of multivariable systems, or natural ecosystems.Therefore, management of ecosystems using biomanipulation requires an understanding of what is the driving variable for switching the system.In the case of a real lake restoration, reducing the abundance of both zooplanktivorous and benthivorous fish is considered to be effective.Our simulations show that reducing the abundance of both phytoplankton (X) and detritus (D) is most effective for lake restoration.

Fig. 2 .
Fig. 2. Local stability of the steady-state solutions of Eq. (1) as a function of I N , the rate of input of nutrient from the environment.N is the amount of dissolved nutrient, X the biomass of phytoplankton, Y the biomass of zooplankton, Z the biomass of zooplanktivorous fish, and D the detrital biomass.Blue and red filled circles are stable steady states, and represent clear (A) and turbid (B) states, respectively.Open circles are unstable steady states.The parameter values used are listed in Table1.

Fig. 3 .
Fig. 3. Temporal changes in the variables N, X, Y, Z, and D when only X or Z were manipulated at time zero for I N = 0.75 mg m -2 day -1 .(a) X was manipulated from a turbid stable state (B) to a clear stable state (A) as indicated by a star and an arrow in the inset of X.The other variables (except for Y) followed this manipulation.The dotted arrows in the insets of N, Z, and D represent these changes.(b) Z was manipulated from a turbid stable state (B) to a clear stable state (A) as indicated by a star and an arrow in the inset of Z.After a short time, all the state variables returned to the initial condition.

Fig. 4 .
Fig. 4. Phase diagrams (quasi-basins of attractions) spanned by two of the five variables in Eq. (1).(a) X and Z manipulation from a turbid state (B), (b) X and D manipulation from a turbid state (B), and (c) X and D manipulation from a clear state (A).All the variables were initially set to either the turbid or the clear stable state (red or blue filled star).Two of the variables were manipulated from the initial state to a certain state within the phase diagram at time zero.The regions represented by red and blue crosses represent basins of attractions of the turbid and the clear states, respectively, i.e., the final state after about 5 years (2000 days) by the two-variable manipulation was the turbid and the clear state, respectively.An open circle represents an unstable state.Local stabilities of the two variables are also shown.Note that the magnitude of the arrows one would expect from the one-variable bifurcation diagram does not coincide with the magnitude in the two-variable phase diagram.