Structural Instability, Multiple Stable States, and Hysteresis in Periphyton Driven by Phosphorus Enrichment in the Everglades

Theoretical Population Biology 61, 1–13 (2002) doi:10.1006/tpbi.2001.1554, available online at http://www.idealibrary.com on

Structural Instability, Multiple Stable States, and Hysteresis in Periphyton Driven by Phosphorus Enrichment in the Everglades

Quan Dong Southeast Environmental Research Center, Florida International University, OE 148 University Park, Miami, Florida 33199

Paul V. McCormick and Fred H. Sklar Everglades Systems Research Division, South Florida Water Management District, 3301 Gun Club Road, West Palm Beach, Florida 33416-4680

and Donald L. DeAngelis Biological Resources Division, U.S.G.S., Department of Biology, University of Miami, Miami, Florida 33124-0421 Received March 10, 1999

Periphyton is a key component of the Everglades ecosystems. It is a major primary producer, providing food and habitat for a variety of organisms, contributing material to the surface soil, and regulating water chemistry. Periphyton is sensitive to the phosphorus (P) supply and P enrichment has caused dramatic changes in the native Everglades periphyton assemblages. Periphyton also affects P availability by removing P from the water column and depositing a refractory portion into sediment. A quantitative understanding of the response of periphyton assemblages to P supply and its effects on P cycling could provide critical supports to decision making in the conservation and restoration of the Everglades. We constructed a model to examine the interaction between periphyton and P dynamics. The model contains two differential equations: P uptake and periphyton growth are assumed to follow the Monod equation and are limited by a modified logistic equation. Equilibrium and stability analyses suggest that P loading is the driving force and determines the system behavior. The position and number of steady states and the stability also depend upon the rate of sloughing, through which periphyton deposits refractory P into sediment. Multiple equilibria may exist, with two stable equilibria separated by an unstable equilibrium. Due to nonlinear interplay of periphyton and P in this model, catastrophe and hysteresis are likely to occur. © 2002 Elsevier Science (USA)

Key Words: catastrophe; hysteresis; threshold; periphyton; phosphorus; Everglades.

⁄

1

0040-5809/02 $35.00 © 2002 Elsevier Science (USA) All rights reserved.

2

INTRODUCTION Phosphorus (P) is the principal nutrient limiting the productivity of many freshwater ecosystems, and P enrichment has become a major problem in these ecosystems worldwide (DeAngelis et al., 1989; Tiessen, 1995; Vymazal, 1995). Algae respond rapidly to increased nutrient loading, and excessive algal growth is one of the common and undesirable consequences of P enrichment in lakes and rivers (Carpenter et al., 1998). Attached algal assemblages, or periphyton, are a ubiquitous and ecologically important feature of wetlands (Goldsborough and Robinson, 1996) and might be expected to respond in a similar manner to increased P loading. However, relatively few studies have investigated the effects of increased P loading on wetland processes in general and, specifically, on periphyton abundance and growth. The Florida Everglades is one of the few wetlands where the ecology of the indigenous periphyton assemblage and its response to P enrichment have been studied in detail. The Everglades is an extensive subtropical wetland that encompasses a large portion of southern Florida. Pristine areas of this wetland consist of extensive stands of the emergent macrophyte sawgrass (Cladium jamaicense) interspersed with more sparsely vegetated areas dominated by low stature emergent (wet prairies) and floating macrophytes (sloughs) (Gunderson, 1994). Attached and floating periphyton mats account for a large fraction of the vegetative biomass and primary productivity in wet prairies and sloughs. These mats are formed by the calcium-precipitating (i.e., calcareous) cyanobacteria Scytonema and Schizothrix and have a three-dimensional, calcified exoskeletonal structure (Gleason and Spackman, 1974). During the wet season, floating mats several centimeters thick can develop, covering most submerged surfaces and much of the water surface, and they provide habitats and food for a variety of animals (Browder et al., 1994; McCormick et al., 1998). The Everglades periphyton assemblage is extremely sensitive to changes in P supply and provides an early signal of eutrophication in this wetland (McCormick and Stevenson, 1998). Increased P loading in the form of agricultural runoff produces dramatic changes in periphyton standing crop, productivity, and species composition (reviewed in Browder et al., 1994 and McCormick and Scinto, 1999). One of the more sensitive and obvious changes caused by enrichment is the loss of the abundant calcareous mats that characterize pristine (low P) areas. The contribution of periphyton photosynthesis to ecosystem primary productivity declines with increased P loading (McCormick et al., 1998).

Dong et al.

Periphyton changes caused by increased P supply may in turn affect P dynamics within the wetland. Periphyton functions as a major transformer and sink for P in the Everglades (Craft and Richardson, 1993; McCormick et al., 1998). Periphyton mats rapidly remove excess P from the water column (Scinto and Reddy, 1997) and may be responsible for maintaining low P availability in pristine areas where periphyton is abundant (McCormick et al., 1998). The importance of periphyton to P removal in enriched areas may be considerably diminished due to low standing crop. Calcareous mats precipitate P in the form of calcium–phosphate complexes, which are relatively insoluble and are deposited as calcitic muds, a common soil type in the pristine Everglades (Gleason and Spackman, 1974). Thus, increased P loading in association with the loss of calcareous periphyton mats may reduce the ability of this wetland to regulate water-column P availability. Many models have been built to simulate P dynamics in wetlands. These models often use a constant (settling) rate to describe P deposition (e.g., Walker, 1995). Estimating the settling rate for such models is often difficult. Uncertainties are probably great but very difficult to quantify and measure due to their numerous uncontrollable sources (see Qian, 1997a, 1997b, and Walker, 1997). Recently, many data have been collected about the relationship between periphyton and P dynamics. The current information about periphyton responses to P supply suggests that the P settling may vary greatly as a consequence of biological processes. Nevertheless, it is not clear if and how the responses alter the settling rate. It is desirable and necessary to quantitatively understand how biotic processes in periphyton affect P fluxes in wetlands in general and in the Everglades in particular. In addition, scientific hypotheses concerning the quantitative relationship between periphyton and P supply have significant implications in the formulation of policies and administrative rules. For example, Florida law requires the establishment of a P criterion that can be used to make rules ‘‘to prevent an ecological imbalance’’ and to restore areas impacted by high P loading. Periphyton has been identified as a biological indicator of the P threshold (Lean et al., 1992). In this study, we used a modeling approach to examine the interaction between periphyton and P supply. We constructed a suite of models to serve as a quantitative framework that can be used to (a) generate hypotheses about the threshold and the retention capacity of P, (b) evaluate importance of various ecological processes and parameters, such as sloughing and mortality of periphyton, and (c) identify critical information missing in current ecological understanding. The models are

