Development and global sensitivity analysis of a closed ecosystem model

Development and global sensitivity analysis of a closed ecosystem model

Ecological Modelling, 30 (1985) 13-47 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands 13 DEVELOPMENT AND GLOBAL SENSITIVITY...

2MB Sizes 0 Downloads 6 Views

Ecological Modelling, 30 (1985) 13-47 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

13

DEVELOPMENT AND GLOBAL SENSITIVITY ANALYSIS OF A CLOSED ECOSYSTEM MODEL

DORIAN LIEPMANN and GREGORY STEPHANOPOULOS

California Institute of Technology, Department of Chemical Engineering, Pasadena. CA 91125 ~U.S.A.) (Accepted 15 April 1985)

ABSTRACT Liepmann, D. and Stephanopoulos, G., 1985. Development and global sensitivity analysis of a closed ecosystem model. Ecol. Modelling, 30: 13-47. A computer model of a materially-closed marine microcosm, under investigation at JPL, has been developed. The model simulates 1-L microcosms, containing biota from anchialine pools in the Hawaiian Islands, that can sustain a constant population of shrimp for years. The stability of the systems depends on the initial number of shrimp and, once stable, the shrimp do not moult or reproduce. The model is an initial step in the investigation of the endogenous control systems responsible for the stability of homeostatic ecosystems and was developed as part of NASA development of controlled ecological life support systems for extended in the model as two separate foodwebs which are coupled by carbon, nitrogen and oxygen cycles. The resultant system of feedback loops creates an endogenous control system, stabilizes the model, and replicates features of a homeostatic ecosystem. The manner and degre of stress required to disrupt the endogenous control and destabilize the model is qualitatively discussed. Finally the Fourier Amplitude Sensitivity Test (FAST) is presented as an extremely efficient and practical alternative to Monte Carlo approaches for global sensitivity analyses and is used to analyze the effects of parametric uncertainty on the response of the model.

INTRODUCTION A t s o m e p o i n t in the d e v e l o p m e n t of s p a c e travel, biological entities o t h e r t h a n h u m a n s will h a v e to b e i n t r o d u c e d i n t o the l i f e - s u p p o r t systems. A t the v e r y least, this will h a v e to b e d o n e for f o o d p r o d u c t i o n . T h e rest of the recycling s y s t e m c o u l d b e a n y w h e r e a l o n g a s p e c t r u m of p h y s i o - c h e m i c a l a n d / o r biological systems. O n e s p e c t r a l e x t r e m e w o u l d b e a s y s t e m with o n l y higher p l a n t s a n d a p u r e l y p h y s i o - c h e m i c a l n u t r i e n t recycle, entirely d e p e n d e n t on e x t e r n a l control. T h e o t h e r e x t r e m e w o u l d involve a c o m p l e x o l i g o t r o p h i c s y s t e m , w h i c h h a s the a d v a n t a g e s of stability a n d r o b u s t n e s s a u t o m a t i c a l l y p r o v i d e d b y an e n d o g e n o u s c o n t r o l system. Less e x t r e m e 0304-3800/85/$03.30

© 1985 Elsevier Science Publishers B.V.

14 would be a more eutrophic system requiring more external control in order to prevent succession and to maintain long-term stability. Such an ecosystem would have the advantage of high production efficiency and low biomass. Any failure of external homeostatic mechanisms or a change in the activity of the support organism would be fatal, since the system would then undergo ecological succession and might no longer support the astronauts. Along this spectrum lies an optimum point for basing the nutrient regeneration system for a Controlled Ecological Life Support System (CELSS). The decision about such optimum and the capacity to engineer biological systems will require a better understanding of ecosystem dynamic behavior than currently exists. The simulation model reported here, created by mathematically describing an existing, materially-closed ecosystem, represents a first approximation of a tool for the investigation of this ecological spectrum. The biological communities, imported from pools on lava flows in Hawaii and Maui, can adapt successfully to enclosure in sealed flasks and so far have survived for about 2 years at the Jet Propulsion Laboratories (JPL). Success depends on the number of shrimp included in the system. Once stabilized this population remains constant at seven or eight shrimp per flask. The model has been developed as an aid in planning the experimental research necessary to understand how these systems recycle nutrients so efficiently, to observe whether or not they reach steady states, and to investigate the endogenous controls which make the enclosed communities so robust. It is also presented as a first step in theoretically examining homeostatic control and changes in the endogenous nature of oligotrophic ecosystems caused by variations in their nutrient composition. The model is made up of nine compartments whose dynamics and interactions are characterized by mass balances over three essential elements: nitrogen, carbon and oxygen. The biota are divided trophically into five compartments which comprise two interdependent foodwebs. These are interconnected by a common top omnivore, the shrimp, and by the nutrient cycles. The elemental cycles are also coupled. By closing these cycles a set of interconnected feedback loops can be formed, thus creating a simple homeostat. Unlike natural ecosystems, simulation models tend to exhibit increasing parametric instability as they become more complex (Margalef, 1978; Tiwari et al., 1978). The use of two foodwebs and the incorporation of an endogenous feedback control system make the ecosystem model relatively stable and parametrically insensitive. Furthermore, it illustrates how nutrients flow through a microcosm and how the resultant feedback loops serve to stabilize and regulate the ecosystem. And finally, although the model represents a highly simplified system, it can demonstrate some mechanisms by which homeostatic control breaks down.

15 THE ECOSYSTEM The biological communities that evolve in oligotrophic ecosystems could prove useful to bioregenerative systems research. As a result of their highly selective environments, these systems have minimized and essentially eliminated any dependence on external nutrient sources. The minerals utilized by the primary producers are recycled within the specific system, as opposed to eutrophic areas where there exists some form of nutrient input from outside the system. Oligotrophic assemblages also tend to be more stable, though less productive. Attributes which make the anchialine assemblage an interesting model for life support systems development are those which optimize complete nutrient recycling: great specific diversity, a complex food web, and a niche structure saturated with highly specialized biotic components. These characteristics reflect an ecologically mature community. The biotic community has evolved system level interactions and characteristics that sustain a constant balance in its biotic and abiotic structures. The functional abilities to complete the essential nutrient cycles and to endure stress are complementary to the structural balance maintained in the ecosystem.

Laboratory system The laboratory ecosystems are 1-L flasks containing 0.862 L of 'instant' seawater, phytoplankton, and some shrimp kept in a water bath at 20°C and exposed to 12 h of light per day. The biota is imported from anchialine ponds in Hawaii. System success apparently depends only on the number of shrimp introduced into the flasks. A handful of plankton, sent to JPL in freeze-dried form, supply the rest of the biota. The systems are viable and visually stable with seven and eight shrimp. If any more were introduced, the systems would crash rapidly, except with fifteen or sixteen shrimp in this case, the system can persist, but the shrimp population eventually declines to seven or eight members. Once the number of shrimp has stabilized, it remains constant and even though the shrimp graze continuously on what appears to be a copious supply of phytoplankton, they do not molt or reproduce. The shrimp (Decapoda, Atyidae) are small, red, atyid prawn about 7-14 m m in length. After their initial discovery early in the century, the shrimp were reclassified by Holthius (1964) as Halocardinia rubra. The species has been found only in shallow, saltwater pools in the lava flows of the Hawaian Islands (Holthius, 1964) and had remained monotypical of the new genus, Halocardinia, until 1975. A more complete ecological description of the anchialine pools is given by Wong (1975) in her master's thesis. David

16 Lawson at JPL desiccated and weighted a sample of 12 shrimp: one of these was singularity large weighing 7.9 rag, while the rest were approximately the same size giving an average weight of 4.66 + 0.15 mg (personal communication, 1982). Plankton is the visually and trophically dominant biotic component in the microcosms and contains a multitude of interdependent species. Because of the complexity of the plankton and the diversity of their biotic structure, each flask is visibly different and the ecosystems do not evolve uniformly. The algal composition of one flask has been analyzed by Frieda Reid at the Institute of Marine Resources, U.C. San Diego. She found that the flank contained mostly small benthic diatoms along with a few blue-green algae and small flagellates; the dominant diatom contained in the flask was a species of Amphora, similar to A. veneta (J. Hanson, JPL, personal communication, 1982). The remarkable feature of this ecosystem are its robustness with respect to initial (non-shrimp) biomass and its practically binary behavior due to the number of shrimp. These omnivores not only graze on the algae but also control the biomasses of the micro-predators. This has s buffering effect on the dynamics of the lower trophic levels and protects the smaller predators from larger carnivores by preferential ingestion. Grazing also selects for the adaptive changes that increase the efficiency of energy transfer to the top organisms. The metabolism of the heterotrophs resupplies the inorganic nutrients for primary production and closes the elemental cycles. The various chemical and biotic interactions form a large set of interconnected feedback loops that create a high degree of interdependence in the ecosystem. Furthermore, these functional links between the organisms couple their dynamic behavior and, thus, provide an endogenous control system.

