CBP/TRS 127/94
December 1994
903R94028 -
Ecosystem Process Modeling of
Submerged Aquatic Vegetation in the
Lower Chesapeake Bay
U.S. EnttkonnwtfaJ Preteciten'
R«glsa ill Inhifiiatiofl Resource
Canter (CPM52)
m Chaslnui Strest
Pr^-Jfiijjhia, PA 13107 /
TD
225
.C54
S915
Chesapeake Bay Program
Printed on
Recycled Paper
-------
Ecosystem Process Modeling of
Submerged Aquatic Vegetation in the
Lower Chesapeake Bay
by
Richard L. Wetzel
and
Mark B. Meyers
School of Marine Science/Virginia Institute of Marine Science
College of William and Mary
Gloucester Point, Virginia 23062
U.S. EfltfiroTKJMmtaJ Protecfen Agency
Region III Infonnalioa Resourca
Center (-PM52) ,
841 Ckssinut Street ",}r ... v,
PA 19183; .-$l
Printed by the U.S. Environmental Protection Agency for the Chesapeake Bay Program
-------
-------
EXECUTIVE SUMMARY
The modeling studies reported here emphasize the polyhaline eelgrass habitats of the lower
Chesapeake Bay, as a part of a larger program of ecosystem modeling sponsored by .the Living
Resources Subcommittee of the Chesapeake Bay Program. This report covers progress made
through April 1993. The primary focus of the polyhaline SAV modeling studies has been on light-
dependent eelgrass productivity and water quality parameters which affect the submarine light
climate. Where eutrophication impacts have been assessed, they have been used to drive increases
in epiphytes, which block light penetration to eelgrass leaves. Eutrophication of the tidal tributaries
diminishes light-dependent eelgrass production and growth and has been implicated as the principal
cause for SAV declines throughout the Chesapeake Bay. The implicit assumptions have been that
eelgrass is not nutrient-limited in its natural setting and that increased nutrient loadings to the
estuary negatively affect SAV growth and survival by promoting the growth of planktonic
microalgae and epiphytes.
The successes of eelgrass modeling to date are illustrated by three major findings. First, the
current model version yields productivities and plant survival characteristics, both short and long
term, which are very similar to those observed in process-oriented studies and in analyses of long
term SAV decline within the Chesapeake Bay. The modeling results suggest that the principal
environmental variable controlling eelgrass growth and survival in the lower Chesapeake Bay is the
availability of light, specifically photosynthetically active radiation. For the water column,
suspended inorganic solids and chlorophyll concentration determine the intensity of light reaching
the plant canopy but not necessarily the amount available for plant photosynthesis which further
depends on epiphytic loads.
Secondly, these modeling studies have demonstrated the importance of epiphytic grazers
in controlling epiphyte density on eelgrass leaves. Epiphytes reduce plant photosynthesis by
attenuating light and interfering with gas (CO2 and O2) exchange at the leaf surface. The absence
of or the reduction in the population density of grazers, as shown in modeling sensitivity studies and
confirmed in mesocosm studies, can hasten the decline of eelgrass under nutrient loadings typical
of present conditions in areas of the lower Chesapeake Bay.
The third finding, particularly suited to elucidation by model simulation analysis, is that
variability in the submarine light climate, determined by variable solar irradiance and attenuation
coefficients, at daily frequencies can have critical consequences with regard to SAV survival. Model
scenarios where light parameters are specified as smoothly varying functions of monthly averages
produced results that differ significantly from model scenarios where the light parameters were
treated as stochastic variables derived from in situ data; the latter produced results that better
reproduced observed conditions. This finding has significance for setting target water quality and
habitat restoration criteria, and for deciding the appropriate means of simulating environmental
.conditions experienced by SAV's in continuing modeling efforts.
Future SAV modeling and simulation efforts will expand along three principal lines. First,
the formulation of the light portion of the model will be revised to reflect a more theoretical basis
rather than the site specific, data-driven empirical formulations used in the current versions. This
will allow better integration with other efforts including the revised 3D hydrodynamic-water quality
-------
-------
model as well as better reflect SAV habitat criteria evaluation of restoration goals. Second, the SAV
model will be expanded to include nutrient cycles, sediment oxygen and nutrient fluxes, and water
column trophic dynamics. This will form the basis of a generic SAV-Littoral Zone model that can
be parameterized for different regions of the Chesapeake Bay and its tributaries. Third, the SAV-
Littoral Zone model will be coupled with the water quality and hydrodynamic modeling efforts in
tributaries of the lower Bay. This will improve the ability to predict the transport of dissolved and
paniculate materials important for water quality and better determine the role of littoral zone
physical and biological processes in the maintenance of water quality.
in
-------
-------
ACKNOWLEDGEMENTS
A modeling project of this nature depends on the input from a variety of sources and
more often than not the sharing of data often before the originator of those data are able to
publish. We particularly wish to thank Robert J. Orth, Kenneth A. Moore, and Hilary A.
Neckles for their generous aid in this regard. To our modeling colleagues in Maryland who
share in this continuing program, Michael Kemp, Chris Madden and Walter Boynton, thanks for
the many ways you have contributed to this effort. We also thank our graduate students who
have been directly involved in these studies, William Seufzer and David Fugate.
Funding for these studies have been provided by several sources. We acknowledge the
funding support provided by the U.S. E.P.A. Chesapeake Bay Program (Contract nos.
CB003909-01 and CB003909-02) and the specific assistance of Carin Bisland and Richard
Batiuk of the Chesapeake Bay Program Office for program management support and technical
assistance. This report has been greatly enhanced by their input and patience. We also
acknowledge the additional funding support provided for this ecosystem modeling effort by the
Chesapeake Corporation and the Commonwealth of Virginia.
Finally, without the support of the Living Resources Subcommittee of the Chesapeake
Bay Program none of this, would have been possible.
This report is also available as Special Report in Applied Marine Science and Ocean
Engineering (SRAMSOE) No. 323, Virginia Institute of Marine Science, Gloucester Point,
Virginia 23062.
IV
-------
-------
TABLE OF CONTENTS
EXECUTIVE SUMMARY ii
ACKNOWLEDGEMENTS iv
INTRODUCTION 1
BACKGROUND '. 2
Ecological Modeling of Complex Systems ; 2
The York River Regional Ecosystem 3
Physical-Chemical Characteristics 3
Biological-Ecological Characteristics 4
Integrative Modeling and Synthesis Program 5
The Approach 5
Ecosystem Process Models 8
CHESAPEAKE BAY SUBMERGED AQUATIC VEGETATION 10
SAV Modeling Goals and Objectives 10
SAV Management Needs 11
SAV Conceptual Models 11
SAV Simulation Models 14
Eelgrass Simulation Model 14
Environmental Variables and Forcing Functions 16
Biological Processes, Interactions and Controls 18
Eelgrass Photosynthesis and Growth 18
Epiphytes 20
Grazers 20
Numerical Computation and Simulation Analysis 21
SIMULATION RESULTS AND DISCUSSION 22
Nominal Simulation Analyses (Version 2: non-variable inputs) 22
Model Responses to Physical Regime and Epiphyte-Grazer Interactions 27
Physical Regimes 27
Epiphyte-Grazer Interactions 27
Effect of Submarine Light Variability on Eelgrass Dynamics (Model Version 4) 32
Solar Irradiance 34
PAR Attenuation 38
Key Simulation Findings 41
FUTURE WORK 44
REFERENCES 47
APPENDIX 51
-------
-------
INTRODUCTION
For the past several decades and in particular since the mid-1970s, the Chesapeake Bay and
its tributaries have been the focus of intense study and directed research to understand the causes
and consequences of reduced water quality and the loss of living resources. A long-term goal of
these varied research efforts has been to provide the causes of and propose the means to reverse
trends in environmental degradation. Research and management programs supported by regional,
state, and federal initiatives have included efforts ranging from the large-scale, three-dimensional
hydrodynamic and water quality modeling to studies focused on specific habitats (e.g., shallow
waters dominated by submerged aquatic vegetation) and populations of economic and/or
recreational importance (e.g., oysters, blue crabs and striped bass). Water quality issues have been
one general focus, examining the phenomena of low dissolved oxygen, nutrient loadings, and upland
land-use practices and management. Living resource issues are another area of research, addressing
particular habitats and specific populations, often out of the context of the larger system of which
they are a part. This leaves a complex and yet unresolved question:
How and to what extent can the results of research on water quality and living
resources, often viewed independently from one another, be integrated to arrive at
sound ecological management strategies of the estuary, its tributaries, and multiple
resources?
To address this question in part, Wetzel and Hopkinson (1991) suggested that ecosystem
modeling and simulation analysis was an appropriate and useful tool for the Chesapeake Bay
Program. Ecosystem modeling and simulation analysis has progressed rapidly over the past two
decades. Advances in computer technology have played a leading role in this development and at
the same time the modeler and model end-user have matured with regard to methods and
expectations. While modeling techniques and the use of mathematical models are not new to marine
science or environmental management, no concerted effort has been given to integrating the various
methods used in the physical, chemical, and biological sciences for studies of the Chesapeake Bay
and it tributaries. An integration that is considered by many as both necessary and appropriate to
deal with the high degree of complexity of natural systems and the environmental management of
large scale systems.
The Living Resources subcommittee of theU.S. E.P.A.'s Chesapeake Bay Program initiated
in 1991 an ecosystem modeling program to develop an integrative modeling framework responsive
to management needs particularly as they related to living resources of the Chesapeake Bay. The
program is a cooperative, integrated modeling effort among the College of William & Mary's
Virginia institute of Marine Science and the University of Maryland laboratories at Solomons (CBL)
and Horn Point (HPEL). The program involves three general modeling approaches: 1) Ecological
Regression Models, which employ statistical models to establish correlative measures between
variables such as nutrient and suspended loads and readily observed environmental variables such
as phytoplankton biomass (chlorophyll a) and dissolved oxygen concentration, 2) Ecosystem Process
Models, which develop mathematical descriptions of ecological processes and simulate the
dynamics in time and space of populations, habitats or entire ecosystems, and 3) Fish Bioenergetics
Models, which relate finfish physiological and behavioral responses to food abundance, temperature,
and dissolved oxygen. Conceptually, these approaches can be linked through the sharing of common
-------
variables and forcing functions. Analyzed as standalone models, they provide insight and some
degree of predictive capability in establishing management criteria and evaluating alternative
management scenarios. Together, they provide a powerful tool for analysis of living marine
resources in the context of environmental variability and quality.
Within the programmatic goals of ecosystem modeling is the conceptualization,
implementation and analysis of three ecosystem process models for spatially and ecologically
distinct regions or habitats of the estuary; 1) submerged aquatic vegetation, 2) emergent intertidal
marshes and 3) the water column-benthos. While these are somewhat arbitrary distinctions, some
environmental management concerns are unique to these areas (e.g. the development of habitat
criteria and restoration goals) but as importantly they represent a logical ecological division for
model development and validation. Also within the programmatic framework, is the coupling of
these models via common forcing functions, environmental water variables and shared state
variables with both ecological regression models and fish bioenergetics models.
This report presents results from ecosystem process models of the polyhaline littoral zone
habitat dominated by eelgrass (Zostem marina) in the lower portions of the Chesapeake Bay. The
emphasis of these modeling studies has been on the role of environmental factors, particularly light,
in governing the distribution, growth and survival of eelgrass.
BACKGROUND
Ecological Modeling of Complex Systems
Ecosystem models are abstractions or simplifications of complex natural systems and as such
are useful tools for scientific analysis; the multitude of interactions is too overwhelming to be
perceived intuitively, but may be penetrable within a quantitative simulation framework.
Deterministic, numerical simulation models of ecosystems synthesize large amounts of information
on individual parts of systems, integrate these data based on a conceptual structure and explore, via
simulation analysis, compartmental behavior under various simulated operating conditions. Models
are now commonly used to plan and guide research, to identify data weaknesses and gaps, to
evaluate management-oriented alternatives, and to provide the basis to formulate hypotheses
regarding a system's structure and function. Wetzel and Hopkinson (1991) have recently reviewed
ecosystem models of coastal systems with reference specifically to their application to the
Chesapeake Bay region.
Models can be classified or characterized according tcTa number of schemes. Two general
categories into which ecological models can be placed are 1) Explanatory Models and 2) Correlative
Models (Gold 1977; Gilchrist 1984). Explanatory models attempt to examine causal relations
between components of the modeled system based on a conceptual structure or hypothesis and
explain cause-effect relationships. Correlative models, on the other hand, statistically describe the
relationships between a dependent variable and one or more independent variables using techniques
of regression and correlation analysis. These models have no explanatory power and cannot
extrapolate relationships beyond the conditions under which they were developed. They can,
however, be used to suggest probable causal relationships that can be explored using other methods.
-------
Often, aquatic ecosystem models combine features of both model types, incorporating statistical
inferences where theoretical constructs are poorly developed or process-oriented data are lacking.
At the Virginia Institute of Marine Science (VIMS), modeling has been used in conjunction
with interdisciplinary research conducted within the Chesapeake Bay and its tributaries. Field and
laboratory studies provide necessary information on processes and distributions. Modeling provides
an ecosystem context into which these results can be placed and allow for experimentation regarding
the sensitivity of an ecosystem to chronic or traumatic perturbations. Poorly understood processes
or quantities uncovered by modeling then feed back into a renewed series of observation and
experimentation. Models then prove very useful for examining environmental, biological, and
ecological relationships and for indicating productive avenues of research (e.g., Wetzel and Neckles
1986). A case study of this process is our research efforts in a tributary of the Chesapeake Bay
which has experienced comparatively low levels of anthropogenic impact, the York River.
The Yoik River Regional Ecosystem
In 1991, the Institute (VIMS) and the School of Marine Science, College of William and
Mary embarked upon a multidisciplinary program designed to address large scale, complex coastal
ecosystems and contemporary management issues. The intent of the program is to provide a
professional infrastructure for the integration of knowledge from various marine science disciplines
relevant to coastal zone management and the acquisition of basic knowledge necessary for
understanding coastal ecosystem dynamics within a large-scale framework. The focus initially is on
the York River estuary and its watershed: the York River Regional Ecosystem. The primary
organizational and systems analysis tools employed will be the development of ecosystem models
coupled with water quality, hydrodynamic and fisheries models. The undertaking is viewed as a
long-term effort and commitment of institutional resources augmented by grants and contracts from
state and federal agencies, local industry, and the private sector. A brief description of the York
River system follows.
Physical-Chemical Characteristics
The York River is a subestuary of the Chesapeake Bay and is formed by the confluence of
the Mattaponi and Pamunkey Rivers at West Point, Virginia. Its drainage basin is 69,000 km2, of
which approximately 70% is forested, 22% is in crop land and pasture, and <2% is classified as
urban (Bender 1987). Gloucester County which makes up much of the north shore of the York is
one of the most rapidly developing Tidewater counties, experiencing a population increase of 49.3%
for the 1980-90 census period. The south shore is somewhat protected for the present from large
scale growth and development due to the large military land holdings, the federal Colonial Historic
National Park and the York River State park.
Total average freshwater inflow to the river is estimated at 70 m3 sec"1. Average annual
rainfall over the watershed is 112 cm which peaks in August but runoff is low due to transpiration
and evaporation. The climate is humid temperate with an annual mean temperature of 15 C.
The salinity structure of the York River is influenced mainly by the interaction of freshwater,
salt water and tidal energy and to a lesser extent by other forcings (e.g., wind). Periods of
stratification/destratification cycle in the York River as a result of energy associated with the spring-
-------
neap tidal cycles. Salinity gradients between the surface and bottom waters tend to be stronger
during periods of neap tides and to disappear during spring tide periods. During low flow periods,
salt water intrudes 20 to 30 km upriver from West Point. At the mouth, periods of bottom water
hypoxia or anoxia generally coincide with stratification events during the summer. The extent to
which hypoxia or anoxia in bottom waters of the York River are changing in either duration or
spatial coverage is not known. The data for examining this question is not extensive enough nor for
sampling reasons designed to address it in terms of an historical data base.
Nutrient concentrations and distributions in the York River show longitudinal and vertical
patterns typical of temperate, coastal plain estuaries. Generally, there is a longitudinal gradient in
nutrient concentrations that fluctuates in magnitude on a seasonal basis. Concentrations increase
upriver reflecting watershed sources and biological processes that modulate in situ concentrations.
Sources of nutrients to the York River estuary include inputs at the fall-line, point sources at West
Point at the head of the estuary, surface runoff, ground-water inputs, and precipitation.
Overall, nutrient and chlorophyll concentrations in the York River are not indicative of a
highly eutrophic system while other indicators suggest some degree of enrichment (e.g., loss of SAV
from historically vegetated areas). Based on data available through the late 1970's, Heinle et al.
(1980) characterized the York as a 'moderately eutrophic' estuary. Because projected growth in the
area (primarily along the northern shoreline and upper watershed) is high for the next decade or
longer and there will be concomitant changes in land use (conversion of forested and agricultural
lands to rural housing developments and perhaps industry), the likelihood that anthropogenic
impacts will increase without proper management is great.
Biological-Ecological Characteristics
Biologically, the York River Regional Ecosystem supports a vast array of species
populations and community types that range from tidal freshwater to estuarine-marine dependent.
Primary production in the York River Regional Ecosystem is determined by several sources.
Emergent wetlands, dominated by various freshwater and salt-tolerant species, are an important
component of the estuary and form extensive marshes in the Pamunkey and Mattaponi Rivers,
smaller tributaries and creeks, and at the mouth of the York River. Submerged aquatic vegetation,
once a dominant autotrophic component in the lower reaches of the York, are now confined to the
lower eight to ten kilometers and are dominated by eelgrass (Zostera marina). Non-vegetated
sediments in the shoal areas (<2 m deep) make up a significant benthic habitat (ca. 38% of the York
River surface area) and contribute to primary production by supporting an active autotrophic
microflora dominated by benthic diatoms. Phytoplankton, because of the large euphotic water mass,
contribute the greatest proportion of primary production to the York River Regional Ecosystem
system and are dominated by larger forms (i.e., diatoms) during the spring bloom (January to May)
and by smaller forms (nanoplankton and picoplankton) during the summer. Primary production
from these sources and allochthonous organic inputs from the watershed support a complex trophic
network that includes organisms ranging in size from bacteria to large finfish.
-------
Integrative Modeling and Synthesis Program
The Approach
The York River Regional Ecosystem is conceptualized as being composed of three,
interacting large-scale components: uplands, wetlands, and aquatic systems (Figure 1). The
interaction between these components is governed primarily by larger-scale hydrologic,
meteorologic and anthropogenic natural forcings and human perturbations. Within each of these
components are smaller-scale units defined or identified by their ecological/biological structure and
organized by the flow of energy, the cycling of essential elements, and the controls imposed by
physical, chemical and biological interactions including those imposed by human activity. It is at
this fundamental level of ecological organization that our modeling efforts are organized. The
development of a single, large-scale, fully-integrated model of the York River Regional Ecosystem,
though a goal to work toward, is beyond the scope of the current efforts. The development and
implementation of the models proposed here should be viewed as the first necessary steps toward
building integrated, large-scale ecosystem models that truly integrate a systems physics, chemistry
and biology over both space and time.
The focus of these first efforts has been on the development, validation and simulation
analysis of ecosystem process models for specific components of the York River Regional
Ecosystem. We have proposed to develop over the first several years conceptual and simulation
models of four principal components of the York River Regional Ecosystem: 1) Submerged Aquatic
Vegetation (SAV), 2) Emergent Intertidal Marshes), 3) Water Column-Benthos, and, 4) Fisheries.
Figure 2 illustrates the general connectivity of the proposed component models and the dominant
factors influencing or controlling their interaction.
The models will be time-dependent and spatially-averaged for five geographic regions of the
York River Regional Ecosystem that are defined by salinity regime for the wetlands and aquatic
components. These are 1) uplands, 2) tidal freshwater (<0.5 ppt salinity), 3) oligohaline (0.5-5 ppt),
4) mesohaline (5-18 ppt), and 5) polyhaline (>18 ppt) or lower estuary/river mouth regions. Each
of these geographic regions can be treated as a hydrodynamic unit characterized by one or more of
the component models. This segmentation scheme has several advantages in that it
• corresponds to the principal hydrodynamic and physical-chemical regimes of the York
River Regional Ecosystem,
• can be adapted to water quality and hydrodynamic modeling efforts,
• includes four sites of the National Estuarine Research Reserve System which assures
long term data acquisition, and,
• follows in general a segmentation scheme proposed for similar modeling studies of the
Patuxent River watershed in Maryland.
-------
Principal Components of York River Regional
Ecosystem
1
I
1
1 •
I *<
1 *
1 -
1 L
1
%&&& VOR
| i \Jr\
^•^
UPLAND |"fr
FOREST _|
AGRICULTURE]
RURAL |
URBAN |
NDUSTRIAL |
K RIVER WATERSH
*
FD U^^Ss^^
dj ^
» WETLAND JH AQUAT|C |
£
*|FOREST *
* SWAMP ^
UMARSH
> TIDAL FRESH
"* OLIGOLHALINE
^* MESOHALINE
NON-TIDAL ^
1 j^j
TIDAL ^
. £S
^
* TIDAL FRESH ^
OLIGOLHALINE g
-------
Conceptual Ecosystem Process Model Interactions
UJPLANDSy
SUBMERGED
AQUATIC
VEGETATION
WETLANDS
WATER
COLUMN
BENTHOS >
ABUNDANCE
PRODUCflVIT^t
MAN
Figure 2. A conceptual diagram illustrating general relationships and couplings between the four
principal ecosystem process models currently under development for the Ecosystem Process
Modeling Program. The diagram illustrates the importance of land use, physical-hydrodynamic
processes and water quality in governing the behavior of Wetlands, Submerged Aquatic Vegeta-
tion and the Water Column-Benthos components and their influence relative to abundance and
productivity on Fisheries and ultimately human use.
-------
Ecosystem Process Models
Ecosystem process models focus on energy and mass transfers among the components of an
ecological system and those physical, chemical and biological variables considered to control these
transfers. The state variables are the masses, or stocks, of the modeled components and can
represent nutrient stocks, population abundance, or the aggregate density of an entire trophic level.
The control of mass or energy transfer among these model components can be a function of
environmental factors (e.g., solar irradiance, temperature, salinity, water depth) and density-
dependent feedbacks. With respect to biotic state variables, environmental factors often act by
influencing the physiologic rate processes, such as growth, respiration, and enzy matically-controlled
nutrient uptake. Density-dependent controls act by imposing limits on, or altering processes related
to, the density or concentration of organisms or biomass. Examples of such controls are the self-
shading of light, and predation.
Ecosystem process models can be used to assess the impacts of organizational complexity
(the number and kinds of modeled components and processes), the temporal cycles and spatial
variabilities with ecosystems, and the impacts of environmental and anthropogenic perturbations.
This latter capability is most often used to address environmental management concerns. In these
respects, models form the basis of simulation experiments that are impractical or too expensive to
perform in nature. Their veracity, however, is limited to the quality of observation and theory used
in their formulation and to the complexity of the processes that they attempt to simulate.
The construction and analysis of ecosystem process models follows the general scheme
illustrated in Figure 3. The modeling effort begins with a clearly defined set of objectives. Once
these are established, the general steps are (1) construct a conceptual model depicting the general
structure and identifying the principal forcing functions and controls, often depicted with a cartoon
or diagrammatically using Odum's symbolic language (Odum 1971), (2) formulate the equations
that mathematically describe the forcing functions, flux equations and feedbacks, (3) develop the
computer code using a programming language (FORTRAN, Pascal, C, etc.) or simulation tool (e.g.,
High Performance Systems' STELLA ) and verify the program, (4) calibrate and validate the
simulation model using independent data sources, (5) perform systems analysis, and (6) re-design
the model to address the results of model studies, incorporate new data or evolve the model to
address new questions.
Presently, ecosystem process models have been developed for two of the four ecosystem
components identified for the York River Regional Ecosystem: submerged aquatic vegetation (S AV)
and wetland salt marshes characteristic of polyhaline and mesohaline regions of the estuary. For
the wetland ecosystem process models (Spartina altemiflora marshes), a conceptual model for marsh
grass productivity and growth, and distribution relative to specific physical and chemical factors has
been developed. For the S AV model, models have been developed, tested and completed simulation
analyses with two versions.
Here, the results of SAV model studies are presented as they relate to water quality and the
restoration goals set for Chesapeake Bay SAV using other methods (Batiuk et al. 1992; Dennison
etal. 1993).
-------
A Generalized Ecosystem Process Modeling Scheme
1. LIST OBJECTIVES
2. MODEL CONCEPTUALIZATION
- DECIDE COMPARTMENTALIZATION SCHEME
(STATE VARIABLES)
- DETERMINE EXTERNAL FORCING FUNCTIONS
- DEVELOP INTERACTION MATRIX
• 3. DERIVE MATHEMATICAL STRUCTURE(S)
- STATISTICAL - EMPIRICAL RELATIONSHIPS
- THEORECTICAL - A PRIORI CONSTRUCTS
- FEEDBACK - INFORMATION FLOWS
- "FORCED" - DATA LOOK-UP TABLES
4. PROGRAMMING
- DEVELOP CODE
- APPLY MODELING PACKAGES (STELLA, MATLAB, ETC)
- VERIFICATION
5. DIGITAL SIMULATION
- VALIDATION
6. SYSTEMS ANALYSIS
- SENSITIVITY ANALYSIS
- STABILITY
7. RE-DESIGN
Figure 3. The generalized scheme followed in development of ecosystem process models of the
principal components of the Chesapeake Bay and its tributaries. Following this scheme leads to
scenario runs that are designed to address management concerns and provide insight to the best
management options.
-------
CHESAPEAKE BAY SUBMERGED AQUATIC VEGETATION
Submerged aquatic vegetation (SAV) in Chesapeake Bay and its tributaries has been a focal
point of basic research, resource management concerns, and public interest and growing awareness
since the late 1970s. Starting in the early 1970s and continuing into the 1980s, SAV began a
chronic, bay-wide decline in distribution and abundance (Orth and Moore 1983). In a previous,
much-publicized decline during the late 1930s, which was pandemic for the North Atlantic, only
eelgrass was greatly affected (Cottam 1935a; 1935b; Cottam and Munro 1954). By contrast, this
recent loss of SAV has involved multiple species and appears to be local to the Chesapeake Bay
where freshwater, mesohaline, and marine species have been adversely affected, suggesting a bay-
wide deterioration in water quality. Meanwhile, other areas of the U.S. East Coast and western
Europe had not appeared to be losing SAV during this time period. In the Chesapeake Bay, a
"down-river, down-Bay" pattern of loss was demonstrated early in the Chesapeake Bay Program
further supporting the hypothesis that changes in water quality were related to SAV decline. At that
time, however, this hypothesis could not be supported scientifically because of a paucity of relevant
data.
In 1980, the first concerted efforts were begun to scientifically investigate Chesapeake Bay
SAV with the ultimate goals of understanding the cause(s) of the declines and the management
changes necessary to conserve and restore SAV in the Chesapeake and its tributaries. Over of the
next decade, SAV were studied at a variety of scales ranging from physiological to community-level
ecological investigations. From the outset, ecosystem process modeling has played a significant role
in the research. The results of this long-term SAV modeling effort are presented here and the
findings compared to the recently published technical synthesis document on SAV habitat
requirements and restoration targets (Batiuk et al. 1992).
SAV Modeling Goals and Objectives
The SAV modeling program goals and objectives fall into two areas: programmatic and
specific. For the programmatic area, there are several points which, in general, can be applied to
all modeling efforts. They are:
• to provide a conceptual framework in support of SAV research;
• to identify data and specific research needs for better understanding of the controls on,
and dynamics of, natural communities;
• to test using simulation analysis the stability characteristics of hypothetical model
structures and relationships; and
• to aid in generating alternative hypotheses that better explain the behavior of natural and
perturbed systems.
The goals and objectives specific to the SAV modeling program focus on those
environmental variables related to water quality that have been shown to be primary factors in
10
-------
controlling SAV distribution and abundance, depth distribution and primary production. For the
studies reported here, they are
• to simulate SAV growth as a function of in situ light, water depth, and temperature as
primary factors controlling photosynthesis and growth,
• to evaluate using simulation experiments physical-chemical controls and epiphyte-grazer
interactions on long-term SAV survival and stability, and,
• to evaluate using simulation experiments the effects of variable physical-chemical
regimes characteristic of in situ conditions on long-term SAV survival and stability.
SAV Management Needs
The 1992 amendments to the Chesapeake Bay Agreement emphasize the need to restore and
enhance living resources and their habitats. Toward this goal, links must be forged between
resource management, habitat restoration, and pollution reduction and prevention, in part by
connecting, scientifically and programmatically, resource management with habitat restoration
priorities. Simulation modeling is an important tool with which such linkages can be made. Access
to an ecologically-based simulation framework can give Bay Program managers a strong foundation
on which to base decisions and establish priorities in the years ahead. The initial support for
ecosystem modeling (funded in summer of 1991) represents a reaffirmation of interest by the Bay
management community in using these tools to investigate the impact of human perturbations on
the production of living resources.
The Ecosystem Process Modeling program is organized around five classes of models: 1)
Plankton-Benthos, 2) Littoral Zone/Submerged Aquatic Vegetation, 3) Emergent Intertidal
Wetlands, 4) Fish Bioenergetics, and 5) Ecological Regression/Visualization models. This report
focuses on submerged aquatic vegetation modeling as a separate component; the future of this
modeling effort, however, is the coupling of SAV with planktonic and benthic components in the
distinctive littoral zone (shallow, well-mixed waters) environment.
A major goal of the Ecosystem Modeling Program is to link model simulation experiments
with definable management endpoints. For instance, modeling studies will address quantitative
living resource restoration goals and/or habitat requirements and restoration targets that are currently
being developed under the direction and leadership of the Chesapeake Bay Program's Living
Resources Subcommittee. In the SAV model presented here, particular attention has been paid to
the relationship between environmental variability and biological interactions with the established
restoration goals for water clarity. The modeling studies have also addressed possibilities for
achieving the Tier n and EQ targets for SAV restoration with the lower, polyhaline portion of
Chesapeake Bay given these water quality standards.
SAV Conceptual Models
Two conceptual and simulation models have been built for SAV communities characteristic
of the lower Chesapeake Bay higher salinity habitats. These SAV communities are generally fco-
dominated by eelgrass, Zostera marina, and widgeongrass, Ruppia maritima. Eelgrass dominates
11
-------
deeper areas to a maximum depth of approximately two meters while widgeongrass occupies the
shallow, near-shore areas. For the results presented here, the information used in constructing and
validating the models reflect available data on the structural and functional ecology of eelgrass-
dominated communities.
The first conceptual and simulation model was an ecosystem model designed primarily as
a means to summarize available data on SAV productivity, to elucidate sensitive processes where
information was lacking, and to help guide a field research effort. The model, illustrated in Figure
4, follows the flow of carbon and includes the major trophic groups characteristic of lower bay
eelgrass communities, including two higher trophic level compartments representing blue crabs and
large predators such as blue fish, trout and sandbar sharks. Extensive simulation analyses were not
performed using this model due to the great uncertainty in many of the parameter estimates used for
simulation. Sensitivity analysis with the model, however, clearly indicated the dependency of many
compartments on the dynamics of the SAV component, whose functioning was greatly simplified
in the model. This led to the second generation conceptual and simulation SAV models designed
to better represent SAV growth and environmental controls.
Concurrent with these first modeling efforts, field and laboratory studies on Chesapeake Bay
SAV were beginning to identify the principal controls on SAV growth, the response of SAV to
various physical, chemical, and biological interactions, and the development of a data base sufficient
for model simulation and validation, data collected independent of model parameterization. During
the 1980s, sufficient data were collected from monitoring programs, field studies, and laboratory
experiments on SAV such that ecosystem process models capable of addressing estuarine water
quality issues and the distribution, abundance, and stability of SAV could be developed and tested.
A second conceptual model was developed which focused on eelgrass growth in relation to
major physical, chemical, and biological factors, most of which operate by modifying the amount
of light available for eelgrass photosynthesis. The conceptual model depicts the factors which
control submarine light intensity, as photosynthetically active radiation (PAR), reaching the eelgrass
canopy. Three PAR intensities are considered in this model:
1) PARD: the daily integrated solar irradiance at the water surface determined by
atmospheric conditions, latitude, time-of-year, and time-of-day;
2) PARz: the instantaneous PAR at a given water depth determined by attenuation within
the water column which, in turn, is primarily a function of water depth and the
suspended particle concentration (both organic and inorganic); and,
3) PARvp: the instantaneous PAR reaching the SAV leaf surface; i.e., PARz that is further
attenuated by epiphytic fouling on the leaf surface and self-shading by the eelgrass
canopy.
Although not explicitly modeled, dissolved inorganic nutrients, both nitrogen and
phosphorus, impact the relationship of light, photosynthesis, and carbon flow by stimulating the
growth of phytoplankton and epiphytes, reducing PAR available at SAV leaf surfaces, and shifting
the balance of inorganic carbon fixation from angiosperm plants (i.e., SAV) to algae. Grazers, in
part, may ameliorate this eutrophication effect by cropping both phytoplankton (not included in this
12
-------
0)
•o
o
a
2
0)
c
o
a.
75
I.
a>
o
o
O
rs a.
13
-------
conceptual model) and epiphyte cover on leaf surfaces, thereby increasing PAR available for SAV
photosynthesis and growth.
Using this conceptual model, simulation models were developed for analysis of specific
factors controlling SAV growth, distribution, and long-term community stability. For the latest
versions of the SAV model (version 4 results reported in a following section), particular attention
was given to evaluating the proposed SAV habitat light attenuation requirements for restoration that
were developed using empirical methods (Batiuk et al. 1992).
SAV Simulation Models
Since 1984, we have developed and analyzed four simulation models that are based on the
aforementioned conceptual model. The first two versions addressed various aspects of physio-
chemical controls on eelgrass photosynthesis, the effects and general role of epiphytic grazing on
plant productivity and survival, and the long-term stability of eelgrass communities under simulated
characteristic environmental conditions for various areas of the lower bay (van Montfrans et al.
1984; Wetzel and Neckles 1986). The site-specific data used for input included information
collected from sites that historically supported SAV but no longer did and sites that had healthy
populations which had been stable for long periods of time (i.e., decades), as determined from aerial
photography.
The third and fourth model versions have been implemented since 1991 as part of the current
Chesapeake Bay Program sponsored ecosystem modeling program. The third version is a
STELLA® implementation of the eelgrass model reported in Wetzel and Neckles (1986).
STELLA® is a commercially available simulation tool well-suited for ecological modeling using
time-dependent systems of ordinary differential equations (High Performance Systems 1992). The
fourth version is a revision of the original SAV model that includes new environmental data for
SAV habitats and addresses the effects of in situ variability relative to specific water quality
parameters. With the exception of version 3, all models were programmed in FORTRAN 77.
Eelgrass Simulation Model
Figure 5 gives the compartmentalization scheme and flow structure for the simulation models
(versions 1-4) using Odum's symbolic language (Odum 1971). The flows shown as solid lines
represent linear, donor-controlled fluxes. The dashed lines represent non-linear, donor and/or
recipient controlled fluxes. The dotted lines represent information flows that control negative
feedbacks which operate on specific fluxes (indicated on the flows by the open arrow symbol).
The flows are characterized by the principal components involved. All abiotic-biotic and
biotic-biotic interactions were modeled with non-linear, feedback-controlled functions. The
feedbacks were derived as density-dependent functions of either substrate (donor) concentrations
or recipient determined spatial constraints. Linear, donor-controlled functions were used to model
processes such as respiration and mortality for the biological components and some physical-
chemical exchanges (e.g., air-water CO2 exchange). Physical factors were modeled with empirical
or statistically derived functions. Wetzel and Wiegert (1983) have outlined the approach and
general techniques followed here. Similar applications have been reported more recently by
Christian and Wetzel (1991).
14
-------
Compartmental Polyhaline SAV Model
TIME OF YEAR
TIME OF DAY
WATER DEPTH
ATTENUATION
ATMOSPHERE
SOLAR
PAR
GRAZERS
EPIPHYTES
SUBMERGED
AQUATIC
VEGETATION
WATER COLUMN
'•" /*.'•'.'*'•'.'•..'•'.**. /•,'','*.'«'.*''•*.'*'•*.'•-*•
Figure 5. Compartmental and flow structure of the carbon-based SAV model for the lower,
polyhaline Chesapeake Bay. The modeled SAV population is eelgrass (Zostera marina),
which is divided into aboveground leaves and belowground roots-rhizome components.
Epiphytic coverage attenuates light and interferes with gas exchange at leaf surfaces. Incident
light follows diurnal and seasonal cycles and is attenuated within the water column by water
itself (depth) and by paniculate inorganic and organic matter and dissolved organic matter.
Losses of carbon from the system include burial, accumulation of detrital carbon, and losses
through grazers to higher trophic levels. Symbols are after Odum( 1971).
15
-------
Environmental Variables and Forcing Functions
The environmental variables and forcing functions modeled were solar and submarine PAR,
PAR attenuation as a function of water column turbidity and epiphytic fouling, tidally variable and
fixed water depths relative to the SAV canopy, and daily water temperature and photoperiod. Two
simulation model versions are referred to in the following sections, model version 2 and version 4.
Model version 2 incorporated fixed inputs for depth and attenuation coefficient designed to represent
annually averaged environmental conditions. Model version 4 incorporated stochastically
fluctuating inputs based on an eight-year data base for York River shoal habitats (see Moore 1992).
Table 1 gives descriptive statistics for those parameters used for version 2 simulations that
were derived from field data for shallow water habitats in the lower Chesapeake Bay.
Table 1. Environmental parameters and forcing functions used for the nominal case of
version 2 of the Polyhaline SAV Model. Annual means and ranges are for
smoothly varying functions, and were derived from observations in the lower
Chesapeake Bay (Wetzel and Neckles 1986). Depth and PAR attenuation
were fixed within model runs, and were changed for different model cases
(see text).
Variable
Temperature
Daily incident PAR
Water depth
PAR attenuation
Photoperiod
Units
°C
Ein nT2 d'1
m
m-1
hours
Annual mean
16.2
28.2
1.0
1.0
11.8
Annual range
2.5-30.0
11.5-45.0
9.5-14.0
For the non-variable model (version 2), the environmental parameters and forcing functions
were simulated using trigonometric functions fit to available in situ data (Table 2). The
downwelling PAR attenuation coefficient (Kd) was treated as a constant (annual mean) for different
simulation scenarios in model version 2 principally because of a lack of spatial and temporal data.
For the variable-forcing model (version 4), the natural variability characteristic of surface incident
solar irradiance (PAR0) and Kd was incorporated by using eight-year (1984-1992) monthly means
and ranges measured in situ as part of a bi-weekly monitoring program in shoal areas of the York
River estuary (Moore 1992). Version 4 model simulations calculated Kd as a variable with a random
component that varied within measured statistical limits. Two temporal scales were explored for
Kd variability, seasonal as defined by temperature and used in establishing habitat requirements
(Batiuk et al. 1992), and monthly as used for solar irradiance estimation. To derive Kd estimates,
three approaches were used: 1) site-specific statistical estimates using seasonal mean values only,
2) site-specific statistical estimates using seasonal means and the standard deviation, and 3) site-
16
-------
specific variation using monthly ranges, the annual minimum and random variability within the
observed range.
Table 2. SAV model version 2 forcing functions, based on environmental data from
the lower Chesapeake Bay.
Water temperature
Incident solar irradiance
Photoperiod
Tidal water depth
Submarine, time-
dependent PAR
irradiance
T = 16.25 - 13.75 cos[2rc(day-25)/365]
PAR0 = 28.25 - 16.75 cos[2rc day/365]
PP = 11.75 - 2.25 cos[27C day/365]
Tlag = 0.842105263 Lday
Z(t) = 1.25 + 0.275 cos[27t(hr-Tlag)/12]
PAR(t, hourly) = PAR0 / (0.63662 PP )
PAR(z) = PAR(t) exp[-Kd z]
cos[7i(hr-12)/PP]+
Note: day = day-of-year (1-365), Tlag = tidal phase lag, Lday = lunar day (0-28).
For the simulation studies, three sites, all in the York River, were used for data input and
parameter estimation for the environmental variables and forcing functions:
1. Guinea Marsh, a shoal monitoring site located at the York River mouth and having stable
eelgrass beds;
2. Gloucester Point, a shoal monitoring site, ca. 10 km upriver from Guinea Marsh, and
having eelgrass beds that are highly variable from year to year; and,
3. Claybank. a shoal monitoring site ca. 20 km upriver from Guinea Marsh which
historically marked the upriver limit for eelgrass, but which has not supported a stable
SAV community since the early 1970s.
These factors were incorporated as controls on various pathways of carbon flow, most of
which also operated within intrinsic biological limits. The light-related parameters controlled
photosynthesis by eelgrass and epiphytes. Temperature controlled the rate coefficients for
photosynthesis and respiration by eelgrass and epiphytes, root/rhizome translocation by eelgrass,
ingestion and respiration by grazers, and the daily ration of higher level predators, further discussed
below.
17
-------
Biological Processes, Interactions and Controls
The biological processes simulated were photosynthesis by eelgrass and epiphytic
microflora, ingestion of eelgrass leaves and epiphytes by grazers, respiration and natural mortality
for all biological components, and seasonally controlled immigration-emigration of grazers and
higher level predators. The mathematical functions derived to simulate these processes follow the
rationale developed originally by Wiegert (1973). Applications for other ecosystem models are
given by Wiegert and Wetzel (1979), Christian and Wetzel (1978), van Montfrans et al. (1984),
Wetzel and Christian (1984), Wetzel and Neckles (1986), Christian et al. (1988), Neckles and
Wetzel (1989) and Christian and Wetzel (1991). Details of the mathematical derivations are given
in the following sections.
Eelgrass Photosynthesis and Growth
The processes modeled that determined eelgrass aboveground (shoots, or leaves) and
belowground (roots and rhizomes) biomass were photosynthesis, respiration, shoot and root/rhizome
translocation, and mortality of shoots and root/rhizome material. Photosynthesis was modeled as
a function of a specific rate coefficient (Py) that was both light- and temperature-dependent, leaf
biomass (Xj), and non-linear, negative feedback controls which reduced photosynthesis when either
CO2 or space (SAV density) became limiting, the terms FBy and FBjj respectively. The general
form of the flux equation, Fy, for CO2 fixation by eelgrass leaves was
F^PjjXjKl-FBgXl-FBjjCij)] (1)
The specific rate coefficient for photosynthesis (Py) was calculated using a rectangular hyperbolic
(Monod) function which is dependent on the light-saturated photosynthesis rate, Pmax at
temperature, T, the PAR intensity reaching the leaf surface, PARyp, and the half-saturating PAR
intensity, Ik':
Py = PmaxCTEMP) [P ARvp/ak'+PARvp)] (2)
Pmax(TEMP) was derived as a linear function of temperature using the data of Wetzel and Penhale
(1983) and Evans (1984); for temperatures greater than 25°C, Pmax(TEMP) was reduced linearly
such that at 30°C, Pmax(30) equaled one-half the rate at 25°C (Equation 3).
Pmax(TEMP) = (0.000162*TEMP + .0041)*[1 - (TEMP-25)/(35-25)] (3)
Ik' then was calculated as a linear function of Pmax(TEMP) based on Penhale (1977) (Equation 4).
Ik' = 1 5220.78*Pmax(TEMP) - 3 1 .86 (4)
The PAR intensity reaching the plant surface (PAR^, ) was derived as a function of daily solar
irradiance (PARD), water depth (Z), the water column PAR attenuation coefficient (Kd), and PAR
absorbance by epiphytes (EPLR). PAR intensity at the top of the plant canopy was calculated as a
negative exponential function of PARD, Z, and Kd. PAR attenuation due to epiphytes was derived
as a hyperbolic function (Equations 5 and 6) of the ratio of epiphyte (X05) to eelgrass leaf biomass
(X03), which approached 75% at maximum epiphytic colonization (Murray 1 983). The relationship
18
-------
was derived from Murray's field and experimental data for leaves with a minimum reported biomass
ratio of 0.5. Since new, unepiphytized leaves are produced continuously within the canopy, direct
application of this relationship would significantly over-estimate epiphyte-induced PAR attenuation.
To account for this, we assumed that 50% of the leaf biomass at any time was composed of young
leaves having no significant epiphyte biomass, and therefore we reduced the predicted epiphyte
attenuation by half.
PARvp = (1 - 0.75*EPLR)*PARZ (5)
EPLR = 0.25*[(X05/X03) - 0.1]1/2, where 0.0 < EPLR < 1.0 (6)
Finally, Sand-Jensen (1977) showed that at high light intensity and high epiphyte biomass,
diffusion of bicarbonate limited eelgrass photosynthesis by 30%. We incorporated this effect in a
manner analogous to epiphytic PAR attenuation as a hyperbolic function of the epiphyte-leaf
biomass ratio (Equation 6).
The final terms in Equation 1, FB^ and FBy, that potentially limit eelgrass photosynthesis
are non-linear feedback functions of CO2 concentration and leaf biomass, respectively. The first
of these terms, FBy, which is the substrate- or donor-control feedback, was derived as a function of
the ambient CO2 concentration (Xj), the CO2 concentration below which carbon became limiting
(Ay), and the CO2 concentration (Gy) at which photosynthesis by eelgrass approached zero
(Equation 7).
y) (7)
The second of these terms, FBy, which is the recipient- or self-control feedback, was derived as a
function of leaf biomass (X:), the leaf biomass (A-) above which space or some density- dependent
factor (e.g., self-shading, crowding, nutrient limitation) limited growth, and the maximum leaf
biomass (Gjj ) that could be maintained metabolically (i.e., analogous to "carrying capacity" in
population models) (Equation 8).
FBjj = (Xj-Ajj)/(Gjj-Ajj) (8)
Both of these density dependent feedbacks are constrained mathematically to assume values in the
range 0.0 to 1.0.
The final term, Cy in Equation 1, is a metabolic correction factor which allows for
maintenance (dXj/dt = 0) at maximum standing stock (X: = GjA It was calculated as
(9)
. where R: equals the specific rate coefficient for leaf respiration and Py is as defined previously. The
correction term was derived such that the instantaneous rate of CO2 uptake exactly equalled the rate
of metabolic loss when no other factors limited growth at the maximum standing stock.
Natural losses of eelgrass leaf biomass occurred through respiration, mortality, and
translocation of organic carbon compounds to roots and rhizomes. Respiration and leaf mortality
19
-------
were derived as linear functions of leaf biomass and a specific rate coefficient. The specific rate
coefficient for respiration was estimated as the sum of a basal rate (Murray 1983; Murray and
Wetzel 1987) and a specific rate operating only during the photoperiod that was linearly dependant
on the realized rate of photosynthesis (B iebl and McRoy 1971; McRoy 1974). Both rate coefficients
were temperature dependent and statistically derived (Biebl and McRoy 1971; Nixon and Oviatt
1972; Murray 1983). The specific rate coefficient for leaf mortality was calculated using a cosine
function and ranged from 0.5% to 3.0% per day with the maximum rates occurring in mid-summer
(Vaughan 1982).
Translocation was perhaps the least documented process affecting eelgrass dynamics.
McRoy (1974) indicated that the maximum potential translocation to roots and rhizomes was 17%
of net leaf organic matter production. We set the maximum specific rate coefficient at 17% of net
productivity and reduced the rate using negative feedback control functions when either
aboveground or belowground compartments became limiting. The derivations are analogous to
Equations 3 and 4.
Epiphytes
Epiphytes were modeled as an autotrophic community dominated by microflora (Murray
1983; Murray and Wetzel 1987; Neckles 1990; Neckles et al. 1993). Processes affecting epiphyte
dynamics were photosynthesis, respiration, natural mortality, and grazing.
Photosynthesis by the epiphytic community was derived mathematically as for eelgrass.
Pmax and temperature were positively and linearly related up to 25°C (Penhale 1977). Above 25°C,
Pmax declined linearly such that at 30°C, Pmax was 75% of the rate estimate at 25°C (Penhale 1977).
Ik' was linearly related to temperature in the range 10 to 30°C and increased from 50 to 150 uEin
m~2 sec"1 (Penhale 1977). PAR intensity reaching the epiphyte community was derived as for
eelgrass except the only reduction in light was due to water column attenuation. Similar to eelgrass
(1), photosynthesis by epiphytes was derived as a function of a specific rate coefficient given
available PAR and temperature, the compartment biomass, and non-linear feedback control
functions. For epiphytes, however, the limiting factor for colonization and growth was leaf surface
area. We reformulated the density- dependent feedback, FBjj, as the ratio of epiphyte (X:) to SAV
leaf (Xj) biomass rather than epiphyte biomass alone. The relationship between leaf surface area
and biomass was statistically derived from our own unpublished data. The feedback function was
formulated as
FB^ = KX/Xj) - Ajj] / (Gjj - AJJ) (10)
where A:: and G^ were biomass ratios. Natural losses due to respiration and mortality were modeled
as for eelgrass except the specific rate coefficient for mortality was fixed at 0.5% per day.
Grazers
Grazing on epiphytes were modeled using data available for isopod and amphipod
populations typical of lower Chesapeake Bay eelgrass communities. Their dynamics were
determined by preferences for, and ingestion of, organic resources, resource-specific assimilation
20
-------
efficiencies, egestion, respiration, natural mortality, seasonally variable immigration-emigration,
and loss via predation by higher trophic levels.
Ingestion of epiphytes by grazers was derived as a function of the preference or selectivity
for a specific resource, the maximum potential specific rate of ingestion, the grazer biomass, and
non-linear feedback control functions dependent on resource availability and grazer density. The
mathematical form of the equation was analogous to Equation 1 for eelgrass photosynthesis.
Ingestion of eelgrass leaves was derived in the same manner. The preference values and the
maximum potential ingestion rates of eelgrass leaves and epiphytes were derived as functions of
temperature and based on data for gammarid amphipods (Zimmerman et al. 1979). As in Equations
3 and 4, the feedback control terms were derived as functions of resource availability and grazer
density.
Immigration and emigration were not modeled explicitly. However, Marsh (1973) found few
epifaunal organisms in Chesapeake Bay eelgrass communities from November to March when water
temperatures are generally below 10°C. Therefore, seasonality in grazer standing stocks was
incorporated by restricting epiphytic grazing to periods when water temperature was greater than
10°C. Mid-summer declines in grazer population densities have also been documented and
correspond to summer increases in predatory fish densities (Orth and Heck 1980; Diaz and Fredette
1982). Loss of grazers to fish predation was incorporated implicitly in the model as a function of
date and temperature and was based on the seasonal pattern of predatory fish abundances in eelgrass
meadows of North Carolina (Adams 1976a 1976b) and their reported food preferences, daily rations,
and assimilation efficiencies for specific prey (Hoss 1974; Adams 1976c; Peters et al. 1976).
Respiration by grazers was modeled as a function of temperature. Grazer mortality was fixed at
1.0% per day.
Numerical Computation and Simulation Analysis
The models were programmed in FORTRAN 77. All time-dependent equations were solved
using simple Euler numerical integration (Wiegert and Wetzel 1974). The time step for integration
was 1.0 hour. Appendix A contains the complete FORTRAN 77 source code listing for the SAV
simulation model.
Simulation analyses consisted of first establishing a "nominal" run by repeated simulations
of the model until limit-cycle behavior (within three years) and long-term stability occurred (over
five years) and the predicted rates, material fluxes, and standing stocks agreed with available field
and/or laboratory data. Various models were then run to explore the effects of selected physical-
chemical and biological perturbations on eelgrass dynamics. The nominal case and test cases all
were started from the same initial conditions. Two complete series of analyses have been completed
to date using two model versions (versions 2 and 4).
Model version 2 simulation analyses addressed physical-chemical controls and biological
interactions without incorporating environmental variability and portions of these studies have been
reported in the literature (Wetzel and Neckles 1986). The results of these past model studies will
be briefly reported here for continuity. The analyses were divided into two series of simulations.
The first series included physical (environmental) controls only. The second series addressed the
potential for epiphyte-grazer interactions to control eelgrass dynamics by varying the grazing
21
-------
pressure on epiphytes. Grazing pressure was varied under both nominal and perturbed physical
conditions in the model.
A final series of simulation studies using version 2 explored implicitly the effects of water
column nutrient enrichment and epiphyte grazing intensity (Neckles 1990). The effects were noted
by varying surrogate variables in model simulation runs. In this case, increases in the rate
coefficient for epiphyte photosynthesis was used as the surrogate for nutrient enrichment.
Version 4 simulation analyses addressed physical-chemical controls on eelgrass growth,
abundance and depth distribution relative to specific water quality parameters. Simulation studies
investigated in situ variability of environmental parameters for shoal areas that encompass SAV
historical or present distribution patterns. With the incorporation of statistically-defined random
variability in the environmental forcings, specifically with regard to the underwater light climate,
modeling scenarios explored whether such fluctuations might be important in the natural cycle of
SAV productivity and in determining water quality criteria for SAV (Batiuk et al. 1992).
SIMULATION RESULTS AND DISCUSSION
Nominal Simulation Analyses (Vetsion 2: non-variable inputs).
Figure 6 illustrates the measured and simulated values for three of the four principal physical
forcing variables in the model, water temperature, solar irradiance, and the downwelling PAR
attenuation coefficient. Tidally varying water depth relative to plant canopy height is not shown.
The high degree of natural variability is evident in all three and this variability is obviously not
completely described by the trigonometric functions used to simulate these variables. Because of
the high in situ variability, general lack of data at the time, and poor agreement between the
simulated and observed values, the attenuation coefficient was fixed at 1.0 m"1, a value
characteristic of all observations from a historically stable eelgrass bed (Moore 1992). Also, for the
version 2 nominal simulation, the depth was fixed at 1.0 m (MLW) for comparative purposes,
because this represented the predominate depth of in situ sampling within existing eelgrass beds.
Therefore, nominal model predictions for the various compartments best represent long term
averages and do not address spatial or temporal variability.
The simulated dynamics for the standing stocks of eelgrass leaves, epiphytes, and grazers under
nominal conditions are given in Figure 7. All components demonstrated limit cycle behavior and
were stable over time (up to ten years of simulated dynamics). Eelgrass biomass maxima occurred
in late April, reaching 143 g C m~2. Epiphyte biomass reached a maximum of 79 g C m~2 in late
June. The grazers feeding on these epiphytes reached their maximum biomass, 0.93 g C m~2, during
the late summer. The epiphyteieelgrass biomass ratio was 0.04 during the early spring growth of
eelgrass and reached 1.5 during the late summer eelgrass decline. These values are consistent with
field estimates (Murray 1983; Neckles 1990; Neckles et al. 1993). At maximum epiphyte biomass,
PAR reaching the plant epidermis was attenuated 22% by epiphytes and eelgrass photosynthesis was
reduced 10% by epiphyte-limited CO2 diffusion. The maximum sustainable standing stock for the
above-ground eelgrass component (leaves) was set to 150 g C m~2 in the density-dependant
feedback control function (Equation 4), which was the maximum reported field density for
22
-------
Physical Forcing Functions
A. Water Temperature
I 30.00
2 20.00
| 10.00
£ 0.00
: ^xr^^
-xX
• » ^/^^^^^^
s? "~v— ^--,
A *-"v_
^^^*^^Mf
^^^ ^^
15 45 75 105 135 165 195 225 255 285 315 345
Day
/**.
I 60.00
1 40.00
& 20.00
5 o.oo
B. Solar
PAR
M/ ^K— — ye— — *— - -NT'
• JIT m^"^ r\ /t* ^^/K~ * —W
X— ^""^""^ ^"^*"X-*-X—
15 45 75 105 135 165 195 225 255 285 315 345
Day
C. Attenuation
1.50j
, 1.00- -
* 0.50 -•
nnn..
*_v- ---.,,-X-y
* *^^K x ^
Xv
_ _ _^ _v v
JK^^* .
15 45 75 105 135 165 195 225 255 285 315 345
Day
Figure 6. Modeled (dashed line) versus observed (*) physical parameters for the polyhaline SAV
model. The in situ values are the monthly means of biweekly observations taken over the period
1985 to .1992 for water temperature and the attenuation coefficient (Kd). The solar PAR data are
the monthly means of daily integrated PAR taken at the VIMS over the period 1988 to 1991.
23
-------
CNI
'E
o
o>
150
Polyhaline SAV Model
Biomasses for the nominal case with
constant Z and Kd, and smoothly-varying PAR
1.00
- 0.75
- 0.50
- 0.25
CVI
'E
o
at
¥
OJ
o
CO
<5
o
0.00
365
730 1095
Model Day
1460
1825
Eelgrass
Epiphyte
- Grazer
Figure 7. Biomasses of eelgrass leaves, epiphytes, and grazers in the polyhaline SAV model nominal
case, with depth fixed at 1.0 m, attenuation coefficient fixed at 1.0 m"1, and using a smoothly
varying annual PAR cycle, designed to mimic healthy habitats in the lower Chesapeake Bay.
24
-------
communities in the lower Chesapeake Bay (Orth and Moore 1983). For the nominal simulation, this
feedback function reduced potential eelgrass growth by 88% during periods of maximum biomass,
suggesting that some density-dependent factor (e.g., self-shading) may limit SAV growth at certain
times of the year.
The physical factors controlling eelgrass production in the simulation are submarine
irradiance and temperature. Simulated eelgrass photosynthesis (Figure 8) is light-limited throughout
the year and particularly during mid-summer. Temperature (Figure 8B) also contributes to the mid-
summer decline as the temperature-dependent respiration rate exceeds the rate of photosynthesis
(Figure 8A). Eelgrass is near its equatorward biogeographic extent and the high summer water
temperatures (>25°C) impose a significant stress on the plant (Wetzel and Penhale 1983; Evans
1984). Together, light limitation and thermal stress contributes to the mid-summer decline in
eelgrass production and biomass. Production (Figure 8A) and biomass (Figure 7) do increase in the
fall, with the increase in biomass continuing through the following spring.
The overall pattern and standing stock predictions agree with field data but the temporal
pattern is out of phase with some data sets; i.e., simulated standing stocks precede in situ
observations by ca. 30 days. This discrepancy could have several explanations. First, given the
temperature dependence of SAV photosynthesis and respiration, unusually warm or cold years
would lead to discrepancies between the temporal pattern of in situ observations versus simulated
results based on multi-annual averaged temperature cycle. Second, the factors that control the
magnitude and pattern of translocation of organic matter to roots and rhizomes are poorly
understood. The assumption of a constant 17% of net carbon fixation is an oversimplification and
requires better information. Third, the model does not account for the energetic costs of non-
vegetative reproduction; flowering and seed production in the spring would lower or delay the
accumulation of above-ground biomass. Fourth, solar irradiance and water column PAR attenuation
are highly variable. For model version 2, a non-variable, continuous function was fitted to local data
to predict daily solar irradiance, and the PAR attenuation coefficient was fixed at an annual mean
value. Both should be modeled as stochastic functions that operate within the range of natural
variability, which are addressed in model version 4 studies (see below).
These nominal simulations suggested that in situ eelgrass photosynthesis and growth are
governed primarily by submarine irradiance and temperature. Maximum standing stock and annual
production appear limited by these physical factors and possibly by density-dependent controls that
operate when submarine irradiance is high and temperatures are optimal. Sediment dissolved
inorganic nutrients, particularly nitrogen, may limit above-ground growth at these times but, in
general, do not limit the current distribution or preclude growth of eelgrass in the lower Chesapeake
Bay (Orth 1977). Although physical and chemical factors limit eelgrass growth and production to
less than its intrinsic maximum, the simulated dynamics of the model's components demonstrate
long term stability and predicted densities agree with field estimates. Qualitatively, the model can
be applied as an analytical tool for simulating the effects of various changes in physical-chemical
and biological interactions on SAV productivity (Wetzel and Wiegert 1983).
25
-------
Polyhaline SAV Model
Eelgrass production and respiration in the nominal model case
0.100
Jan
Figure 8. Polyhaline SAV model nominal case with fixed depth and light attenuation coefficient.
A. Rates of photosynthesis and respiration in the nominal model case. B. The annual temperature
cycle used in the model. While photosynthesis is light-limited much of the year, temperature drives
respiration rates higher than photosynthesis during the summer months.
26
-------
E 150
O
100 -
g
in
fO en
0)
III
Polyhaline SAV Model
Sensitivity to fixed values of depth and light attenuation
A. Various depths, Kd = 1.0 m"
0.5m
1.0 m 1.5 m
2.0 m 2.5 m
B. Various Kd's, depth = 1.0 m
730 1095
Model Day
1460
1825
1.0m'1 1.5m'1
" 1.75m"1 2.0m"1 2.25m'1
Figure 9. Sensitivity of eelgrass leaf biomass predicted by the polyhaline SAV model to changes
in depth and light attenuation coefficient (Kd). A. Various depths with Kd fixed at 1.0 m"1. B.
Various Kd's with depth fixed at 1.0 m.
28
-------
Model Responses to Physical Regime and Epiphyte-Grazer Interactions
Physical Regimes
The sensitivity of the model to different physical regimes was evaluated by varying water
depths and average annual water column PAR attenuation. Changes in water depth were simulated
by fixing the depth at 0.5,1.0 (nominal), 1.5,2.0, and 2.5 m. For 0.5 m, the upper temperature limit
was also increased to 33°C (Wetzel, unpublished data). The PAR attenuation coefficient was varied
through fixed values 1.0,1.5, 1.75,2.0, and 2.25 m"1, based on field data for existing or historical
eelgrass habitats in the lower Chesapeake Bay (Moore, unpublished data).
Changing the mean water depth from 1.0 to 2.0 m resulted in stable but lower eelgrass annual
peak biomass estimates, from 143 (1.0 m) to 110 g C m"2 (2.0 m) (Figure 9A) and a decrease in the
importance of density-dependent feedback control. At 1.0 m fixed depth, density-dependent
controls limited eelgrass leaf growth by a maximum 93%, whereas at 2.0 m, growth was limited by
50%. This suggests that the relative importance of density-dependent factors changes not only
temporally but also as a function of depth within the community. At fixed depths less than 1.0 m
and greater than 2.0 m, a predicted loss of above-ground eelgrass biomass occurred as a result of
temperature effects in very shallow water and PAR limitation at depth. Additionally, eelgrass leaf
biomass was sensitive to changes in the average water column PAR attenuation (Figure 9B).
Simulations indicate that eelgrass will not survive over long time periods (i.e., years) with
attenuation coefficients greater than 2.0 m"1. Modeled eelgrass population survival was somewhat
less sensitive to Kd than has been suggested for natural population (Wetzel and Penhale 1983).
Orth and Moore (1983) surveyed the relative density and depth distribution of SAV
communities throughout lower Chesapeake Bay and reported that monospecific Eelgrass
communities typically occurred to a maximum depth of 1.2-1.6 m relative to mean sea level. Long
term studies (surveying, transplantations, and experiments) across a continuum of potential SAV
sites, from presently non-vegetated to well-established, has shown that eelgrass will not survive
where the water column PAR attenuation coefficient averages >1.5 m"1 (Moore 1992; Batiuk et al.
1992). The simulation analyses are consistent with these data and support the hypothesis that
submarine irradiance is a principal factor determining productivity and long term survival of
eelgrass in Chesapeake Bay.
Epiphyte-Grazer Interactions
The correspondence of eelgrass biomass and productivity with the annual variation of
physical factors suggested that the epiphyte and grazer components had little effect on eelgrass
productivity. To test this, additional simulations were run without epiphytes and their grazers. In
these simulations, predicted dynamics and standing stocks of eelgrass leaves were nearly identical
to those in the nominal simulation. Predicted photosynthetic rates also were similar most of the
year, but increased slightly during mid-summer relative to nominal simulations. It can be inferred,
then, that in the nominal case, grazing maintains epiphyte density below levels that would limit
eelgrass photosynthesis and growth. As a corollary, any perturbation that would allow epiphytes
to outgrow their grazers might be expected to have an adverse effect on eelgrass productivity. In
the absence of impacts on the epiphyte-grazer components, any perturbation which would limit
eelgrass photosynthesis would also adversely effect eelgrass productivity.
27
-------
The simulated response of eelgrass to changes in epiphytic grazing pressures combined with
various levels of water column PAR attenuation is shown in Figure 10. Under nominal physical
conditions (Kd =1.0 m"1), eelgrass biomass declined with decreased epiphyte grazing pressure.
However, even at very low grazing intensity, the community persisted. The effect of reduced
grazing and the concomitant increase in epiphyte biomass became more dramatic as PAR
attenuation increased. For example, nominal grazing with Kd < 1.75 nT1 resulted in long term SAV
stability, whereas grazing pressures less than 25% of those presently observed in the environment
caused a predicted diminution or, under the more severe cases, a complete demise of eelgrass leaf
biomass where light attenuation coefficients were equal to or greater than 1.5 m"1 (Figure 10). Of
note, at Kd's of 1.75 and 2.0 nT1, where modeled eelgrass populations survived under nominal
grazing pressures (Figure 9), they failed to survive in combination with severely reduced grazing
(Figure 10).
These simulations indicate that under reduced light conditions, additional factors, such as
grazing, which control epiphytic growth can impact eelgrass production and community stability.
For eelgrass communities typical of the lower Chesapeake Bay which are stressed naturally by sub-
optimal submarine irradiance and high, late-summer water temperatures (>25°C), relatively small
changes in the factors controlling epiphyte growth may greatly affect community stability, recalling
that epiphytes interfere with both light transmission to the leaf surface as well as gas diffusion across
the leaf surface. For example, increased water column nutrients and reduced grazing favor epiphyte
growth and reduce SAV photosynthesis (Murray 1983; van Montfrans et al. 1984; Borum 1985;
Twilley et al. 1985). These simulation experiments suggest that under conditions of reduced grazing
intensity and increased PAR attenuation, the SAV community would not survive over longer time
frames. One might also expect that conditions which accelerate intrinsic epiphyte growth would
also have an adverse impact on eelgrass communities.
To more closely examine the impact of the epiphyte assemblage on eelgrass, Neckles (1990)
extended the SAV simulation model using the results of mesocosm experiments. These laboratory
experiments investigated plant-epiphyte-grazer interactions relative to eelgrass photosynthesis,
growth and community stability. Grazer effects (presence/absence) on epiphyte density and SAV
growth were examined under both ambient and nutrient-enriched (three-fold ambient concentration)
conditions in mesocosms.
However, because the model does not explicitly include nutrients, a surrogate variable was
used for these simulation studies. It has been shown that a primary effect of increased water column
nutrient concentration is an increase in both the realized growth rate of microalgae and the leaf-area-
specific density of epiphytes (Murray 1983). To simulate this, Neckles (1990) adjusted the model's
epiphyte Pmax upward until predicted epiphyte densities agreed with observed densities from the
experimental mesocosms. Increasing Pmax by two- and three-fold appeared to best represent
increased nutrient effects relative to lower bay SAV. The simulations were run for ten-year periods
to investigate the effects of eutrophication, epiphytes, and grazers on community stability.
Figure 11 shows model predictions (solid and dashed lines) and mesocosm observations (x
± SE) of the mass ratios of epiphyte to SAV leaf. The upper panels (A & B) in Figure 11 give the
results for ambient nutrient conditions with and without grazers present. The observed and model
predicted values are remarkably similar for this condition. The outlying observation (panel A during
the fall) resulted from an increase in ambient nutrient concentrations during the mesocosm studies
29
-------
150
Polyhaline SAV Model
Annual Eelgrass biomass maximum and the
interaction of light attenuation and epiphytic grazing
1.00
Grazing Intensity
100%
50%
25%
I I 10%
(Moribund/No survival)
1.25
1.50
1.75
2.00
PAR Attenuation Coefficient (rrf1)
Figure 10. Predicted eel grass leaf biomass with combinations of attenuation coefficients and grazing
pressure (relative to the nominal case) using polyhaline SAV model version 2.
30
-------
Epiphyte:SAV Biomass Ratios - Grazing and Eutrophication
Nominal Epiphytic Growth
Nominal Grazing No Grazing
3 -
.22-
co
CC
1 -
0
Jan
3 -
02-
03
CC
1 -
- 3
- 2 .0
(0
CC
- 1
Apr Jul Oct Jan Apr Jul Oct Jan
Eutrophication Cases
Nominal Grazing No Grazing
Jan Apr Jul Oct Jan Apr Jul Oct Jan
- 3
-20
CO
CC
- 1
Month
Month
Figure 11. Predicted biomass maximum of eelgrass leaves during model year 5, under combinations
of PAR attenuation coefficient (m"1) and relative grazing intensity (percent of the nominal case) in
the Polyhaline SAV Model using a fixed depth of 1.0 m.
31
-------
which the model did not simulate (Neckles 1990). The lower panels (C & D) in Figure 11 give the
same comparison except under nutrient enriched conditions. With grazers absent (panel D), model
simulations and experimental observations agreed very well except during summer. The nutrient
enriched mesocosms during this time developed dense macroalgal populations which the model
cannot simulate. With grazers present (panel C), the agreement between model predictions and
mesocosm results was not as good. The late summer observations agreed with a two-fold Pmax
increase while the fall mesocosm experiments agreed with a three-fold Pmax increase. Neckles
(1990) noted that the mesocosms during the fall experiment had nutrient concentrations much higher
than at other times. The poor agreement during the summer resulted in part from the dense growth
of macroalgae in the mesocosms as noted before.
To investigate the long term effects of increased nutrients and grazer interactions on
community stability, simulations were run using nominal and a three-fold increase in epiphyte Pmax
both with and without grazers present for periods of five years. Relative to natural grazer densities
and in situ nutrient concentrations in lower Chesapeake Bay tributaries where eelgrass did or
currently does exist, these conditions would represent extremes. All other parameters of the model
were held at nominal conditions typical of contemporary lower bay eelgrass habitats.
Under ambient nutrient conditions (Figure 12A), eelgrass biomass remained stable over
five-year simulations both with and without grazers present. The model indicated a ca. 20-25%
reduction in peak standing stock, however, without grazers. These model results are consistent with
all previous simulation studies and the available field data and literature information.
For enriched nutrient conditions (Figure 12B), a more dramatic pattern emerged. With
grazers present, eelgrass persisted over the five-year simulation period, but suffered a 30 to 40%
reduction in abundance relative to nominal simulations. With grazers absent, the eelgrass population
was not stable over time and demonstrated a protracted decline in abundance. In nature, for
conditions such as these, once eelgrass populations reach a threshold minimum density, they would
cease to be viable and would be replaced with a different community structure (e.g., macroalgal or
planktonic).
These results indicate the importance of the interactive effects of physical, chemical and
biological controls for determining the distribution, abundance and long-term stability of eelgrass
communities in the lower bay and, by analogy, SAV communities throughout the bay. In addition
to the direct influence of physical-chemical (i.e., temperature, salinity, hydrodynamic regime)
properties on SAV productivity and community composition, a major control appears to be the
intensity of submarine PAR reaching the plant leaf to support photosynthesis. Both grazing and
nutrient enrichment appear to act as controls via modification of the submarine light environment.
Effect of Submarine light Variability on Eelgrass Dynamics (Model Version 4)
Field, laboratory, and simulation studies all indicate that the submarine light environment
is fundamental to the depth distribution, productivity, and stability of SAV communities.
Understanding the factors attenuating photosynthetically available radiation (PAR) is critical for
developing effective SAV management policies and implementing strategies for SAV conservation
and enhancement.
32
-------
150
Polyhaline SAV Model
Impact of grazing pressure and eutrophjcation effects
on epiphyte cover as reflected in SAV biomass
A. Nominal nutrient conditions
B. Eutrophic conditions
365
730 1095
Model Day
1460
1825
Nominal grazing No grazing
Figure 12. Epiphyte:Eelgrass leaf biomass ratios predicted by the polyhaline SAV model (with
fixed depth of 1.0 m) under combinations of grazing and eutrophication conditions. Dots are mean
and standard error of ratios from manipulated microcosms (Heckles 1990; Neckles et al. 1993). A.
Nominal grazing and epiphytic growth rates. B. Nominal epiphytic growth rate without grazing.
C. Nominal grazing rate with enhanced epiphytic growth rate, as a proxy for eutrophication. D.
No grazing with enhanced epiphytic growth rate. Figure modified from Neckles 1990.
33
-------
The simulation model as configured in version 2 treated PAR and the factors attenuating
PAR in the water column simply and without regard for environmental variability. The simulation
results presented for version 2 modeled solar irradiance as a simple cosine function fitted to local
observations. Water depth and the attenuation coefficient were fixed or held constant for a given
simulation.
Variability in solar irradiance, the downwelling PAR attenuation coefficient, and variable
water depth driven by tides, as well as absolute depth relative to MLW, determine the intensity and
temporal variability of PAR reaching the SAV canopy. To gather data to better model subm'arine
PAR, an intensive water quality monitoring program was began in 1984 in the York River estuary
to characterize, among other things, solar irradiance, submarine PAR, and PAR attenuation at
selected shallow water sites that either currently (Guinea Marsh) or historically (Claybank)
supported SAV or represented marginal SAV habitat (Gloucester Point). For the purposes of model
revision and simulation studies, the analyses concentrated on using these three site-specific data sets.
These areas represent SAV habitats that grade from healthy and stable to unsuitable for SAV under
present water quality conditions. In the following, the results of simulation studies with model
version 4 are presented on environmental variability in solar PAR, PAR attenuation, and tidally
varying water depth relative to eelgrass depth distribution, abundance (biomass) and long term
community stability. All other model characteristics were the same as in version 2 (e.g., grazer-
epiphyte interactions).
Solar Irradiance
The first revisions to the model were made to incorporate the natural variability of solar
irradiance and determining the effect on selected rates and state variables. Figure 13 shows the
daily record of solar irradiance collected at VIMS for the period 1988 to 1991. These data indicate
the high day to day variability as a result of local weather conditions and seasonal patterns for the
area. Pooling all observations by month indicated that the greatest increase occurred between March
and April, the highest average daily irradiance occurred during May and the greatest decrease
occurred between September and October. The degree of variability ranged from ca. 50%
(coefficient of variation) during January, February, and March to 25-30% (coefficient of variation)
from May to August. This differs greatly from the simple cosine model used for version 2.
Two stochastic functions were developed to better simulate this behavior in solar irradiance
compared to known variability, specifically based on the VIMS data set. The first description
(RANDOM-1) for daily irradiance was based on monthly averages and standard deviations
(Equation 11). This daily solar PAR (PARD; Bin m~2 d'1) was derived as
PARD = [AVG - SD]+ 2*RAND*SD (11)
where AVG and SD are the monthly means and standard deviations calculated from the VIMS data
set and RAND is a random number in the range 0.0 to 1.0. A second description (RANDOM-2) was
based on monthly ranges, in place of standard deviations, of the pooled data and PARD derived as
PARD = RAND*RANGE -I- MIN (12)
34
-------
:s *•%'••••• • •
. / • . • •. . • •
«/«««. .• •*• • ,
•• --Z" i •••:.•••
*w A • •
—* •'••«.•• "V '
T- •..•*<*•.
5> vr .-. . . •• • ••
? •#• • • •'.
T iVV; •: • •
00 /V •*. '"••
o> •*^J»''«* *. *i • " •• •
T- . •. V * • •
«*^ ., .j • C
1 .. r-fe •.-':"
s> •'
5 s ' • * .'.•• • ' ' . .. .1 .
* •"-^•V. ••••• •:•
o
CL .. .,•.. . •
. «Y. r f •».
•• « • p
411 .*!••••«
M * • 1 . • « • ,
tfl • f. • ••'• •
in " f/* • • * *
§ j /•-/'. .-.-;. ' .-.-•
•^ •^*»*«* "
0 ' 'i. '.< ' • '• •'
O« j. . . • • ••
•.«•„.*/• . * _
i-i •«•/*•."
^2 • • •„ • •• . ••
CB **j*«€ r •
*** • • • ^. • • •
«••*•• «^
g • :.'•••• •; /-• • •».:. •
c • ••«.../ •*•... \ j.
CB •«••%. /• •
•—• • -/ 1 t • .••••«
•o ' . ' » ...>•.••
co ' ^ •.. ; .*. -
I "'*•.;••• '-. /•*•:
& &>'. '«.
& X -. • • JL •
S **X'*-*
•a ..'\ V. ' • • • . ". •
{ ,.:-''!V".----. V
S? " ^t •*••
ffi ^ fc.. •••••• .
z '• •> • *. •
_ • » • • •
«••«•••.••
».'..•• • - • - ••
t« • • . * . •
. . u -. . 1 1 1 |'"^''l~
8 S S 8 8 2
! 1
"3
^
A
"S 9
- o> -o
»— . „
- 1
— ' 1
OO
S
•»-N
S— ^
<4
• v4
g,
i^
_ o
S «
_cn •-
_ ft-
t^
O
— *J
v>
I
- O
__ O
ta
UJ OT
i- ^
< 7
Q -o
- \
- en .S
00 UJ
"en 11
_ »^ **-s
U
O
(4
1
S
• V*
ptf
^<
ft-
>>
^
"O
u
:. I
00 3
- OT J2
_ *~ ^
- &
IM
§
n.
0 2
S
.
/4>'
U
.S
"3
1/5
8
'i
2
<«->
o
s
1
(A
5
,_-
ra
•v«
taJl
a*
>'
i
35
-------
where MIN is the minimum monthly irradiance. Both derivations allowed for intra-annual and
interannual variability in daily solar irradiance that is characteristic of lower bay climatic conditions.
Table 3 gives the descriptive statistics for the eight-year VIMS data base and for one-year
simulations using RANDOM-1 and RANDOM-2. RANDOM-1 mean monthly PAR predictions
compared well with the observed means but the standard deviations and coefficients of variation
were significantly underestimated. RANDOM-2 underestimated the observed mean monthly PAR
during parts of the year but was a much better predictor of in situ variability overall. Figure 14
gives the observed daily integrated PAR compared to the smooth model used in version 2 and the
stochastic model (RANDOM-2) used in version 4. It is apparent that the stochastic model captures
the variability and annual pattern iitsolar PAR very well. However, changing the solar PAR model
did not greatly affect the simulated dynamics of eelgrass under otherwise nominal conditions, i.e.,
with fixed attenuation and water depth (Wetzel, unpublished data).
Table 3. Comparison of solar irradiance (Ein m ! d l) observed at Gloucester Point,
Virginia (1984-1992) and simulated for a single model year using two
stochastic functions (see text). X = mean; SD = standard deviation; CV =
coefficient of variation.
Month
Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec
Observed
X SD CV
14.8 7.4 50.0
19.7 9.7 49.2
27.7 13.5 " 48.9
38.1 14.9 39.1
44.4 11.8 26.6
43.2 12.8 29.6
43.0 11.2 26.0
39.4 9.4 23.9
33.2 10.9 32.8
23.6 8.9 37.7
18.2 5.8 31.9
13.8 5.8 42.0
Stochastic on SD
X SD CV
14.0 3.8 27.1
19.0 6.0 31.6
29.9 8.0 26.8
41.4 7.8 18.8
44.8 6.8 15.2
42.3 5.6 13.2
42.1 5.7 13.5
38.6 5.7 14.8
34.6 5.9 17.1
24.5 4.9 20.0
' 18.8 3.3 17.6
14.1 2.4 17.0
Stochastic on range
X SD CV
12.8 6.4 50.0
18.4 9.5 51.6
25.7 12.0 46.7
29.9 13.4 44.8
38.6 10.8 28.0
35.4 13.6 38.4
33.9 14.1 41.6
29.7 11.5 38.7
26.0 11.8 45.4
22.6 8.0 35.4
14.0 6.1 43.6
12.8 5.9 46.1
36
-------
Daily Solar Irradiance
Comparison of measured and modeled PAR used in SAV simulations
T3
CVI
'E
•^
LLI
CC
0)
2
'o
10 -
92
274
365
Day
Figure 15. A comparison of daily surface irradiance from a long-term data base at Gloucester Point,
Virginia (open circles) with modeled irradiance. The solid line shows the smoothly varying
function, PAR = 28.25 - 16.75 cos(2rcDAY/365), fitted from the observations. The dashed line
shows a stochastically varying function derived from the daily integrated minimum and range from
the same data base.
37
-------
PAR Attenualion
The second component determining submarine PAR intensity in model version 4, and
perhaps still the least predictable, is the downwelling PAR attenuation coefficient, Kd. S AV model
version 2 treated Kd as a constant for a given simulation. However, these earlier simulation results
indicated that the depth distribution and long term stability of S AV were sensitive to changes in Kd.
Attenuation coefficients of ca. 2.0 irf1 and depths of 2.0 m were the maximum values allowable for
long term eelgrass survival in the simulations. Data accumulated from field studies and the shoal
monitoring program indicated that these values were higher than maxima for Kd values and water
column depths observed in extant eelgrass beds. We proposed that these overestimates could be
explained by lack of'simulated variability in Kd particularly when coupled with solar PAR
variability. That is, inclusion of higher frequency fluctuations in the modeled light environment
(both incident and submarine) would result in a greater sensitivity of modeled eelgrass productivity
to parameters that determine the in situ light environment.
The impact of time-varying Kd was investigated in two ways. First, temporal averages based
on the temperature-delineated growth phases of eelgrass in the lower Chesapeake Bay (Moore 1992)
were constructed for three sites in the Shoal Monitoring Program in the York River, Virginia.
Guinea Marsh, Gloucester Point, and Claybank represent the spatial continuum of SAV habitat
quality, from healthy extant eelgrass meadows of the first site to a presently denuded site upriver
at the third site. The seasons are defined as follows; see Moore 1992 for a detailed description of
this delineation:
Winter 13°C -> 0°C -> 9°C
Spring: 9°C -» 23°C
Summer: 23°C -» 30°C -» 25°C
Fall: 25°C -» 13°C
The observed seasonal variability within and between these sites Js evident in Table 4. In
comparison with the earlier simulations (version 2) where Kd was fixed, it can be seen that the
Guinea Marsh site has Kd values similar to the SAV model's nominal case value of 1.0 m"1 during
the spring and fall eelgrass growth periods, while Claybank has Kd's higher than the predicted
maximum value for a sustainable eelgrass community (Batiuk et al. 1992). Seasonal attenuation
coefficients at the Gloucester Point site are marginally high for predicted eelgrass survival.
The impact of these differences in site-specific and temporally varying Kd's is evident in the
simulated annual cycle of eelgrass shoot biomass (Figure 15). Simulated eelgrass populations
exhibited long term stability under variable Kd forcing for Guinea Marsh and Gloucester Point site-
specific data. Despite the introduced variability, the Guinea Marsh case supported a biomass within
90% of that sustained using the fixed conditions and smoothly varying PARD of the nominal case.
Similarly, the Gloucester Point site simulation exhibited somewhat reduced, but stable eelgrass
population density. Using seasonal patterns at Claybank, however, a simulated eelgrass population
rapidly (over three years) declined. Of course, presently there is no eelgrass population at Claybank.
Temporal variability in Kd was further enhanced in another set of simulations, whereby
stochastic variations were introduced on a daily basis to the seasonally averaged data (see Table 4).
The method for introducing this variability was identical to that used for introducing stochastic
38
-------
(4
I
•a
o
~*—' CM
•SS
IS.
60 s» O
v-i 4> -rt
!T 3 ta
> 13 '?
^ > 5
»> -O "°
Sll
•a 1"8
o 2 *
JH t/3 "S
5 BZ
« «'x
c3 .X •s.*-
CM
1
U
•S
Q
oo
£
a
o
CS
•
2 8 "»
S «.g
M
,0
§ 60 "~
so C .0
a Hi
0>
1
H
o
cs
o
g
S
V)
ON
«—<
00
I
CO
39
-------
Polyhaline SAV Model
Using seasonally-varying Kd for three sites (York River, Va.)
150
365
I ' '
730 1095
Model Day
1460
1825
Claybank
Gloucester Pt.
- Guinea Marsh
Figure 16. Predicted eelgrass leaf biomass using seasonally varying light attenuation coefficients
from three sites in the York River, Virginia. In all cases, depth was fixed at 1.0 m.
40
-------
variations in the incident PAR data, as described above. Using the Guinea Marsh case as an
example, daily random variations about seasonal means imposed little further limitation on
simulated eelgrass populations; biomass levels were similar under the imposition of both frequencies
of variation in Kd (compare the Guinea Marsh case in Figure 15 with the 1.0 m case in Figure 16).
After examination of the effects of variable versus annual mean Kd on the prediction of
eelgrass population survival, it is useful to ask, given these site-specific and temporally varying
results, can these simulations reproduce the depth-restricted distributions presently observed, and
what will the predicted depth limits be using the restoration criteria for Kd established for SAV
habitats (see Batiuk et al. 1992)? Examining the impact of variable Kd for different isobaths
provides a test for these seasonally defined restoration criteria. Such tests were performed using
Kd's, including the addition of daily stochastic variations, characteristic of those at the extant Guinea
Marsh eelgrass beds (see Table 4). Predicted stable populations were maintained for plant canopies
1.5 m from the surface (Figure 16). Assuming a canopy height of ca. 0.30 m, this corresponds to
isobaths of ca. 1.80 m, somewhat shallower than the Tier IE target (2.0 m; see Batiuk et al. 1992)
for restoration. That is, even though this site has water clarity that is within the restoration criteria,
the model does not predict that stable populations can be maintained in waters 2.0 m deep.
Presently, eelgrass populations at this site, while dense and apparently healthy, remain restricted to
depths <1.5 m (MLW) (e.g., Orth and Moore 1988; Orth et al. 1992, 1993).
The final issue addressed by SAV modeling during this work period dealt with tidal
variations in depth as a third source of high frequency variability associated with eelgrass habitat.
In conjunction with the light attenuation coefficient, this is another factor influencing the light
available at the leaf surfaces of eelgrass. To examine the potential impacts, relative to the previous
model cases with constant depth, a series of simulations were run using the conditions at Guinea
Marsh, which, of the three sites used for simulation analysis, has the most productive extant eelgrass
population. Each case was run using a stochastically varying Kd based on seasonally grouped
statistics as described above (Table 4).
In the face of stochastically varying PAR and Kd and tidally varying depth, eelgrass
populations were not stable for plant canopies deeper than 1.0 m below the surface (Figure 17). This
is 0.5 m shallower than in simulations lacking tidal variability (Figure 16), a significant difference
with respect to both the restoration targets and the potential area! extent of SAV habitat. Again,
assuming a ca. 0.30-m canopy height for eelgrass populations in the lower Chesapeake Bay,
conditions in this model scenario would meet the Tier n target for maintaining or restoring eelgrass
to a depth of 1.0 m, but not the Tier HI target of 2.0 m. Thus, for the same water quality conditions,
the addition of tidal variability to the model led to a greater light sensitivity for eelgrass.
Furthermore, these scenarios suggest that the adopted habitat restoration criterion for Kd may not
be sufficient to achieve the Tier in goal.
Key Simulation Findings
The results of several phases of polyhaline SAV modeling, focusing on eelgrass as the
dominant species, have provided strong support for the efficacy of such models in addressing
management oriented questions regarding water quality and habitat restoration. While future work
will enhance and integrate SAV models into larger ecosystem models capable of addressing a
broader scope of living resource questions (e.g., relating water quality and habitat restoration goals
41
-------
Polyhaline SAV Model
Using stochastic PAR and Kd with fixed depth
365
730 1095
Model Day
1460
1.0m
1.25 m 1.5 m
•-- 1.75m
1825
2.0m
Figure 17. Predicted eelgrass leaf biomass with increasing water depth (mean low water) using
stochastically varying light (PAR) and light attenuation coefficients (Kd), the latter derived from
data at Guinea Marsh, York River, Virginia.
42
-------
Polyhaline SAV Model
Using stochastic PAR and Kd with tidally-varying depth
150
365
730 1095
Model Day
1460
1825
•
0.5 rn
f\ ~7C r-ri
1A n-i
1 .25 m
1 .5 m
Figure 18. Predicted eelgrass leaf biomass with tidal fluctuations at five mean low water depths
using stochastically varying light (PAR) and light attenuation coefficients (Kd), the latter derived
from data at Guinea Marsh, York River, Virginia.
43
-------
to the productivity of living resources), several key findings have arisen from the standalone SAV
models presented here. These are summarized here:
1. In a nominal model case, with depth fixed at 1.0 m and Kd fixed at 1.0 m"1, both
indicative of conditions within healthy, extant eelgrass meadows, simulations yielded
biomasses and the annual cycle of eelgrass production similar to that observed in situ.
Further cases with Kd values reflective of degraded habitats showed decline and failure
of modeled eelgrass populations, also in accord with field observations. These results
indicate that the extent of eelgrass populations in lower Chesapeake Bay is primarily
controlled by light available for photosynthesis; nutrients do not appear to be limiting
production.
2. While the nominal model captured the sensitivity of eelgrass to light attenuation, the
results were less sensitive than is suggested by observations and the distribution patterns
of eelgrass in the lower Chesapeake Bay.
3. The inclusion of daily scales of variability for incident light and for attenuation, based
on long term field data in the region and its various habitats, increased slightly the
sensitivity of the model predictions to light attenuation. However, the addition of tidally
varying depth greatly increased the sensitive of modeled eelgrass biomass to light
attenuation. While this may be non-intuitive, it should be recalled that light is attenuated
exponentially with depth and that eelgrass populations, while light-limited for much of
the year and experience severe thermal stress during the summer, such that the additional
stress of low light when high tide is coincident with periods of daily photosynthetic
maxima, eelgrass production is sufficiently depressed to lead to loss in biomass. This
sensitivity, detected in the modeling effort, may indicate that habitat restoration targets
with respect to water clarity may not be sufficient to allow the restoration of eelgrass
beds to the 2.0 m Tier ffl target.
4. Modeling the interactions of eutrophication on epiphyte density, which in turn affects
the photosynthesis of SAV, and the control of these epiphytes by grazing activity,
indicates that, at present, grazers are sufficiently cropping epiphytes within, though the
impact of eutrophication on phytoplankton and its light attenuation was not explicitly
addressed. However, any perturbations which reduce the grazer populations or increase
algal (planktonic and epiphytic) growth rates will have a negative impact on SAV by
allowing phytoplankton and epiphytes to accumulate to levels which would severely
reduce SAV photosynthesis.
FUTURE WORK
In our simulation analyses, in corroboration of continuing studies of SAV communities and
water quality, the nature of the highly variable light environment is a critical factor in SAV
productivity. Often, even in flourishing SAV communities, light criteria can be near critical values
for SAV productivity (Moore 1992). The submarine light climate is determined by depth, suspended
particulates (including phytoplankton, non-living organic detritus, and inorganic sediments,
44 . .
-------
measured collectively as TSS), and dissolved organic matter (DOM). In terms of anthropogenic
impacts, alteration in the TSS from upstream perturbations or increased nutrient loadings which lead
to increased phytoplankton and epiphytic growth will reduce the light available to SAV.
Our previous simulation analyses have examined the relationship of light availability and
SAV productivity by incorporating, directly, data on incident irradiance and downwelling
attenuation into simulation studies. Obviously, this approach is largely dependent site-specific data
and thus is limited in its ability to explain or predict features in other area of the Chesapeake Bay
and its tributaries. Also, the data-driven approach limits our ability to use the model to examine
complex relations between water column (microalgal) and benthic (microalgal and macrophytic)
productivity. Competition for light and nutrients between water column and benthic autotrophs
determines the nature of the food web for the community, whether it is based in the benthos or in
the pelagial.
In our future work, we will extend our simulation efforts in three main directions. The first
direction will be to extend the theoretical basis of our light description. This will be done by
calculating incoming solar irradiance and downwelling attenuation from a more theoretical basis
(e.g., Kirk 1983). The total diffuse downwelling attenuation coefficient can be described as the
linear sum of individual coefficients for water, detritus, inorganic suspended matter, phytoplankton
pigments, and dissolved organic matter. Particularly for inorganic sediments and phytoplankton,
this would allow for more complex examinations of interactions of these readily-measured
parameters, light availability and ecosystem processes in relation to water quality and living
resources restoration goals.
The second area of fruitful expansion of our modeling efforts will be the inclusion of nutrient
cycling, sediment oxygen and nutrient fluxes, and water column trophic dynamics (including
phytoplankton, grazing zooplankton, and predators). Of course, a description of phytoplankton
biomass will be necessary in the aforementioned description of light attenuation. Additionally,
inclusion of this aspect of the SAV-littoral ecosystem will allow for complex examinations of
eutrophication, specifically, the cycles of nitrogen and phosphorus. Also included will be the
dynamics of dissolved oxygen, which is intimately couple in cycles of organic production,
degradation, and in the fluxes between the sediment and water column.
Thirdly, we are current in the initial stages of integrating a SAV ecosystem process model
with the water quality and hydrodynamic modeling efforts of colleagues at VIMS. This was
discussed in the Background section, and is focused on the York River in a program called the York
River Regional Ecosystem project. The coupling of physical and biological processes two (x, y,
depth-integrated) and three dimensions (x, y, z) improves the ability to describe and predict the mass
transport of materials important in determining water quality (nitrogen, phosphorus, suspended
matter, oxygen), as well as examining the spatial scales associated with biological processes
important for water quality and habitat restoration. It can also indicate potentially sensitive
locations, based on water flow and sediment suspension characteristics, within the tributaries where
improved water quality is crucial for the re-establishment of SAV in ways in which spatially
averaged models cannot predict.
While the last of these improvements is a long term effort, because of the mathematically,
theoretical and computational complexities, all three of these areas of improvement will have
45
-------
significant impacts on our understanding of the complex biological interactions within shallow water
aquatic SAV habitats. In concert with long term survey data and with process-oriented field and
laboratory experiments, the modeling efforts can aid the identification of critical and sensitive
aspects of the ecosystem (Figure 4). In addition, the coupling of complex biological processes with
water quality and hydrodynamic models should refine and improve predictive and management
capabilities.
46
-------
REFERENCES
Adams, S.M. 1976a. The ecology of eelgrass, Zostera marina (L.) fish communities. I. Structural
analysis. J. exp. mar. Biol. Ecol. 22: 269-291.
Adams, S.M. 1976a. The ecology of eelgrass, Zostera marina (L.) fish communities. II. Functional
analysis. J. exp. mar. Biol. Ecol. 22: 293-311.
Adams, S.M. 1976c. Feeding ecology of eelgrass fish communities. Trans. Am. Fish. Soc. 4:514-
519.
Batiuk, R.A. and others. 1992. Chesapeake Bay submerged aquatic vegetation habitat requirements
and restoration targets: a technical synthesis. U.S. EPA, Chesapeake Bay Program contract no. 68-
WO-0043, Annapolis, MD
Bender, M. E. 1987. The York River: A Brief Review of its Physical, Chemical and Biological
Characteristics. Report to American Petroleum Institute, Health and Environmental Sciences Dept.,
Washington, D.C.
Biebl, R. and C.P. McRoy. 1971. Plasmatic resistance and rate of respiration and photosynthesis of
Zostera marina at different salinities and temperatures. Mar. Biol. 8: 48-56.
Borum, J. 1985. Development of epiphytic communities on eelgrass (Zostera marina) along a
salinity gradient in a Danish estuary. Mar. Biol. 87: 211-218.
Christian, R. R. and R. L. Wetzel. 1978. Interactions between substrate, microbes and consumers
of Spartina detritus in estuaries. In: Wiley, M. (ed). Estuarine Interactions. Academic Press, New
York, pp 93-114.
Christian, R. R. and R. L. Wetzel. 1991. Synergism between research and simulation models of
estuarine microbial food webs. Microb. Ecol. 21:111-125.
Cottam, C. 1935a. Further notes on past periods of eelgrass scarcity. Rhodora 37:269-271.
Cottam, C. 1935b. Wasting disease of Zostera marina. Nature 135:306.
Cottam, C. and D. A. Munro. 1954. Eelgrass status and environmental relations. J. Wildl. Mgt. 18:
449-460.
Diaz, R.J. and T. Fredette. 1982. Secondary production of some dominant macroinvertebrate species
inhabiting a bed of submerged vegetation in the lower Chesapeake Bay. In: R. J. Orth and J. van
Montfrans (eds.), Structural and Functional Aspects of the Biology of Submerged Aquatic
Macrophyte Communities in the Lower Chesapeake Bay, Vol. ffl. Interactions of Resident
Consumers in a Temperate Estuarine Seagrass Community: Vaucluse Shores, Virginia. SRAMSOE
267, Virginia Inst. Mar. Sci., Gloucester Point, VA, pp. 95-123.
47
-------
Evans, A.S. 1984. Temperature adaptation in seagrasses. M.A. Thesis, College of William and
Mary, Williamsburg, VA 76 pp.
Gilchrist, W. 1984. Statistical Modelling. John Wiley & Sons, New York. 339 pp.
Gold, H. J. 1977. Mathematical Modeling of Biological Systems; An Introductory Guidebook. John
Wiley & Sons, New York. 357 pp.
Heinle, D. R. and others. 1980. Historical Review of Water Quality and Climatic Data from
Chesapeake Bay with Emphasis on Effects of Enrichment. Chesapeake Research Consortium, Inc.,
Publication no. 84, Solomons, MIX
High Performance Systems, Inc. 1992. STELLA n An Introduction to Systems Thinking. Hanover,
NH. 176 pp.
Hoss, D.E. 1974. Energy requirements of a population of pinfish Lagodon rhomboides (Linnaeus).
Ecology 55: 848-855.
Kirk, J.T.0.1983. Light and photosynthesis in aquatic ecosystems. Cambridge Univ. Press, New
York. 399pp.
Marsh, G.A. 1973. The Zostera epifaunal community in York River, Virginia. Ches. Sci. 14:87-97.
McRoy, C.P. 1974. Seagrass productivity: carbon uptake experiments in eelgrass, Zostera marina.
Aquaculture 4:131-137.
Moore, K.A. 1992 Regional SAV study area finding—York River, in Chesapeake Bay submerged
aquatic vegetation habitat requirements and restoration targets: a technical synthesis. U.S. EPA,
Chesapeake Bay Program contract no. 68-WO-0043, Annapolis, MD
Murray, L. 1983. Metabolic and structural studies of several temperate seagrass communities, with
emphasis on microalgal components. Ph.D. Dissertation, College of William and Mary,
Williamsburg, VA 90 pp.
Murray L. and R.L. Wetzel. 1987. Oxygen production and consumption associated with the major
autotrophic components in two temperate seagrass communities. Mar. Ecol. Prog. Ser. 38:231-239.
Neckles, H.A. 1990. Relative effects of nutrient enrichment and grazing on epiphyton-macrophyte
(Zostera marina L.) dynamics. Ph.D. Dissertation, College of William & Mary, School of Marine
Science, Gloucester Point, VA.
Neckles, H.A., R.L. Wetzel and R.J. Orth. 1993. Relative effects of nutrient enrichment and grazing
on epiphyte-macrophyte (Zostera marina L.) dynamics. Oecologia 93: 285-295.
Nixon, S. W. and C. A. Oviatt. 1972. Preliminary measurements of midsummer metabolism in beds
of eelgrass, Zostera marina. Ecology 53: 150-153.*
48
-------
Odum, H.T. 1971. Environment, Power and Society. J. Wiley, New York. 331pp.
Orth, R.J. 1977. Effect of nutrient enrichment on growth of the eelgrass, Zostera marina in the
Chesapeake Bay, Virginia, U.S.A. Mar. Biol. 44: 187-194.
Orth, RJ. and K.L. Heck, Jr. 1980. Structural components of eelgrass (Zostera marina) meadows
in the lower Chesapeake Bay—fishes. Estuaries 3: 278-288.
Orth, R.J. and K.A. Moore. 1983. Chesapeake Bay: An unprecedented decline in submerged aquatic
vegetation. Science 222: 51-53.
Orth, RJ. and K.A. Moore. 1988. Distribution of Zostera marina L. and Ruppia maritima L. sensu
lato along depth gradients in the lower Chesapeake Bay, U.S.A. Aquat. Bot. 32: 291-305.
Orth, R.J., J.F. Nowak, G.F. Anderson, K.P. Kiley, and J.R. Whiting. 1992. Distribution of
submerged aquatic vegetation in the Chesapeake Bay and tributaries and Chincoteague Bay—1991.
A Final Report submitted to the U.S. E.P.A., Chesapeake Bay Program Office, Annapolis, MD. 268
pp.
Orth, R.J., J.F. Nowak, G.F. Anderson, and J.R. Whiting. 1993. Distribution of submerged aquatic
vegetation in the Chesapeake Bay and tributaries and Chincoteague Bay—1992. A Final Report
submitted to the U.S. E.P.A., Chesapeake Bay Program Office, Annapolis, MD. 268 pp.
Penhale, P. 1977. Macrophyte-epiphyte biomass and productivity in an eelgrass (Zostera marina)
community. J. Exp. Mar. Bio. Ecol. 26: 211-224.
Peters, D.S., M.T. Boyd and J.C. DeVane, Jr. 1976. The effects of temperature, salinity and food
availability on the growth and food-conversion efficiency of postlarval pinfish. In: G. W. Esch and
R. W. MacFarlane (eds.) Thermal Ecology n. National Technical Information Service, CONF-
750425, Springfield, VI, pp. 106-112.
Twilley, R.R., W.M. Kemp, K.W. Staver, J.C. Stevenson and W.R. Boynton. 1985. Nutrient
enrichment of estuarine submersed vascular plant communities. 1. Algal growth and effects on
production of plants and associated communities. Mar. Ecol. Prog. Ser. 23: 179-191
van Montfrans, J., R.L. Wetzel and RJ. Orth. 1984. Epiphyte-grazer relationships in seagrass
meadows: consequences for seagrass growth and production. Estuaries 7: 289-309.
Vaughan, D.E. 1982. Production ecology of eelgrass (Zostera marina L.) in Little Egg Harbor, New
Jersey. Ph.D. Dissertation, Rutgers University, New Brunswick, New Jersey, 126 pp.
Wetzel, R.L. and R.R. Christian. 1984. Model studies on the interactions among carbon substrates,
bacteria and consumers in a salt marsh estuary. Bull. Mar. Sci. 35:601-614.
Wetzel, R. L. and C. S. Hopkinson, Jr. 1990. Coastal ecosystem models and the Chesapeake Bay
Program: philosophy, background and status. In: Haire, M. and E. C. Krome (eds.). Perspectives on
49
-------
the Chesapeake Bay, 1990, Advances in Estuarine Sciences. CBP/TRS41/90, Chesapeake Research
Consortium, Inc., Solomons, MD. pp 7-23.
Wetzel, R. L. and H. A. Neckles. 1986. A model of Zostera marina L. photosynthesis and growth:
simulated effects of selected physical-chemical variables and biological interactions. Aquatic Bot.
26: 307-323.
Wetzel, R.L. and P.A. Penhale. 1983. Production ecology of seagrass communities in the lower
Chesapeake Bay. Mar. Tech. Soc. J. 17: 22-31.
Wetzel, R.L. and R.G. Wiegert. 1983. Ecosystem simulation models: tools for the investigation of
nitrogen dynamics coastal and marine ecosystems. In: E. J. Carpenter and D. G. Capone (eds.)
Nitrogen in the Marine Environment. Academic Press, New York, pp. 869-892.
Wiegert, R.G. 1973. A general ecological model and its use in simulating algal-fly energetics in a
thermal spring community. In: P. W. Geier, L. R. Clark, D. J. Anderson and H. A. Nix (eds.) Studies
in Population Management. Vol. 1. Occasional Papers, Ecol. Soc. Australia, Canberra, pp. 85-102.
Wiegert, R.G. and R.L. Wetzel. 1974. The effect of numerical integration technique on the
simulation or carbon flow in a Georgia salt marsh. Proc. Summer Computer Simulation Conf., Vol.
2, Houston, TX, pp. 275-277.
Wiegert, R. G. and R. L. Wetzel. 1979. Simulation experiments with a 14-compartment salt marsh
model. In: Dame, R.(ed). Marsh-Estuarine Systems Simulation, Univ. South Carolina Press,
Columbia, pp. 7-39.
Zimmerman, R., R. Gibson and J. Harrington. 1979. Herbivory and detritivory among gammaridean
amphipods from a Florida seagrass community. Mar. Biol. 54:41-47.
50
-------
APPENDIX
PROGRAM SEAGRASS
C
C A SIMULATION MODEL OF ZOSTERA MARINA PHOTOSYN-
C THESIS, RESPIRATION, MORTALITY AND LOSS DUE TO
C GRAZING; EPIPHYTE COLONIZATION, PHOTOSYNTHESIS,
C RESPIRATION, MORTALITY AND LOSS DUE TO GRAZING
C AND THE ROLE OF EPIPHYTIC GRAZERS (MODELED AFTER
C THE ISOPOD, IDOTEA BALTICA) . THE MODEL IS DRIVEN
C BY LIGHT, WATER TEMPERATURE, PAR ATTENUATION AND
C WATER DEPTH (A FUNCTION OF A PROGRESSIVE TIDAL
C CYCLE WITH A LUNAR PERIOD OF 30 DAYS. THE MODEL RUNS
C ON THREE TIME SCALES: HOURS, DAYS AND YEARS.
Q**********************************************************A*******************
C DEFINE PROGRAM ENVIRONMENT (ASSIGN NAMES TO DATA STRUCTURES, ALLOCATE
C STORAGE SPACE) AND FORMATS FOR MODEL OUTPUT
(2**** **************************************************************************
IMPLICIT DOUBLE PRECISION (A-H.O-Z)
IMPLICIT INTEGER (I-N)
DOUBLE PRECISION MLW, MSL
CHARACTER*! PARFCT, TIDE
CHARACTER* 3 2 FNAME
CHARACTER RUNHDR*79, SITE*2
C$LARGE : FLUX, AFLUX, PHYSICAL, TRATES IMS-Fortran-specific compiler instruction
C
DIMENSION FUMMY(16) ,X(7) ,DX(7) ,FB(12,500) ,SFB(12)
DIMENSION FLUX(16, 500) ,SIFLUX(16) ,STFLUX(16) ,DFLUX(16) ,TFLUX(16)
DIMENSION AFLUX ( 16 , 500 )
DIMENSION DAFLUXO)
DIMENSION DARATE(2)
DIMENSION PHYSICAL (7, 500) , PHYSDAT (7)
DIMENSION TRATES ( 14 , 500 ), SRATES ( 14 )
EQUIVALENCE (FUMMY(l) , A0203) , (X(l) ,X01) , (SFB(l) , FB0203)
EQUIVALENCE (SIFLUX(l) ,F0203)
EQUIVALENCE ( PHYSDAT ( 1 ) , TEMP )
EQUIVALENCE ( SRATES ( 1 ) , EPLR)
COMMON/FIXED/A0203,G0203,A0303,G0303,A0404,G0404,A0205,G0205,
+A0505,G0505,A0306,G0306,A0606,G0606,A0506,G0506
C
COMMON/ STATE /X01 , X02 , X03 , X04 , X05 , X06 , X07 , DT , NDPRT, JYRMX
C
COMMON/TIME/ ITIME( 500)
C
COMMON/FLUX/F0203 , F0302 , F0304 , F0402 , F0307 , F0205 , F0502 , F0507 ,
+F0306,F0506,F0602,F0607,F0600,F0006,F0102,F0201
C
COMMON/FB/FB0203 , FB0303 , FB0404 , TF0203 , FB0205 , FB0505 , TF0205 ,
+FB0306,TF0306,FB0506,TF0506,FB0606
C
COMMON/PHYSICAL/TEMP, Z, PP.PARD, PARSEC , PARZ , AKZ
C
COMMON /RATES /EPLR, PARVP, VPIK, PM03 , T0203 , R0302 , EPIK, PM05 ,
+T0205,R0502,DIFFL,T0306,T0506,R0602
write(6. '
WRITE(6,*)' Welcome to the Seagrass Ecosystem...'
write(G,*)' PAR, T, Kd, Z — > Zostera, Epiphyta, Grazer'
write(6,*)' — > C02, Fish'
write(6, ' (///) ')
write(6,1000)
read(5, ' (a79) ') RUNHDR
1000 formate ' , 'Give a model run descriptive header (<80 chars)...'
& , ' ', '>')
c WRITE(6,1001)
C 1001 FORMAT (5X, 'ENTER FILENAME FOR DAILY INTEGRATED FLUX OUTPUT ' \)
C READ (5, 1002) FNAME
51
-------
C 1002 FORMAT(A)
NFLX = 20
FNAME='DFLX.PRN'
OPEN(NFLX,FILE=FNAME,STATUS='NEW' )
C WRITE(6,1003)
C 1003 FORMAT(5X,'ENTER FILENAME FOR DAILY INTEGRATED RATE OUTPUT '\)
C READ(*,1002)FNAME
NRAT = 21
FNAME='DRAT.PRN'
OPEN(NRAT,FILE=FNAME,STATUS='NEW1 )
C WRITE(*,1004)
C 1004 FORMAT('5X, 'ENTER FILENAME FOR STANDING STOCK OUTPUT '\)
C READ(*,1002)FNAME
FNAME= ' SS. PRN'
NSST =22
OPEN(NSST,FILE=FNAME,STATUS='NEW1)
C WRITE(*,1005)
C 1005 FORMAT(5X,'ENTER FILENAME FOR PHYSICAL DATA OUTPUT '\)
C READ(*,1002)FNAME
FNAME='PHYS.PRN'
NPHY = 23
OPEN(NPHY,FILE=FNAME,STATUS='NEW')
C WRITE(*,1006)
C 1006 FORMAT(5X,'ENTER FILENAME FOR FLUXES AND MASS BALANCE OUTPUT '\)
C READ(*,1002JFNAME
FNAME='MB.PRN'
NMBL = 24
OPEN(NMBL,FILE=FNAME,STATUS='NEW')
C
C WRITE FILE HEADERS TO STD STOCK FILE
WRITE(NSST,*) RUNHDR
WRITE(NSST,200)
C
C TABULAR OUTPUT FORMATS FOR COMPARTMENT BIOMASS, DAILY AND
C SUMMED FLUXES AND FEEDBACK TERMS.
C
C OUTPUT FORMATS FOR STATE VARIABLES
200 FORMAT(3X, 'DAY1,5X, 'CO2-AIR',5X, 'C02-H20',5X, 'PLT-LVS',5X,
+'PLT-R/R',5X,'EPIPHY1,6X,'GRAZER',5X,'DETRITUS',5X,'EPI/PLT-LVS1)
210 FORMAT(IX,15,2X,7F12.4,5X,F7.4)
C
C OUTPUT FORMAT FOR DAILY INTEGRATED FLUXES
599 FORMAT(I5,9F10.4)
598 FORMAT(I5,2F10.5)
C
C OUTPUT FORMATS FOR DAILY AND SUMMED FLUXES
C
229 FORMAT(3X,'MODEL PREDICTIONS: SPECIFIC FLUX RATE ',
+'(GM-C / M2 / HOUR) AT 12:00 HOURS ON GIVEN DAY1)
230 FORMAT(3X,' DAY',2X,' F0203 F0302 F03041,
+' F0402 F0307 F0205 F0502 F0507•)
231 FORMAT(3X,' DAY',2X,• F0306 F0506 F0602•,
+ ' F0607 F0600 F0006 F0102 F0201')
232 FORMAT(3X,I5,2X,8F10.4)
233 FORMAT(3X,'MODEL PREDICTIONS: CUMULATIVE FLUX ',
-t-'BY SPECIFIC PATHWAY1)
C
C OUTPUT FORMATS FOR FEEDBACK TERMS
C
239 FORMAT(3X,'MODEL PREDICTIONS: FEEDBACK CONTROL ',
+1VALUES AT 12:00 HOURS FOR THE GIVEN DAY1)
240 FORMAT(3X,' DAY',2X,' FB0203 FB0303 FB04041,
+' TF0203 FB0205 FB05051)
242 FORMAT(3X,' DAY',2X,' TF0205 FB0306 TF03061,
+ 1 FB0506 TF0506 FB06061)
241 FORMAT(3X,I5,2X,6F10.4)
C
C OUTPUT FORMATS FOR COMPARTMENT MASS BALANCE
C
249 FORMAT(10X,'MODEL PREDICTIONS: ANNUAL COMPART1,
•»-'MENTAL MASS BALANCE')
250 FORMAT(1OX,'COMPARTMENT',1OX,'SUM INPUT1,10X,'SUM OUTPUT',
+10X,'MASS BALANCE1)
52
-------
251 FORMATdOX. 'X01: C02-AIR ' , 6X,F10 .4.10X,F10 . 4,10X.F10 .4)
252 FORMATdOX,'X02: C02-H20 ' ,6X.F10.4,10X,F10.4,10X,F10.4)
253 FORMATdOX, 'XO3: PLANT-LVS1 , 6X,F10 .4, 10X,F10 . 4, 10X,F10 .4)
254 FORMATdOX, 'X04: PLANT-R/R', 6X,F10 .4, 10X, F10 .4,10X,F10 .4) •
255 FORMATdOX, 'X05: EPIPHYTES ' , 6X.F10.4, 10X, F10 . 4,10X,F10. 4)
256 FORMATdOX,'X06: GRAZERS ' , 6X.F10 .4,10X, F10 . 4,10X, F10 . 4)
257 FORMATdOX, 'X07: DETRITUS • , 6X, F10 .4,10X,F10 . 4,10X, F10 .4)
C
C OUTPUT FORMATS FOR PHYSICAL DATA (ANNUAL CYCLES)
C
299 FORMAT(3X, 'MODEL PREDICTIONS: PHYSICAL DATA USED ',
+ 'FOR FORCING FUNCTIONS')
300 FORMAT(3X, ' DAY',IX, ' TEMP ',' DEPTH ',
+ ' DAYLT ',' PAR/DAY ',' PAR-INC ',
+' PAR-SUB ',' ATT-COEF ')
301 FORMAT(9X,' (C) ',' (M) ',' (HRS)
+ ' (E/M2/DY)',' (E/M2/S) ',' (UE/M2/S) ' , ' (M-2) ')
302 FORMAT(3X,I5,7F10.3)
C
C OUTPUT FORMATS FOR CALCULATED RATE COEFFICIENTS
C UNITS ARE MG-C/MG-C/HOUR
C
399 FORMAT(3X, 'MODEL PREDICTIONS: PHOTOSYNTHESIS ',
+ ' PARAMETERS AND RATE COEFFICIENTS FOR ZOSTERA MARINA AND ' ,
+'EPIPHYTES')
400 FORMAT(3X,' DAY1,' ZOSTERA ',' ZOSTERA ',
+ ' ZOSTERA ',' ZOSTERA ',' ZOSTERA ',' ZOSTERA ' , 5X,
-»•' EPIPHYTE ' , ' EPIPHYTE ' , ' EPIPHYTE ' , ' EPIPHYTE ' )
401 FORMAT(8X, ' EPLR ',' PARVP ',' VPIK
+' P-MAX ', ' GPP ' ,' RESP ',5X, ' EPIK
+' P-MAX ', ' GPP ',' RESP ' , ' DIFFL ')
402 FORMAT(3X,I5,6F10.5,5X,5F10.5)
403 FORMAT('MODEL PREDICTIONS: RATE COEFFICIENTS ',
+'FOR THE GRAZER COMPARTMENT')
404 FORMAT(3X,' DAY',' GRAZER ',' GRAZER ',
+' GRAZER ' )
405 FORMAT(8X,' ING:PLTS ',' ING:EPI ',' RESP ')
406 FORMAT(3X,I5,3F10.5)
C
Q******************************»***********************************************
C GENERAL MODEL STRUCTURE
C
C COMPARTMENTS:
C X01 C02-AIR
C X02 C02-WATER
C X03 ZOSTERA MARINA LEAVES
C X04 ZOSTERA MARINA ROOTS AND RHIZOMES
C X05 ALGAL EPIPHYTES
C X06 GRAZERS
C X07 DETRITUS
C
C FLOWS:
C FIJ=FLOW FROM COMPARTMENT I TO J; UNITS GRAMS CARBON/M2/H
C
C TRANSFER COEFFICIENTS:
C TIJ, UIJ, OR RIJ=SPECIFIC RATE OF CARGON TRANSFER FROM COMPARTMENT
C I TO J; UNITS GRAM/GRAM/HOUR
C »»NOTE THAT SOME TERMS REMAIN CONSTANT THROUGHOUT A GIVEN SIMULATION
C WHEREAS OTHERS ARE CALCULATED WITH EACH ITERATION. EXOGENOUS
C VARIABLES (E.G. TEMPERATURE, LIGHT) EXERT CONTROL BY AFFECTING
C THESE VARIABLE RATE COEFFICIENTS.
C
C FEEDBACK CONTROL:
C FBIJ=FEEDBACK FUNCTION WHICH REDUCES FLUX FROM COMPARTMENT I TO J WHEN
C CONCENTRATION OF RESOURCE I BECOMES LIMITING.
C FBJJ=DENSITY DEPENDENT FEEDBACK FUNCTION WHICH REDUCES FLUX TO
C COMPARTMENT J WHEN SPACE (NUTRIENTS, ETC.) BECOMES LIMITING;
C CONSTRAINS DENSITY OF COMPARTMENT J BELOW A MAXIMUM MAINTAINABLE LIMIT
C
£******************************************************************************
C
C INITIALIZE FIXED PARAMETERS
C
53
-------
PI = ACOS(-l.OO)
LCB = 0
LGP = 0
LGM = 0
C MBM: LAHEY F77L3 FCT TO SEED/PICK A RANDOM NO.
C (MIMICKS MS-FTN, WHICH USES AN INITIAL SEED OF 1.0)
RAND = RANDS(1.0)
WRITE(6,3005)
READ(5,'(F6.2)') AKD
AKZ = AKD
IF (AKD .EQ. 0.0) THEN
WRITE(6,4001)
READ(5,'(A2)') SITE .
IF (SITE .EQ. 'Cb1 .or. SITE .EQ. 'CB') LCB = 1
IF (SITE .EQ. 'gp' .or^SITE .EQ. 'GP') LGP = 1
IF (SITE .EQ. 'gm' .or. SITE .EQ. 'GM') LGM = 1
ENDIF
WRITE(6,3006)
WRITE(6,3007)
READ(5, '(F4.2) ') MLW
WRITE(6,3008)
READ(5,'(F4.2)') TAMP
WRITE(6,3009)
READ(5,'(Al)') TIDE
MSL = MLW + 0.5*TAMP
WRITE(6,3011)
READ(5,'(Al)') PARFCT
WRITE(6,3012)
READ(5,'(F5.2)') EUTROF
. WRITE(6,3013)
READ(5,'(F5.2)') GRAZF
IF (AKD .EQ. 0) WRITE(6,4010) SITE
WRITE(6,3010) AKD,MLW,TAMP,TIDE
WRITE(6,3015) PARFCT,EUTROF,GRAZF
3005 FORMAT(' ','Enter Kd (x.xx, 0. For computed values ) : ')
3006 FORMAT(' ','Enter nominal depth (m MLW), tidal range (m) and',
&' whether tidal variations',/
&, ' ',T5,'should be calculated (MSL = MLW -t- 0. 5*range) ... ')
3007 FORMATC ',T15,'Z (m MLW) : ')
3008 FORMATC ',T15,' Range (m) : ')
3009 FORMATC ',T15,'Tide (Y|N)_: ')
3010 FORMATC • , 'Kd = ',F6.2,' 1/m Z(mlw) = ',F4.2,' Range = ',F4.2
&,' m Tide calc''d= ',A1)
3011 FORMAT(5X,'Pick a (S)mooth or (N)oisy PAR signal : ')
3012 FORMAT(5X,'Eutrophication in this model impacts epiphyte growth.
&/,' ','Assign a factor (0.0-99.9, 1.0=nominal, xx.xx) : ')
3013 FORMAT(5X,'Grazing pressure impacts epiphyte biomass.',
&/,' ','Assign a factor (0.0-99.9, 1.0=nominal, xx.xx) : ')
3015 FORMATC ','PAR = ',A1,' EUTROF = ',F5.2,' GRAZF = ',F5.2)
4001 FORMATC ','Computed Kd''s are available for 3 sites in the',/
& ' ','York River, Virginia. Enter your choice using',/
& ' ','the initials as indicated:',//
& ' ','Claybank (cb) Gloucester Point (gp) '
& ,'Guinea Marsh (gm)',/
& ' ', ' : ')
4010 FORMAT(/,' ','Site = ',a2)
C
C ENVIRONMENTAL EXCHANGE PARAMETERS
-C
C > T0102=0.0
C > T0201=0.0
C > T0600=0.0
C > T0006=0.0
C
C ROOT/RHIZOME TRANSFER COEFFICIENTS
C
T0304=0.17
R0402=0.00042
C
C FIXED MORTALITY TERMS (FRACTION OF STANDING STOCK PER DAY)
54
-------
00507=0.0002083
U0607=0.000417
C
C GRAZER PREFERENCE VALUES AND ASSIMILATION EFFICIENCIES
C
P0306=0.05
P0506=0.95
AE0306=0.05
AE0506=0.518
C
C THE FOLLOWING CONTINUE STATEMENT IS THE PROGRAM REFERENCE
C POINT FOR MULITPLE YEAR SIMULATIONS.
C
1 CONTINUE
C
C SET YEARLY SIMULATION CONTROL VARIABLES (FLAGS)
C
DT=1.0
NDPRT=5
JYRMX=5
IYEAR=1
C
C WHERE
C DT IS THE ITERATION INTERVAL (DAYS)
C NDPRT IS THE PRINTING INTERVAL (DAYS)
C • JYRMX IS THE NUMBER OF YEARS SIMULATED
C IYEAR IS THE YEARLY COUNTER
C
C SET INITIAL CONDITIONS FOR STATE VARIABLES AND FIXED
C PARAMETER VALUES FOR FEEDBACK TERMS
C
X01=25.
X02=25.
X03=15.0
X04=7.5
X05=3.0
X06=0.01
X07=0.0
A0203=15.
G0203=5.
A0303=75.
G0303=150.
A0404=100.
G0404=125.
A0205=15.
00205=5.0
A0505=1.0
G0505=2.0
A0306=30.
00306=15.
A0606=0.5
G0606=1.0
A0506=0.25
G0506=0.1
C
C SET MAXIMUM STORAGE AND DERIVE THE NUMBER OF ITERATIONS PER DAY
C
ISTMAX=365
ITER=IFIX(24/DT+O.001)
C
C INITIALIZE COUNTERS AND LOOP CONTROLS
C
IEND=365
IDAY=0
LDAY=1
IPRT=1
ISTORE=1
IHR=12
JCDAY=1
C
C WHERE
C IEND IS THE TOTAL NUMBER OF DAYS SIMULATED
55
-------
C IDAY IS THE DAY COUNTER
C LDAY IS THE LUNAR DAY COUNTER
C IPLT IS THE COUNTER FOR PLOTING INTERVAL
C IPRT IS THE COUNTER FOR PRINTING INTERVAL
C ISTORE IS THE COUNTER FOR STORAGE USED
C IHR IS THE HOUR COUNTER FOR EACH DAILY ITERATION
C JCDAY IS THE JULIAN DAY COUNTER
C
C INITIALIZE FLUX STORAGES TO ZERO
C
DO 10 1=1,16
STFLUX(I)=0.0
10 CONTINUE
DO 11 1=1,16
SIFLUX(I)=0.0
11 CONTINUE
C
C INITIALIZE DAILY INTEGRATION STORAGE ARRAYS TO ZERO
C
DO 50 1=1,9
DAFLUX (I)=0.0
50 CONTINUE
DO 60 1=1,2
DARATE(I)=0.0
60 CONTINUE
C
C COMPUTE FEEDBACK DENOMINATOR TERMS FOR FIXED VALUES OF
C ALPHA AND GAMMA FOR BOTH DONOR AND RECIPIENT CONTROLS.
C
C NOTE: THEY ARE CALCULATED AS THE INVERSE
C NOTE: FB0505 (I.E. EPIPHYTE SPATIAL LIMITATION) IS A
C VARIABLE FEEDBACK DEPENDENT ON PLANT LEAF BIOMASS
C (IMPLICITLY, SURFACE AREA AVAILABLE FOR COLONIZA-
C TION) AND THE DATA ARE INPUT AS BIOMASS RATIOS.
C
C CALCULATE FEEDBACK DENOMINATOR TERMS THAT REMAIN CONSTANT
C THROUGHOUT A GIVEN RUN
C
00203=1./((A0203-G0203))+.1E-15
D0303=1./((G0303-A0303))+.1E-15
D0404=l./((G0404-A0404))+.lE-15
D0205=l./((A0205-G0205))+.1E-15
D0505=l./((G0505-A0505))+.1E-15
D0306=l./((A0306-G0306))+.1E-15
D0506=1./((A0506-G0506))+.1E-15
D0606=l./((G0606-A0606))-i-.lE-15
C
C PRINT OUT INITIAL CONDITIONS AND RESET IDAY COUNTER
C CALCULATE INITIAL EPIPHYTE:PLANT LVS RATIO
RA0503=(X05/(X03+0.1E-15))
WRITE(NSST,210)IDAY,X.RA0503
IDAY=1
C
C MODEL SIMULATION STARTS HERE
C
C THREE DO LOOPS CONTROL SIMULATION FOR TIME OF DAY, DAY
C OF YEAR AND NUMBER OF YEARS SIMULATED.
C
DO 120 JYEAR=1,JYRMX ' -
DO 110 JYDAY=1,IEND
C
C MS-Fortran random number generator:
C CALL RANDOM( RAND )
C MBM: LAHEY F77L3 FCT REF TO GET A PRE-SEEDED RANDOM NO. ONCE A DAY:
RAND = RNDO
C
£*********************************************************************
C***************** FITTED OR ASSUMED ***************
c***************** FUNCTIONAL RELATIONSHIPS ***************
c***************** THAT VARY DAILY (NOT W/IN A DAY) ***************
£*********************************************************************
C
C COMPUTE THE FOLLOWING ENVIRONMENTAL VARIABLES TO 'DRIVE1
56
-------
c
c
c
c
c
c
c
c
c
VARIOUS COMPARTMENTAL FLOWS
1. TEMP: WATER TEMPERATURE BY DAY (R; 0 TO 30 C)
2. PARD: DAILY INCIDENT PAR
(R: 11.5-45 EINSTEINS/M-2/CALENDAR DAY)
3. PP : DAILY PHOTOPERIOD IN HRS BY CALENDAR DAY
(R; 9.5 TO 14 HOURS)
6. AKZ : VERTICAL DOWNWELLING PAR ATTENUATION COEFFICIENT
IF (MLW .LE. 0.50) THEN
TEMP=17.75-(15.25*COS((2.*PI*(IDAY-25))/365.»
ELSE
TEMP=16.25-(13.75*COS((2.*PI*(IDAY-25))/365.))
ENDIF
PP=11.75-(2.25*COS((2.*PI*IDAY)/365.))
£******
£*****
£*****
£*****
c*****
£*****
c*****
Q* * * * *
£*****
Q*****
£*****
£*****
CHOICE OF PAR FCTS
1. SMOOTHLY VARYING OVER AN ANNUAL CYCLE
PARD ASSUMES A COSINE FUNCTION WHICH
PRODUCES A CONTINUOUS FUNCTION.
2. RANDOMLY NOISY DAILY VARIABLITY
PARD IS CALCULATED USING THE RANGE AND LOWER LIMIT OF
MEASURED PARD (VIMS DATA) AND A RANDOM NUMBER FUNCTION
THE GENERAL FORM OF THE EQUATION IS:
PARD = (RANDOM* * RANGE) + LOWER LIMIT
WHERE RANGE = 4*SD
******
******
******
******
******
******
******
******
******
******
******
******
IF (PARFCT .EQ. 'S1 .OR. PARFCT .EQ. 'S') THEN
PARD=28.25-(16.75*COS((2.*PI*(IDAY))/365.))
ELSE
IF (JCDAY.LE.31) THEN
C MONTH=1 !JANUARY
PARD = RAND*22. + 3.
ELSE IF (JCDAY.LE.59) THEN
C . MONTH=2 [FEBRUARY
PARD = RAND*31. + 3.
ELSE IF (JCDAY.LE.90) THEN
C MONTH=3 !MARCH
PARD = RAND*45. + 3.
ELSE IF (JCDAY.LE.120) THEN
C .MONTH=4 !APRIL
PARD = RAND*46. + 9.
ELSE IF (JCDAY.LE.151) THEN
C MONTH=5 !MAY
PARD = RAND*40. + 19.
ELSE IF (JCDAY.LE.181) THEN
C MONTH=6 1JUNE
PARD = RAND*47. + 13.
ELSE IF (JCDAY.LE.212) THEN
C MONTH=7 1JULY
PARD = RAND*48. + 9.
ELSE IF (JCDAY.LE.243) THEN
C MONTH=8 !AUGUST
PARD = RAND*39. -t- 11.
ELSE IF (JCDAY.LE.273) THEN
C MONTH=9 !SEPTEMBER
PARD = RAND*42. + 6.
ELSE IF (JCDAY.LE.304) THEN
C MONTH=10 '.OCTOBER
PARD = RAND*33. + 4.
ELSE IF (JCDAY.LE.334) THEN
C MONTH=11 !NOVEMBER
PARD = RAND*24. + 3.
ELSE
C MONTH=12 [DECEMBER
PARD = RAND*20. + 2.
END IF
END IF
£***********************************4
t************************
57
-------
C*********
C*********
C*********
C*********
C*********
C*********
C*********
C*********
C*********
C*********
C*********
c*********
C*********
c*********
C**********
AKZ, THE VERTICAL PAR ATTENUATION COEFFICIENT WAS *****
TREATED AS A CONSTANT (ANNUAL AVERAGES) FOR DIFFERENT *****
SIMULATION SCENARIOS IN EARLY MODEL VERSIONS. BASED
ON YORK RIVER SHOAL DATA (1984-92), AKZ IS NOW SET
SEASONALLY (TEMP BASED USING MOORE'S DEFINITIONS)
FOR DIFFERENT TIMES OF THE YEAR AND REPRESENTATIVE OF
DIFFERENT HABITATS WITHIN THE LOWER BAY (GUINEA, VIMS
AND CLAYBANK DATA BASES).THREE SCENARIOS ARE POSSIBLE
FOR EACH AREA: 1) SIMULATIONS USING SEASONAL MEANS;
2) SEASONAL MEANS + OR - 1.0 S.D.; 3) RANDOM VARIATION
WITHIN A SEASON USING THE RANGE AND LOWER LIMIT AS FOR
P'ARD (I.E. ((RANDOM* * RANGE) + LOWER LIMIT). ALL ARE *
POSSIBLE IN THE PROGRAM SECTION THAT FOLLOWS BY *
"COMMENTING" THE APPROPRIATE EXECUTABLE STATEMENTS.
*************************************************************
C
C<
C'
C'
C"
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C'
C'
c
IF (AKD .EQ. 0.0) THEN
********
********
* *******
********
THE FOLLOWING SET OF STATEMENTS ARE USED TO ***********
SET THE ATTEN.COEFF. TO EITHER SEASONAL MEANS ***********
OR TO SEASONAL RANGES DETERMINED BY THE LOWER ***********
LIMIT, MEASURED RANGE (+- ONE STAND. DEV.) AND ***********
AND A RANDOM NO.
***********
IF(JCDAY.LE.85 .OR. JCDAY.GE. 315) THEN
SEASON=WINTER
AKZ=1.596
AKZ=( (RAND* (0.616) )+ (1.288) )
ELSE IF (JCDAY.GE. 86 .AND. JCDAY.LE. 145) THEN
SEASON=SPRING
AKZ=2.247
AKZ=RAND*0.910 + 1.792
ELSE IF (JCDAY.GE. 146 .AND. JCDAY.LE. 260) THEN
SEASON= SUMMER
AKZ=2.069
AKZ=RAND*0.804 + 1.951
F*L SE
SEASON=FALL
AKZ=1.604
AKZ=RAND*0.514 + 1.347
END IF
! CLAYBANK
! CLAYBANK
! CLAYBANK
'.CLAYBANK
******** THE LAST METHOD FOR ESTIMATING THE ATTN. COEFF. **********
******** IS DERIVED AS ABOVE FOR PARD USING MONTHLY MEANS *********
******** THE RANGE IS ESTIMATED AS 2X THE STD. *********
IF (JCDAY.LE.31) THEN
AKZ = LGM*(RAND*0.50-l-0.78)
& + LGP*(RAND*0.50+0.88)
& + LCB*(RAND*0.96+1.17)
ELSE IF (JCDAY.LE.59) THEN
AKZ = LGM*(RAND*0.50-i-0.79)
& + LGP*(RAND*0.56+0.85)
& -t- LCB*(RAND*1.38+1.13)
ELSE IF (JCDAY.LE.90) THEN
AKZ = LGM*(RAND*0.80-f0.49)
& + LGP*(RAND*0.68+0.77)
& + LCB*(RAND*1.16-t-l. 09)
ELSE IF (JCDAY.LE.120) THEN
AKZ = LGM*(RAND*0.50-i-0.53)
& + LGP*(RAND*0.50+0.64)
& + LCB*(RAND*1.46+1.17)
ELSE IF (JCDAY.LE.151) THEN
AKZ = LGM*(RAND*0.66+0.58)
& + LGP*(RAND*0.88+0.95)
& -t- LCB*(RAND*1.56-H.78)
'ELSE IF (JCDAY.LE.181) THEN
AKZ = LGM*(RAND*1.12+0.62)
& + LGP*(RAND*0.96-t-1.40)
& + LCB*(RAND*3.02-l-l.ll)
ELSE IF (JCDAY.LE.212) THEN
AKZ = LGM*(RAND*1.04-i-0.69)
& + LGP*(RAND*0.78-1-1.39)
58
-------
& + LCB*(RAND*0.72+1.49)
ELSE IF (JCDAY.LE.243) THEN
AKZ = LGM*(RAND*0.72+0.78)
& + LGP*(RAND*0.70-t-1.12)
& + LCB*(RAND*0.76^1.53)
ELSE IF (JCDAY.LE.273) THEN
AKZ = LGM*(RAND*0.82+1.00)
& + LGP*(RAND*0.82 + 1.03)
& + LCB*(RAND*1. 56-1-1.34)
ELSE IF (JCDAY.LE.304) THEN
AKZ = LGM*(RAND*0.90+0.81)
& + LGP*(RAND*0.48-t-0.96)
& -t- LCB*(RAND*1.18+1. 08)
ELSE IF (JCDAY.LE.334) THEN
AKZ = LGM*(RAND*0.74+0.66)
& + LGP*(RAND*0.32+0.93)
& + LCB*(RAND*0.68-t-1.06)
ELSE
AKZ = LGM*(RAND*0.82+0.62)
& -t- LGP*(RAND*0.56+0.72)
& + LCB*(RAND*0.78+0.92)
END IF
C
ENDIF
Q******************************************+***************************
c*************** SET DO LOOP FOR WITHIN DAY ITERATIONS ************
0**********************************************************************
C
DO 100 NITER=1,ITER
C
C FEEDBACK CONTROL CALCULATIONS
C RESOURCE CONTROL FB TERMS
C
FB0203=DIM(1.,(DIM(X02,G0203)*D0203))
FB0205=DIM(1., (DIM(X02,G0205)*D0205))
FB0306=DIM(1.,(DIM(X03,G0306)*D0306))
FB0506=DIM(1.,(DIM((X05/X03),G0506)*D0506) )
C
C NOTE: DIM FUNCTION PRODUCES THE POSITIVE DIFERENCE OF TWO VALUES;
C IF X>Y, DIM(X,Y)=X-Y; ELSE DIM(X,Y)=0.
C
C SPATIAL CONTROL FB TERMS
C
FB0303=(1.-DIM(1.,DIM(X03,A0303)*D03 03))
FB0404=(1.-DIM(1.,DIM(X04,A0404)*D0404))
FB0505=(1.-DIM(1.,DIM((X05/(X03+.1E-15)),A0505)*D0505))
FB0606=(1.-DIM(1.,DIM(X06,A0606)*D0606»
C
C TOTAL FEEDBACK CALCULATIONS
C CALCULATE FIJ AND FJJ PRIME VALUES FOR MULTIPLICATIVE
C INTERACTION
C
FBP0203=DIM(1.,FB0203)
FBP0205=DIM(1.,FB0205)
FBP0306=DIM(1.,FB0306)
FBP0506=DIM(1.,FB0506)
C—> FBP0303=DIM(1.,FB0303)
C--> FBP0505=DIM(1.,FB0505)
C—> FBP0606=DIM(1.,FB0606)
C
C CALCULATE REALIZED FLUX PREFERENCE VALUES
C
C COMPUTE THE DENOMINATOR TERM
C AS THE INVERSE
C
PD06=1./((P0306*(1.-FB0306)+P0506*(1.-FB0506))+.IE-IS)
C
C COMPUTE FLUX PREFERENCES
C
TP0306=P0306*(1.-FB0306)*PD06
TP0506=P0506*(1.-FB0506)*PD06
C
C*******************************************.**************************
59
-------
c***************** FITTED OR ASSUMED ***************
C***************** FUNCTIONAL RELATIONSHIPS ***************
c***************** THAT VARy WITHIN A DAY ***************
C*********************************************************************
C
C COMPUTE THE FOLLOWING ENVIRONMENTAL VARIABLES TO 'DRIVE'
C VARIOUS COMPARTMENTAL FLOWS
C
C 4. PARHR: HOURLY INCIDENT PAR (EINSTEINS/M-2/HOUR-l)
C 5. PARSEC:INSTANTANEOUS INCIDENT PAR (EINSTEINS/M-2/SEC-l)
C 7. TLAG: TIME LAG FOR TEMPORAL PROGRESSION OF SUCCESSIVE
C TIDES OVER A LUNAR CYCLE (PERIOD =30 DAYS)
C 8. Z : WATER DEPTH IN METERS (R; 1.0 TO 1.55) BASED
C ON DEPTHS TYPICAL OF Z. MARINA IN THE LOWER
C CHESPEAKE BAY.
C 9. PARZ: SUBMARINE IRRADIANCE AT DEPTH Z
C (MICROEINSTEINS/M-2/SEC-l)
C
TLAG=LDAY*0.842105263
Z=MLW
IF (TIDE .EQ. 'Y' .OR. TIDE .EQ. 'y')
& Z = MSL + 0.5*TAMP*COS(2.*PI*(IHR-TLAG)/12.)
C
C EVALUATE PARSEC AND PARZ ONLY DURING DAYLIGHT PERIOD
C
DAYLT=DIM(COS((2.*PI*(IHR-12))/(PP*2.)) , 0.0)
IF(DAYLT.GT.O.O) THEN
PARHR=(PARD/(0.63662*PP))*DAYLT
PARSEC=277.78*PARHR
PARZ=EXP((-AKZ*Z)+LOG(PARSEC))
ELSE
PARHR=0.0
PARSEC=0.0
PARZ=0.0
ENDIF
C
C*********************************************************************
C COMPUTE VARIABLE FLUX COEFFICIENTS AS
C DRIVEN BY THE ABOVE FITTED OR ASSUMED
C FUNCTIONAL RELATIONSHIPS.
C
C VASCULAR PLANT PHOTOSYNTHESIS (T0203) AS A FUNCTION OF
C TEMP AND PARZ WHICH INCLUDES THE EFFECTS OF WATER COLUMN
C PAR ATTENUATION AND EPIPHYTE GROWTH
C
EPLR=(1.-DIM(1.,(SQRT(DIM((X05/X03),0.1)/(3.0-0.1)))))*0.5
PARVP=DIM((PARZ-(PARZ*0.75*EPLR)),0.0)
PM03=((.000162*TEMP)+.0041)*
+ (l.-(DIM(TEMP,25.)/(35.-25.»)
VPIK=((15220.78*PM03)-31.86)
DIFFL=DIM(1.0,(0.3*<(DIM((X05/(X03+.1E-15)),.l))/(3.-.l))))
T0203=DIM( (PM03MPARVP/ (VPIK+PARVP+l.E-15) )) ,0.0) *DIFFL
C
C PARVP= PAR REACHING THE PLANT CANOPY
C VPIK= MICHAELIS-MENTEN HALF SATURATION COEFF.
C DERIVED AS 'FUNCTION OF PMAX (PENHALE 1977)
C PM03= THE PREDICTED P-MAX AT GIVEN TEMP AND PARZ.
C INCREASES LINEARLY WITH TEMPERATURE TO 25 C;
C ABOVE 25C, PMAX DECLINES LINEARLY, SUCH THAT AT
C 30 C, PMAX=.5(MAXIMUM).
• C EPLR= PAR ATTENUATION DUE TO EPIPHYTE GROWTH, ESTIMATED
C AS A FUNCTION OF THE SPATIAL FEEDBACK TERM, FB0505.
C AS THE RATIO OF EPIPHYTE TO LEAF BIOMASS APPROACHES 3,
C THE LIGHT REDUCTION DUE TO EPIPHYTE ATTENUATION
C APPROACHES
C 75%. ABSORBANCE AT ANY TIME IS REDUCED BY 50% TO
C ACCOUNT FOR THE PRESENCE OF NEW, UNEPIPHYTIZED LEAVES
C IN THE CANOPY.
C DIFFL= REDUCTION IN PHOTOSYNTHETIC RATE DUE TO LIMIT TO
C BICARBONATE DIFFUCION CAUSED BY EPIPHYTE LAYER;
C APPROACHES 30% AS RATIO OF EPIPHYTE TO LEAF BIOMASS
C APPROACHES MAXIMUM.
60
-------
c
£********************»************************************************
C
C VASCULAR PLANT RESPIRATION (R0302) AS A FUNCTION OF TEMPERATURE
C AND PHOTOSYNTHESIS.
R0302= «0.305*DIM«( (.0104*TEMP)+. 3432) *T0203), 0.0))+
+EXP((.1370*TEMP)-10.09))
C
C VASCULAR PLANT LEAF MORTALITY AS A FUNCTION OF TEMPERATURE
C AND DERIVED SUCH THAT AT 0 DEGRESS, LEAF MORTALITY IS 0.5% PER DAY
C AND AT 32 DEGRESS, LEAF MORTALITY IS 3.0% PER DAY.
C
U0307=(0.0175 - 0.0125*COS(2.0*PI*IDAY/365.))/24.0
C
£*********************************************************************
C
C EPIPHYTE PHOTOSYNTHESIS (T0205) USING THE SAME FUNCTIONAL
C RELATIONSHIPS AS FOR VASCULAR PLANT PHOTOSYNTHESIS.
C
EPIK=50.+(100.*(DIM(TEMP,10.)/(30.-10.)))
PM05=EUTROF*((.0003801*TEMP)*(1.0-(DIM((TEMP-25.),0.0)/
+(45.-25.))))
T0205=DIM((PM05*(PARZ/(PARZ+EPIK+.1E-15))) , 0.0)
C
C EPIPHYTE RESPIRATION
C
C NOTE: ASSUMES THE SAME FUNCTIONAL RELATIONSHIPS AS FOR
C THE VASCULAR PLAN EXCEPT IT IS 10% THE SPECIFIC RATE
C
R0502=0.5*((0.5*DIM((((.0104*TEMP)+.3432)*T0205),0.0))+
+EXP((.1370*TEMP)-10.09))
C
p**********************************************************************
C INGESTION AND RESPIRATION OF GRAZERS
C**********************************************************************
C
T0306=GRAZF*(.00325*EXP(.0921*TEMP)/24.)*DIM(TEMP,10.)/(30.-10.)
T0506=GRAZF*0 . 805* ( (DIM(TEMP, 1. 0) / (-30. -10 .) ) /24 . 0)
R0602= (.0001046*TEMP)+.0008009
C
£*********************************************************************
Q***************** ***********************
C****************** END CALCULATIONS FOR ***********************
C****************** VARIABLE COEFFICIENTS ***********************
C****************** ***********************
C* *********************************************************************
C
C COMPUTE METABOLIC CORRECTION TERMS FOR
C SELF-CONTROLLED FEEDBACK TERMS
C
C0203=DIM((1.-(R0302/(T0203+.1E-15))),0.0)
C0205=DIM((1.-(R0502/(T0205+.1E-15))),0.0)
C0306=DIM((1.-(R0602/((T0306*AE0306)+.1E-15))),0.0)
C0506=DIM((1.-(R0602/((T0506*AE0506)+.1E-15))),0.0)
C
C CALCULATE TOTAL MULTIPLICATIVE FEEDBACK TERMS
C
TF0203=DIM(1.,(FBP0203*(1.-(FB0303*C0203))))
TF0205=DIM(1.,(FBP0205*(1.-(FB0505*C0205))))
TF0306=DIM(1.,(FBP0306*(1.-(FB0606*C0306))))
TF0506=DIM(1.,(FBP0506*(1.-(FB0606*C0506))))
C
C CALCULATE HOURLY FLUXES
C
C VASCULAR PLANT: ROOT/RHIZOME EXCHANGES
C
F0203=T0203*X03*(1,-TF0203)
F0302=R0302*X03
F0304=T0304*DIM(F0203,F0302)*(1.-FB0404)
F0402=R0402*X04
F0307=U0307*X03'* (DIM( (TEMP-10 .) , 0. 0) / (30.-10 .) )
C
C F0307=U0307*X03
61
-------
c
C EPIPHYTE EXCHANGES
C
F0205=T0205*X05*(1.-TF0205)
F0502=R0502*X05
F0507=(U0507*X05)+(DIM((X05/X03+.1E-15),0.0)*F0307)
C
C GRAZER EXCHANGES
C
F0306=TP0306*T0306*X06*(1.-TF0306)
F0506=TP0506*T0506*X06*(1.-TF0506)
F0.602=R0602*X06
F0607=(U0607*X06)-i-(F0306*(l.-AE0306) ) + (F0506* (1.-AE0506) )
C
C ENVIRONMENTAL EXCHANGES
C
F0102=0.0
F0201=0.0
C
C LOSS OF GRAZERS TO FISH PREDATION INCORPORATED USING
C PRED=SEASONAL PATTERN OF FISH ABUNDANCE IN ZOSTERA MEADOWS OF
C NORTH CAROLINA (R: 0 TO 1 GRAMS CARBON/M2), AND
C DR=DAILY RATION (%/DAY) OF PREDATORY FISH BASED ON TEMPERATURE.
C
PRED=0.5-(0.5*COS((2*PI*IDAY)/365))
DR=0.3*(DIM((-13.1+(2.29*TEMP)-(0.032*(TEMP**2))) , 0 . 0)/100./24 .) *
+(DIM(1.0,(DIM(0.2,X06)/(0.2-0.1))))
F0600=DR*PRED
F0006=0.0
C
C ASSIGN FLUX VALUES FOR EACH ITERATION
C
DFLUX(1)=F0203
DFLUX(2)=F0302
DFLUX(3)=F0304
DFLUX(4)=F0402
DFLUX(5)=F0307
DFLUX(6)=F0205
DFLUX(7)=F0502
DFLUX(8)=F0507
DFLUX(9)=F0306
DFLUX(10)=F0506
DFLUX(11)=F0602
DFLUX(12)=F0607
DFLUX(13)=F0600
DFLUX(14)=F0006
DFLUX(15)=F0102
DFLUX(16)=F0201
C
C SUM AND STORE FLUX VALUES IN STFLUX(I)
C
DO 12 1=1,16
TFLUX(I)=(DFLUX(I)*DT)+STFLUX(I)
12 CONTINUE
DO 13 1=1,16
STFLUX(I)=TFLUX(I)
13 CONTINUE
C
C SUM HOURLY FLUXES IN DAFLUX(I)
C
DO 51 1=1,9
DAFLUX(I)=DAFLUX(I)+(DFLUX(I)*DT)
51 CONTINUE
C
C SUM HOURLY RATES IN DARATE(I)
C
DARATE(1)=DARATE(1)+(T0203*(1.-TF0203))
DARATE(2)=DARATE(2)+R03 02
C
C COMPUTE ITERATIVE COMPARTMENTAL CHANGES
C
DX(1)=F0201-F0102
DX(2)=0.0
62
-------
DX(3)=F0203-(F0302+F0304+F0307+F0306)
DX(4)=F0304-F0402
DX(5)=F0205-(F0502+F0507+F0506)
DX(6)=(F0306+F0506+F0006)-(F0602+F0607+F0600)
DX(7)=(F0307+F0507+F0607)
C
C CALCULATE STATE VARIABLE CHANGE
C
DO 14 1=1,7
X(I)=DIM((X(I)+DX(I)*DT),0.0)
14 CONTINUE
C
C CALCULATE THE RATIO BETWEEN EPIPHYTE AND PLANT BIOMASS
C
RA0503=(X(5)/(X(3)+0.1E-15))
C
C INCREMENT DAY-HOUR COUNTER (IHR)
C
IHR=IHR+DT
IF(IHR.GT.24)THEN
IHR=1
ELSE
IHR=IHR
ENDIF
C
100 CONTINUE
C
C
£*****************************************************************
Q*************** END DAILY LOOP ********************************
c******************************************************************
C
C INCREMENT PRINTER AND DAY COUNTERS
C
IPRT=IPRT+1
IDAY=IDAY+1
LDAY=LDAY+1
JCDAY=JCDAY-i-l
C
C RESET DAY-HOUR COUNTER
C
IHR=12
C
C CHECK LUNAR DAY COUNTER FOR COMPLETION OF LUNAR CYCLE
C
IF(LDAY.GT.30)THEN
LDAY=0
ELSE
LDAY=LDAY
ENDIF
C
C CHECK THE PRINT/STORAGE INTERVAL COUNTER
C
IF(ISTORE.GT.ISTMAX) THEN
ISTORE=ISTMAX
GO TO 21
ELSE
ISTORE=ISTORE
ENDIF
C
IF(IPRT.LT.NDPRT)GO TO 21
C
C IF IPRT = NDPRT PRINT OUT STATE VARIBLE STANDING STOCKS
WRITE(NSST,210)IDAY,X,RA0503
C
C WRITE DAILY INTEGRATED FLUXES TO FILE —• UNITS G C / M2 / D
WRITE(NFLX,599)IDAY,DAFLUX
C
C WRITE DAILY INTEGRATED RATES TO FILE — UNITS G C / G C / DAY
WRITE(NRAT,598)IDAY,DARATE
C
C RESET THE PRINTING INTERVAL TO ZERO
C '
63
-------
IPRT=0
C
20 CONTINUE
C
C INCREMENT THE STORAGE COUNTER FOR ALL OTHER OUTPUTS
C
ITIME(ISTORE)=IDAY
C
C STORE THE FEEDBACK VALUES
C
DO 16 1=1,12
FB(I,ISTORE)=SFB(I)
16 CONTINUE
C
C STORE THE DAILY FLUX VALUES
C
DO 17 1=1,16
FLUX(I,ISTORE)=SIFLUX (I)
17 CONTINUE
C
C STORE THE ACCUMULATED (SUMMED) FLUX VALUES
C
DO 18 1=1,16
AFLUX(I,ISTORE)=STFLUX(I)
18 CONTINUE
C
C STORE PHYSICAL DATA FOR TABLUAR OUTPUT
C
DO 19 1=1,7
PHYSICAL(I,ISTORE)=PHYSDAT(I)
19 CONTINUE
C
C STORE CALCULATED RATE COEFFICIENTS FOR TABULAR OUTPUT
C
DO 199 1=1,14
TRATES(I,ISTORE)=SRATES(I)
199 CONTINUE
C
C ADVANCE THE STORAGE COUNTER AND RESET THE PLOTTING
C COUNTER
C
ISTORE=ISTORE+1
C
C
21 CONTINUE
C
C RESET THE DAILY FLUX STORAGE TO ZERO
C
DO 53 1=1,9
DAFLUX(I)=0.0
53 CONTINUE
C
C RESET THE DAILY RATE STORAGE TO ZERO
C
DO 61 1=1,2
DARATE(I)=0.0
61 CONTINUE
C
110 CONTINUE
C
C**********************************************************************
£******************** END ONE YEAR LOOP *********************
Q**********************************************************************
C
C CALCULATE ANNUAL MASS BALANCE FOR EACH COMPARTMENT
C • WHERE: SX I = SUM INPUTS
C SX O = SUM OUTPUTS
C SDX = MASS BALANCE
C
SX01I=STFLUX(16)
SX01O=STFLUX(15)
SDX01=SX01I-SX010
SXO 21=STFLUX (15) + STFLUX (2) + STFLUX (4) + STFLUX (7) + STFLUX (11)
64
-------
SX020=STFLUX(16)
SDX02=SX02I-SX02O
SX03I=STFLUX(1)
SX03O=STFLUX(2)+STFLUX(3)+STFLUX(5)+STFLUX(9)
SDX03=SX03I-SX030
SX04I=STFLUX(3)
SX040=STFLUX(4)
SDX04=SX04I-SX04O
SX05I=STFLUX(6)
SX050=STFLUX(7) +STFLUX (8) -t-STFLUX (10)
SDXO 5=SXO 51-SXO 5O
SX06I=STFLUX(9)+STFLUX(10)+STFLUX(14)
SX06O=STFLUX(11)+STFLUX(12)+STFLUX(13)
SDX06=SX06I-SX06O
SX07I=STFLUX(5)+STFLUX(8)+STFLUX(12)
SX07O=0
SDX07=SX07I-SX070
C
C BEGIN NEXT YEAR SIMULATION AND INCREMENT YEAR COUNTER
C AND RESET JULIAN DAY COUNTER (JCDAY)
C
IYEAR=IYEAR+1
JCDAY=1
C
120 CONTINUE
22 CONTINUE
23 CONTINUE
C
C OUTPUT ANNUAL FLOWS, FEEDBACK AND COMPARTMENTAL DYNAMICS
C AND PHYSICAL DATA IN TABULAR FORMATS
C
C OUTPUT TABLE: PHYSICAL DATA
WRITE(NPHY,299)
WRITE(NPHY,300)
WRITE(NPHY,301)
WRITE(NPHY,302)(ITIME(J),(PHYSICAL)I,J),1=1,7),J=1,ISTORE)
C
C OUTPUT TABLE: CALCULATED RATE COEFFICIENTS
WRITE(NMBL,399)
WRITE(NMBL,400)
WRITE(NMBL,401)
WRITE(NMBL,402)(ITIME(J),(TRATES(I,J),1=1,11),J=l,ISTORE)
WRITE(NMBL,403)
WRITE(NMBL,404)
WRITE(NMBL,405)
WRITE(NMBL,406)(ITIME(J),(TRATES(I,J),1=12,14),J=1,ISTORE)
C
C OUTPUT TABLE: FLUXES INCREMENTED BY DAY
WRITE(NMBL,229)
WRITE. (NMBL, 230)
WRITE(NMBL,232)(ITIME(J),(FLUX(I,J),1=1,8),J=1,ISTORE)
WRITE(NMBL,229)
WRITE(NMBL,231)
WRITE(NMBL,232)(ITIME(J),(FLUX(I,J),1=9,16),J=1,ISTORE)
C
C OUTPUT TABLE: FLUXES SUMMED OVER THE SIMULATION PERIOD
WRITE(NMBL,233)
WRITE(NMBL,230)
WRITE(NMBL,232)(ITIME(J),(AFLUX(I,J),1=1,8),J=1,ISTORE)
WRITE(NMBL,233)
WRITE(NMBL,231)
WRITE(NMBL,232)(ITIME(J),(AFLUX(I.J),1=9,16),J=1,ISTORE)
C
C OUTPUT TABLE: FEEDBACK TERMS
WRITE(NMBL,239)
WRITE(NMBL,240)
WRITE(NMBL,241)(ITIME(J),(FB(I,J),1=1,6),J=1,ISTORE)
WRITE(NMBL,239)
WRITE(NMBL,242)
WRITE(NMBL,241) (ITIME(J), (FB(I,J),1=7,12),J=1,ISTORE)
C
C OUTPUT TABLE: COMPARTMENTAL ANNUAL MASS BALANCE
WRITE(NMBL,249)
65
-------
WRITE(NMBL,250)
WRITE(NMBL,251)SX01I,SX01O,SDX01
WRITE(NMBL,252)SX02I,SX02O,SDX02
WRITE(NMBL,253)SX03I,SX03O,SDX03
WRITE(NMBL,254)SX04I,SX04O,SDX04
WRITE(NMBL,255)SX05I,SX05O,SDX05
WRITE(NMBL,256)SX06I,SX06O,SDX06
WRITE(NMBL,257)SX07I,SX07O,SDX07
STOP
END
66
-------
------- |