3

Catastrophe and Hysteresis in Periphyton

designed to address two questions: (a) how do periphyton communities respond to P supply, and (b) how does the response of periphyton communities in turn affect P dynamics. The present paper is largely restricted to description of one particular model with limited structural complexity. With this model, we used equilibrium analysis, stability analysis, and bifurcation analysis to reveal general features of interactions between periphyton and P and the system dynamics.

MODEL STRUCTURE Our model contains two dependent variables: (1) the amount of dissolved, available P per unit volume in the water, P, and (2) the amount of P bound in living periphyton per unit volume, A. A system of two equations describes the dynamics of periphyton and P: dA/dt=F(A, P),

(1)

dP/dt=G(A, P).

(2)

and

The forms of F(A, P) and G(A, P) specify the way periphyton and P behave and interact. These forms are based on empirical information. Observations from field, laboratory, and mesocosm studies suggest that P enrichment in a previously oligotrophic marsh could have two major effects on periphyton (Steward and Ornes, 1975; Swift and Nicholas, 1987; Gleason and Spackman, 1974; Flora et al., 1988; McCormick and O’Dell, 1996; McCormick et al., 1996, 1998; Richardson et al., 1996; McCormick and Scinto, 1999; McCormick and Stevenson, 1998). First, on periphytometers (artificial substrata), the biomass accumulation increases with P additions on the short term (a couple of weeks), in field dosing experiments. This indicates that the growth rate per unit biomass of periphyton increases with P, particularly in the low concentration range of water column P. However, the biomass accumulation decelerates and soon stops in several weeks. Second, increases in P supply exert a detrimental effect on calcareous periphyton mats. At high P concentrations, the mat matrices lose their architectural integrity. Mats collapse and dissolve. Without structural support of calcareous mats, the standing crop of periphyton remains low. We incorporated these two effects into a simple model. A plausible representation of the net growth rate per unit biomass of periphyton is F(A, P)/A=umax (P/(ku +P))(1 − A/Amax (P)) − m, (3)

where, umax is the maximum growth rate of periphyton (g/g · day); Amax (P) is the maximum standing crop of periphyton, incorporating self-limiting effects (g/m 3); ku is the half-saturation coefficient for growth rate of periphyton as a function of P(g/m 3); and m is the mortality (g/g · day). In Eq. (3), the term umax P/(ku +P) is the classic Monod function (Monod, 1950) for periphyton growth and P uptake, and it has been used widely to represent the P limitation on algal growth in freshwater systems in the absence of toxic and prohibitive effects (Grover, 1990, 1991; Borchardt et al., 1994). This type of model has been used previously in the empirical studies of interactions between periphyton and P in South Florida wetlands (Scinto and Reddy, 1997; Hwang et al., 1998). The Monod function assumes that the growth rate of periphyton increases rapidly with P at low concentrations of P, while the increase per unit increment of P diminishes and becomes insignificant at high P concentrations. The next factor in Eq. (3), (1 − A/Amax (P)) is in the form of a modified logistic function. The logistic equation also has been used widely in algal studies (Rodriguez, 1987; Momo, 1995). In this model, the factor represents a biomass-dependent, self-limiting effect on growth, resulting from a combination of competition for substrates, light, or nutrients other than P. The limiting effect of P and depletion of P are explicitly included in the equations. It is assumed here that the architecturally supportable maximum of standing crop, Amax , is a decreasing function of P, Amax (P). This can be represented by the nonlinear curve shown in Fig. 1. Field and experimental observations (Flora et al., 1988; Richardson et al., 1996; McCormick et al., 1998; McCormick et al., 2001; McCormick and Stevenson, 1998; McCormick and Scinto, 1999) indicate that phosphorus effects on the periphyton mat matrices are slight at low P concentrations, but become substantial as concentrations increase above background levels. For simplicity, the Amax (P) curve can be approximately represented in a piecewise linear form

˛

a1 , (P < P1 ), a1 +(P − P1 )(a1 − a2 )/(P1 − P2 ), Amax (P)= (P1 < P < P2 ), a2 , (P, P2 ).

(4)

(see Fig. 1) At any given point on the Amax (P) curve, the derivative, d(Amax (P))/dP, measures the effect of changes in P on the maximum supportable standing crop of periphyton assemblage at a particular P concentration. Note that

4

Dong et al.

A, Amax , P, m, and v \ 0,

and

Pin , umax , ku , a1 , a2 , P1 , and P2 > 0. When the uptake of P by periphyton is taken into account, the equations describing changes in periphyton and in water P through time are, dA/dt=A[umax (1 − A/Amax (P))(P/(ku +P)) − m], (6) and dP/dt=v(Pin − P) − A[umax (1 − A/Amax (P))(P/(ku +P)) − m] − fmA. FIG. 1. A hypothetical curve depicting the effect of P on the maximum supportable standing crop of periphyton. The horizontal axis is water-column P concentration. The vertical axis indicates the maximum supportable standing crop of periphyton in the unit of P. A piecewise linear representation is used to approximate the smooth curve. Roughly, there are three regions along the P axis, divided by points P1 and P2 . The region between P1 and P2 is the ‘‘transitional region’’ for periphyton, in which the periphyton response to P is most sensitive and conspicuous architectural changes occur. In the other two regions, periphyton standing crop is relatively stable to changes in P.

Amax (P) is a monotonically decreasing curve (Fig. 1), i.e., d(Amax (P))/dP is either negative or zero. Thus, −d(Amax (P))/dP measures the strength of the P effect. As we show below, this measure is an important parameter. Finally, in Eq. (3) the mortality, m, represents the loss rate per unit biomass of periphyton through sloughing. During sloughing, dead periphyton tissues deposit into sediment. We assume that m is a constant. We assume a constant flow of P into the system and a loss rate due to washout that is dependent on watercolumn P concentration. In the absence of periphyton, the equation for P can be written, dP/dt=(Pin − P) v,

(5)