Attractive features of the laboratory system This sytem is of interest due to its physical characteristics and unique behavior. It essentially has only two phases; the light source is a controlled step function; and the water-bath makes it isothermal. Thus the simplifications ordinarily used in the development of simulation models are inherently consistent with the actual system. Its size, economy and robustness make it attractive from a laboratory point of view. The systems are judged stable in reference to their ability to adapt asymptotically to their new environment after a severe perturbation and they can serve as good models of confinement stresses developed in closed systems. The ecological consequences of predator pressure are of importance to CELSS not only because astronauts will be trophic components, but because this type of pressure controls the behavior of an ecosystem. Herbivore

17

grazing appears to increase the rate of primary production in aquatic microcosms, since the resultant decrease in the primary standing crop does not result in less food material for consumers (Hargrave, 1970: Cooper, 1973). Harte et al. (1980) found that their 50-L mesotrophic microcosm was unable to support a small, 1-g fish and estimated that the microcosm would have to hold approximately 2.5 x 103-10 s L in order to satisfy the fish's natural requirements. In comparison, the anchialine ecosystem is far more productive. It seems reasonable to infer that either the rate of primary production in these flasks is substantially greater than is normally found experimentally, a n d / o r Halocardinia has extremely modest and flexible nutritional requirements. Similarly, the presence of phagotrophs has similarly been shown to increase the growth efficiency and ingestion rate of bacteria as well as improve their ability to decompose refractory waste material (Bick, 1973; Stout, 1973; Barsdate et al., 1974; Baath et al., 1981: Sherr et al., 1982). Barsdate et al. (1974) found that protozoan predation increased the rate and efficiency of phosphorus mineralization in a pond. Bacteria, having the lowest carbon to nitrogen ( C : N ) ratio in the food web, release little or no inorganic nutrients, for which they have a high demand. Instead, they concentrate nutrients for the next trophic level in the cycle (cf. Sorokin, 1978) and nutrient regeneration is then performed by the phagotrophs in the system. Their food source has a higher nitrogen concentration than they need for growth so the excess nutrient is respired and excreted as ammonia or nitrate. Legner (1973) has reported evidence indicating that microphagous ciliates release metabolites which stimulate microbial activity, creating a nutrient subcycle with its own feedback mechanism. Legner also suggests that this may be one of the causes of the extreme effectiveness of aquatic ecosystems. Protozoa further serve energy transfer along the food chain by feeding on particles, such as bacteria, which are inefficiently grazed by larger zoopplankton and provide a source of food that is assimilated much more readily by the shrimp (Porter et al., 1979; Crisman et al., 1981). While it is recognized that a regeneration of the nutrients from detritus is essential to the overall function of the ecosystem, it also has been found that there exists much faster nutrient exchange extrinsic to the organic substrate. To a large extent, this mineral cycling controls decomposition processes and the structure of communities of organisms belonging to decomposer food chains (Barsdate et al., 1974). TIlE MODEL

In order to simplify the mathematical description of the system we ignore any spatial variations in the flasks and assume them to be isobaric. With these simplications and the isothermal laboratory conditions, the model

18 consists of only mass balances which, under the lumped parameter assumption, take the form of ordinary differential equations. The system is further simplified by selecting a manageable number of compartments that approximate the behavior of the original in ways relevant to the questions at hand. Only three elements have been considered, and carbon and nitrogen mass balances have been used to represent the material cycles in the closed ecosystem as the basis for mathematically characterizing trophic-dynamic behavior and its effect on nutrient regeneration. Carbon relates to the relative size of an organic compartment while the C : N ratio reflects its caloric content. The respiratory quotient can be correlated with the net energy gain from catabolizing organic matter, interconnecting all three elemental cycles. The model also includes constitutive relationships which quantitatively characterize observed biotic rates. Figure 1 diagrams the interactions between the compartments (Table 1) includes in the model. The microtrophic foodweb has been included in the system to further couple trophic-dynamic feedback and nutrient regeneration. The more traditional components, algae and herbivorous protozoa, make up a separate web. The omnivorous shrimp interconnect the two biotic systems. One advantage of this architecture is that the microtrophs can be modeled more accurately, since the herbivores and bactivorous protozoa are dynamically quite different. The shrimp are assumed to have constant biomass. The state variables and differential equations which make up the model are listed in Tables 1 and 2 (names of the variables and parameters are unconventional, but make the equations less cryptic).

GAS PHASE q

SHRIMP

Ii

DISSOLVED

DISSOLVED PROT~

/

INORGANIC MATTER

PROT2 : I1

ALGAE

1

II, BACT~

PARTICULATEORGANICMATTER

Fig. 1. Process flows in the model.

ORGANIC

MATTER

19 TABLE 1 State variables Varable (1) Protozoa 1, herbivore (2) Protozoa 2, bactivore (3) Algae (4) Bacteria (5) Dissolved organic carbon (6) Dissolved inorganic carbon (7) Atmospheric carbon dioxide (8) Dissolved organic nitrogen (9) Dissolved inorganic nitrogen (10) intracellular algal carbon (11) Intracellular algal nitrogen (12) Particulate organic carbon (13) Particulate organic nitrogen (14) Dissolved oxygen (15) Atmospheric oxygen

Symbol Pl P2 A B DOC DIC CO 2 DON DIN AC AN POC PON DO 02

Unit ~gCL i ggCL ~gL i txgCL i ~gCL ~g C L I ~g CO2-C L ~ gas ~gN L p~gNH 4 N L g C g ~algae g N g- ~ algae }xg N C 1 ~gN L ~g O: L i ~g 0 2 L i gas

Heterotrophic growth T h e g r o w t h rates of the bacteria and p r o t o z o a are m o d e l e d by multiple s a t u r a t i o n kinetic m o d e l s expressions. F o r these species the ingestion a n d r e s p i r a t i o n rates b o t h d e p e n d on the availability of food. T h e r e f o r e if the o r g a n i c m a t t e r in the m i c r o c o s m decreases, then the rates of n u t r i e n t m i n e r a l i z a t i o n a n d p r i m a r y p r o d u c t i o n will increase and thus stabilize the system. T h e p r e d a t o r - p r e y i n t e r a c t i o n s of the p r o t o z o a are stabilized b y the resultant decrease in p r e d a t o r efficiency at low p r e y c o n c e n t r a t i o n s ( H a m ilton and Preslan, 1970; D a n s o and A l e x a n d e r , 1975: Bader et al., 1976). Using c a r b o n or biomass c o n c e n t r a t i o n , x/, in M i c h a e l i s - M e n t e n kinetic expressions, the specific rate expressions for c o n s u m e r i are the following.

- - Carbon uptake

L

s,i

Xj

--- Growth efficiency EFF,.= ~

K,,,+xj

- - Biotic" production GRTH,=KiY

, VMXi

xj K~ i + x l

xz Kt,i+xl

20 TABLE 2 Differential equations ~f'P1 = X p I ( K p 1 EFFm XAC CUPA - D p 1 ) - S u p , P l

Xp2 = Xp2 ( Kp2 EFFp2 C U P s - Dp2) - Sup.p 2 fi(A = AGRTH -- X A D A -- Xp1 CUPA -- SUP,A

J(s = XB (KB E F F s C U P o o c - D . ) -

Xp2 CUP s - Sup,s + SBX(Scu e - SCEx )

XDOC = Xp1 (1 -- Kpl ) XAc CUPA + Xp2 (1 - Kp2)CUPB + (1 - K A ) XACAGRTH

-- X B K B CUPDo c + SOX SCEx + Xpo c ( K D O C + Y D O C XDOc C U P ~ / X D o c ) Z~'DiC = Xp1 (1 -- EFFp~ ) Kp~ XAC CUP A + Xp, (1 - EFFp2 ) K 2 C U P s - A cuP

+ X B ( 1 - E F F B ) K B CUPDo c + (1 -- SOX) SCEX + D L A R E A ( X c o 2 H c - XDIC) X co2 = - DE A R E A ( Xco2H¢ - XDIC ) / V B U F XDOh = Xpl(1 - Kp1) XAN CUPA + Xp2(1 - Kp2 ) N : C B C U P B + (1 - KA)XANAGRTH -- X B K B ( X D o N / X D o C ) C U P D o C + S O X SNE x

+ Xpo N ( K D O N + Y D O N X s C U P o o c / X D o c ) "~DIN = XPa NEXp1 4" Xp2 NEXp2 - ANU P + X B NEX s + (1 - SOX)SNEx

XAC = ( ~ c u P -

XAcAG~TH)/XA

XAN = ( A N u P -- X A N A G R T H ) / / X A X p o c = Xp, Dp1 + Xp2Dp: + XADAXAc + X s D . + (1 - s a x ) ( S c u p

- SEE×)

- Xpo c ( K D O C + YDOC X B C U P D o c / X D o c ) Xpo N = X m D p l N : C m + Xp2Dp2 N : Cp + XADAXAN + XBD B N : C~ + (SNuP -- SNE×)--SBX N : C . ( S c u p - SCE×) -- Xpo N ( K D O N + Y D O N X s CUPDo C / X D o C ) )f'DO = RQ[ Xpl (1 - EFFp1 ) Kp1XAC C U P A + Xp2 (1 - EFFp2 ) Kp2 C U P B - A c u p

+ X B ( 1 - E F F s ) K s C U P D o c ] + DE A R E A ( X o 2 H o - XI)o) X02 = - DE A R E A ( X o H 0 - X D o ) / V B U F where CUPDo c = VMX B Ks,B + XDOC

0.04Ks,B + XDON ]

S I Z E / m i n ( 0 , X i - STH,i } SpREF =

4 i=1

SIZE i rain(0, X i - STH,i }

SuP,i = SSGRzSpREF,i m i n ( 0 , X, - STH,i } SCU P = Sup,p I -1- Sup,p 2 -1- SUP,AXAc -[- SUP,B SNU p = Sup.p ] C : N p l -J- Sup,p 2 C : N p 2 + SUP,AXAN -}- SUP,B C : N B

SNEx = rain{ SSREsP, SNUV } SCEx = min ( SNEx ( S c u v / S N u P), Scu P }

21 where K i is the (dimensionless) assimilation efficiency, and Y, the maximum metabolic efficiency (also dimensionless); VMX i - ~i.ma×/Ki Yi represents the maximum specific ingestion rate (h-l), determined from published specific growth rate data; Ks. ~ is the substrate concentration at which CUP, = V M X J 2 , or the half-saturation constant for ingestion, and Kt. ~ the halfsaturation constant for metabolic efficiency. Organic excretion is proportional to ingestion. - - Nitrogen mineralization

NEX, =

CUP, ( N ' C j - N .C, EFF,)

For the herbivorous protozoa and the bacteria, N : C j is a variable and equals A N / A C or D O N / D O C . Other microbial rates, which depend on the energy content of the dissolved organic matter, also reflect any variations in the latter term (Fenchel and Blackburn, 1979). If nitrogen ingestion is insufficient to meet metabolic needs, then NEX becomes negative and N H ] is assimilated rather than excreted. The bacteria have an additional dependence on the C : N ratio of the DOM which reduces microbial ingestion as the availability of organic nitrogen decreases; this term discourages bacterial depletion of NH 4 and inhibits competition between the bacteria and alga for the DIN. A major simplification has been made by grouping all heterotrophic bacteria in one compartment. The resulting kinetic models of mixed bacterial cultures require curve-fitting in order to be accurate. Parametric estimates for the kinetic expressions are complicated by the effects of grazing pressure (Barsdate et al., 1974); mixed cultures (Vaccaro and Jannasch, 1967; Jannasch, 1967, 1968); mixed substrates (Williams, 1970; Wright, 1974; Robinson et al., 1982); succession (see Fukami et al., 1981); and surface growth, especially that on the slime film of the shrimp feces (Zobell, 1943; Paerl, 1974). Finally, the lumped parameter assumption makes it impossible to consider anaerobic bacteria as a separate biotic entity, since this would require enlarging the scope of the model to include spatial variations, at least, in the concentration of dissolved oxygen. The shrimp are not a variable in the model. The assumption that their biomass remains constant assures a constant predator population, representative of a mature community, and includes the stabilizing effects of an omnivore in the simulation. The shrimp's low, constant grazing and excretion rates protect the producers at the bottom of both foodwebs by dampening any violent fluctuations among the micropredator populations and by providing a constant source of abiotic substrates. Because the shrimp never molt or reproduce in the stable microcosms, one can safely assume that they are bordering on starvation. Fecal losses and basal metabolism completely

22 balance ingestion. The shrimp, therefore, must eat everything they can possibly find, which would be proportional to the availability of their prey. The expression for predation by the shrimp includes a normalized grazing efficiency based on the prey's size and a threshold concentration, below which the prey are inaccessible. The former term allows greater flexibility in modeling omnivorous grazing in a system where size is considered to be a major factor in prey selection (Griffith, 1975). In the simulations performed so far, the parameter SIZE B has been set to zero because zooplankton cannot efficiently ingest free-floating bacteria (Marshall, 1973; Paerl, 1974; Porter et al., 1979). The ecological validity of these factors is reviewed in more detail by Conover (1978) and Steele and Mullin (1977). The expressions and parameters in the ecosystem model have been adapted in part from Canale and co-workers' (1976) model of a plankton based foodweb. The respiration of a shrimp is given by the mass of nitrogen it must catabolize per hour for its basal metabolism. Therefore, the NH~- and DON excreted by the shrimp remain constant while the CO 2 and DOC are proportional to the overall N : C ratio in its ration. If the shrimp excrete less nitrogen than their basal requirement, they are starving. About 15% of the nitrogen excreted by tropical meso-zooplankton is organic (Johannes and Satomi, 1966; Mayzaud, 1973; Parsons et al., 1977). Primary production

Algae is represented by three state variables, biomass and the intracellular concentrations of carbon and nitrogen. The approach used by Lehman (1975) and J0rgensen (1976) has been used to characterize the dynamic behavior of the algae, assuming that benthic diatoms are predominate in the compartment. Separate rate expressions are used in the model for algal growth and nutrient assimilation. Growth depends on the internal nutrient concentrations (Droop, 1973). The use of three state variables to define one compartment increases the complexity of the model but it results in a more realistic description of the system (J0rgensen, 1976). - - Algal nitrogen uptake

ANUP = AVMXAN AN~max--~A~Nmin K1.AN+ DIN - - Algal growth rate

[ AN - ANmin] [ AC - ACmin] AGRTH= AAmax AN AC Because the nutrient factors are always less than 1, the parameter Areax is adjusted so that, under optimal conditions, the algae can grow at its maximum rate. Organic excretion is proportional to the growth rate.

23 The model assumes that the phytoplankton does not excrete nitrogen while respiring or, if so, that the NH~- is reabsorbed immediately. Algae respires inorganic carbon at a rate that depends on a maximum specic rate, YAC, and the ratio of AC to ACma × raised to 0.67 power, in accord with the empirical 'surface law' of respiration (Lehmann et al., 1975; cf. Platt et al. 1977). The net photosynthesis, then, is the CO 2 absorbed minus the CO 2 respired. - - Primary production