where Pin is the P concentration in water that flows into the system from external sources (g/m 3); and v is the daily water flux through the system, measured as a proportion of the total water volume of the system. v=vŒ/w. vŒ is the velocity of water through the system (m 3/day). w is the total volume of the system (m 3). In Eq. (5), vPin is the per space unit P loading flux from external sources (g/day/m 3). All of the state variables and parameters in the above equations are nonnegative and most parameters are positive:

(7)

In Eq. (7), we have made the assumption that a constant portion of the dead tissue material, f, is recalcitrant and does not decompose before going into the sediment. The remainder of the dead periphyton biomass, which dies at the rate mA, decomposes and releases its P back to the water. The deposition of refractory fraction of tissue P into sediment after death removes P from the pools of periphyton and water. The importance of periphyton on P cycling lies in the deposition. The term, fmA, is the rate of P transport from water to sediment by the periphyton. It can be regarded as the system’s biological assimilative capacity of P.

EQUILIBRIA AND STABILITY Equilibria and their stability are two useful attributes, delineating the long-term behavior of a system. We evaluated equilibria and their stability analytically and graphically, and examined how different processes can affect them. We were interested only in those biologically reasonable equilibria for which A and P are positive (A > 0, and P > 0). 1. Steady State Equilibria We evaluate the equilibria first. The equilibria, (A g, P g), were determined by setting the right-hand sides of Eqs. (6) and (7) to zero: dA/dt=0 and dP/dt=0. Solving for the A in terms of P in both equations, we first noted that the equilibria are easily found as the intersections of the dA/dt zero isocline A=Amax (P)(1 − m(1+ku /P)/umax )

(8)

5

Catastrophe and Hysteresis in Periphyton

with the line A=v(Pin − P)/fm.

(9)

Figures 2, 3, and 4 show possible configurations of the zero isoclines, dA/dt=0 and dP/dt=0. We will discuss different configurations below. In general, the intersection points of these two isoclines represent equilibria. The straight line defined by the Eq. (9) will always intersect the equilibria. The line defined by Eq. (9) intersects the P axis at the point Pin , therefore illustrating the relation between Pin and the equilibrium point. Equation (9) suggests and the straight line in Figs. 2, 3, and 4 shows that P concentration of the in-flow, Pin , is the dominant factor determining both the number and position of equilibria. Nevertheless, the number and position of equilibria also depend on several other parameters, particularly, the slope of Eq. (9), −v/fm, and the steepest downward slope of the dA/dt=0 curve, min(dA/dP). The pattern of A g and P g and their responses to Pin can be examined under the two opposite conditions: min(dA g/dP) \ −v/fm,

and

min(dA g/dP) < −v/fm. First, consider the case that min(dA g/dP) \ −v/fm. Only one positive equilibrium point exists, at the intersection of the two isoclines (Figs. 2a and 2b). Roughly four regions can be seen, in each of which A g and P g respond differently, with respect to the Pin . These four regions are divided by three levels of P; P0 (P0 =mku /(umax − m)), Pa , and Pb (Fig. 2a). When Pin is extremely low, Pin < P0 , the periphyton, A, cannot maintain positive growth, and thus disappears. If Pin decreases from a higher level, negative growth may occur. If the Pin level is between P0 and Pa (P0 < Pin < Pa ), A g will increase rapidly with Pin , while P g does not show much change. The Monod growth function plays a large role in shaping the equilibrium pattern in the low Pin environment. In the third region, Pa < Pin < Pb , the algal biomass and water column P respond sensitively to changes in Pin . A dramatic but continuous transition of the steady state equilibrium occurs. When Pin is high, Pin > Pb , P g is high and A g is low. A g is insensitive to Pin (Fig. 2). P g is almost linearly increasing with Pin in this range. It is interesting to see the influence of P effect on periphyton mats, d(Amax (P))/dP, on the relationship of A g and P g. By taking the derivative of A in Eq. (8) with respect to P, we obtain dA/dP=(1 − m(1+ku /P)/umax ) d(Amax (P))/dP +Amax (P) mku /(umax P 2).

FIG. 2. Equilibrium and its stability, and the effect of changes in Pin on the equilibria and stability. The solid curves are isoclines, where dP/dt=0 and dA/dt=0. Their intersection point is the nonnegative equilibrium. The straight broken line representing the equation, A=v(Pin − P)/fm, intersects at equilibrium. Sloughing, measured by fm, is one of the major determinant of the positions and stability of equilibria. In this case, the value of fm is relatively low, thus one equilibrium exists. Arrows show the moving direction of system around that point. Inward arrows indicate a tendency to return to the equilibrium, and thus a stable equilibrium. Outward arrows indicate unstable equilibrium. Instability may occur only when Pin is between Pa and Pb . In (a) the slope of Amax (P) curve is moderate, and the equilibrium is stable. In (b) the slope of Amax (P) curve is steep and the equilibrium is unstable. See text for more details.

This implies that the downward slope of dA/dP becomes steeper as d(Amax (P))/dP becomes more negative. The sloughing rate, fm, is another important parameter. Equation (9) (and Fig. 2) shows how fm affect the system attributes. Note that the difference between Pin

6

Dong et al.

FIG. 4. Equilibria and their stability. The solid curves are isoclines, where dP/dt=0 and dA/dt=0. Their intersection points are nonnegative equilibria. The straight broken line representing the equation, A=v(Pin − P)/fm, lies across the equilibria. In this case, the value of fm is relatively high. The steepest slope of Amax (P) curve is steeper than the slope of the broken line, −v/fm. Thus, multiple equilibria may exist. Stable equilibria are surrounded by inward arrows from all directions. Outward arrows from an equilibrium to any direction indicate instability. The equilibria E1 and E3 are both stable nodes, while the equilibrium E2 is a saddle point.

dA/dt=A[umax (P/(ku +P))(1 − A/Amax (P))], dP/dt=v(Pin − P) − A[umax (P/(ku +P))(1 − A/Amax (P))].

FIG. 3. Equilibrium and its stability; outcome from the model that does not include the effect of sloughing. The solid lines are isoclines, where dP/dt=0 or dA/dt=0. The intersection point of two nonlinear isoclines is the nonnegative equilibrium. At equilibrium, P g=Pin . One equilibrium exists for a given Pin . Inward arrows indicate a tendency to return to the equilibrium, and thus a stable equilibrium. Outward arrows indicate unstable equilibrium. Instability may occur only when Pin is between P1 and P2 . (a) The equilibrium in the transitional region is stable, if the slope of Amax (P) curve is moderate, and (b) the equilibrium is unstable when the slope is steep.