A c u p = AVMXAc ACma x __ ACmin ][ K s , A C ~ DIC

LT(t, FLT)

- A YAc [ A---Cma l When the light is assumed constant, the rate expression for algal assimilation of carbon is analogous to that for nitrogen. The additional term, LT(t, FLT), provides diurnal cycles. Since the available light is assumed not to limit carbon assimilation, LT(t, FLT) is a 'box' function: LT = 1 for FLT hours starting at the beginning of the day and LT = 0 at 'night' when photosynthesis stops. DISCUSSION OF SIMULATION RESULTS The 15 ordinary differential equations (Table 3) were integrated numerically using EPISODE, an improved version of GEAR's method developed by Hindemarsh and Byrne (1975) at the Lawrence Livermore Laboratories. Because EPISODE uses a truly variable step size, it is more stable than its predecessor and can integrate extremely stiff sets of equations. The model uses single precision and allows a maximum relative error of 0.001 per integration step. Total elemental mass balances are calculated to check global error terms. Each graph presented in Figs. 2 through 13 has been generated from over 1000 data points. The simulations presented in the following are the first trials performed and use the initial estimates made for the parameter values (Table 4) from available literature. Kinetic parameters reported for species with high growth rates were chosen to characterize the algae, bacteria, and protozoa. The a priori reasons for these selections were to account for the influences of predation, mixed cultures, and wall-growth. The first four sets of graphs (Figs. 2-5) demonstrate the response of the model with and without diurnal forcing. Runs 1 and 2 are the nominal responses of the model from the FAST sensitivity analysis discussed below.

24 TABLE 3 Parameters Symbol

Units

AREA ACmax

cm2 L-1 water g C g- 1 algae g C g- ~ algae g N g-1 algae g N g 1 algae cm2

ACmi m

ANmaX ANmi. DL Di

Hc

Ho KDOC KDON N:C RQ S S~Ez

SRESP STH SBX SIZE SOX VBUF YDOC YDON

Description

Volume specific area of the gas-liquid interface Maximum internal algal carbon concentration Minimum internal algal carbon concentration Maximum internal algal nitrogen concentration Minimum internal algal nitrogen concentration Molecular diffusivity of O 2 and CO 2 in sea water h- 1 Specific death rate of species i Distribution (mass law) constant for carbon [CO2(aq)]/[CO2 (g)] dioxide at 1 atm; the pressure is assumed to remain constant in the flask Distribution (mass law) constant for oxygen [02 (aq)]/[O2 (g)l First order rate constant for the dissolution of POC First order rate constant for the dissolution of PON gNg -1C Biomass ratio of nitrogen to carbon Average respiratory quotient assumed for the biota g O g -I C Number of shrimp in the ecosystem Constant filtering rate of the shrimp L h- 1 per individual i~g N h 1 per individual Rate of nitrogen excretion due to shrimp respiration I~g prey L-1 Threshold prey concentrations for grazing by the shrimp Fraction of bacteria contained in fecal pellets excreted by the shrimp (Johannes and Satomi, 1966) Shrimp electivity coefficients based on size (volume) of prey (Canale et al., 1976) Fraction of the shrimps' metabolic losses excreted in organic form (Corner and Davies, 1971) Volumetric ratio of gas to liquid in the flask Proportionality constant for the dissolution of POC due to microbial activity (Tiwari et al., 1978) Proportionality constant for the dissolution of PON due to microbial activity.

M o s t of the p a r a m e t e r s have the values listed in T a b l e 4. T h e four p a r a m e ters a n a l y z e d in Section 5 are set to the m e a n values o f their global u n c e r t a i n t y . T h e o t h e r p a r a m e t e r s have the following values: the m a x i m u m rate o f p h o t o s y n t h e s i s equals 0.10 ~g C O 2 - C h -1 per m g A; the m a x i m u m ingestion rate o f the bacteria is 0.75 ~g D O C h - 1 per m g B; the respiration rate of the s h r i m p is 0.17 ~g N h -1 per shrimp. I n b o t h cases the m o d e l is able to rapidly reach a stable state a n d d a m p out the large fluctions c a u s e d b y the initial c o n d i t i o n s . F o r c i n g p r i m a r y p r o d u c t i o n generates oscillations t h r o u g h o u t the system (Figs. 4 a n d 5). Because the light f u n c t i o n is d i s c o n t i n u o u s , the diurnal

0.20

(x0)

Data:

Inttia[ condttton

Data:

Integration parameters

Data:

M a s s transfer across g a s - l i q u i d interface

Data:

P O M dissolution constants

Data:

AIgae

Data:

Shrimp

STH

9.0

DON 50.0

PI 69.0

days 250

Kt 060

0.008

K DOC

A max 0.55

S

40.0

0,156 0.50

N:C SIZE

D

0.40 0.002

0.85

150 100

P1

Parameters

Y

Kt K

Ks

VMX

Data

Block d a t a m~xtule for e c o s y s t e m m~udel

TABLE 4

0.7

28.0

68.0 DIN

P2

hours 4.0

155.0

tt t

KDON 0.008

0.65

A C , .....

SRISP 0.15

0.45 35.0

0.156

0.002

(I.45

50.0 0.9

150

P2

AC 0.55

A 825

0.0

To

0.0268

Ht)

Y DOC 0.05

A C mM 0.35

S(3RZ 0.025

450

0.05

0.01 0.004

0.85

0.02 25.0

AN

500 AN 0.076

B

12.0

FLT

VBUF 0.25

0.05

YDON

0.12

A N .....

0.15

SOX

0.0

0.0

0.26

0.001

1.0 0.60

250 100

0.95

AC

PON POC 1350

175

DI(" 2.28E4 DOC 374

AREA 0.06

ANmm 0.03

0.05

SBX

0.01

0.10 500

CO: 147 DO 7.2t.3

02 85E5

t~ Gel

26 BIOTIC

CARBON 1

:1 E, LIJO

01

~D ~O

4

8

12

16

20

24

TIME

28

32

36

40

44

:32

36

40

44

48

(DAYS)

O

"~ xl

t

0

4

8

12

iS

20 TIME

2:4

28

48

(DAYS) 0 ID

o Lr~I0

S

OO&~t

0~0 i~_ rib

a

O

4'

I~,

12

16

20 TIME

2'4

28

(DAYS)

Fig. 2. Temporal behavior of the biotic variables computed without diurnal cycles.

changes generated throughout the system are much more intense and abrupt than normally observed in models of natural systems. Because the algae can assimilate NH~- constantly, it becomes deficient in carbon during the 12 h of darkness. During the next 12 h, unlike a natural system, algal assimilation of

27 DISSOLVED

NUTRIENTS 1

8

-

Z

~0

4

B

12

~6

20

2~,

TIME

2B

32

36

40

44

48

a-2

36

4o

4:4

48

(DAYS)

Z

Z 0

,~

8

12

i-6

20

24

TIME

28

(DAYS)

c,J

co

0

too

0'~" a

/

¢xl

?;

~,

13

1"2

' I-6-

20

24

28

32

36

40

44

4B

TIME (DAYS)

Fig. 3. Temporal behavior of the nutrient variables computed without diurnal cycles.

C() 2 is not limited by light. The rapid changes in algal biomass, carbon, and growth result in increased DOM excretion, death, and grazing rates by the shrimp and herbivorous protozoa. Fluctuations in detrital concentrations and C : N ratios generate oscillations in the bacterial and bactivorous protozoan populations. The system responds just as quickly when the light turns off, but the detrital foodweb dominates the cycles.

28 BIOTIC

CARBON 2

O0

0 <

~:cJ Wc~ U <

m~

cO

0

4

B

12

16

20 24 28 TIME (DAYS)

32

36

40

44

48

~2

16

20

32

36

40

,i4

48

0

~o ~/ O

r ~0

i ,~

B

54

28

TIME (DAYS)

O.

tO~_

~u3 o tLOj tO

o

n-

I

~O

21

B

i2

i6

20

24

28

32

36

210

214

4B

TIME ( D A Y S )

Fig. 4. Temporal behavior of the biotic variables computed with diurnal cycles.

Simulations with constant light The following section focuses on the manner and degree of stress which must be applied to the model in order to disrupt the feedback control loops

29

DISSOLVED

NUTRIENTS 2

~o z

®

o3

z:L o

o C

4

8

12

16

20 24 28 TIME (DAYS)

32

36

40

,44

48

4

8

12

16

20 24 28 TIME (DAYS)

32

36

40

44

48

"T d

z~ 2. 0

LO

0

2, ¸ o

Lo. ~J

2. o

/

~q

vv, VvvvvvvwWwvvwvvW

. "4~" I Oa - J

(.9 ® v

'-'S £q

o

4

8

12

16

20 24 28 TIME (DAYS)

32

36

40

44

48

Fig 5. Temporal behavior of the nutrient variables computed with diurnal cycles.

and how these stresses manifest themselves in the dynamic behavior of the model. The descriptions of ecosystem dynamics refer only to computer-generated simulations and not to physical systems unless explicitly stated. Simulations run without diurnal cycles are generally more stable, but less realistic, than those generated with diurnal forcing; however, under certain

30 parametric and initial conditions the model's behavior becomes highly oscillatory. The superimposed oscillations can obscure aspects of the system's behavior. Sjoberg (1977), for example, noted that instability could not be discerned in his model when time varying forcing functions were included. The overall ratio of carbon to nitrogen is a key factor in the performance of the ecosystem. When the total amount of nitrogen in the system is increased, there is a corresponding increase in the biomass. Furthermore, if the overall C : N ratio remains greater than that of the biota, then nitrogen will limit the rate of production in the ecosystem. Because algae can have the lowest intracellular nitrogen content, a major fraction of the new biomass will be phytoplankton. Therefore, an increase in the relative availability of nitrogen can dramatically shift the balance between production and respiration and, in order for the system to remain stable, a balance must be maintained between biotic nitrogen and DIN. This requirement, on a smaller scale, expresses the complementary existence of chemo- and homeostatis. The availability of nitrogen also influences mineralization by the bacteria and so provides an extremely effective lever with which to stress nutrient-regenerative abilities of the ecosystem. The overall C : N ratio of Run 3 and Run 1 differ by only 4%: 18.24 versus 17.60. Yet, this difference and the parameter changes, on the average have increased the algal carbon by over 60% or about 200 ~g C L - l ; have tripled the bacterial population; and have depleted the NH~ concentration by a factor of 8 (Figs. 6 and 7). The biota in Shriner's Pond reacted similarly to nutrient enrichment in experiments done by Kerr and co-workers (1972). The drop in the system's DIN is due to the rise in the biotic nitrogen and is reflected by an increase in the A C : A N ratio from 5.86 to an average value of 10.2. Nitrogen mineralization by the bacteria is always negative, fluctuating between - 0.45 and - 0.82 ~g N H ~ - N L - 1 h - 1 (Table 5). Except for the rate of primary production, the simulation that generated Fig. 8 is parametricaly identical to Run 3. The maximum algal growth rate has been decreased by about 18% to 0.15 h -1 (Area x = 0.45 h - l ) . The model responds similarly to the earlier simulation, but the envelopes of the output variable are shaped differently. Not only does the variation in Area× change the behavior of the phytoplankton variables, but the resulting shift in the system's oscillations also modifies the harmonic interferences and modulations in the model's output. The decrease in Area × has two basic effects: (a) by lowering the rate of primary production, it decreases the total biomass in the system; and (b) by making growth relatively less dependent on intracellular nitrogen,,it decreases the C : N ratio of the algae. This, then, decreases the C : N ratio of the detritus and the range within which BNE x oscillates (Table 5).

31

BIOTIC

O

3

CARBON

b

~ 8 j

20

40

60

20

40

60

80

100

120 140 TIME (DAYS

160

180

200

220

240

o IJJ {,q <~ Lq

~t

o "~'0

80

100

120 140 TIME (DAYS)

160

180

200

220

!40 ,,q,. o

~,,l~i~l ~ ~ l'ilti(lti !li lit Ill ,v~t

~- ~

~0

It' i I~ !l',r ~U~

2'O

40

,~

80

~ 160

,

[o

HlvfIltvt lVlIIl~Vll~ttvH,"ql.

120 140 TiME (DAYS)

160

1gO

260

~

250

240

Fig. 6. Response of the biotic variables in Run 3 (without diurnal forcing): the system's total C:N ratio equals 17.6.

More nitrogen is introduced into the system by increasing the initial PON concentration which lowers the overall C : N ratio to 16.6 and further increases the biotic content of the ecosystem. The oscillations, however, have constant amplitudes, but are quadruple the size of the previous runs (Fig. 9).

32

o0

LJ z~

<

~o to

0

20

4b

6b

80

160

120

TIME

+40

160

180

260

220

240

160

180

200

220

240

(DAYS)

+I ~0 r+_j Z

~ :::k 0~.

~o

20

40

60

80

1 O0

120

140

TIME (DAYS) oo ai ('M

0 CO' O ,,~..

aJ

0 O ~LI3

% ~u

O :i ~0

aJ U £3

tq.

0

m'O

20

40

60

80

1OO

120

140

160

180

200

220

240

TIME (DAYS)

Fig. 7. Response of the nutrient variables in Run 3 (without diurnal forcing): the system's total C : N ratio equals 17.6.

The dissolved inorganic, organic, and particulate nitrogen increase; the latter increases to the point where the organic matter in the system has almost the same C : N ratio as the algae, so the bacteria can alternate between assimilating and regenerating NH~-.

33 INORGANIC

NITROGEN 4

to

z v

100

120

140

160

180

200

220 240 260 TIME (DAYS)

280

300

320

340

360

Fig. 8. DIN concentration (fig NH~--N L 1) in Run 4: Am~X has been decreased by 18% over that in Run 3.

A similar disruption of homeostatic control occurs if one shrimp is removed from Run 3 (Fig. 10). The reduction in predatory pressure produces the same imbalance as the increase in production, even though the overall C : N ratio is higher (Table 5). Compared to the previous response, the fluctuations in this sytem are greater for AC and BNEX, while reduced for the bacteria, protozoa 2, and DIN. The C : N ratios of the detritus and algae are larger than any of the previous runs and, as a result, the uptake of NH + by the bacteria is also larger (0.31-1.70 ~tg NH~- N L -1 h - l ) . The overall C : N ratio is increased by the loss of the nitrogen contained in a 4.5-mg shrimp.

TABLE 5 Effects of the C : N ratio on nutrient regeneration in the ecosystem model Unforced simulations Run

C : N ratios Total

ae

18.2 17.6 17.6 16.6 18.7

5.86 10.2_+0.8 9.1_+0.5 8.5_+1.2 11.0_+1.8

NL-'h~q~

1 3 4 5 6

B~Ex

DIN

Algal C

DOM

POM

(~xg NH~

(~xg NH~

(~g C

5.78 9.15 8.40 7.70 9.55

5.92 9.93 8.99 8.22 10.28

0.54 -0.63_+0.18 -0.24_+0.05 0.02_+0.40 -l.01_+0.70

96.6 10.5_+ 4.0 13.0_+ 3.0 24.0_+15.0 10.0_+ 6.0

324.0 532.0_+28.0 505.0_+12.0 557.0_+38.0 565.0_+55.0

NL ~')

YL-')

34 INORGANIC NITROGEN 5

~J rd

Z

100

120

140

160

180

200

220 240 260 TIME (DAYS)

280

300

320

340

360

Fig. 9. DIN concentration (fig N H ~ - N L -1) in Run 5: the C : N ratio has been decreased to 16.6.

INORGANIC

:t

NITROGEN 6

i I1~ 1tI J!!lII1tittlIt!tljl~!~ 111It1!!!1t I I,~IJI1~

zO

©

O

20

40

6"0

80

100

120 140 160 TIME (DAYS)

180

200

220

2,40

260

Fig. 10. DIN concentration (~g NH4~ - N L -1) in Run 6: same parametric conditions as in Run 3 except with one less shrimp in the system, raising the C : N ratio to 18.70.

Simulations with diurnal cycles The forced simulations demonstrate essentially the same response to the overall C : N ratio, becoming oscillatory as the ratio decreases and the biomass increases. However, the dark phases of the simulation allow the predators in the system to 'catch up' with the increased rate of primary production and dampen the system's oscillations more effectively than without forcing. Large changes in the parametric and/or initial conditions of the model will disrupt its stability, but light forcing makes the simulations less sensitive to shifts in the system's C : N ratio. Although the total produc-

35

tion increases with the amount of nitrogen in the ecosystem, the diurnal light cycles, in effect, limit the carbon available to the phytoplankton. Figure 11 shows the response of a system (the accumulation of numerical error causes the overall decrease of DIN in simulations 7 and 8; diurnal forcing increases the number of functional and jacobian evaluations by an order of magnitude) similar to the Run 3, with an overall C : N ratio of 17.6, to the addition of 12-h light cycles. The diurnal forcing reduces the biomass and biotic nitrogen in the ecosystem (approximately to the concentrations in Run 1) and increases the nutrient availability. The higher maintenance energy of the predators due to the lower and oscillatory prey populations is an additional factor in the higher NH~ concentration. Because of the increase in DIN and the decrease in photosynthesis, the algal C : N ratio varies within 3.1-6.0, versus 9.4-11.0 without forcing. The extrema of the diurnal NH~ cycles indicate that periodic changes in the system's behavior still occur, but their magnitude has been diminished by diurnal forcing. Lowering the C : N ratio of the system to 16.6 finally disrupts the feedback control (Fig. 12), Although in comparison to Run 5 the system's response appears more chaotic, the oscillations, the diurnal cycles seem to reduce the impact on the system. The algal C : N ratio varies only from 3.1 to 6.4, while the C : N ratios of the DOM and POM remain around 5.4 and 5.5, respectively. The bacteria never ingest NH~-. Therefore, unstable behavior in this simulation is not caused by competition between the algae and bacteria for NH 4 , but only by the unbalanced increase in production.

INORGANIC NITROGEN 7

l','~'t~,!#,'~g¢~',v~l~,~'l,~vtl,~t,~'~,W,~,,~, ' ~,a~y~,tt~,ttt,,l,t~~,,,~w~,,¢~,~,~ o

z

ID

t

-

o.q

,D-

0

2o

40

60

8b

100

120

140

[email protected]

leo

2o0

a20

240

260

TIME (DAYS)

Fig. l l . DIN concentration ( ~ g N H ~ - N L i) in Run 7: the effect of diurnal forcing on the third simulation run.

36 INORGANIC

O

NITROGEN 8

ij Z

100

120

140

160

180

200

222 240 260 TIME (DAYS)

280

300

320

340

36:)0

Fig. 12. DIN concentration (~g NH~ ~N L -a) in Run 6.

Unstable behavior

Competition between the algae and bacteria for DIN is common to all of the unforced simulations which were oscillatory. Even though the microbial ingestion rate approaches zero as the DON is depleted, BNEX can still overshoot this curve and become negative. This breaks down the stabilizing effect of the nutrient cycle by creating two parallel negative feedback loops. In the third simulation, the amplitudes of the biotic oscillations are initially less than those in Run 1 (Fig. 2). As the total biomass increases, which forces the DIN concentration to decrease, the simulation becomes unstable. Once the bacteria are unable to obtain enough nitrogen from the dissolved organic matter and start to assimilate NH~-, the availability of inorganic nitrogen diminishes, which limits primary production in the system, and increases the algal C : N ratio. The nitrogen content of particulate and dissolved organic matter, much of which is supplied to the system by the algae, also drops, resulting in further microbial depletion of the available DIN. Because of the high assimilation rate and low C : N ratio of the bacteria, This cycle would ultimately crash the system if not for the influence of microtrophs. The rising bacterial population is checked by the bactivorous protozoa, which lag the microbial oscillations slightly. This starts a sequence of responses as a nitrogen 'shock wave' travels through the detrital foodweb: the bacterial population blooms and depletes the DIN concentration; the bactivorous protozoa bloom slightly after the bacteria; the shrimp ingest relatively more protozoa 2, because of their increased availability. Nitrogen regeneration increases during the second two steps. As predation decreases the concentration of the bacteria, the nitrogen concentrated in their biomass

37 is returned by the predators to the ecosystem. The reduced grazing pressure and increased nutrient availability allows the phytoplankton population to surge, which sets off another sequence: during their rapid growth, the algae release DOM at a proportionately greater rate; the system's population of herbivorous micro-zooplankton peak shortly after the algal bloom; the shrimp ingest more algae due to greater availability, and consequently, excrete more carbon. These events increase the relative fraction of particulate and dissolved waste in the system that is of algal origin. Subsequently, the bacteria face a system with a high substrate concentration and very few bactivores. And the cycle restarts. In the simulations performed so far, this cyclic behavior shows no indication of dying out, nor does there appear to be a reason for it doing so. Once instability of this type develops, the oscillations are conveyed by continuity through the nutrient fluxes and are spread throughout the system. Major factors in modulating the output variable envelopes are fluctuations in the C : N ratios of the algae, DOM, and POM; highly complex harmonic interferences are introduced by the carbon and nitrogen composition of these compartments, which are distinct variables in the model. Therefore, they can oscillate independently of each other with different frequencies and amplitudes. The C : N ratio depends on the dynamic behavior of the two subcomponents: the shape of the envelope depends on the difference between the periods of oscillation of the two nutrient concentrations; their anaplitudes and phase lag determine the range within which the C : N ratio varies.

Effect of grazing pressure on the ecosystem A sufficient heterotrophic population, relative to primary producers, is a preprequisite for homeostatic control. Decreasing the C : N ratio in the ecosystem raises the rates of production at the primary trophic levels of the foodwebs to the point where predation can no longer balance production. The average size of the algal and microbial populations and the amplitude of their oscillations become so large that the grazing pressure of the shrimp is not able to limit the response of the protozoa to the availability of their prey. Once the ability of the protozoa to find food becomes the key factor in limiting their growth rate, ' b o o m and bust' cycles dominate their predator-prey interactions. Because the nutrient cycle dynamically couples the two foodwebs, the fluctuations in the biotic and nutrient composition of the ecosystem reinforces the unstable behavior of both subsystems. Homeostatis and chemotaxis mutually self-destruct. Pearson (1982) observed a similar response to organic matter input to a benthic ecosystem in Loch Eil. The ecosystem was able to remain stable

38

when the organic carbon input was increased by a factor of ten. Up to this point the assemblage was able to transfer the additional carbon up the foodweb, but when the rate was doubled the total biomass of the assemblage decreased by half. Furthermore, after the addition of this much DOC, Pearson noted that the macrofauna did not appear to control the levels of microbial activity as they had previous to the final enrichment.

FAST SENSITIVITY/UNCERTAINTY ANALYSIS

Introduction to the method

An essential consideration in the development and use of a model is the effect of uncertainties, which are inherent in any mathematic characterization of a complex system, on its performance. Two factors are important in such an assessment: the level of uncertainty in each parameter; and the sensitivity of the model prediction to the parameter. For example, the value of a parameter may be uncertain to a large degree but have little influence on performance of the model, or the situation may be reversed. Furthermore, the effects of several parameters varying simultaneously may have a much greater influence combined than individually. A sensitivity analysis which accounts for all feasible combinations of parameter values is specified as global. The Monte Carlo method describes a random process in an explicit manner that provides extreme flexibility and robustness even though it is inherently inefficient. The same information can be obtained much more efficiently using a global sensitivity/uncertainty analysis developed and refined by Cukier, Shuler and coworkers (Cukier et al., 1978). The method was further improved and incorporated in a general purpose FORTRAN program (Koda et al., 1979; McRae et al., 1982). The mathematical basis of FAST and its application to the ecosystem model are described by Liepmann (1983). The Fourier Amplitude Sensitivity Analysis (FAST) provides an almost automatic global sensitivity/uncertainty analysis through the implementation of a general purpose program that: readily accommodates arbitrarily large variations in up to 50 parameters; provides greater statistical accuracy than the Monte Carlo method given the same number of trial solutions; calculates the partial variance for each parameter, a normalized measure of each parameter's relative contribution to the total variance; makes detection of unexpected parametric coupling of multivariate effects easy; provides a means to determine which parameters are not independent with reasonable ease and an economical amount of computer time.

39

FAST analysis of the ecosystem model Both the local and global sensitivity of the ecosystem model to key parameters have been analyzed using FAST via the program written by Greg McRae. The following analysis uses three indices (referred to as state variables) to judge the response of the ecosystem: (1) the total biotic production, the sum of all net growth rates in the system (~g C L ~ h-~): ratio of biotic carbon to total carbon; and (3) the dissolved inorganic nitrogen (DIN) concentration (~tg N-NH~- L l). The Fourier Amplitude Sensitivity Test Application Program (FASTAP) has been used to investigate the sensitivity of the model to key parameters with and without forcing b) diurnal light cycles. Fifteen parameters were considered. As the investigation progressed, the most influential parameters in the system were identified and the scope ol' the analysis was narrowed to focus on these: the maximum algal growth rate: 0.06 > pt.... > 0.23 with diurnal forcing and 0.03 > ~ .... > 0.19 without forcing the grazing rate of the shrimp; 0.020 > S~;~ > 0.035 the maximum predation rate of the microtrophic protozoa: 0.40 >~ VMXp2 >~0.85 (Banse, 1982) the maximum metabolic efficiency of the bacteria: 0.30 > YB > 0.66 The model is very sensitive to the grazing rate of the shrimp, for small perturbations, especially the second state variable. Because the shrimp do not grow but merely decrease the other biota and add to the organic matter in the system, they effect the total biotic production rate less than either the ratio of biotic carbon to total carbon or the DIN (Fig. 14); of these two, the former is most directly influenced by the omnivores' behavior. The maximum growth efficiency of the bacteria, which grow faster than their predator, ranks third except for state variable 3, indicating that the protozoa mineralize slight more NH~-. When these parameters are varied globally, the algal growth rate becomes paramount to the behavior of the system, due to its degree of uncertainty and its (local) influence on the performance of the model (Fig. 15). Because algal dynamics are such a major factor of the system, even for small perturbations, the primary production rate has a greater effect on the NH 4 concentration and total production than the ingestion rate of the shrimp. Earlier results have shown that these two variables are also effected by algal assimilation rates; since the algae depend on the two uptake parameters for A c u e and ANUP,the maximum rate of each of these terms should have about -

-

-

-

-

-

-

-

40

the same impact as Amax on the biotic carbon and NH~- concentrations, respectively. These parameters, in other words, are coupled; i.e. if two were considered random, as both terms move between their minimum to their maximum values, the change in the model's behavior would exaggerate the actual influence of the algal growth rate. When light cycles are included in the analysis, the partial variances of all of the parameters decrease, in part due to the non-linear response of the model to a discontinuous forcing function. The parametric sensitivity of the model does not change however and most of the observations made previously about the unforced system remain applicable. For example, the partial variances tend to stabilize as the initial oscillations of the biota damp out and global variations of the algal growth rate having the greatest effect on the response of the system. In both sets of FASTAP runs, the nominal output values of the biotic production rate and the NH~- concentration differ from the mean response, low in the first case and much too high in the second (Fig. 13). The nominal and mean responses do not agree at very small values of t and persist through the whole simulation run. Large discrepancies between the nominal and expected output values imply that the model's predictive ability will be poor as a result of parametric uncertainty (Cukier et al., 1978). Furthermore, the low production rate and high NH~- concentration are caused by the relatively low rate of algal growth, as the limiting factor in the system; this suggests that the maximum growth a n d / o r uptake rate constants of algae could be too low.

INORGANIC

NITROGEN

i A ,,o z

2 v~

0

80

160

240

320

400

480 .560 640 TIME (HOURS)

720

800

880

9,60

Fig. 13. Nominal output and mean values (*) of the DIN concentration from the global uncertainty analysis of the model (the error bars have a total length of one standard deviation).

41

Multivariate effects If all of the system's parameters were completely independent then the total sum of their partial variances would equal one. Without forcing, the sums of the partial variances returned by FASTAP have all been greater than 0.65: the value suggested by McRae (personal communication, 1982) as a 'rule of thumb' borderline, below which indicates the existence of significant coupling or extreme non-linearity in the problem. The first state variable was affected most by parametric coupling, as one would expect, because the system's predator prey interactions make the rate at which one component of the model produces biomass a function of the biomass of another. During the initial stages of the (unforced) analysis, while the biotic population are oscillatory, the multivariate effects are most pronounced: the partial variance sums (PVS) increase from an initial value of 0.80-0.92 after 336 h when they stabilize. The initial sensitivity to the bacterial and microtrophic rates also results from excess organic matter 'introduced' into the system. The sums of the partial variances for the third state variable (Table 6) exhibit a similar increase during the first 14 days of the simulation (0.76-0.85). The PVS of all three analyses fluctuate within +0.03 after the system stabilizes. All three state variables for the forced system have partial variance sums which drop below the 0.65 threshold. Total biotic production and the biotic fraction of carbon are consistently lower, with the PVS of the former fluctuating by as much as 0.35 and the latter stabilizing at 0.44_+ 0.03. Diurnal increases multivariate effects. Because of the overall reduction in photosynthesis, the combined uncertainties of the other biotic process have a 0 O.

S L A T E VARIABLE 3

PAR

Z

>4 o d

Ci

< c,j 0 0 C)

1 4 0

80

160

240

320

400

480

560 TIME

640

720

800

880

960

1040

MAX UPTAKE RATE FOR PROT 2 MAX ALGAL GROWTH RATE + MRX SHRIMP ORRZINO RATE × 8RCTERIAL GROWTH EFFIC]ENCT

Fig. 14. Time-varyingpartial variancesof the major parameters affectingthe DIN concentration (small parameter variations).

42

STATE VARIABLE 3

0 O ,=

j

<

P A R.

Ci

go

0

80 0 A + X

160

240

1320

400

480

¢'J6O TIME

640

720

800

8~30

96Q

1040

MAX UPTRKE RATE FOR PBOT 2 MRX RLGRL GROHTH RATE MRX 5HR]MP GRAZING RRTE BACTERIAL OROHTH EFFICIENCT

Fig. 15. Time-varying partial variances of the major parameters affecting the DIN concentration (large parameter variations).

TABLE 6 Time-varying statistical reponse of the DIN concentration, with and without forcing, to global variations in the values of four key parameters: A . . . . S6Rz, VMXp2, and YB Dissolved inorganic nitrogen Time

With diurnal light cycles

With constant light

(h)

Mean value

Variance (×10 -s)

Coeff. of var.

Partial var. sum

Mean value

Variance ( × 1 0 -3 )

Coeff. of var.

Partial var. sum

96.0 144.0 192.0 240.0 288.0 336.0 384.0 432.0 480.0 528.0 576.0 624.0 672.0 720.0 768.0 816.0 864.0 912.0 960.0

47.32 61.05 69.10 70.59 74.58 84.11 87.80 85.81 83.62 80.70 80.08 80.77 81,68 79.72 79.13 80.20 81.31 80.44 79.41

1.06 1.27 1.35 1.42 1.34 1.36 1.61 1.54 1.41 1.54 1.67 1.60 1.60 1.69 1.72 1.88 2.07 2.04 2.07

0.687 0.584 0.532 0.533 0.490 0.439 0.457 0.457 0.448 0.486 0.510 0.495 0.490 0.516 0.524 0.541 0.559 0.561 0.572

0.83 0.87 0.85 0.87 0.91 0.69 0.55 0.56 0.64 0.74 0.73 0.73 0.79 0.78 0.76 0.70 0.62 0.61 0.61

40.39 50.99 59.69 59.05 62.99 72.05 71.78 70.89 70.08 70.03 71.25 71.66 70.70 70.38 69.40 69.39 69.43 68.65 67.31

1.38 1.89 2.25 2.34 2.30 2.76 2.97 3.00 3.04 3.11 3.09 3.07 3.0? 3.13 3.00 2.92 2.85 2.91 2.96

0.919 0.852 0.821 0.820 0.761 0.729 0.759 0.773 0.787 0.796 0.781 0.773 0.784 0.795 0.789 0.778 0.769 0.786 0.808

0.76 0.82 0.81 0.81 0.87 0.82 0.87 0.85 0,84 0.85 0.86 0.85 0.85 0.86 0.87 0.85 0.86 0.86 0.83

43 greater effect on the balance between production and respiration. Even when the parameters are varied by only 5%, the PVS of the two analyses are as low as those from global variations, implying that coupling rather than non-linearity is responsible. Cukier and co-workers (1978) mention that in general for chemical systems, as the product concentrations increase with time, so do the effects of coupling. State variables 1 and 2 would be analogous to the reaction rate and total production of such a system and their partial variance sums, in fact, do decrease. This is not the case with DIN since the process flow rates of dissolved inorganic nitrogen tend to vary inversely with NH~- concentration. Furthermore, the third output variable provdes the clearest picture of the situation because LT(t) only indirectly influences its behavior. Although the local and global sensitivity of the model to the algal, bacterial, and protozoan parameters are relatively similar (Fig. 14 and 15), the PVS decreases from 0.97 to about 0.6 (Table 6). The partial variance curves of the bactivore and bacteria become more similar. This behavior, like that of the algal uptake and growth rates mentioned previously, suggests that these two parameters are also coupled. Because the algal partial variance also drops and because the algae depends on N H ~ , it is also possible that Am~Xcontributes to the low PVS. DISCUSSION The importance of the balance between production and respiration to the performance of the model has been demonstrated by the FAST sensitivity/ uncertainty analysis. Parametrically induced changes in this ratio effect the inorganic nutrient cycle more than the bio-carbon dynamics in terms of either its relative abundance of its rate of production. This is mainly due to limitations placed on photosynthesis by NH~ and CO 2 availability; the algae is the predominant biotic variable and the original source of all organic matter produced in the ecosystem. As a result, the nutrients initially introduced into the system determine its total biomass to a large extent. The endogenous control of the model works because there is negative feedback between predators and their prey, e.g. a decrease in prey availability reduced the predation rate. The prey populations, at lowest trophic levels of both foodwebs, are also dependent on substrate concentrations. In a stable system the bacteria and algae do not utilize the same nutrient sources, and thus they avoid potential extinction due to the combined effects of predation and competition. Feedback control is responsible for the stability of the simulation model. Only two nutrient cycles are used to control on the performance of the model. The feedback loops, incorporated in the carbon and nitrogen cycles,

44

are able to damp out oscillations and stabilize the system. Compared to the model of a tundra pond developed by Tiwari and co-workers (1978), in which a + 20% change in the maximum photosynthetic rate constant causes the algal biomass rapidly to become almost zero, the present closed-ecosystem model displays a remarkable degree of stability and parametric insensitivity. The endogenous control system built into the model demonstrates that even an extremely simplified approach can behave similarly to a homeostatic ecosystem and show the complementary nature of chemotaxis and homeostasis in ecological stability. The endogenous control which stabilized the ecosystem model is similar to the homeostatic control responsible for the stability of the anchialine microcosm. Both systems exhibit a critical dependence on the number of shrimp and are destabilized when this number is decreased by one. Other conditions which disrupt the stability of the model are qualitatively similar to those which would disrupt a homeostatic ecosystem; e.g. the observations made by Pearson (1982) at Loch Eil and by Kerr and coworkers (1972) at Shriner's Pond. Qualitatively, the simulation runs have indicated the following: The disruption of homeostatic and chemotaxic control and resultant inability to dampen out the oscillatory behavior caused by an insufficient population of shrimp (heterotrophs). - - How omnivorous grazing serves to stabilize the community by dampening any predator prey oscillations and prevents the extinction of primary producers by controlling the populations of the secondary consumers. The potential importance of microtrophs to nutrient regeneration in a closed ecosystem. - - How diurnal light cycles and intercellular nutrient storage may help phytoplankton to succeed in unstable nutrient environments and successfully compete with microbes for inorganic nutrients. - - A possible mechanism for the failure of an unbalanced microcosm. -

-

-

-

The model is of course a highly simplified imitation of the actual system; a static compartment represents a community whose biotic structure shifts constantly and kinetic expressions approximate the gross activity of the community. The parameters for these expressions must be adjusted to provide an accurate quantitative representation of the microcosm. ACKNOWLEDGEMENTS

The assistance of Dr. Greg McRae, Carnegie-Mellon University, Department of Chemical Engineering, who also supplied the FAST application program used to analyze the model, is greatly appreciated; the contribution by Mr. Joe Hanson and Dr. David Lawson of JPL in the form of helpful

45 d i s c u s s i o n s a n d s u g g e s t i o n s is also a c k n o w l e d g e d w i t h sincere a p p r e c i a t i o n . This work was supported by NASA, Ames Research Center, Agreement No. NCC2-122. REFERENCES Baath, E., Lohm, U., Lungren, B., Rosswall, T., S6nderstr6m, B. and Sohlenins, B., 1981. Impact of microbial-feeding animals on total soil activity and nitrogen dynamics: a soil microcosm experiment. Oikos, 37: 257-264. Bader, F.G., Fredrickson, A.G. and Tsuchina, H.M., 1976. Dynamics of an algal-protozoan grzing interaction. In: R.P. Canale (Editor), Modeling Biochemical Processes in Aquatic Ecosystems. Ann Arbor Science, Ann Arbor, MI, pp. 257-250. Banse, K., 1982. Cell volumes, maximal growth rates of unicellular algae and ciliates, and the role of ciliates in the marine pelagial. Limnol. Oceanogr., 27: 1059-1071. Barsdate, R.J., Prendtki, R.T. and Fenchel T., 1974. Phosphorus cycle of model eco-systems: significance for decomposer food chains and effect of bacterial grazers. Oikos, 25: 234-251. Bick, H., 1973. Population dynamics of Protozoa associated with the decay of organic materials in fresh water. Am. Zool., 13: 149-160. Canale, R.P., DePalma, L.M. and Vogel, A.H., 1976. A plankton based food web model for Lake Michigan. In: R.P. Canale (Editor), Modeling Biochemical Processes in Aquatic Ecosystems. Ann Arbor Science, Ann Arbor, MI, pp. 33-74. Conover, R.J., 1978. Transformation of organic matter. In: O. Kinne (Editor), Marine Ecology. Vol. IV, Dynamics. Wiley, Chichester, pp. 221-499. Cooper, D.C., 1973. Enhancement of net primary production by herbivore grazing in aquatic laboratory microcosms. Limnol. Oceanogr., 18: 31-37. Corner, E.D.S. and Davis, A.G., 1971. Plankton as a factor in the nitrogen and phosphorus cycles in the sea. Adv. Mar. Biol., 9: 101-204. Crisman, T.L., Beaver, J.R. and Bays, J.S., 1981. Examination of the relative impact of microzooplankton and macrozooplankton on bacteria in Florida lakes. Int. Ver. Theor. Angew. Limnol. Verh., 21: 359-362. Cukier, R.I., Levine, H.B. and Shuler, K.E., 1978. Nonlinear sensitivity of multiparameter model systems. J. Comp. Phys. 6: 1-42. Danso, S.K.A. and Alexander, M., 1975. Regulation of predation of prey density: the P r o t o z o a n - R h i z o b i u m relationship. Appl. Microbiol., 29: 515-521. DiToro, D.M., 1980. Application of cellular equilibrium and Monod theory to phytoplankton growth kinetics. Ecol. Modelling, 8: 201-218. Droop, M.R., 1973. Nutrient limitation on asmotrophic protista. Am. Zool., 13:209 214. Fenchel, T. and Blackburn, T.H., 1979. Bacteria and Mineral Recycling. Academic Press, London, 225 pp. Fukami, K., Simidu, U., and Taga, N., 1981. Fluctuations of the communities of heterotrophic bateria during the decomposition process of phytoplankton J. Mar. Res., 55: 171-184. Griffith, D., 1975. Prey availability and the food of predators. Ecology, 56: 1209-1214. Hamilton, R.D. and Preslan, J.E., 1970. Observations on the continuous culture of a planktonic phagotrophic protozoa. J. Exp. Mar. Biol. Ecol., 5: 94-104. Hargrave, B.T., 1970. The effects of a deposit-feeding amphipod on the metabolism of benthhic microflora. Limnol. Oceanogr., 15: 21-30.

46 Harte, J., Levy, D., Rees, J. and Saegebarth, E., 1980. Making microcosms an effective assessment tool. In: J.P. Giesy, Jr. (Editor), Microcosms in Ecological Research. Proc. Symp., 8-10 November 1979, Augusta, GA. CONF-781101, Technical Information Center, U.S. Department of Energy, Washington DC, pp. 105-137. Hindemarsh, A.C. and Byrne, G.D., 1975. EPISODE: An experimental package for the integration of ordinary differential equations. Rep. UCID-30112, Lawrence Livermore Laboratory. Holthius, L.B., 1963. On red coloured shrimps (Decapoda, Caridea) from tropical land-locked saltwater pools. Zool. Meded. Leiden, 38: 261-279. Jannasch, H.W., 1967. Growth of marine bacteria at limiting concentrations of organic carbon in seawater. Limnol. Oceanogr., 12:264-271 Jannasch, H.W., 1968. Growth characteristics of heterotrophic bacteria in seawater. J. Bacteriol., 95: 722-723. Johannes, R.E. and Satomi, M., 1966. Composition and nutritive value of fecal pellets of a marine crustacean. Limnol. Oceanogr., 11:191-197. Jorgensen, S.E., 1976. A eutrophication model for a lake. Ecol. Modelling, 2: 147-165. Kerr, P.C., Brockway, D.L., Paris, D.F.~and Barnett, J.T., 1972. The interrelation of carbon and phosphorus in regulating heterotrophic and autotrophic populations in an aquatic ecosystem. In: G.E. Likens (Editor), Proc. Symp. Nutrients and Eutrophication, the Limiting-Nutrient Controversy. Spec. Symp. 1, American Society of Limnology and Oceanography, 11-12 February 1971, Michigan State University. Allen Press, Lawrence, KS, pp. 41-62. Koda, M., McRae, G.J. and Seinfeld, J.H., 1979. Automatic sensitivity analysis of kinetic mechanisms. Int. J. Chem. Kinetics, 11: 427-444. Legner, M., 1973. Experimental appraoch to the role of protozoa in aquatic ecosystems. Am. Zool., 13: 177-192. Lehman, J.T., Botkin, D.B. and Likens, G.E., 1975. The assumptions and rationals of a computer model of phytoplankton population dynamics. Limnol. Oceanogr., 20: 343-364. Liepmann, D., 1983. A closed ecosystem model: development and global sensitivity analysis via the Fourier Amplitude Sensitivity Test. M.S. Thesis, Department of Chemical Engineering, California Institute of Technology, Pasadema, CA, 84 pp. Margalef, R., 1978. General concepts of population dynamics and food links. In: O. Kinne (Editor), Marine Ecology. Vol. IV, Dynamics. Wiley, Chichester, pp. 617-704. Marshall, S.M., 1973. Respiration and feeding in copepods. Adv. Mar. Biol., 11: 57-120. Mayzaud, P., 1973. Respiration and nitrogen excretion of zooplankton. If. Studies of the metabolic characteristics of starved animals. Marine Biol., 21: 19-28. McRae G.J,, Tilden, J.W. and Seinfeld, J.H., 1982. Global sensitivity analysis - - a computational implementation of the Fourier amplitude sensitivity test (FAST) Computers Chem. Eng., 6: 15-25. Paerl, H.W., 1974. Bacterial uptake of-dissolved organic matter in relation to detrital aggregation in marine and freshwater systems. Limnol. Oceanogr., 19: 966-972. Parsons T.R., Takahashi, M. and Hargrave, B., 1977. Biological Oceanographic Processes (2rid Edition). Pergamon Press, Oxford, 332 pp. Pearson, T.H., 1982. The Loch Eil Project: assessment and synthesis with a discussion of certain biological questions arising from a study of organic pollution of sediments. J. Exp. Mar. Biol. Ecol., 57: 93-124. Platt, T., Denman, K.L. and Jassby, A.D., 1977. Modeling the productivity of phytoplankton. In" E.D. Goldberg (Editor), The sea. Ideas and Observations on Progress in the Study of the Seas, 6. Wiley-Interscience, New York, NY, pp. 807-856.

47 Porter, K.G., Pace, M.L. and Battey, J.F., 1979. Ciliate Protozoans as links in freshwater planktonic food chains, nature, 277:563 563. Robinson, J.D., Mann, H.H., and Novitsky, J.A.~ 1982. Conversion of the particulate fraction of seaweed ditritus to bacterial biomass. Limnol. Oceanogr., 27: 1072-1079. Sherr, B.F., Sherr, E.B. and Berman, T., 1982. Decomposition of organic detritus: a selective role for microflagellate protozoa. Limnol. Oceanogr., 27:765 769. Sjoberg, S., 1977. Are pelagic systems inherently unstable? A model study. Ecol. Modelling, 3: 17 37. Sorokin, Yu. I., 1978. Decomposition of organic matter and nutrient regeneration. In: O. Kinne (Editor) Marine Ecology. Vol. IV, Dynamics. J. Wiley, Chichester, pp. 501 616. Steele, J.H. and Mullin, M.M., 1977. Zooplankton dynamics. In: E.D. Goldberg (Editor), The Sea. Ideas and Observations on Progress in the Study of the Seas, 6. Wiley-Interscience, New York, NY, pp. 857-890. Stout, J.D., 1973. The relationship between protozoan populations and biological activity in soils. Am. Zool. 13: 193-201. Stumm, W. and Morgan, J.J., 1970. Aquatic Chemistry, Wiley-Interscience, New York, NY, 583 pp. Tiwari, J.L., Hobbie, J.E., Reed, J.P., Stanley, D.W. and Miller, M.C., 1978. Some stochastic differential equation models of an aquatic ecosystem. Ecol. Modelling, 4: 3-28. Vaccaro, R.F. and Jannasch, H.W., 1967. Variations in uptake kinetics for glucose by natural populations in seawater. Limnol. Oceanogr., 12: 540-542. Williams, P.J. LeB., 1970. Heterotrophic utilization of dissolved organic compounds in the sea. I. Size distribution of populations and relationship between respiration and incorporation of growth substrates. J. Mar. Biol. Assoc. U.K., 50: 859-870. Wong, D.C.L., 1975. Algae of the Anchialine pools at Cape Kinau, Maui, and aspects of the trophic ecology of halocaridina rubra Holthius (Decapoda, Atyidae). Masters Thesis. Botanical Sciences, University of Hawaii, HI, 63 pp. Wright, R.T., 1974. Minerization of organic solutes by heterotrophic bacteria. In: R . R Colwell and R.Y. Morita (Editor), Effect of the Ocean Environment on Microbial Activities. Proc. 2nd United States Japan Conference on Marine Microbiology, 25-30 August 1972, University of Maryland. University Park Press, Baltimore, MD, pp. 546-565. Zobell, C.E., 1943. The effect of solid surfaces on bacterial activity. J. Bacteriol., 4 6 : 3 8 - 5 9