and P g measures the removal of P by the system at steady state. The difference between P g and Pin decreases, as −v/fm decreases with fm. Periphyton removes a smaller portion of P from the water column. The special case, m=0, means that sloughing does not occur in the model. Many earlier models that simulated interactions between algae and nutrients ignored P sedimentation through sloughing (e.g., Edelstein-Keshet, 1988), and thus were constructed in this way. In this case, the model is equivalent to

The equilibrium solution of this model is P g=Pin , and A g=Amax (Pin ). Thus, first the P concentration will reach Pin at steady state, without sloughing. The periphyton biomass does not affect the equilibrium value of P. Second, P limitation on growth is unimportant in determining the steady states. Periphyton biomass is high in the low Pin region (Pin < P1 , in Figs. 3a and 3b). Thus, the position of equilibrium (A g, P g) and its response to Pin differ significantly from the case in which fm > 0, in the low Pin regions (Figs. 3a and 3b). In the high Pin region, periphyton biomass is low (P2 < Pin , in Figs. 3a and 3b). The effect of changes in Pin , on periphyton is small within both low and high P regions. With intermediate levels of Pin , P1 < Pin < P2 , the significant transition occurs. The periphyton biomass responds sensitively to a change in Pin and declines approximately with a magnitude of (a1 − a2 )(P1 − P2 ) for a unit increase in Pin . There are only three regions with respect to the driving force Pin (divided by the points P1 and P2 ). Third, the water flow is also insignificant in determining the equilibrium. The comparison between results of models with fm > 0, and with fm=0 reveals the contribution of sloughing to system behavior and demonstrates the significance of fm

7

Catastrophe and Hysteresis in Periphyton

and its explicit representation in model. It is interesting to note that the width between Pa and Pb could be narrower than that between P1 and P2 (Figs. 2a, and 1, or, 3a). Thus the transition of A g may appear more dramatic than Amax (P), due to the interactive feedback of biomass change. Second, consider the opposite case, where min(dA g/dP) < −v/fm. Multiple equilibria may now exist. There are still roughly four regions with respect to the Pin , in which the positions and the number of equilibria differ (Fig. 4). The pattern of A g and P g is similar in the first, second, and fourth regions, as before in the case in which min(dA g/dP) \ −v/fm. One equilibrium exists. A g either stays at zero, rises, or is insensitive as Pin increases, in the first, second, or fourth regions, respectively. In the third region, where Pin is between Pa and Pb , the isoclines may cross more than once. Up to three equilibria with positive A g and P g could exist, one at high A g and low P g (E1 ), one at low A g and high P g (E3 ), and one at intermediate A g and P g (E2 ) (Fig. 4). As discussed in more detail below, the equilibrium at intermediate A g and P g is unstable. Thus, transition from one equilibrium to another is an abrupt change, in contrast to the continuous change in the case, min(dA g/dP) \ −v/fm. The abrupt change can occur even when the response of mat integrity to P is gradual and continuous, due to the interplay of the periphyton response to P and the consequent changes in the P removal by periphyton. With high A g and relatively low Pin , periphyton removes a relatively large portion of P, and thus maintains P g at low levels. With low A g and relatively high Pin , periphyton only removes a small portion of P from the water column. P g remains high. Thus, a positive feedback operates in certain circumstances in the interaction between periphyton and P. The architectural response of the periphyton mat seems to be an important mechanism that shapes the equilibrium pattern. Overall, multiple steady states may occur. P loading (Pin v) is driving the system. The sloughing process that removes P from the system at the rate fmA and the P effect on periphyton mat architecture (Amax (P)) are major processes determining the system states. Periphyton affects wetland P dynamics mainly through P uptake from the water and P deposition in sediment after death. Thus, growth, death, decomposition, and deposition are all critical processes affecting the capacity of P retention and assimilation in periphyton-dominated marshes. 2. Stability We examined the local stability of equilibria and the determinants of the stability with the Jacobian of the system of equations; i.e., the first-order equations linearized

about the equilibrium points (e.g., Edelstein-Keshet, 1988). Such an analysis helped to identify potential stabilizing and destabilizing processes. The real parts of eigenvalues (l) of the Jacobian are measures of stability. l1, 2 =(b ± `b 2 − 4c)/2; where we define, b=(dF/dA+dG/dP)A*, P* , and, c=[(dF/dA)(dG/dP)− (dF/dP)(dG/dA)]A*, P* . If one or both of the eigenvalues are positive, the system is unstable (Lyapunov instability). The sufficient condition for Lyapunov instability of the equilibrium point is either, b > 0, or, c < 0. The rate at which the system returns to equilibrium is termed ‘‘relative stability.’’ A system with negative real eigenvalues has a relative stability in proportion to the absolute values of eigenvalues. In a loose sense, the largest l increases as b increases and as c decreases. Calculation of b, c, and l suggests that the general condition for an equilibrium to be unstable is either b=m − umax P g/(ku +P g) − v − (dF/dP)A*, P* > 0, (10) or c=fm(dF/dP)A*, P* − v(m − umax P/(ku +P)) < 0. (11) For A g to be positive, it must be true that m − umax P/ (ku +P) < 0. Thus, Eqs. (10) and (11) imply that unstable equilibria could occur only when (dF/dP)A*, P* is negative. (dF/dP)A*, P* is the derivative of Eq. (6) with respect to P at the equilibrium point (dF/dP)A*, P* =Aumax [(ku /(ku +P) 2)(1 − A/Amax (P)) +(A(P/(ku +P))/Amax (P) 2) × d(Amax (P))/dP]. This suggests, the instability occurs mainly in the transitional region of Pin with a steep slope of Amax (P) (Figs. 2b, 3b, and 4). A more negatives d(Amax (P))/dP with a larger fm tends to destabilize the system.

8

Dong et al.

Inequality (10) also shows that b (and thus the real parts of the eigenvalues and the instability of the system) increases with an increase in −d(Amax (P))/dP. d(Amax (P))/dP is the slope of the tangent lines of Amax (P) curve in Fig. 1. It stands for the strength of the deteriorating effect of P on periphyton mats, and is a critical parameter to stability. If Amax (P) were an increasing function of P, so that d(Amax (P))/dP > 0, then the equilibrium would always be stable. However, field and experimental observations suggest, and we are assuming, the opposite. The inequalities (10) and (11) suggest that an increase of v, water flux, is stabilizing. Among other parameters, umax and ku tend to play opposite roles in stabilization. Inequality (11) shows that the equilibrium is stable in both the high and low Pin regions. Since (dF/dA)A*, P* =m − umax P/(ku +P) < 0, we can rewrite Eq. (11) as (dA/dP)A*, P* < −v/fm.

(12)

Equation (12) means that the equilibrium is unstable if the slope of the straight line is less steep than the slope of periphyton isocline at the intersection in Fig. 4. This suggests that the equilibrium in the middle is unstable when three steady state equilibria exist. Figure 4 reveals the mathematical nature of the three equilibria. Points E1 and E3 are both stable nodes, while E2 is a saddle point. Thus, the plane is divided into two basins of attraction, one containing all trajectories that move toward the point E1 and one containing all trajectories that move towards E3 (also see Fig. 5). The high value A g has a high removal capacity of P and maintains low P, and vice versa. This operates as a self-reinforcing mechanism, keeping A g high and P g low, or A g low and P g high. This mechanism creates the local stability, but global instability. With the special case, m=0, A g=Amax (P g), and g P =Pin , the instability condition is simplified. First, c= vumax P/(ku +P) > 0. Thus, only functions and parameters composing the expression for b determine the stability of the equilibrium, and b=(−d(Amax (P))/dP|P=Pin − 1) umax Pin /(ku +Pin ) − v. We have b < 0 if −d(Amax (P))/dP|P=Pin < v(ku +Pin )/umax Pin +1.

(13)

FIG. 5. Bifurcation. Panels (a) and (b) show how the steady state periphyton and P respond to Pin changes, respectively. Broken arrows indicate the direction of change as a consequence of the increase of Pin and solid arrows indicate the direction of change with the reduction of Pin . A gradual increase in Pin tends to result in a stable equilibrium with increasing values of A g until Pin goes beyond Pp . Then a dramatic drop of A g occurs, because the mat structure cannot support the high biomass. The state plane suddenly becomes a basin of attraction for the stable equilibrium at low level of A g and high level of P g. When Pin decreases from a high level, the system remains at low A g, even after the level of Pin decreases below Pp . One must decrease Pin below Pr to restore a high A g equilibrium. The state plane then becomes a basin of attraction of the equilibrium of low P g and relatively high A g.

Inequality (13) is the stability criterion for the single equilibrium, when fm=0. In summary, in this model system the in-flow P concentration, Pin , is the major determinant of the

Catastrophe and Hysteresis in Periphyton

equilibrium. The strength of the detrimental effect of P increase on periphyton mat is a major determinant of the stability. The sloughing rate, fm, is another important determinant of the equilibrium and its stability.

PERTURBATION, HYSTERESIS, AND NONEQUILIBRIUM BEHAVIOR The P input, Pin , is a driving force in this type of system. A question of much management interest for ecosystems such as the Everglades, parts of which have been enriched by increased loads of P in recent decades, is what will happen if (a) loading continues and expands into the areas that are currently not impacted, and (b) loading decreases in areas that have already been enriched. These two scenarios were explored by increasing or decreasing Pin in the model. Figure 5 explicitly shows the outcome of these two scenarios. First, when Pin increases from a low level gradually, a single stable equilibrium exists. At this phase of equilibrium state, P serves as a limiting factor and A g increases with Pin (Fig. 5a). A large portion of P is taken up by periphyton from water column and then deposited into the sediment. Thus, P g shows a little change and remains at a very low level (Fig. 5b). As A g reaches the peak value, a further small increase of Pin beyond a threshold level, Pp , would lead the equilibrium to an abrupt switch, from a stable state with high A g and low P g to another stable state at low A g and high P g. This drop in A g represents the dissolution of periphyton mats. P is no longer limiting. Such a large change in the equilibrium, as a response to a small difference in Pin around the threshold, is a result of a structural instability or a ‘‘fold catastrophe’’ (see Zeeman, 1977). The abrupt change occurs even when the response of mat integrity to P is gradual and continuous, due to the interplay of the periphyton response to P and the consequent changes in the P removal by periphyton. Pp is a threshold for the equilibrium at high A g and low P g to exist. We may call this threshold the ‘‘protection threshold.’’ The equilibrium with high A g and low P g represents a thick periphyton mat in an oligotrophic environment, which is a feature of the slough communities dominated by native periphyton in the Everglades. Now, suppose that efforts are made to reduce P loading from a high level gradually. As Pin is reduced, the system remains at low A g and high P g, even after the level of Pin decreases below the protection threshold. In fact, Pin may have to be decreased to a much lower level, Pr , to restore the high A g, low P g equilibrium (Fig. 5). Pa indicates a threshold for the system to recover from a crash in A g due to high Pin . Below this Pin level, the entire state

9 plane becomes a basin of attraction for high A g, and the system returns to that attractor (Fig. 5, also see Fig. 6). We may call Pr the ‘‘restoration threshold.’’ If the protection threshold differs from the restoration threshold, hysteresis exists. Hysteresis means that the system going through a change in a control variable (e.g., the change of Pin from a to b in Fig. 5) does not return to its initial state when the control variable is returned back to its initial point (e.g., the change in Pin from b back to a in Fig. 5). Instead, the control variable must be returned beyond its starting point before recovery can occur. Hysteresis occurs in this system because of a strong, nonlinear interplay of periphyton and water column P. Periphyton removes P from the water. The removal capacity depends on the standing crop. When the P loading rate is low, periphyton standing crop and the capacity for periphyton to remove P are high. Under these conditions, it may be possible for the marsh to assimilate increased P loads and maintain low P availability. In contrast, if P inputs are reduced from previously high loading rates, conditions under which periphyton standing crop and removal capacity are low, then the capacity for P removal may remain low. Therefore, even under a scenario of reduced P loading, excess P may remain in the water column and inhibit the formation of calcareous mats that support high standing crops. The occurrence of this type of hysteresis is likely to be an important feature of this system. Environmental perturbations and stochastic events often push the system into nonequilibrium states. When multiple stable states exist, it is interesting to examine how new nonequilibrium states caused by perturbation determine the asymptotic behaviors. We ran simulations with the model to investigate the nonequilibrium behaviors. The parameter values used in the simulation are based on empirical studies, and thus they are realistic. As earlier analyses indicated, with low and high Pin , the system approaches to one attractor of high A and low A, respectively. With an intermediate Pin level, there are two attractors. A small difference in the initial condition may lead to very different stable states (Fig. 6). This implies that below a certain level, perturbations may not affect the periphyton and P much, as they tend to approach one asymptotic state. A small increase in the perturbation at a certain magnitude may push a large change and the system approaches to another asymptotic state. In this model, we assumed that the self-limiting effect operates through growth. The self-limiting effects may also operate via mortality instead of growth. We modified the model with this assumption and conducted equilibrium and stability analyses. The analysis of the new model also showed possible catastrophe and hysteresis. Therefore, general patterns presented above seem to be robust, regardless of some details in the model.

10

FIG. 6. Multiple attractors and nonequilibrium simulation. Simulation trajectories are shown starting from various initial conditions. For this value of Pin , the plane is largely divided into two basins of attraction. One contains all trajectories that move toward the equilibrium point E1 with high A g and low P g, and one contains all trajectories that move towards the equilibrium E3 with low A g and high P g. Two stable equilibria are likely to occur and both equilibria are locally stable. The parameter values are: m=0.05, f=0.5, Umax =2, ku =25, v=0.4, a1 =1400, a2 =400, P1 =10, P2 =30.

DISCUSSION We constructed a two-variable model to simulate the interaction between periphyton and P dynamics. Despite the simplicity of the model structure, our model generated some new insights and identified some important ecological processes and missing links. For example, a major prediction of our model is the likely occurrence of multiple stable states, and structural instability. One of the characteristics of this kind of system is that there can be a large change in the state of the system resulting from a small change in the control variable, which in this case is the P loading. This result was supported by two lines of evidence. First, field observations have demonstrated that changes in periphyton biomass and taxonomic composition occur rather abruptly along P gradients in the marsh. Second, mesocosm dosing experiments showed that the changes in periphyton are associated with rather modest increases in water-column P concentrations (McCormick and O’Dell, 1996; Richardson et al., 1996; McCormick et al., 1998; McCormick and Scinto, 1999). Our analyses also show the potential for hysteresis in this type of system. The likely existence of hysteresis implies that two thresholds may exist, one for protection and one for restoration. Studies of lake restoration have shown this could occur. In some lakes, algal biomass and

Dong et al.

trophic status do not always recover in a desired direction, and the symptoms of eutrophication remain after a significant reduction of P loading (Sas, 1989; Tilzer et al., 1991). We predicted that increased loading rates do not result in proportional increases in water column P, at least at low P loading levels. This has been supported by the studies in the Everglades (McCormick et al., 1998) and other ecosystems (e.g., Pomeroy, 1960). McCormick et al. (1998) found that added P is scavenged and retained efficiently in the periphyton due to the P-limited nature of the Everglades marsh, and they concluded that periphyton may represent the major sink for this excess P in open-water habitats, which occupy large areas of the marsh interior. We believe this model captures some realism and generality in the Everglades and similar ecosystems. Structural instability, multiple stable states, and hysteresis in this system are mainly a consequence of nonlinear interactions between periphyton and P, and particularly two opposite effects of P on periphyton. First, periphyton growth is limited by P at low P levels. Second, the physical structure of periphyton mats breaks and the maximum standing crop declines at high P level. This type of structural instability and hysteresis has been found in several algal and plant studies (Kempf et al., 1984; Gatto and Rinaldi, 1987; Momo, 1995; Scheffer et al., 1997; see Loehle 1989 for review of catastrophe in various ecological systems). These models are similar to ours in that they are process-oriented, contain the Monod growth response to nutrients and the logistic density limitation, and often include two or more opposite effects occurring in a relationship of two variables. This kind of model seems more realistic and offers more mechanistic insights than Rosenzweig’s model (1971), which leads to the hypothesis of enrichment paradox. We found that even simple modifications to the basic form of the Monod or logistic equations are usually sufficient to produce a fold catastrophe and hysteresis. Structural instability and hysteresis might be common features in ecological systems. Further, these and our study suggest that structural instability and hysteresis characterizes responses of natural systems to anthropogenic influence, due to the nonlinear nature of the responses. Structural instability and hysteresis also can be sources of variability in system states and presents a challenge to empirical studies. Data analyses that are based on simple, static, linear statistical models or the assumption of a single stable equilibrium state usually fail to deal with this kind of variability and instability. As important, potential features of this system, structural instability and hysteresis may have significant implications for the Everglades restoration. Structural instability may entail difficulties in preserving the

11

Catastrophe and Hysteresis in Periphyton

Everglades marshes, and necessitate cautious conservation strategies and ‘‘threshold’’ studies. The potential for hysteresis implies that two thresholds could exist, one for protection and one for restoration. This study predicts that the threshold for recovery is likely to be lower than the threshold for protection, with respect to P loading. Some ongoing P dosing experiments with periphyton in mesocosms will soon terminate. Continuous monitoring of these mesocosms may provide data to test and quantify this prediction. Structural instability, multiple equilibria, and hysteresis also imply a likely falsification of a hypothesis suggested by Hopkinson et al. (1997), which is important for rule-making in the restoration and conservation of the Everglades. The hypothesis states that small increases in P inputs over a long time will have the same end result as larger increases in P inputs over a short time. This is not necessarily true. Within certain ranges, a steady P input may lead to a stable equilibrium with a high standing crop and a fluctuating P input may lead to a switch between equilibria of high and low standing crops. We found that the mortality, m, and the refractory fraction of dead tissue, f, are critical demographic parameters, and the detrimental effect of P on periphyton, Amax (P), is a critical biological function. They determine the number and positions of equilibria, the stability of equilibria, and dynamic behaviors of the system. Unfortunately, estimates, in situ or in vitro, of many these rates, such as mortality and sloughing (overall consequence of death and deposition), are rarely made (Reynolds, 1984; McCormick and Stevenson, 1991). The architectural response of periphyton to P is poorly understood. These parameters and processes deserve more attention in periphyton studies and need to be explicitly represented in model. The mortality and refractory fraction of dead tissue determine the P flux from water through periphyton to sediment. This flux is an important measure of the assimilative capacity of marshes. Our study suggests that the efficiency of P removal is high at low P concentrations and low with high P loading in periphyton-dominated slough areas. The P flux to sediment could even decline when P concentration increases beyond a certain level. This proposition is supported by observations of Craft and Richardson (1993), but deviates from several widely accepted theories and model assumptions. First, this means the P flux is not directly proportional to watercolumn P concentrations, which is assumed by many water quality models that ignore biological processes (e.g., Walker, 1995). Second, this also means that the periphyton systems do not follow the theory that P recycling is high in low P environments (see DeAngelis

et al., 1989). According to this well-accepted theory, a high level of recycling and reuse of nutrients in oligotrophic environment are evolutionarily adaptive. The marshes in the Everglades appear to offer an example of an exception to this generalization. The adaptive argument needs to be further examined. In the Everglades, P enrichment leads green algae to replace abundant filamentous cyanobacteria. In other freshwater systems, P enrichment usually leads to blooms of filamentous cyanobacteria and the replacement of green algae (Carpenter et al., 1998). These deviations suggest a need to improve theories in this area and call for more studies. Periphyton is an important component of the food web and ecosystem of Everglades marshes. Our periphyton model is a part of our food web model and ecosystem model in our meta-frame modeling approach. The simple model described here ignored many factors and processes that may influence the dynamics of periphyton. The mathematical analysis of this model and the hypotheses of structural instability and multiple stable states help us to identify the important components in model development. For example, we need to look at the processes that cause the collapse of mat architecture. Periphyton community shifts from calcareous to noncalcareous when periphyton standing crop changes as a response to P enrichment. Periphyton community succession also could play a significant role in P cycling in the Everglades wetland. P cycling between water column and periphyton mats may differ significantly from P cycling within mats. The mucilage production and microbial activity may be critical processes determining mat architecture. Currently, we are increasing the model complexity, including the number of community types, other environmental variables and the nonlinearity and discontinuity in the functional relationship. The purpose is to capture more realism and to improve our understanding of the interaction among those important processes and variables in the Everglades marshes. We are developing a suite of models progressively in a metaframe modeling approach to answer different questions, along a spectrum of temporal and spatial scales and degree of details and specifications. In this approach, simulations will be performed to link ecological processes at different temporal and spatial scales, reconstruct system dynamics, and examine the scenarios of management or engineering approaches more specifically.

ACKNOWLEDGMENTS This paper benefited from comments by S. Dailey, T. Fontaine, E. Gaiser, P. Geddes, A. Hastings, K. Havens, G. Huxel, J. Richards, L. Scinto, and an anonymous reviewer on various drafts of

12 the manuscript. The sawgrass leaf blades, mosquito songs, heat, and thunderstorms in the Everglades marshes also sharpened and stimulated our thoughts. Financial support was provided by the South Florida Water Management District, the U.S. Department of the Interior South Florida Ecosystem Restoration Program ‘‘Critical Ecosystems Studies Initiative’’ (administrated through the National Park Service), and the U.S. Geological Survey, Florida Caribbean Science Center.

REFERENCES Borchardt, M. A., Hoffmann, J. P., and Cook, P. W. 1994. Phosphorus uptake kinetics of Spirogyra fluviatilis (Charophyceae) in flowing water, J. Phycol. 30, 403–417. Browder, J. A., Gleason, P. J., and Swift, D. R. 1994. Periphyton in the Everglades: Spatial variation, environmental correlates, and ecological implications, in ‘‘Everglades: The Ecosystem and Its Restoration’’ (S. M. Davis and J. C. Odgen, Eds.), St. Lucie Press, Delray Beach, FL. Carpenter, S., Caraco, N. F., Correll, D. L., Howarth, R. W., Sharpley, A. N., and Smith, V. H. 1998. Nonpoint pollution of surface waters with phosphorus and nitrogen, Ecol. Appl. 8, 559–568. Craft, C. B., and Richardson, C. J. 1993. Peat accretion and N, P, and organic C accumulation in nutrient-enriched and unenriched Everglades peatlands, Ecol. Appl. 3, 446–458. DeAngelis, D. L., Mulholland, P. J., Palumbo, A. V., Steinman, A. D., Huston, M. A., and Elwood, J. W. 1989. Nutrient dynamics and food-web stability, Ann. Rev. Ecol. Systemat. 20, 71–95. Edelstein-Keshet, L. 1988. ‘‘Mathematical Models in Biology,’’ Random House, New York. Flora, M. D., Walker, D. R., Scheidt, D. J., Rice, R. G., and Landers, D. H. 1988. ‘‘The Response of the Everglades Marsh to Increased Nitrogen and Phosphorus Loading. Part I. Nutrient Dosing, Water Chemistry, and Periphyton Productivity,’’ unpublished report Everglades National Park, Homestead. Gatto, M., and Rinaldi, S. 1987. Some models of catastrophic behavior in exploited forests, Vegetatio 69, 213–222. Gleason, P. J., and Spackman, W., Jr. 1974. Calcareous periphyton and water chemistry in the Everglades, in ‘‘Environments of South Florida: Present and Past’’ (P. J. Gleasoon, Ed.), Memoir No. 2, pp. 146–181, Miami Geological Society, Coral Gables, FL. Goldsborough, L. G., and Robinson, G. G. C. 1996. Patterns in wetland, in ‘‘Algal Ecology: Freshwater Benthic Ecosystems’’ (R. J. Stevenson, M. L. Bothwell, and R. L. Lowe, Eds.), Academic press, San Diego. Grover, J. P. 1990. Resource competition in a variable environment: Phytoplankton growing according to Monod’s model, Am. Nat. 136(6), 771–789. Grover, J. P. 1991. Non-steady state dynamics of algal population growth: Experiments with two chlorophytes, J. Phycol. 27, 27–79. Gunderson, L. H. 1994. Vegetation of the Everglades: Determinants of community composition, in ‘‘Everglades: the Ecosystem and Its Restoration’’ (S. M. Davis and J. C. Odgen, Eds.), pp. 323–340, St. Lucie Press, Delray Beach. Hopkinson, C. S., Jr., Mulholland, P. J., Pomeroy, L. R., Twilley, R. R., and Whigham, D. F. 1997. ‘‘External Panel Report to the Florida Department of Environmental Protection, Overview and Evaluation of Everglades Nutrient Threshold Research for the period October 1996 to May 1997.’’

Dong et al. Hwang, S., Havens, K. E., and Steinman, A. D. 1998. Phosphorus kinetics of planktonic and benthic assemblages in a shallow subtropical lake, Freshwater Biol. 40, 729–745. Kempf, J., Duckstein, L., and Casti, J. 1984. Relaxation oscillations and other non-Michaelian behavior in a slow-fast phytoplankton growth model, Ecol. Model. 23, 67–90. Lean, D., Rechkow, K., Walker, W., and Wetzel, R. 1992. ‘‘Everglades Nutrient Threshold Research Plan,’’ Report of the Research and Monitoring Subcommittee, Everglades Technical Oversight Committee. Loehle, C. 1989. Catastrophe theory in ecology: A critical review and an example of the butterfly catastrophe, Ecol. Model. 49, 125–152. McCormick, P. V., and O’Dell, M. B. 1996. Quantifying periphyton responses to phosphorus in the Florida Everglades: A synopticexperimental approach, J. N. Am. Benthol. Soc. 15, 450–468. McCormick, P. V., Rawlik, P. S., Lurding, K., Smith, E. P., and Sklar, F. H. 1996. Periphyton water quality relationships along a nutrient gradient in the northern Everglades, J. N. Am. Benthol. Soc. 15, 433–449. McCormick, P. V., and Scinto, L. J. 1999. Influence of phosphorus loading on wetlands periphyton assemblages: A case study from the Everglades, in ‘‘Phosphorus Biogeochemistry in Subtropical Ecosystems’’ (K. R. Reddy, G. A. O’Connor, and C. L. Schelske, Eds.), pp. 301–320, Lewis Publishers, Boca Raton. McCormick, P. V., and Stevenson, R. J. 1991. Mechanisms of benthic algal succession in lotic environments, Ecology 72(5), 1835–1848. McCormick, P. V., and Stevenson, R. J. 1998. Periphyton as a tool for ecological assessment and management in the Florida Everglades, J. Phycol. 34, 726–733. McCormick, P. V., O’Dell, M. B., Shuford, R. B. E., III, Backus, J. G., and Kennedy, W. C. 2001. Periphyton responses to experimental phosphorus enrichment in a subtropical wetland. Aquatic Botany 71, 119–139. McCormick, P. V., Shuford, R. B. E., III, Backus, J. G., and Kennedy, W. C. 1998. Spatial and seasonal patterns of periphyton biomass and productivity in the northern Everglades, Florida, U.S.A., Hydrobiologia 362, 185–208. Momo, F. R. 1995. A new model for periphyton growth in running waters, Hydrobiologia 299(3), 215–218. Monod, J. 1950. La technique de culture continue, theorie et applications, Ann. Inst. Pasteur (Paris) 79, 390–410. Pomeroy, L. R. 1960. Residence time of dissolved phosphate in natural waters, Science 131, 1731–1732. Qian, S. S. 1997a. An illustration of model structure identification, J. Am. Water Res. Assoc. 33(4), 811–824. Qian, S. S. 1997b. Reply to discussion by William W. Walker, Jr., J. Am. Water Es. Assoc. 33(6), 1403–1404. Rader, R. B. 1994. Macroinvertebrates of the northern Everglades: Species composition and trophic structure, Florida Sc. 57, 22–33. Reynolds, C. S. 1984. ‘‘The Ecology of Freshwater Phytoplankton,’’ Cambridge University Press, Cambridge, UK . Richardson, C. J., Bush, M., Cooper, S., Craft, C. B., Qualls, R. G., Rader, R., Ramanowicz, E., Stevenson, J., Vaithiyanathen, P., and Vymazal, J. 1996. ‘‘Effects of Phosphorus and Hydroperiod Alterations on Ecosystem Structure and Function in the Everglades,’’ 1996 Annual Report to Everglades Agricultural Area Environmental Protection District, Duke Wetland Center. Rodriguez, M. A. 1987. Estimating periphyton growth parameters using simple models, Limnol. Oceanogr. 32(2), 458–464. Rosenzweig, M. L. 1971. The paradox of enrichment: Destabilization of exploitation ecosystems in ecological time, Science 171, 385–387.

Catastrophe and Hysteresis in Periphyton Sas, H., (Coord.) 1989. ‘‘Lake Restoration by Reduction of Nutrient Loading: Expectations, Experiences, Extrapolations,’’ Academic Verlag, Richarz. Scheffer, M., Rinaldi, S., Gragnani, A., Mur, L. R., and Nes, E. H. 1997. On the dominance of filamentous cyanobacteria in shallow turbid lakes, Ecology 78(1), 272–282. Scinto, L. J., and Reddy, K. R. 1997. ‘‘Phosphorus Retention by Periphyton,’’ final report to South Florida Water Management District, C-6630, University of Florida, Gainsville, Fl. Steward, K. K., and Ornes, W. H. 1975. Assessing a marsh environment for wastewater renovation, J. Water Pollut. Contr. Fed. 47, 1880–1891. Swift, D. R., and Nicholas, R. B. 1987. Periphyton and water quality relationships in the Everglades water conservation areas, 1978–1982, Technical Publication 87-2, South Florida Water Management District, West Palm Beach, Florida.

13 Tienssen, H. (Ed.) 1995. ‘‘Phosphorus in the Global Environment: Transfers, Cycles and Management,’’ Wiley, Chinchester/New York. Tilzer, M. M., Gaedke, U., Schweizer, A., Beese, B., and Wieser, T. 1991. Interannual variability of phytoplankton productivity and related parameters in Lake Constance: No response to decreased phosphorus loading? J. Plankton Res. 13(4), 755–777. Vymazal, J. 1995. ‘‘Algae and Element Cycling in Wetland,’’ Lewis, Boca Raton, FL. Walker, W. W., Jr. 1995. Design basis for Everglades stormwater treatment areas, Water Res. Bull. 31, 671–685. Walker, W. W., Jr. 1997. Discussion, Am. Water Res. Assoc. 33(6), 1125–1127. Zeeman, E. C. 1977. ‘‘Catastrophe Theory: Selected Papers, 1972– 1977,’’ Addison–Wesley, Reading, MA.