EPA Report Collection
Regional Getter for Environmental Information
U.S. EPA legion HI
Philadelphia, PA 19103
-------
-------
ReyionaK t'im-r lot 1 in ironmunr.il Inlnini
I'SH'X RcL'i.m III
IOSO \KhSl
PhiUlelpliu I'A I'MIH
U.3. EPA Region HI
Regional Center for Environmental
Information
1650 Arch Street (3PM52)
Philadelphia, PA 19103
Chesapeake Bay Ecosystem Modeling Program
The Chesapeake Bay Ecosystem Modeling Program explores how
water quality, the growth of plants and animals, and the physical and
chemical forces of Chesapeake Bay affect each other. Model
simulations help predict how things may change over time or under
different conditions. The Bay Program's ecosystem models help clarify
how the Bay's plant and animal life interact with the environment.
Ecosystem models emphasize nutrient and organic matter sources and
cycles, interactions among food web connections, and habitat
structures. These state-of-the-art models help explain how and why the
things we observe in Chesapeake Bay happen.
The Strategy for the Restoration and Protection of Ecologically
Valuable Species directs Bay Program partners to pursue development
of simulation models of the Chesapeake Bay ecosystem. Simulation
models are part of a bigger package designed to restore and protect
Bay species, at all trophic levels. Meeting the 1987 Chesapeake Bay
Agreement goal to "provide for the restoration and protection of the
living resources, their habitats and ecological relationships" requires
understanding the physical, chemical, and biological processes at work
in the Bay. The Ecosystem Modeling effort is developing a series of
interlinked models that address relationships in the Bay by simulating
critical habitats of Chesapeake Bay. These simulations will be used for
management decisions concerning land use, nutrient loadings, and fish
production.
-------
ECOSYSTEM MODELS
OF
CHESAPEAKE BAY
1994 - 1996
Progress reports
prepared for the
Living Resources Subcommittee
Chesapeake Bay Program
August 1997
Printed on recycled paper by the U.S. Environmental Protection Agency
for the Chesapeake Bay Program
-------
TABLE OF CONTENTS
Summary 1
Empirical Ecosystem Models 6
Wetlands
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands:
Ecosystem Processes and Habitat Linkages 21
Submerged Aquatic Vegetation
Trophic Control of SAV Responses to Nutrient Enrichment: Simulation
Modeling of Mesocosms Experiments 52
A Model of SAV Feedback Effects on Water Quality 58
Marsh
Simulation Model of Biogeochemical Processes in Marsh Mesocosms 72
Zooplankton
Zooplankton Process Model for Chesapeake Bay 78
A Stage-Structured Model of Zooplankton Dynamics in Chesapeake Bay 86
Plankton-benthos
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay: Model
Simulations of Responses to External Per terbations 95
Fish
Fish Bioenergetics: Relating Nutrient Loading to Production of Selected Fish
Populations 106
-------
Ecosystem Models of Chesapeake Bay
Summary
Chesapeake Bay is a very complex and diversified estuarine ecosystem. Freshwater pours
into the Bay from five major river systems. The ocean penetrates the Bay's mouth, with saltwater
sliding under the freshwater. Salinity, which affects everything that lives in the Bay, varies from
source to mouth, top to bottom, and year to year. Although much of the Bay is fairly shallow, a
deep trench stretches down the mainstem where the ancestral Susquehanna flowed. Water
temperatures fluctuate, depending on depth and season. Plants and animals crucial to the
ecosystem are often small or hidden. Microscopic phytoplankton and zooplankton inhabit the
water column, providing food for fish and utilizing nutrients. Bottom dwellers are important links
in the food web, not only as predators and prey, but because they help recycle nitrogen and
carbon. Submerged aquatic vegetation (SAV) or Bay grasses growing along the edges of the Bay
supply oxygen to the water and capture suspended sediments. SAV also provides habitat for
young fish and crabs. The Bay's hydrological and biochemical processes are complex and
interrelated with each other as well as the biota. Models use mathematical relationships to
summarize and analyze these relationships. Once created, models can help project future biotic
responses, based on environmental factors.
In addition to natural processes, human effects can't be ignored when examining the Bay's
ecosystem. Activities on land affect water temperatures, oxygen levels, and the amount of
sediment and nutrients in the Bay. If water quality improvements are to be sustained and
restoration of plants and animals successful, the Bay's natural processes and how they affect each
other must be understood. Human activities can then be factored into those processes. The
Chesapeake Bay Program depends on the scientists who gather information about the Bay, but
resource managers need to quickly apply what is learned to the Bay's pressing problems.
Ecosystem models can help managers and decision-makers direct restoration of Chesapeake Bay.
In 1984, the Bay Program initiated several long-term monitoring programs which provide
information on water quality, hydrology, and living resources in the Bay. The Chesapeake Bay
Ecosystem Modeling Program was a natural extension of these monitoring efforts. Ecosystem
models explore how water quality, the growth of plants and animals, and the physical and
chemical forces of Chesapeake Bay affect each other. Models provide a picture of the Bay's
relationships and simulations by the models help predict how things may change over time or
under different conditions.
A key component of the Ecosystem Modeling Program is effective integration of modeling
approaches. The ecosystem process models use nutrient loading and other information to predict
the quality and quantity of food and habitat available to fish populations. The fish bioenergetics
models use this food and habitat information to predict the potential production of striped bass,
bluefish, weakfish, Bay anchovy, menhaden, spot and white perch. The models will be combined
to incorporate ecological feedbacks associated with top-down control by fish of their prey and
ecosystem components, such as SAV and water quality. The two modeling frameworks for
Chesapeake Bay fish bioenergetics and ecosystem process models ~ will also be merged with
hydrodynamic models to evaluate changes in habitat conditions, flow, and nutrient loading on
living resources.
The Chesapeake Bay Ecosystem Modeling Program for 1994-1996 included models for:
submerged aquatic vegetation (SAV); littoral areas, marshes and wetlands; plankton-benthos;
zooplankton; and ecosystem and fish bioenergetics. "Box models" describe ecosystem processes
Summary 1
-------
as they change through time and space. Empirical modeling helps to determine responses of the
ecosystem to nutrient loading on various time-scales. Elements of the ecosystem models are also
being combined with the Chesapeake Bay Water Quality Model.
Model Specifics
Empirical Ecosystem Models explore the large-scale driving forces that influence plant
and animal growth. The models link growth of plants and algae, oxygen conditions of the water,
nitrogen concentrations, deposits from decaying plants, and the amount of water flowing from
rivers into the Bay. Important features of the empirical ecosystem models are linked with
vegetation and fish models to help explain how the growth of these living resources changes
through space and time.
To demonstrate the ecological connections and interactions, the ecosystem is broken
down by important components into boxes. Collectively, these box models portray regions of the
estuary as a network with patterns of physical exchange. Model simulations make it possible to
determine how residence time of biochemical, chemical and biological components vary
seasonally, with river flow, and spatially within the estuary. For example, empirical model results
indicated that 88% of the variability in freshwater residence time for Patuxent River could be
predicted by the freshwater input rate to the Patuxent. Previous reports described a strong
positive relationship between river flow and hypoxia in the Bay's mainstem. The new box models
showed that the stratification associated with river flow and the oxygen demand associated with
nutrient and organic matter loading created this relationship between flow and hypoxia.
These empirical ecosystem models offer some surprising results. They confirmed earlier
conclusions that phytoplankton biomass decreased in the upper estuary as river flow increased.
However, in the upper estuary, model results contradicted the conventional nutrient enrichment
concepts that phytoplankton biomass would be highest when river flow and nutrient loading was
lowest. According to the models, in dry years when nutrient loading was lower, water quality
tended to degrade rather than improve. In middle estuary simulations, phytoplankton biomass and
production did conform to traditional theories.
The Littoral Zone and Wetland Model describes the dynamic exchange of primary
production, paniculate and dissolved substances, and smaller biological components among the
mosaic of different habitats in shallow-water areas. Littoral environments of Chesapeake Bay
exhibit patterns of aquatic productivity, sediment processes, and biogeochemical cycling distinct
from those of adjacent channel waters.
Model results indicated that the biomass of diatoms and other plankton in each of the
habitats modeled was greatly influenced by inter-habitat exchanges. The magnitude of exchange
among habitats was much greater than exchanges associated with production and loss of
components such as carbon, chlorophyll a, and nitrogen with habitats. Estimates of annual
primary production by phytoplankton, sediment microalgae, and plants, however, is essential to
understanding shallow-water ecosystems. Phytoplankton contributed about 15% of littoral
zone/fringing wetlands production. Over one-third of annual littoral zone production was from
sediment microalgae. Production rates for microalgae were greater in nonvegetated sites than in
vegetated habitats. Eelgrass productivity accounted for 15% of total production, with
productivity from shoots four times that from epiphytes. Cordgrass productivity accounted for
almost 36% of the littoral zone ecosystem's production. Among the four different habitats
2 Ecosystem Models of Chesapeake Bay: 1994 -1996
-------
modeled, the vegetated intertidal habitat, which was the smallest habitat, accounted for 43% of
total productivity from the entire littoral zone and wetlands modeled.
The Marsh Model is currently in the early stages of development. Twelve artificial
marshes, called mesocosms, were created to determine the effects of time and biological
community complexity on marsh functions. The model and mesocosm experiments are
interactive; feedback from each is helping to refine both the experiments and the marsh model. To
date, the marsh model reasonably duplicates east coast salt marsh behavior. The model is proving
useful as an interactive tool for designing mesocosm experiments. Preliminary model results
indicated that permanent nitrogen loss through denitrification was dependent on productivity of
the marsh plants, but the importance of nitrogen removal vs. nitrogen storage in plant tissues has
not been determined.
Submerged Aquatic Vegetation (SAV) Models describe the complex feedbacks that affect
SAV growth. One model explores how the effects of nitrogen on SAV depends on both the
amount of nitrogen and the level of Bay grasses present. The model showed that the presence of
vegetation significantly affected nutrient levels in grass beds. Bay grasses reduced nitrogen levels
in water by 50%, compared with unvegetated sites. However, at high nitrogen loads, beds were
eventually eliminated and acted like unvegetated sites. The effects of nutrient loading were
increased because impaired SAV had less ability to clear sediments from the water. This impaired
self-cleaning resulted in more suspended sediment and an even lower light environment.
Another SAV model examines the consequences to grasses of interactions among fish and
invertebrate grazers at different nutrient levels. Preliminary results showed that grazing
amphipods could not control epiphytic blue-green algae. Fish grazing produced complicated
results: adding fish increased nitrogen cycling in moderate nutrient conditions and decreased
cycling in high nutrient conditions.
Zooplankton Models explore the interactions among different species and sizes of
phytoplankton and zooplankton. Zooplankton feed on phytoplankton and smaller zooplankton, as
well as cannibalize juveniles. In turn, they provide food for smaller fish. Two models will be
linked to help understand how nutrient loading interacts with this trophic level. One model uses
size-frequency distribution and stage-specific zooplankton abundances to estimate food
availability for planktivorous fish. A second model examines the environmental factors that
influence zooplankton developmental stages. Both models are in the early stages of development.
To date, patterns of species abundance and progression through developmental stages have been
successfully modeled for mesohaline regions of the Bay.
Plankton-Benthos Process Models examine the ways that plankton and benthos are
affected by environmental factors. Eutrophication and associated bottom water hypoxia can
cause dramatic shifts in microbial metabolism, as well as seasonal mortality of benthic animals.
Hypoxic conditions also inhibit the nitrification-denitrification activities of the benthos. The
plankton-benthos model helps predict complex interactions among fish, plankton, nutrients and
bottom oxygen.
Model results indicated that most of the nitrogen for production came from river flow in
spring and from regeneration in summer. Most of that regenerated nitrogen came from
organisms in the water column. The physical activities of benthic deposit feeders stimulated
Summary 3
-------
nitrification-denitrification. Fish grazing on phytoplankton created mixed results, depending on
season. When fish increased, the subsequent increased nutrient release stimulated phytoplankton
production. Model results suggested that fish feeding, whether on benthos or plankton, had little
effect on primary production or sedimentation. However, benthic suspension feeders had a
potentially significant effect.
Fish Bioenergetics Models look at the growth rate of important fish under specific
environmental conditions. They describe habitat suitability and use by specific fish species.
Bioenergetics models are complete for striped bass, bluefish, and weakfish, and two of their
preferred foods, bay anchovy and Atlantic menhaden. Bioenergetics models estimate how much
food fish will consume, based on growth rates, temperature history, and type of food. Preliminary
models also describe growth of two bottom-feeders, white perch and spot. Traditional models of
fish population production are usually based on average conditions over large areas and assume
constancy of environment. Because habitat quality in Chesapeake Bay is fragmented, three-
dimensional (3-D) modeling of the Bay is necessary. Spatial model results of carrying capacity for
menhaden showed a large degree of patchiness, depending on season.
Fish success is gauged by the number of individuals in a population, over time, and the
growth rates of individuals. Growth rate is species- and size-specific, depending on temperature
and food availability. Growth rate is also closely linked to female reproductive capabilities.
Completed models indicated that smaller striped bass generally had fewer areas in the Bay that
supported positive growth, compared with larger striped bass, except during the winter. As
striped bass aged, they increased their use of open-water food sources. Weakfish and bluefish
tended to utilize food resources on the bottom and their growth was affected by dissolved oxygen
and temperature. It appeared that migration of older striped bass to coastal waters and the
restricted use of Chesapeake Bay by older weakfish and bluefish may be due to conditions created
by certain combinations of lower temperatures and limited food supply.
In another modeled scenario, striped bass growth rate potential was modeled as a function
of the fish's horizontal position and water column depth in Chesapeake Bay. Simulations
indicated a dramatic seasonal/spatial difference in growth rate potential. During May, the best
growth potential stretched from the western shore to midchannel areas and in shallow regions
along the eastern shore of the Bay. During July, growth rate potential was good only in a small
region where the salty and freshwater layers meet and near the bottom on the eastern shore.
According to the model, striped bass growth was limited in July by two things: higher
temperatures and anoxia.
Models for Atlantic menhaden, bay anchovy and spot highlight the importance of
understanding each aspect of Chesapeake Bay's food web. Atlantic menhaden and bay anchovy
consume different types of plankton; spot feed on small bottom-dwellers. Bay anchovy numbers
do not appear to be limited by one particular food type.
Atlantic menhaden spend only a short part of their life in the Bay. Menhaden are an
important food for larger, commercially valuable fish. However, model results showed menhaden
impact on the bay, represented by how much plankton they consumed, was small. In addition to
being an important food source for larger fish, Atlantic menhaden are also commercially
harvested. The Fish Bioenergetics Models could be an important link to understanding how
commercial fishing affects the sportfish that feed on menhaden.
Ecosystem Models of Chesapeake Bay: 1994 -1996
-------
The Big Picture
Bay researchers are linking ecological models to relate fish production and water quality in
Chesapeake Bay and its tributaries to changes in nutrient and organic material inputs. This
involves developing the individual models mentioned above and analytical tools, like spatial
modeling and Geographic Information Systems (GIS). Connections among model types and with
the analytical tools will make it possible to view Chesapeake Bay as the complex ecosystem it
really is.
Ecosystem processes, habitat, and plankton, benthos and fish models will be linked with
the Bay Program's Time Variable Water Quality Model. This will create a dynamic, three-
dimensional, simulation package capable of showing the effects of spatial variation in water
quality and living resources . A user-friendly model interface, complete with user-guide, is slated
for the final stages of the project (1997-1998). With it, scientists and managers will be able to
map habitat and living resources information and simulate management scenarios. Display
possibilities include: maps of the surface later of the Bay for any factor on any date; plots of shore
outline; specific user-defined locations; and animation of output for any variable, over time, at the
surface layer, bottom layer, or 3-D view of the Bay.
Summary
-------
Empirical Ecosystem Models
J.D. Hagy and W. R. Boynton
Chesapeake Biological Laboratory
University of Maryland
Solomons, MD
Background
Reports in past years (e.g. Kemp et al. 1992, Brandt et al. 1995) demonstrated the
potential to describe relationships between nutrient loading and responses of estuarine ecosystems
in a relatively simple, direct, and empirical way. Some ecosystem responses were localized in
particular regions of Patuxent River estuary (e.g. winter, spring and summer phytoplankton,
summer hypoxia), suggesting that a careful investigation of spatial patterning would improve the
empirical ecosystem models developed earlier. The data interpolation and visualization software
developed for spatial patterning did slightly improve the models. The software also expanded the
scope of the models beyond individual stations to regions of the estuary. Moreover, the
interpolator also made it possible to calculate volume-based parameters, such as the hypoxic
volume-days parameter used for hypoxia, and included these parameters in the empirical
ecosystem models.
The broad empirical relationships indicate that, on seasonal to annual time scales,
Patuxent River estuarine ecosystem responds to nutrient loading predictably. In particular,
during years with high river flow, chlorophyll a concentrations tend to be higher in the
middle estuary and summer hypoxia tends to be more intense and longer in duration. The
responses are substantial with respect to ecological reference points, such as the habitat
criteria for Bay grasses. Water quality in high-nutrient load years tends to violate submerged
aquatic vegetation (SAY) habitat requrements more often and over a larger area of the estuary
than in low-load years. Despite these initial successes, several questions remained, including:
(1) how and why do ecosystem responses to nutrient loading vary regionally within estuaries,
among estuaries, and throughout the seasons of the year; (2) what is the relative impact of
freshwater input rates and the highly correlated nutrient loading rates on estuarine ecosystems;
and (3) what is the usefulness of the empirical models in a predictive mode and how can it be
improved?
The ecosystem models developed most recently (Hagy 1996) utilized a range of
empirical approaches in addition to regression. "Empirical ecosystem models," as opposed to
"ecosystem regression models" (the term used earlier within the ecosystem modeling component
of the Chesapeake Bay Program) best describe current ecosystem modeling work. Additionally,
box models have been added to the range of models used. These primarily descriptive models are
related to ecosystem process models in an interesting way: process models use an understanding
of ecosystem processes to investigate how both state variables and ecosystem processes change
through time and space; box models use the direct measurements of state variables through time
and space to calculate how rates of ecosystem processes changed through time and space. Thus,
box models "hindcast" important ecosystem rates into the past, permitting a whole new set of
predictive empirical ecosystem models, including regression models, to be developed. This is a
substantial extension from past work which primarily investigated responses of state variables
6 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
such as chlorophyll a or hypoxia (an exception is sediment nutrient fluxes) to nutrient loading and
other independent variables. These new box models permit much more interpretation of
mechanisms responsible for the empirical relationships we observe.
Box Models for Patuxent River
The box models for Patuxent River were developed for the following purposes, each
of which was successfully addressed with the models:
(1) Estimate residence times for several different regions of the estuary- during each month of the
year, under a range of river-flow conditions.
(2) Estimate the transport of materials of interest (i.e. organic carbon, inorganic nitrogen,
dissolved oxygen) among different regions of the estuary and between the surface and bottom
layer.
(3) Estimate the net production or consumption of the same materials in different regions of the
estuary at different times of the year and under different river flow/nutrient loading conditions.
Essentially, a box model is a system of linear equations describing the balance of salt
and water within each of several regions, or boxes. Collectively, the boxes portray a network of
regions of the estuary with an assumed pattern of physical exchange (i.e. two-layer estuarine
circulation). The equations can be solved to estimate coefficients describing the physical
exchanges. Transport of other materials can be estimated as the product of concentration times
water flow, and net production or consumption of the same material can then be calculated by
mass balance. Box models were first described by Pritchard (1969) and the technique was
extended substantially by Officer (1980). Taft et al. (1978) described a notable application of box
modeling to Chesapeake Bay. The box models developed for Patuxent River and a large portion
of the results obtained from them were described in detail in Hagy (1996); however, essential
features of the models will be described here along with selected results.
Model Structure
The Patuxent River box models divided the estuary into six regions along the
longitudinal axis of the estuary. Five of the regions were subsequently divided at the pycnocline
into two layers (Figure 1). For any one box, the possible exchanges included: lateral advective
and non-advective exchanges in two directions; vertical advective and non-advective exchanges;
and freshwater input. The salt balance is described as:
ds
m
V - -=Q .5 +Q s' -Q s +E (s'-s)+[E , (s -s)+E As -s ) (1)
"m-1 m-1 ^vm in ^jn m vm v m m' L m-l,m* m-1 m m,m+l v m-i-1 m' *
,.
and the water balance is:
QnrQ^+Qvm+Qfm (2)
Empirical Ecosystem Models 7
-------
where the terms are defined as follows:
Qm Advective transport to the down-estuary box
Qm,, Advective transport from the up-estuary box
Qvn Vertical advective transport into the box
Qfa, Freshwater input into the box
Em.lm Non-advective exchange with the up-estuary box
Emm+1 Non-advective exchange with the down-estuary box
Em Vertical non-advective exchange
sm Salinity in the box
sm_j Salinity in the up-estuary box
sm+l Salinity in the down-estuary box
s'm Salinity in the vertically adjacent box
₯ Volume of the box
m
For the model to be solvable when the water column was stratified, it was necessary to
eliminate two unknown coefficients by assuming that salt transport by horizontal dispersion was
negligible compared with horizontal advection. This was a common assertion in this type of
model (Officer 1980) and appeared reasonable for Patuxent River much of the time. Each of the
exchange coefficients shown in Figure 1 was estimated for each month of the 1985-1992 average
year and for each month between January 1985 and December 1992.
Residence Times
Estimates of hydraulic residence times were one of the objectives of the box modeling
effort. Hydraulic residence times were the average amount of time that a parcel of water spent in
the estuary. This was equivalent to the time required to reduce the concentration of a
conservative dissolved material by one exponential factor. The residence time for water entering
at the head of the estuary was the freshwater residence time. Residence times were useful for
interpreting both the time available for biogeochemical transformations to take place in the
estuary and the rates of dilution of the reactants and products of these reactions (related
concepts). Residence times were estimated from the estimated physical transport by a simple
simulation in STELLA.
Using the simulations, is was possible to determine how residence time varied
seasonally, with river flow, and spatially within the estuary (Figure 2). For example, residence
time was shorter if a material entered the outflowing lower estuary surface water than if it entered
with the inflowing bottom water or upstream at the head of the estuary. Note that these models
were not hydrodynamic simulations models.
The model results indicated that 88% of the variability in freshwater residence time for
Patuxent River could be predicted by the freshwater input rate to Patuxent River (Figure 2). The
relationship was described by the hyperbolic function:
T= - - - , r2=0.88
0. 009+0.
where T = the freshwater residence time in days and Qf= the freshwater input rate in units
of m3 s"1. Other residence times for Patuxent River did not relate as well to freshwater
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
input rates to Patuxent River (Figure 2). In fact, for the lower estuary, flushing was better
predicted by the rate of inflowing bottom water, gravitationally driven flow, and tidally
driven flow. These flows were determined by the density structure of the water column at
the mouth of Patuxent River (Figure 3). Therefore, the middle and lower reaches of
Patuxent estuary showed no increase in flushing rate during high river flow to mitigate the
effects of the high-flow induced, high nutrient loading from Patuxent River.
Nonconservative Box Models
In addition to the residence time estimates, box models estimated a range of
nonconservative nutrient and organic matter fluxes for the period between 1985 and 1992. The
net nonconservative flux of materials in any box was calculated from the mass balance equation
for that material in the box. Many of these conclusions are discussed in detail in Hagy (1996).
The box models were used to estimate the net production of total organic carbon
(TOC) and dissolved inorganic nitrogen (DIN; which included ammonium, nitrate, and nitrite) in a
10 km reach of Patuxent estuary located between St. Leonard's Creek and Broomes Island. The
box model-based estimate of annual TOC production was amended with estimates of organic
matter deposition to sediments and transfer of organic matter to fish. These estimates were then
converted to nitrogen to estimate the inorganic nitrogen demand associated with net organic
matter production (Figure 4). An estimate of inorganic nitrogen uptake made it possible to
estimate, by difference, the amount denitrified. Denitrification in the middle Patuxent was
estimated at approximately 736 mmol N m"2 y"1. This amount was within the range of
denitrification values reported in the literature, but was three times greater than the denitrification
ra'te estimated for Patuxent River by Twilley and Kemp (1985) suggests that denitrification may
be a greater sink for nitrogen in Patuxent River than previously believed. However, other sinks
for nitrogen, such as marshes, were also significant and may have accountted for some of the
difference.
Eleven years of high quality monitoring data are now available. Therefore, it may be
possible to make this calculation for years with different levels of nutrient loading and model how
nutrient loading has affected denitrification rates in Patuxent River. If in fact denitrification
becomes a "free" source of nitrogen reduction in Chesapeake Bay as oxygen conditions improve,
this positive feedback may already be apparent in the water quality record for Patuxent River.
New Predictive Models
New empirical models relating external forcing of the Patuxent River ecosystems to
ecosystem responses were possible as a result of the flux estimates generated by the box models.
Special attention was given to: (1) bottom layer oxygen demand and its relation to river flow and
hypoxia; and (2) net phytoplankton production and its relationship with nitrogen loading and
phytoplankton biomass.
Hypoxia
Previous reports described a strong positive relationship between river flow and annual
hypoxia (hypoxic volume-days). Box model-generated bottom layer oxygen demand estimates
showed that the relative role of stratification effects associated with river flow and oxygen
demand effects associated with nutrient and organic matter loading created the positive
Empirical Ecosystem Models 9
-------
relationship that was observed. Moreover, the functional form of the hypoxia vs. river flow
relationship was reevaluated (Hagy 1996; Chapter 3).
During 1985-1992, annual hypoxia (the annual integral of hypoxic volume) varied
between 1,497 106 m3 d y'Mn 1986 and 6,923 106 m3 d y'1 in 1987. To put these values in
perspective, we estimated that maximum likely value that would be expected, based on
assumptions determined from experience with the Chesapeake Bay Water Quality Monitoring
Program Data, and called it the "maximum likely hypoxia." These assumptions were that: (1)
hypoxia was not likely to occur when water temperatures were too cold for rapid metabolism to
occur or when the water column was not well stratified; (2) with limited exceptions, only waters
below the pycnocline would be persistently hypoxic or anoxic, thereby providing a maximum
volume; and (3) hypoxia would not affect the entire maximum volume immediately, but would
increase in extent steadily to the maximum, then decrease. The first decrease in dissolved oxygen
concentrations from saturation levels were usually observed in May. Thus, significant hypoxia
would not be expected prior to June 1. Because a breakdown of stratification usually occurred
toward the end of August, hypoxia would be expected to persist only rarely beyond August 3 1 .
Assuming a sinusoidal increase to the maximum hypoxic volume beginning June 1, then a similar
decline concluding on September 1, the maximum likely hypoxia for Patuxent River would be
9,869 10* m3 d y"1. The maximum observed annual hypoxia during 1985-1992 was 70% of this
maximum. Thus, hypoxia in Patuxent River, in the worst years, was severe but could have been
worse.
Calculating "maximum likely hypoxia" indicated that the relationship between nutrient
loading and hypoxia could not be linear unless it was assumed that the relationship broke down at
high levels of nutrient loading. Although the relationship between river flow and annual hypoxia
had a linear component, the best relationship was obtained using a hyperbolic functional form
(Figure 5).
The best-fit hyperbola (r^O.92) was:
Hypoxia=9562 - ( 4 ,
Flow
where Hypoxia = annual hypoxia (103 m3 d y"1) and Flow = water-year mean fall line river flow
(m3 s'1)- The estimated asymptotic value of the equation (9,562 103 m3 d y"1) was very close to
the estimated maximum likely annual hypoxia (9,869 103 m3 d y"1). The flow level at which no
hypoxia would be expected was 5.25 m3 s'1, barely lower than the average flow that occurred in
1986. To investigate the relative role of water column stratification and bottom layer oxygen
demand in generating the above relationship, both factors were individually related to river flow.
There was no significant relationship between hypoxia and river flow averaged over any period
and stratification. However, bottom layer oxygen demand increased with river flow and
approached an asymptotic value (Figure 6). The relationship (1^=0.88) was:
00=1492--- (5)
Flow
10 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
where OD = annual bottom layer oxygen demand (g O2 m"2 y"1) in equation 4 and Flow =
water-year mean fall line river flow (m3 s"1). The proper functional form for the relationship
between bottom layer dissolved oxygen demand and annual hypoxia was not statistically clear
from the data. However, only the linear relationship had conceptually valid predicted values.
The best fit linear line relating oxygen demand to hypoxia was:
Hypoxia= -43 94 +9.48(Q£>) {6)
where the terms of the equation are as described above. The composite function obtained
by combining eq. (5) and eq. (6) was similar to eq. (4). These results indicated that, in
Patuxent River, oxygen demand was related to river flow and appeared to drive the
formation of hypoxia. Because there was no direct conceptual link between river flow and
oxygen demand, the nutrient and organic matter loading associated with river flow must
drive the oxygen demand. Therefore, a relationship between nutrient loading and organic
carbon production in the estuary was expected. The interannual variability was significant
in terms of its capacity to generate hypoxia in Patuxent River.
Phytoplankton Biomass and Net Production
In previous reports, a series of regression models related chlorophyll a concentrations
to fall line nutrient loading and river flow (Brandt et al. 1995; Kempetal. 1992). A negative
'relationship between river flow and chlorophyll a was noted in the upper estuary and a positive
relationship in the middle estuary. It was hypothesized that the upper estuary relationship was
primarily driven by the dilution-effects of river flow; whereas, the middle estuary relationship was
driven by nutrient enrichment. To test this hypothesis, the regional difference in the effect of
Patuxent River discharge rate on residence times was quantified. The results did not contradict
the hypothesis.
To further examine the hypothesis, net phytoplankton production was estimated for
the upper estuary and middle estuary using the box models. The new series of box models used
individual cruises as observations, rather than annual means used by the early regression models,
and a non-parametric approach as opposed to regression. In general, the results from the box
models were much more conclusive and the statistical approach more robust (Hagy 1996).
The box models confirmed earlier conclusions from the regression models.
Phytoplankton biomass decreased in the upper estuary as river flow increased (Figure 7).
Conversely, phytoplankton biomass was highest when river flow and nutrient loading was lowest,
contradicting conventional nutrient enrichment concepts.. During much of the year, and
especially the winter months, all nutrient concentrations were many times greater than limiting
concentrations, eliminating nutrient limitation as an important factor. Given the high turbidity at
all times of the year, light limitation was likely to occur and self-shading was probably important.
Thus, the dilution effects related to high flow would be expected to both decrease phytoplankton
biomass and increase net phytoplankton production in the upper estuary. This is exactly the result
that was obtained by the model (Figure 8). Net phytoplankton production, especially during
summer, increased with river flow. From a management perspective, these results may be
disconcerting. In dry years when nature offers a temporary reprieve from some diffuse source
nutrient loading, water quality tends to degrade rather than improve. However, limiting nutrient
Empirical Ecosystem Models 11
-------
concentrations are occasionally observed during summer, and nutrient loading reductions will
likely increase the frequency of these occurrences.
Although patterns of phytoplankton biomass and net production in the upper estuary
did not conform to the "agricultural paradigm," middle estuary phytoplankton biomass and
production did conform. Phytoplankton biomass tended to increase with nutrient loading (Figure
9) and the increase was related to in situ phytoplankton production (Figure 10), rather than
advection of phytoplankton from the upper estuary into the middle estuary. Because nitrogen
concentrations in the middle estuary were frequently below limiting concentrations during summer
and there was no discernible residence time effect of river flow, the causal chain of high nutrient
load -* nutrient enrichment -> increased phytoplankton production -4 increased phytoplankton
biomass was plausible for the middle estuary. However, a notable effect of increased nutrient
loading to the middle estuary was to increase the variability in phytoplankton biomass in addition
to the mean biomass.
Mainstem Box Models
The larger focus of the Chesapeake Bay Program modeling efforts dictated moving the
focus of analysis from Patuxent River into the mainstem Chesapeake Bay. This required
developing a box model for the mainstem Bay. The main goals of this effort are similar to those
for Patuxent River; however, the possibilities for synthesis of research approaches may be greater
due to the larger number of studies conducted on Chesapeake Bay. The effort will be more
difficult, however, because of the greater complexity of Chesapeake Bay.
The Chesapeake Bay box models are branched, two-dimensional, time-variable
models. The only aspect of the models that is fundamentally different from the Patuxent box
models is the branching. The mainstem box models simulate exchanges between the mainstem
Bay and the James, York, Rappahannock, Potomac, Patuxent, and Choptank Rivers. The
structure of the box models is shown in schematic form in Figure 11. At present the models
results seem reasonable; however, the models are still at a developmental stage.
12 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
References
Brandt, S. B., W. R. Boynton, W. M. Kemp, R. Wetzel, R. Bartleson, J. D. Hagy, K. J. Hartman,
J. Luo, C. J. Madden, M. Meyers, T. Rippetoe, D. D. Fugate, and R. Batiuk. 1995.
Chesapeake Bay Ecosystem Modeling Program. 1993-94 Technical Synthesis Report to
the Chesapeake Bay Programs's Living Resources Subcommittee and Modeling
Subcommittee. Chesapeake Bay Program. Annapolis, MD.
Hagy, J.D. 1996. Residence times and net ecosystem processes in Patuxent River estuary. M.S.
Thesis. University of Maryland at College Park.
Kemp, W. M., S. B. Brandt, W. R. Boynton, R. Bartleson, J. Hagy, J. Luo, and C. J. Madden.
1992. Ecosystem models of the Patuxent River estuary relating nutrient loading to
production of selected fish populations. Final Report, Year 1. Maryland Department of
Natural Resources. Tidewater Administration. Chesapeake Bay Research and Monitoring
Division. Annapolis, MD.
Officer, C. B. 1980. Box models revisited. In P. Hamilton and R. B. Macdonald (eds.).
Estuarine and Wetland Processes. Marine Sciences Series. Vol. 11. New York:
Plenum Press.
Pritchard, D. W. 1969. Dispersion and flushing of pollutants in estuaries. American Society of
Civil Engineers. Journal of the Hydraulics Division. 95(HYI): 115-124.
Taft, J. L., A. J. Elliott, and W. R. Taylor. 1978. Box model analysis of Chesapeake Bay
ammonium and nitrate fluxes, pp. 115-130. In M. L. Wiley (ed). Estuarine Interactions.
New York: Academic Press.
Twilley, R. R., W.M. Kemp, K.W. Staver, J.C. Stevenson, 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
Empirical Ecosystem Models 13
-------
Distance from Estuary Mouth (km)
40 Upper Estuary
Middle Estuary
Lower Estuary
Figure 1. Schematic diagrams of the box model structure showing the box
boundaries, the exchange coefficients that were estimated, and model input s.
The estimated exchanges include seaward advection (QJ landward advection
(Q'.p), vertical advection (CXJ, vertical diffusive exchange (E vm), and
horizontal dispersion (Qmm+i)- Inputs included the volume in each box (Vmand
V'J concentration of salt'for each box (Sm or S'J, the input of freshwater to
each box (QfJ, and the salinity at the seaward boundary (Scb). The sill noted
at the boundary of box 1 and box 2 effectively terminates the landward flow in
the bottom layer, forcing the water into the surface layer.
14
-------
10 20 30 40
Freshwater Input Rate (m3 s"')
50
Figure 2. Residence times for surface layer boxes 1 -6.
This is the time required to flush from the estuary 63% of a
material introduced to a single box in an instantaneous
pulse RTi is equivalent to the residence time of
freshwater. In addition to RTor the average year, RiT
for each February and August during 1985-1991
Exponential least-squares regression lines are shown for
RTi, RTzandRB
200 300 400
Seawater Inflow (m3 s'1)
500
Figure 3. Residence times (RT) for surface layer boxes
4-6 related to seawater influx (Q'cb). RT is defined as in
Figure 2. Exponential regression lines are shown for
RTs and RTe, but not for RT4
15
-------
401 NEM
103 Fish
819
Denitrification
357
Burial
units: mmol N rf? y"1
Figure 4. The balance of nitrogen for the middle Patuxent estuary (box 4). The
inward arrows indicate the combined sources of inorganic nitrogen as determined
by the box models. NEM is the nitrogen export associated with the net production
or organic matter that is exported. Loss of nitrogen to fish is from Hagy (1996, p.
93), nitrogen burial is from Boynton et al. (1995), and denitrification is calculated
by difference.
>
a
>.
a
i 10 15
Flow (m s )
20
Figure 5. The relationships between water-year mean fall line
river discharge and annual hypoxia (upper panel) and between
summer (June-August) mean fall line river discharge and annual
hypoxia (lower panel). Correlation coefficients were calculated
without the 1987 observation, and the curves were hand drawn.
16
-------
o
3
C
O
O
So
1200
1100
1000
900
800
700
600
500
198Z
5 7 9 11 13 15
Water-Year River Flow
(m3s-1)
Figure 6. The relationship between water-year river
flow and annual bottom layer oxygen consumption
in the middle estuary for the period 1985-1991. The
correlation coefficient was calculated without the
1987 observation. The curve is hand drawn to sho
the functional form only.
f
V
3
Increasing Algal Biomass (m "2 )
Figure 7 The relationship between phytoplanWon biomass
and nutrient loading m the upper Patuxent River estuary
The nutrient loading and algal biomass observations were
divided into quartiles by rank and the graph depicts the
frequency in each of the 16 resulting cells of a contingency
table The expected row percentage for each column is
therefore 25% and departures from this value represent
non-random association
17
-------
1?
t?A
3-5
txpeaea=oir/o
Increasing Net
Phytoplankton Production (m ~2)
Figure 8. Upper estuary net phytoplankton production vs.
river flow and nutrient loading. The ranked nutrient foadin
and net phytoplankton production observations were
divided into two groups at the median and the graph
depicts the frequency in each of the four resulting cells of
contingency table. The expected row percentage for each
column is therefore 50% and departures from this value
represent non-random association.
CO
_c
"5
CO
o
c
g>
'_
"5
CO
CO
-------
O)
_c
'a
CO
o
o>
c
CO
eo
£
u
c
Increasing Net Production (m2-)
Figure 10. The relationship between middle Patuxent estuary
net phytoplankton production and total nitrogen (TN) loading.
The nutrient loading and net production observations were
divided into quartiles by rank and the graph depicts the
frequency in each of the 16 resulting cells of a contingency table.
The expected row percentage for each column is therefore 25%
and departures from this value represent non-random
association.
19
-------
usquehanna R.
Bottom Layer Advection --
Surface Layer Advection &
Surface Layer Dispersiono-o
Bottom Layer Dispersion .*-».
Vertical Advection ©
Vertical Diffusion o
Eastern Bay
Choptank R.
Patuxent R.
Tangier Sound
Potomac R.
Rappahannock R.
Pocomoke Sound
Figure 11. A schematic diagram of the box model being developed for
Chesapeake Bay. The arrows pointing toward the outside of the model
boxes indicate estimated freshwater inputs which are inputs to the box model
calculation. The model has 9 segments in the mainstem Bay, of which 8 are
divided into a surface and bottom layers. The mainstem segments exchange
with the James, York, Rappahannock, Potomac, Patuxent, and Choptank
Rivers. Pocomoke Sound, Tangier Sound, and Eastern Bay are included in
the time-variable salinity calculations, but exchanges with the mainstem are
not estimated by the models.
20
-------
Modeling the Lower Chesapeake Bay Littoral Zone & Fringing Wetlands:
Ecosystem Processes and Habitat Linkages
Christopher P. Buzzelli, Richard L. Wetzel and Mark B. Meyers
School of Marine Science / Virginia Institute of Marine Science
College of William & Mary
Gloucester Point, VA
Background
The estuarine littoral zone is a mosaic of different habitat types which are connected by
the dynamic exchange of primary production, particulate and dissolved substances, and animal
populations (Correll et al. 1992; Cbilders et al. 1993; Kneib and Wagner 1994; Rozas 1995). A
number of coastal studies have focused on subsystem interactions within coastal marsh and
shallow near-shore ecosystems (Wolaver et al. 1983; Stevenson et al. 1988; Dame et al. 1991;
Correll et al. 1992; Vorosmarty and Loder 1994). These studies quantified material production
and exchange in habitats fringing channel and upland environments. Although approximately 40%
of the subtidal area of Chesapeake Bay is < 2.0 m below mean low water (MLW), littoral zone
ecosystems have not been included in efforts to simulate baywide environmental processes (Kuo
and Park 1995; Spinner 1969). Biogeochemical processes in the fringing environments differ
from those of the adjacent channel, but the two estuarine zones are linked on daily, seasonal, and
annual time-scales (Malone et al. 1986; Kuo and Park 1995). Watershed factors, such as river
discharge and nutrient runoff, can influence the annual patterns of production and nutrient cycling
in the estuarine littoral zone (Correll et al. 1992). Understanding of the ecosystem processes and
habitat patterns that occur within these fringing estuarine environments will help assess the role of
wetlands and the littoral zone in Chesapeake Bay ecosystem dynamics.
The littoral zone environments of Chesapeake Bay exhibit patterns of aquatic
productivity, sediment processes, and biogeochemical cycling distinct from those of adjacent
channel environments (Kuo and Park 1995; Malone et al. 1986). Few published studies have
employed mechanistic models to analyze habitat interactions among coastal ecosystem
components to identify the probable linkages to other areas of the landscape (Boumans and Sklar
1990; Childers 1992; Costanza 1990). Understanding the synergistic interactions among littoral
zone habitats would provide an essential link between the preservation of environmental quality
and the protection of living resources, such as plant communities and fishery species (Dennison et
al. 1993; Heck and Thoman 1984; Kneib and Wagner 1994). Mechanistic models can help
address issues related to environmental change in coastal environments (Costanza et al. 1990;
Wetzel and Hopkinson 1990).
Each of the different littoral zone habitats have a suite of primary producers including
water-column phytoplankton, sediment microalgae, submerged aquatic vegetation (SAV) with
attached epiphytic communities, and marsh grasses. Submarine irradiance, along with other
meteorologic (seasonally in temperature and rainfall) and hydrodynamic factors (nutrient runoff
and river flow), influence estuarine phytoplankton processes (Mallin 1994). Sediment microalgae
significantly contribute to primary production in nonvegetated and vegetated subtidal
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 21
-------
environments in many ecosystems, including those in Massachusetts (Gould and Gallagher 1990),
South Carolina (Pinckney and Zingmark 1993a; 1993b) Mississippi (Moncreiffer al. 1992;
Sullivan and Moncreiff 1988), and Denmark (de Jonge and Colijn 1994; Sand-Jensen and Borum
1991). Sediment microalgae also play an important role in the fluxes of oxygen and nutrients
across the sediment-water interface (Rizzo et al. 1992; Sundback and Graneli 1988). SAV
meadows are complex and productive littoral zone ecosystem components (Moncreiffe? al.
1992; Murray and Wetzel 1987; Roman et al. 1990; Sand-Jensen and Borum 1991). SAV
meadows can serve as indicators of water quality because the plants and attached epiphytes are
sensitive to submarine light attenuation and the concentrations of chlorophyll a, suspended
sediments, and inorganic nutrients (Bach 1993; Dennison et al. 1993; Wetzel and Neckles 1986).
Salt marshes are areas of increased rates of productivity and nutrient cycling (Childers et
al. 1993; Dame and Kenny 1986; Pinckney and Zingmark 1993 a). Few studies, however, have
focused on the estuarine fringing wetlands of Chesapeake Bay (Drake and Read 1981; Gross et
al 1991; Wolaver et al. 1983). Ecosystem modeling offers the opportunity to include all of the
principal autotrophs in an analysis of primary production and water quality over multiple habitats
in the littoral zone.
Previous work focused on the development and simulation analysis of SAV models and
conceptual modeling of emergent intertidal marsh communities. The SAV models clearly showed
the importance of environmental factors (submarine light, temperature) and biological factors
(epiphytic fouling, grazing) for controlling SAV growth, distribution, and long-term population
survival (Wetzel and Neckles 1986). The SAV stand-alone model proved an accurate predictor
of water quality-SAV response and habitat criteria for SAV survival (Wetzel and Meyers 1994).
Over this past year, the SAV model was revised and expanded to include other components of the
littoral zone. This effort will make it easier to relate littoral processes - which includes the
benthos, SAV, and pelagic habitats - to models of hydrodynamics and water quality extant for
Chesapeake Bay and its major tributaries.
The Lower Chesapeake Bay Littoral and Fringing Wetlands Model was developed,
calibrated, and validated for habitats characteristic of mesohaline and polyhaline regions of the
Bay. Simulation analyses of these ecosystem process models was conducted for specific, highly-
distributed components of the estuary which emphasized intertidal wetlands, SAV habitats, and
other principal components of the littoral zone. The conceptual models of the principal habitats of
the littoral zone were refined and implemented into numerical simulation models. Incorporating
spatially varying information such as salinity, nutrient concentration, and bathymetry as forcing
functions can suggest how SAV-driven, phytoplankton-driven, and detrital and benthic
microflora-driven food webs function along the tributaries and Chesapeake Bay. One of the goals
was to formulate both spatially and temporally varying forcings in ways that will enable the incor-
poration of biological productivity and biologically driven elemental cycling (e.g., for carbon,
oxygen, and nitrogen) into large-scale, water quality and hydrodynamic models.
Model Description
Reference Site Description
Data for model development, calibration and verification were taken from the literature
and to the extent possible from a single reference site in the lower Bay: the Goodwin Islands
National Estuarine Research Reserve (GINERR). The GINERR is an 800 hectare (ha) littoral
22 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
zone ecosystem at the mouth of York River, Virginia, in the lower Chesapeake Bay (37° 12' 46"
N, 76° 23' 46" W). The islands are owned by the College of William and Mary and are managed
by the Chesapeake Bay National Estuarine Research Reserve System in Virginia (CBNERRS-VA)
of the National Oceanic and Atmospheric Administration (NOAA). The research reserve includes
the islands and a buffer zone that extends seaward to the -2.0 m depth contour (MLW). The -2 m
depth contour is operationally defined as the littoral zone depth limit and the corresponding lateral
boundary condition for the tributary water quality-hydrodynamic model. The GINERR is an
oblong island system with a large subtidal shoal extending between the shoreline and the -2.0 m
depth contour. Between -1.0 and -0.5 m (MLW) is approximately 120 hectares of subtidal SAV,
mostly comprised of eelgrass (Zostera marina L). There is approximately 100 hectares of
nonvegetated intertidal habitats, with fine sands and silty sediments that surround 90 hectares of
intertidal marsh vegetated primarily by smooth cordgrass (Spartina alterniflord). Some marsh
regions are also vegetated by meadow cordgrass (Spartina patens), spikegrass (Distichlis spicatd)
and needlerush (Juncus roemerianus). The intertidal marsh grades into a salt bush habitat that
includes marsh elder (Iva frutescens) and groundselbush (Baccharis halimifolia). The largest
island has a small amount of maritime forest and upland vegetated by loblolly pine (Pinus iaedd)
and several hardwood species. Intertidal and subtidal habitat patterns vary over time (seasonally,
interannually) and space (10's to 100's ha). Historical aerial photography depicts long-term
stability in the GI NERR eelgrass meadows, but overall erosion and some horizontal migration for
intertidal marshes.
Conceptual Model Structure
Four, hydrodynamically linked submodels were developed to represent the principal
littoral zone habitats: (1) nonvegetated subtidal habitats (NVST) composed primarily of coarse
sand; (2) vegetated subtidal habitats (VST) dominated by eelgrass; (3) nonvegetated intertidal
(NVIT) typical of sand and mudflats; and (4) vegetated intertidal (VIT) dominated by smooth
cordgrass. These four habitats and, thus, the models were selected based on abiotic and biotic
characteristics relative to the depth gradient along which they were located. Figure 1 depicts the
conceptual model structures for the four littoral zone habitat models based on the four habitat
types.
The global forcing functions were tidally varying water level, solar and submarine
irradiance, and water temperature. The subtidal and intertidal nonvegetated models had seven
state variables each, including: large and small phytoplankton size classes (diatoms and other
plankton, respectively); labile and refractory paniculate organic carbon (LPOC and RPOC);
dissolved organic carbon (DOC); total dissolved inorganic nitrogen (TDIN); and sediment
microalgae (SM). In addition to these seven state variables, the vegetated subtidal and intertidal
habitat models included additional state variables representing epiphyte carbon (ZepiC), and shoot
and root-rhizome carbon and nitrogen for eelgrass (ZSC, ZSN, ZRRC, ZRRN) and smooth
cordgrass (SSC, SSN, SRRC, SRRN). Table 1 gives a list of the model state variables, their
mathematical abbreviation and modeled units.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 23
-------
A
111
Q
H
Q
O
O
24
OJ
SJ
3)
-------
Table 1. List of state variables for habitat models. Each habitat model includes the first 7 state variables
listed. In addition to the basic seven the vegetated subtidal habitat model (VST) includes those related to
Zostera marina; whereas, the vegetated intertidal habitat model (VIT) has those related to Spartina
alterniflora.
ABBREV.
DIA
OP
LPOC
RPOC
DOC
TDIN
SM
ZSC
ZSN
ZRRC
ZRRN
ZepiC
SSC
SSN
SRRC
SRRN
DESCRIPTION
Diatom Carbon Mass
Other Plankton Carbon Mass
Labile Particulate Organic Carbon
Refractory Particulate Organic Carbon
Dissolved Organic Carbon
Total Dissolved Inorganic Nitrogen
Sediment Microalgae
Zostera marina Shoot Carbon
Zostera marina Shoot Nitrogen
Zostera marina Root-Rhizome Carbon
Zostera marina Root-Rhizome Nitrogen
Zostera marina Epiphytic Biomass
Spartina alterniflora Shoot Carbon
Spartina alterniflora Shoot Nitrogen
Spartina alterniflora Root-Rhizome Carbon
Spartina alterniflora Root-Rhizome Nitrogen
UNITS
gc
gc
gc
gc
gc
/uM
gCnr2
gCm'2
gNm"2
gCm'2
gNm'2
gCm'2
gCrn'2
gNm-2
gCm-2
gNm'2
Hydrodynamic Model Design
The four ecosystem process models were hydrodynamically linked by tidal exchange
across the boundaries of a sequence of modeled cells representing the NVST, VST, NVIT, and
VIT habitats. The cells filled and drained in sequence relative to the depth gradient, with the
output from one cell providing the input for the next in the sequence dependent on tidal stage (i.e.
ebb or flood periods). The nonvegetated subtidal habitat was bounded laterally by an infinite
source/sink representing the offshore channel whereas the vegetated marsh was bordered by the
upland with no exchange across the marsh-upland boundary. Watershed exchanges were assumed
to be zero because the Goodwin Islands have little upland area and are isolated from the mainland.
The cell (habitat) volume changed with each model time or integration interval (dt). Tidally
driven exchanges for water column constituents (i.e. phytoplankton, DOC, LPOC, RPOC, and
TDIN) were derived using finite difference solutions to equations for mass exchanges between a
channel and an adjacent control volume in both flood and ebb directions (Kuo and Park 1995).
This approach assumed no diffusion, no advection, and that the water within each cell was
instantly mixed and homogeneous during each time-step. The change in tidal height with each
time step was multiplied by habitat wet area to derive the changes in habitat volumes used in the
simulation of water column processes. By definition, subtidal habitat wet areas were constant.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 25
-------
Intertidal habitat wet areas varied as a function of tidal inundation and were derived using a
hypsometric relationship. Hypsometry was used because it provided a concise method for
representing the cumulative characteristics of basin morphology (Friedrichs and Aubrey 1994;
Strahler 1952). The area-height relationship of a hypsometric curve provided a better
approximation for basin inundation regimes than did a linear 2-D profile, because it included the
effects of shoreline curvature (Boon and Byrne 1981; Friedrichs and Aubrey 1994). Also,
hypsometric determination of inundation can be useful in the analysis of wetland biogeochemical
cycling (Childers et al. 1993; Eiser and Kjerve 1986).
The tidal exchange equation for a constituent (e.g., chlorophyll a) of mass M; where the
subscript I = {!,..,4} represented each of four habitats is given in equation 1 below.
AA4;
A/
Note that M; is the total mass of a water borne constituent contained in the cell or habitat volume
and can be calculated as C; * h(t)* A(t), where Q is the water column concentration, h(t) is time
varying water depth and A(t) is the time varying wetted area of the habitat. A(t) was constant for
the subtidal habitat, but variable for the intertidal habitats. The tidally varying water height, h, was
referenced to mean sea level, and its change from one model time-step to the next was
represented as Ah. Other processes affecting state variable masses were growth or biochemical
production (a), losses from biological uptake or mortality and/or grazing (m); and exchanges
with the benthos (b). In the present model, 1=0 represented the channel boundary condition.
Ecosystem Processes Model Design
Ecological processes represented in the models for the various habitats included the
principal factors controlling the uptake and loss of organic carbon and nitrogen as the primary
limiting nutrient in these habitats. Detailed treatment of the mathematical structure, interaction
coefficients and data sources for the governing equations are given in Buzzelli et al. (1995) and
Buzzelli and Wetzel (1996) and are beyond the scope of this report. Given below is a general
description of the ecological processes modeled and the state variables affected.
Primary production (gC m"2 or m"3 d"1) was modeled using the rates of gross production,
respiration, and loss through mortality or grazing. Phytoplankton (diatoms, DIA; and other
plankton, OP) were also influenced by exudation, sedimentation, and transport to adjacent
habitats. The mathematical representations of the basic metabolic rate processes in diatoms, other
plankton, sediment microalgae, and Spartina alterniflora were all similar. Gross production was
affected by temperature, irradiance, and dissolved inorganic nitrogen. Respiration followed an
exponential relationship with temperature (Cerco and Cole 1994). Production arid mortality were
represented by Gaussian functions with temperature. Phytoplankton exudation and sedimentation
were also modeled according to Cerco and Cole (1994). Sediment microalgae were lost through
resuspension and grazing by higher trophic levels. The formulations for carbon productivity by
Zostera marina and its epiphytes were taken from Wetzel and Neckles (1986) and Wetzel and
Meyers (1994). Nitrogen uptake by the shoots and root-rhizomes of Zostera marina were
modeled using Michaelis-Menten kinetics, limited by feedback functions based on the maximum
and minimum nitrogen contents of the tissues. Zostera marina shoots and root-rhizomes
maintained C:N ratios through the proportional nitrogen loss terms. Nitrogen was only
translocated from root-rhizomes to shoots in order to meet shoot nitrogen demand. Nitrogen
26 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
translocation was also limited by feedback functions based on the maximum and minimum
nitrogen contents of the tissues. The formulations for nitrogen state variables ofSpartina
alterniflora were similar to those ofZostera marina, except there was no shoot uptake of
nitrogen in Spartina alterniflora.
Water column paniculate organic carbon (POC; gC m"3) was influenced by production,
hydrolysis, settling, and exchange between adjacent habitats. POC was produced from
phytoplankton and a fractional loss term added to that gained through resuspended sediment
microalgae. POC was divided into labile and refractory fractions and rates of hydrolysis were
calculated using an exponential relationship with temperature (Cerco and Cole 1994). LPOC and
RPOC both settled from the water column and were exchanged laterally. DOC was influenced by
production, remineralization, and exchange with adjacent habitats. Hydrolyzed POC provided the
DOC production rate and the remineralization rate was controlled by a temperature function and
the refractory DOC fraction (Cerco and Cole 1994). Water column TDIN (mmoles m"3) was
influenced by production, autotrophic uptake, sediment-water fluxes, and exchange with adjacent
habitats. Production was calculated using the DOC remineralization rate and the C:N ratio of
dissolved organic matter. TDIN was removed from the water column through uptake by
phytoplankton in all habitat models and by Zostera marina in the vegetated subtidal habitat
model. TDIN was exchanged vertically between the sediment and the overlying water column,
based on rates determined from core incubations (Buzzelli 1996).
Model Simulation and Analysis
Sensitivity Analysis: A number of factors potentially influence the simulated behavior of
the state variables (Buzzelli et al. 1995). The sensitivities of the model state variables to the
integration interval (dt), integration routine (Euler vs Runge-Kutta), boundary conditions, and
various parameters were investigated using a systematic series of model trial runs. The interval
for Eulerian integration (dt) was initially set at 0.0625 day (1.5 hours). Analyses included
determining variation in a particular state variable over successive years of the same model run, as
well as the comparison of Year 2 results among a series of different sensitivity runs. The
integration interval (dt) was halved during successive calibration runs to check the effects of dt on
the water column concentrations in each habitat model. Euler and Runge-Kutta integration
routines were compared at similar values of dt. Four to six individual parameters (i.e. rate
coefficients) were selected based on their ecological significance for each state variable to analyze
their effects on the model behavior. Year 2 simulations results were used in the analyses to avoid
initial condition effects. Each parameter was varied by ±10% in separate model runs and the root
mean square deviation (RMS) between the stable, nominal model case and the sensitivity run was
calculated (Cerco 1993):
(2)
where P; = model nominal run; O;= sensitivity run; and n = number of dt in Year 2 simulation
(n=5840). The RMS was compared with the average mean concentration of the nominal run. In
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 27
-------
the cases of the carbon state variables ofZostera marina and Spartina alterniflora, the potential
interactions among two or three varied parameters were investigated for Year 2 output.
Validation: Validation data specific to the Goodwin Islands reference site were available
for only several model state variables. As for the sensitivity analysis, model validation was
addressed by comparing Model Year 2 output of water column chlorophyll a, total paniculate
organic carbon (TPOC), and TDIN from the nonvegetated and vegetated subtidal habitat models
with monitoring data available for the lower York River. Shoot, root-rhizome, and epiphyte
carbon state variables ofZostera marina were compared with data collected at the Goodwin
Islands by Moore et al. (1994). Spartina altemiflora shoot and root-rhizome carbon biomass
were validated using data taken from the literature (including Mendelssohn 1973; Smith et al.
1979; Ornes and Kaplan 1989; Gross et al. 1991). There were no data available to validate model
output for simulated patterns of littoral zone water column DOC, sediment microalgal production
and biomass, and habitat-specific and inter-habitat variations in sediment-water and lateral
(channel boundary condition) material exchanges.
Model Application
The modeled processes were integrated for Year 3 simulation results to derive annual
rates for various, ecologically significant processes and to provide estimates for comparison with
information derived through our own field-laboratory studies and the published literature.
Processes included were: phototrophic net production, phototrophic nitrogen demand and uptake,
and the exchanges of total phytoplankton, TPOC, DOC and TDIN. Annual carbon and nitrogen
dynamics were compared among the autotrophs. Annual carbon production and nitrogen demand
were compared among the four habitats. The annual material exchange that each habitat had with
each of its two adjacent boundary habitats were then summed to derive an annual import or
export estimate for the individual habitat. The annual net carbon production and suspended
material budgets for the entire littoral zone of the Goodwin Islands NERR were then calculated
using the summed process estimates for each habitat. These estimates were compared with those
derived from studies on other mid-Atlantic coast marsh ecosystems.
Model Simulation Results
Model Sensitivity Analysis and Validation
The model was tested extensively for the effects of numerical integration routine and
integration interval, initial conditions and rate coefficient values both singularly and in various
logical combinations. Potential sources of error in simulation analyses and sensitive or controlling
parameters that might require higher degrees of confidence in their estimation were identified.
The results of these multiple runs were analyzed by sensitivity analysis and are discussed in detail
in Buzzelli et al. (1995) and Buzzelli (1996). It is beyond the scope of this report to present these
results but in summary the model was optimized for integration routine and interval. The most
sensitive parameters with regard to effecting relatively large changes in model output (i.e.
predictions of state variable behavior and concentrations) were rate coefficients associated with
the primary producers and, in particular, those coefficients governing the interactions of light and
photosynthesis, translocation and aboveground mortality. For the parameters identified as
28 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
sensitive, the values used for simulation were those derived from the best available data and
values that resulted in model behavior consistent with field and/or laboratory observations.
Following these extensive model sensitivity tests, a model version was decided as the test
case and validated against extant data available from a variety of sources. Model validation was
carried out by comparing simulation output with observed seasonal behavior and concentrations
of the principal state variables. Each is discussed separately below.
Subtidal Water Column
Concentrations: The
modeled
concentrations of
chlorophyll a, total
POC (labile +
refractory), and TDIN
in the water column of
the nonvegetated and
vegetated subtidal
habitats were validated
using data collected
during intensive field
studies at the Goodwin
Islands NERR (Fig. 2;
Moore etal. 1994).
The intensive field
studies were conducted
7- 17 June 1993.
Figure 2A, 2B, and 2C
depict the relationships
between the field data
and concentrations
output from the VST
model. VST model
chlorophyll a was
approximately 5
mg/m3; whereas, the
field data were varied
and ranged
0 Goodinn Islands Water Column Data (Juno 1993)
Goodvin Islands Modal Output (June 1993)
V egota* d Subtidal Habitat NonVe getated Subtidal Habitat
(Stirit]oji2) (Stift MX 4}
CO*""25
op 20
£>
^ 1S|
§. 10
§
a 5
n
r fAVfc 5"r fD\ "I'S
a (AT* I*'/
>. 4
s, , * *: «. i t°0 ;0 ?
"*' " - "*. -30yi^q^I "82:
. »;.**;.. ; ' LoPvWi^Unft,d^^
*0
MAAwwt/vw»AMMA/WW
i i i i i i i i i i i i «
"vvwvuyivyvvv §
0 dSr "
- 1
1 1 1 1 1 1 1 1 1 1 1 L 0
8 9 10 11 12 13 14 15 16 17 1$ 19 89 10 11 12 13 14 15 If 17 18 19
7 (B\ fr (E)
(*> 6
* S
&<
<^
0 3
CU
1 *<
F-. 1
0
0
6
0 5
4
« B ^
°o° .9
1 AAi .oiAtitol iilAI
Wtoto^^ *
tvrt~w# j4rj«^pvuwy mynf1* v v u» j
i i r i i i i i 1 i i 1 0
1, t
0 ' »
0
0 "*," *
'On. 4 O
*° .
0* 1) » t
5 » " * °1 *
^^^^^^^j^^t
vvw^OvwVVvAAMfVV
i i i i i i i i i I i i
8 9 10 11 12 13 14 IS 16 17 18 19 89 10 11 12 13 14 IS 16 17 18 19
6
5
|<
2 *
£3 3
Q
3 2
,0
^ i
Q
(C) "I" (?) .,
15 rj, ihUfllilBlll
- ' '. "WtM^Mim
. <* s> Q
-« « 9
./"P- '"'
riruri t BB nJO 0 nfnf
iTO\0 ' ^°.B« wMr *
^,^v^j^y,
» «
o
« , "'
. $ *H o*.0,,'^
*^*B1 BqflJ»*».« «B0DD%>
T *l T«n»| | | | °?i dwP I I
8 9 10 U 12 13 14 15 16 17 18 19 89 10 11 12 13 14 IS 16 17 18 19
Figure 2
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 29
-------
Shoot BJoroass
between 5 and 25 mg m"3 (Fig. 2A). Vegetated subtidal model concentrations of TPOC ranged
between 1 and 3 gC m"3 and were within the range of values recorded in the field (Fig. 2B).
Water column TDIN computed by the VST model were within the range of field data during the
first few days of simulation, but declined to very low values beginning around 11 June (Fig. 2C).
There was some variability in the concentrations measured in the field (0-5 //M). The
comparisons between model output and field data for the nonvegetated subtidal habitat are shown
in Figure 2D, 2E, and 2F. Model chlorophyll a was approximately an order of magnitude lower
than field determinations (note different y-axis scales in Fig. 2D). NVST model chlorophyll a
concentrations were slightly lower than those output from the VST model; whereas,
concentrations measured at the
Goodwin Islands were slightly
increased in the offshore nonvegetated
habitat (Fig. 2A and 2D). NVST
model TPOC concentrations were
slightly lower than those determined in
the field (Fig. 2E); however, both
model output and field data were
similar among vegetated and
nonvegetated habitats (Fig. 2B and
2E). NVST model TDIN
concentrations ranged from 10-16 gC
m"3; field observations ranged from 0
to 8 gC m'3 (Fig. 2F). The NVST
model concentrations were
considerably greater than those output
from the VST model, although field
concentrations between the two
subtidal habitats were similar (Fig. 2C
and 2D).
<$) Zostera marina Root-Rhizame Biomass
[ | Goadvria tkuds ]>»u
Zost«rt.Mgdtl Output
(C) Zosten marina Epiphytic Biomass
Figure 3
Zostem marina Biomass:
Comparison of modeled Zostera
marina shoot, root-rhizome, and
epiphytic biomass with field
observations are shown in Figure 3.
The validation data were collected at
the Goodwin Islands NERR (Moore et
al. 1994). The model sufficiently
represented the annual patterns in the
biomass of these three state variables.
Although the model predicted summer
shoot biomass of approximately 30 gC
m"2, actual shoot biomass was below 20
gC m'2 (Fig. 3A). Predicted root-
rhizome biomass was consistent with
field data, except for the large peak in
30
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
biomass recorded at the Goodwin Islands NERR in April 1993 (Fig.SB; Orth and Moore 1986).
Although there were few data collected for epiphytic biomass at the Goodwin Islands NERR,
model output was within the range reported and agree with other data collected in the York
River, Virginia (Murray 1983; Neckles 1990).
OO
'e
o
00
250 -
200 -
150--
100--
(&)Sjpsrtu2&
Shoots
fj = Me*a«iss(*» Thesis 1973
^ = Omes.fr Kaplan 198 S>
= Spaitiiia Moid
6 7
Month
8 9 10 11 12
Root-Rhizomes
800
7001
300-
200-
100-
0-
1
Q = Smitk et d. 1979
A = Gross «ti. 199i
= GI SpartinaModrf
i i i i i
i i
2345678
Month
l i i i
9 10 11 12
Figure 4
with latitudes similar to Chesapeake Bay (Smith et al. 1979; Gross et al. 1991).
Spartina alterniflora Biomass: The
model was calibrated using field data
collected at the Goodwin Islands NERR.
The annual patterns in shoot and root
rhizome biomass of Spartina alterniflora
generated by the model were validated with
data assembled from the literature (Fig. 4).
Shoot biomass was compared with data
from the York River, Virginia, and from
South Carolina (Mendelssohn 1973; Ornes
and Kaplan 1989). The root-rhizome output
was validated with data collected in New
Jersey and Delaware (Smith et al. 1979;
Gross et al. 1991). Shoot carbon biomass
was initialized at 3 gC m"2 and stayed low
until the spring pulse of carbon translocated
from below-ground stocks (Fig. 4A). Root-
rhizome carbon biomass was initialized at
635 gC m"2 and decreased in April due to
the upward carbon translocation to shoots
(Fig. 4B). Shoot and root-rhizome carbon
biomass increased through May and June.
Shoot biomass continued to increase,
reaching a maximum of 160 gC m"2 by early
September. The root-rhizome biomass
declined during the summer due to increased
below-ground respiration, which resulted
from higher temperatures (Fig. 4). Shoot
carbon biomass showed a precipitous decline
in the fall, because carbon was translocated
below-ground to the root-rhizome pool as
both state variables returned to their initial
values. Shoot carbon biomass predicted
from the model agreed with field data from
South Carolina (Ornes and Kaplan 1989).
Root-rhizome carbon biomass was within
the range of data reported for other marshes
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 31
-------
Model Predictions of Littoral Zone Processes and Exchanges
Productivity: Annual production by phytoplankton state variables of the Goodwin Islands
NERR habitat models was estimated at 66.0 gC m"2 (Table 2).
Table 2. Estimates of annual net production and contribution to ecosystem production in the littoral zone of
the Goodwin Islands NERR using the four habitat models. Phytoplankton productivity was summed over all 4
habitats and intertidal habitat size used in this summation is the average areal inundation during model
simulation time (m2). The habitats are nonvegetated subtidal (NVST), vegetated subtidal (VST), nonvegetated
intertidal (NVIT), and vegetated intertidal (VIT).
riiutuauiuu oumw
Component
Phytoplankton
Sed. Microalgae
NVST
VST
NVIT
VIT
Zostera marina
Epiphytes
Shoot
RR
Spartina alterniflora
Shoot
RR
TOTAL
Annual Net
Production
gC m-2 yr1
66.0
127.6
101.2
169.0
162.5
55.9
241.3
54.2
830.8
319.7
Habitat Size
104m-2
671
420
120
100
85
120
120
120
85
85
Annual
Net Production
lO'gCyr1
442.7
535.9
121.4
169.0
138.1
67.1
289.6
65.0
706.2
271.7
2806.7
Percent of Total
Ecosystem
%
15.8
19.1
4.3
6.0
4.9
2.3
10.3
2.3
25.2
9.7
99.9
The nonvegetated and vegetated subtidal areas were added to the average inundated area of
each of the two intertidal habitats in order to calculate the total ecosystem size for phytoplankton
production (671 m2). Phytoplankton production rates are listed in Table 2.
The Zostera marina community included productivity from the shoots, attached epiphytes,
and the root-rhizomes. Zostera marina epiphytes and root-rhizomes produced at a similar rate of
approximately 55 gC m"2 yr"1 (Table 2). The Zostera marina community of the Goodwin Islands
NERR produced approximately 421.7 x 106 gC yr"1. The shoots of Spartina alterniflora had the
greatest annual net productivity of any of the model autotrophs.
32
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Nitrogen demand: Nitrogen demand of each phototroph was calculated using the net carbon
production rate and the optimal C:N ratio. Nitrogen uptake was calculated for only the
phytoplankton and the two plants, Zostera marina and Spartina alterniflora. Nitrogen uptake
was calculated for the plants and phytoplankton state variables of each habitat model using
Michaelis-Menten kinetics. There were no formulations to represent nitrogen uptake by sediment
microalgae, although dissolved inorganic nitrogen was exchanged vertically within each habitat
model, based on empirical data (Buzzelli 1996). Table 3 summarizes the annual nitrogen demand
and uptake by each of the phototrophic components of the Goodwin Islands NERR habitat
models. The annual carbon production and nitrogen demand of each of the autotrophs in the
habitat models was calculated to compare the four different littoral zone habitats (Table 4).
Flux: The four habitat models were used to estimate the annual net material fluxes for each
habitat as well as the entire littoral zone of the Goodwin Islands NERR (Table 5). The water
column constituents included total phytoplankton (gC yr'1), TPOC (gC yr'1), DOC (gC yr'1), and
TDIN (mmoles N yr"1). Net import was designated as a negative flux and net export was shown
as a positive flux. Values are provided in Table 5. To assess the interactions between the
Goodwin Islands littoral zone and the surrounding estuary, the annual total exchanges were
summed among the habitats. The totals that were calculated using the four habitat models
provided annual imports of phytoplankton C (-5.9 x 107gC), TPOC (-2.7 x 108 gC), and TDIN (-
1.4 x 109 mmoles N) and an annual export of DOC (1.5 x 108) for the littoral zone of the
Goodwin Islands NERR.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 33
-------
Table 3. Estimates of annual nitrogen demand and uptake for estuarine autotrophs using the Goodwin
Islands habitat models. Demand is calculated using the net carbon production and the optimal C:N ratio.
Uptake is calculated using a Michaelis-Menten relationship based on external nitrogen concentration, a half-
saturation value, and the maximum uptake rate. Phytoplankton nitrogen processes were summed over the
four separate habitat models. The habitats are nonvegetated subtidal (NVST), vegetated subtidal (VST),
nonvegetated intertidal (NVTT), and vegetated intertidal (VIT).
Photoautotrophic
Component
Phytoplankton
Sediment Microalgae
NVST
VST
NVIT
VIT
Zostera marina
shoots
root-rhizomes
total
Spartina alterniflora
shoots
root-rhizomes
total
Annual Nitrogen Demand
gN m"2 yr"1
11.5
22.4
17.8
29.6
28.5
15.1
0.89
16
26
1.53
27.5
Annual Nitrogen Uptake
gN m"2 yr"1
15.7
na
na
na
na
2.09
3.86
5.95
na
11.5
11.5
34 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Table 4. Estimates of net annual carbon production and nitrogen demand of each of the four littoral zone habitats of
the Goodwin Islands NERR using the four habitat simulation models. The habitats are nonvegetated subtidal (NVST),
vegetated subtidal (VST), nonvegetated intertidal (NVIT), and vegetated intertidal (VIT). Each habitat model includes
diatoms, other plankton, and sediment microalgae. In addition to algae the vegetated subtidal and intertidal habitat
models include the net shoot and root-rhizome production by Zostera marina and Spartina alterniflora, respectively.
Habitat
NVST
VST
NVIT
VIT
Size
(ha)
420
120
100
85
Percent of
Total
Size
51.9%
18.5%
12.3%
11.1%
Annual C
Production
gC
740 xlO6
562x1 0s
170 xlO6
1116xlOs
Percent of
Total C
Production
28.6%
21.7%
6.6%
43.1%
Annual N
Demand
gN
130xl06
44 x 106
30 x 10s
47xl06
Percent of
Total N
Demand
51.7%
17.4%
11.9%
19.0%
Table 5. Estimates of annual material exchanges for the four littoral zone habitats of the Goodwin Islands NERR
using the four habitat simulation models. The habitats are nonvegetated subtidal (NVST), vegetated subtidal (VST),
nonvegetated intertidal (NVIT), and vegetated intertidal (VIT). The exchanges of phytoplankton carbon, total
particulate organic carbon (TPOC), dissolved organic carbon (DOC), and total dissolved inorganic nitrogen (TDIN)
between a habitat and its two adjacent boundaries were integrated annually and summed to calculate net import (-)
or export (+).
NVST
VST
NVIT
VIT
TOTALS
Phytoplankton
(gCyr1)
-3.9 xlO7
-1.4 xlO7
-4.5x1 0s
-1.4 xlO6
-5.9 xlO7
TPOC
(gCyr-1)
-4.7 xlO7
-1.7x10"
-4.7 xlO7
-1.4xl07
-2.7 xlO8
DOC
(gCyr-1)
1.4 xlO8
2.4 xlO7
-l.OxlO7
-l.OxlO7
1.5 xlO8
TDIN
(mmolesNyr1)
-l.lxlO9
-2.2 xlO8
-4.7 xlO7
-1.5 xlO7
-1.4 xlO9
Discussion
This ecosystem processes modeling approach was developed to simulate ecosystem dynamics
of littoral zone habitats characteristic of mesohaline and polyhaline shallow water areas of the
lower Chesapeake Bay (Buzzelli 1996). The hydrodynamically coupled models simulate primary
production and water column processes, integrate field and geographic data, and link distinct
aquatic habitats with water quality and living resources of the estuary. The models also provided
a framework to assemble available data, identify missing information, estimate ecosystem and
habitat productivity, and investigate the potential impacts of altered environmental factors on
ecosystem dynamics in the Chesapeake Bay littoral zone.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 35
-------
Model Sensitivity and Validation
Although not presented here, the extensive sensitivity analyses were carried out with the
littoral zone model (see Buzzelli et al. 1996; Buzzelli, 1996). Computed -concentrations of the
various water column constituents for the intertidal habitat models (VIT and NVIT) proved very
sensitive to changes in the integration interval (dt) and were caused by the interactions between
the choice of dt and tidal inundation (i.e. water volume exchange between adjacent habitats).
Water column concentrations (gC or mmole N m"3) were calculated using the change in volume
that resulted from exchanges between adjacent habitats. Because the marsh was not inundated
over the entire tidal cycle, a large dt caused very large and sudden changes in flooded area and
tidal prism volume. The numerical effects were mitigated when dt was reduced to time scales
consistent with those that regulated changes in tidal height (i.e. minutes). A small dt created
smoother hypsometric and volume relationships used to calculate marsh inundation and tidal
volume. Based on considerations of model complexity and output versus computation time, an
integration interval of 11.25 minutes (0.0078125 d) was chosen as the time step for the intertidal
habitat models.
The computed concentrations for diatoms, labile-POC, and sediment microalgae during Year
2 of simulation in the vegetated subtidal model were very robust with respect to 10% deviations
in key controlling parameters. The majority of mathematical expressions governing these state
variables have been calibrated and utilized for a number of years (Cerco and Cole 1994; Kuo and
Park 1995). In most cases, the computed concentrations of water column chlorophyll a, TPOC,
and TDIN by the nonvegetated and vegetated subtidal models were consistent with data recorded
at the Goodwin Islands NERR (Moore et al. 1994) and are within the range of longer-term
measurements made in the lower York River (Batuik et al. 1992). The primary exceptions were
for predicted chlorophyll a concentrations. Model output was an order of magnitude less than
field data for the nonvegetated subtidal habitat. The comparatively low values of chlorophyll a
generated by the NVST model resulted from the very large volume used to calculate the
concentrations. Model chlorophyll a concentrations were lower than those predicted for the
surface waters of the mainstem Chesapeake Bay (10-20 mg m"3; Cerco 1993). The mass of
diatoms and other phytoplankton in each of the habitat models were greatly influenced by the
inter-habitat exchanges, because the magnitudes of the exchange rates were much greater than
those associated with the associated production and loss terms. The TPOC concentrations from
the Goodwin Islands subtidal habitat models were similar to those reported in Cerco (1993). The
TDIN concentrations computed by the subtidal models were within the range of the surface and
bottom values given in Cerco (1993).
Model simulation ofZostera marina shoot, root-rhizome, and epiphytic biornasses were also
fairly robust relative to sensitivity analysis, although epiphytic biomass could change by 40% if its
basal metabolic rate was increased or decreased by 10%. The model replicated the annual
changes in Zostera marina biomass and was used to estimate net annual primary production for
eelgrass meadows of lower Chesapeake Bay. The equations that represented Spartina
alterniflora were highly parameterized. The shoot and root-rhizome carbon biornass was
sensitive to changes in shoot maximum photosynthetic rate (P,^, the root-rhizome basal
respiration rate, and the carbon translocation potential. The connectivity between above- and
below-ground carbon pools was demonstrated by the effects of these three parameters on both
shoot and root-rhizome carbon state variables. Net production was translocated downward, a
pulse of carbon was translocated upwards in the spring, and a large fraction of shoot carbon
remaining in the fall was translocated to the root-rhizomes. P^ appeared to be the most sensitive
36 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
parameter. Values calculated from the literature varied with methods, geographic locations, and
conversion units and ranged from 0.01-0.36 d"1. The maximum rate of 0.15 d"1 used in this study
was the average value calculated from the other studies ofSpartina alterniflora primary
production which are summarized in Table 6.
Table 6. Comparison ofSpartina alterniflora maximum photosynthetic rates (d"1) calculated from
literature sources. The research method referenced in the literature source is provided. A 12-hour day
was used to convert between hourly and daily rates.
METHOD RATE (d"1) SOURCE
Gas flux chambers
Gas flux chambers
Gas flux chambers
Curve fit from growth study
Gas flux chambers
Gas flux chambers
Nitrogen uptake experiments
Goodwin Islands model
o.o r
0.13"
0.04°
0.26d
0.36'
0.06f
0.36s
0.1 5h
Blum et al, 1978
Giurgevich and Dunn, 1 979
Drake and Read, 1981
Morris, 1982
Morris et al., 1984
Pezeshkie/a/., 1987
Morris and Bradley, 1990
This study
"Estimated using 0.4 gC gdw"1 and 1045 gdw m'2.
'Estimated empirically from data provided.
"Estimated using 0.4 gC gdw"1 and 500 gdw m"2 for a Spartina patens community.
""Estimated assuming 30 °C
'Estimated using 0.43 gC gdw1
"Estimated using 0.4 gC gdw"1 and 900 gdw m"2
*Estimated using 0.006 gN gdw"1 root-rhizome tissue
""Average calculated from other studies listed for use in Goodwin Islands model
Model Simulation of Littoral Zone Processes
Primary Production: One principal objective of this modeling exercise was to estimate the
annual rates of net primary production by phytoplankton, sediment microalgae, Zostera marina,
and Spartina alterniflora in the Goodwin Islands NERR (Table 2). The net annual rate of
phytoplankton production (66.0 gC m"2 yr'1) accounted for 15.8% of total annual ecosystem
production (Fig. 5A) and was within the range of values reported in the literature (Table 7). The
annual patterns of chlorophyll a biomass computed using the subtidal habitat models were similar
to long-term patterns evident in data collected in the lower York River, Virginia (Batuik et al.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 37
-------
(A) Goodwin Island NHRR Ecosystem Primary Production
Phylopltnklon
IS !*
Zoftira martiv
U 9*
(B) Nwlh Mel. SC Eosystcn Pin «ry Production
Toilwra nur
13 5%
Sedeicnl
Micrralga
1992). For comparison, annual net
productivity rates for other water bodies
are given in Table 7.
The net annual productivity of
sediment microalgae predicted by the four
habitat models of the Goodwin Islands
NERRranged from 101-169 gC in2 yr'1
(Table 2) and accounted for 34.3% of the
total annual littoral zone production
(Figure 5A). The rate in the nonvegetated
intertidal habitat (NVIT) was greater than
that of the other three habitats. This
resulted from the combined effects of
reduced light attenuation, due to the
depth of the overlying water column in
the NVST and VST habitats and sediment
shading by the plant canopy in the VST
and VIT habitats. Light attenuation due
to overlying water was reduced in the
NVIT habitat, because it was inundated
only 46% of the time over a simulated
annual cycle (21,505 of 46,720 time
steps). The effects of canopy shading are
particularly evident in the differences
between the productivity in the deeper
sand habitat (NVST; 127.6 gC m"2 yr'1)
relative to the shallower SAV habitat
(VST; 101.2 gC m'2 yr1). Although
sediment microalgal productivity
estimates vary with geographic location and habitat, the rates estimated using the Goodwin
Islands habitat models were in overall agreement with those calculated from the biomass data
collected as well as literature values (Table 7).
Zostera marina shoot net annual productivity predicted by the VST model was 241.3 gC m"2
yr"1, approximately four times that computed for the epiphytes (55.9 gC m"2 yr'1) or root-rhizomes
(54.2 gC m"2 yr"1; Table 2). Zostera marina community productivity accounted for about 15% of
the total production in the littoral zone of the Goodwin Islands NERR (Figure 5A). The annual
biomass curves for the three carbon state variables related to Zostera marina were similar to field
data collected in the Goodwin Islands SAV meadow and were within the range of long-term data
for the lower York River, Virginia (Orth and Moore 1986).
Spartina alterniflora annual shoot and root-rhizome biomass changes predicted using the
model agreed with literature values. Model estimates of primary production were similar to those
calculated using Goodwin Islands biomass data. Spartina alterniflora shoot and root-rhizome
productivity were estimated at 830.8 and 319.7 gC m"2 yr"1, respectively, and these rates were
similar to the short form shoot and root-rhizome annual productivity predicted by Dai and
Wiegert (in press) using a canopy model (749 and 397 gC m'2 yr"1; Table 7). The similarities
among the model of the Goodwin Islands Spartina alterniflora and those estimated for short form
Figure 5
38
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Table 7. Summary of annual net production rates (gC m'2 yr"1) taken from published literature. Estimated using
linear regression equation provided,2 Averaged from values provided. R/R in table denotes below ground
production by roots and rhizomes.
Phototroph/Location
Phytoplankton : Chesapeake Bay
Narragansett Bay
Narragansett Bay
Neuse River, NC
North Carolina Estuaries
Goodwin Islands Models
Sediment Microalgae: Mudflat in England
Subtidal in Denmark
Marsh in Mississippi
Mudflat in Massachusetts
Seagrass Mississippi
Marsh in South Carolina
Goodwin Islands Models
Zostera marina : Shoots in Massachusetts
Shoots in Netherlands
Modeled Shoots
R/R in Netherlands
R/R in North Carolina
Modeled R/R
Spartina aherniflora : Shoots in South Carolina
Shoots in Georgia
Modeled Shoots
R/R in South Carolina
R/R in Georgia
R/R in Virginia
R/R in New Jersey
Modeled R/R
Annual Rate
20.26'
101.61
91.252
373.4
52-500
66.0
143.0
89.0
57.4
250.0
339.0
55-234
101-169
155-345
160-412
241.3
53-132
55-102
54.2
289-875
749-1421
830.8
945-2178
397-872
270-857
880.0
319.7
Literature Source
Maloneefa/. 1986
Keller 1989
Keller 1988
Boyererai 1993
Mallin 1994
This Study
Joint 1978
Colijn and DeJong 1984
Sullivan and Moncreiff 1988
Gould and Gallagher 1990
Daehnick et al 1992
Pinckney and Zingmark 1993
This Study
Roman and Able 1988
van Lent and Verschuure 1994
This Study
van Lent and Verschuure 1994
Kenworthy and Thayer 1984
This Study
Dame and Kenny 1986
Dai and Wiegert in press
This Study
Dame and Kenny 1986
Dai and Wiegert in press
Blum 1993
Smith etal. 1979
This Study
Spartina aherniflora from Georgia (Dai and Wiegert, in press) resulted primarily from the
inclusion of seasonal cycles of internal carbon translocation in both. Spartina aherniflora whole
plant production accounted for almost 36% of the total ecosystem production in the Goodwin
Islands littoral zone (Fig. 5A).
Nitrogen demand: The annual Goodwin Islands phytoplankton nitrogen demand was
estimated to be 11.5 gN m"2, based on a C:N weight ratio of 5.7 (Table 3). The annual
phytoplankton nitrogen uptake rate was estimated to be in excess of nitrogen demand, at 15.7 gN
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 39
-------
m"2. This disparity resulted because, unlike Zostera marina and Spartina alterniflora, there were
no mechanisms in the model that limited nitrogen uptake as a function of internal C:N ratio. This
difference may reflect potential luxury nitrogen uptake by phytoplankton. The differences in the
nitrogen requirement of sediment microalgae among the four habitat models resulted from the
differences in the net annual carbon productivity (Tables 2 and 3). Although nitrogen uptake by
sediment microalgae was not explicitly modeled, the models included a vertical exchange of TDIN
between the water column and sediment, based on empirical data collected in subtidal and
intertidal habitats (Neikirk 1996). These field studies measured sediment-water column
exchanges only during the daytime. Other studies are being conducted to determine the diel
variability of sediment-water biogeochemcial fluxes in littoral zone environments of lower
Chesapeake Bay (K. A. Moore and I. C. Anderson, Virginia Institute of Marine Science,
Blacksburg, VA).
Nitrogen was taken up from the water column by the shoots and from the sediments by
the root-rhizomes of Zostera marina. Other studies determined that the sediment was the primary
source of nitrogen for eelgrass (lizumi and Hattori 1982; Short and McRoy 1984). Nitrogen was
translocated from root-rhizomes to the shoots in order to meet the shoot nitrogen requirement for
growth in the Goodwin Islands model (Buzzelli, 1991). Nitrogen uptake by the shoots and root-
rhizomes was influenced both by the external concentration and by feedback limitation terms
based on the maximum and minimum C:N ratios of the tissues. The difference between the annual
nitrogen demand of Zostera marina (16.0 gN m"2 yr'1) and the annual nitrogen uptake (5.95 gN m"
2 yr"1) was attributed to the role of translocation and internal recycling (Table 4). Based on the
Goodwin Islands model, approximately 63% of the plant nitrogen requirement was met through
internal recycling. This value was within the range of annual estimates made by Borum et al.
(1989; 64%), but was approximately twice the short-term rates of translocation measured by
Buzzelli (1991; 34%). Later refinements to this model will include bi-directional nitrogen
translocation within individual plants, as well as carbon and nitrogen translocated from adjacent
root-rhizomes connected below-ground.
Whole plant nitrogen demand and root-rhizome uptake of Spartina alterniflora was
calculated using the vegetated intertidal marsh model. A similar approach to that used to model
the nitrogen relationships of Zostera marina was adopted for the shoot and root-rhizome nitrogen
state variables of Spartina alterniflora, except that there was no nitrogen uptake by the shoots .
As for eelgrass, the whole-plant nitrogen requirement for growth of Spartina alterniflora (27.5
gN m"2 yr'1) was in excess of nitrogen taken up by the plant (11.5 gN m"2 yr"1; Table 3).
Approximately 58% of the plant nitrogen requirement was met through internal recycling, which
agrees with the 54% estimated from studies in a Georgia marsh (Hopkinson and Schubauer
1984). Further field and laboratory studies should include the determination of the actual short-
and long-term rates of carbon and nitrogen uptake and translocation in Spartina alterniflora. A
refinement to the model was the inclusion of bi-directional translocation of nitrogen to
synchronize with seasonal carbon translocation.
Habitat Relationships: Despite the fact that the VIT was the smallest habitat, the annual
production by phytoplankton, sediment microalgae, and Spartina alterniflora (1116x 106gC)
accounted for 43.1% of total in the littoral zone of the Goodwin Islands NERR (Table 4). Over
80% of the intertidal primary production and 34.1% of the total for the littoral zone was
attributable to Spartina alterniflora. This characteristic explained the comparatively low fraction
40 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
of the total ecosystem nitrogen demand required by the VTT (Table 4), because the C:N ratio of
Spartina alterniflora shoots and root-rhizomes was 7-10 times that of the phytoplankton or
sediment microalgae. Conversely, phytoplankton and sediment microalgae primary production in
the NVST was only 28.6% of the total production in the littoral zone of the Goodwin Islands
NERR, although it was the largest of the four habitats (Table 4). The NVST did require 51.7%
of the total littoral zone nitrogen demand, due to the low C:N ratio, compared with the habitats
that include plants. The annual C production by the vegetated subtidal habitat (VST; 562 x 106
gC) was approximately half that of the vegetated intertidal habitat (1116 x 10s gC), but the annual
nitrogen demand and fraction of total ecosystem nitrogen requirement were similar (44 x 106 vs
47 x 106 gN). The nonvegetated intertidal habitat had the least influence on the annual ecosystem
carbon production (6.6%) and nitrogen requirement (11.9%) of the four littoral zone habitats.
Figure 6 A-D depicts the annual net exchanges for each habitat and water column
constituents. An arrow into the habitat denotes a net annual import into the habitat from the
adjacent habitats; an arrow out of a habitat represents a net export of the constituent across its
two boundaries. The nonvegetated and vegetated subtidal models (NVST and VST) both
predicted net annual imports of phytoplankton, particulate organic carbon, and dissolved
inorganic nitrogen; both predicted net annual exports of dissolved organic carbon (Table 5). The
nonvegetated and vegetated intertidal models (NVIT and VIT) predicted net annual imports of all
four water column constituents including dissolved organic carbon (Table 5). The subtidal net
DOC production and export were caused by the increased exudation of the comparatively large
phytoplankton population that was imported (Fig. 6A and 6C). The intertidal net DOC imports
resulted from the decreased exudation and import of phytoplankton, compared with the subtidal
habitat models (Fig. 6A and 6C). Over an annual cycle, the nonvegetated intertidal habitat was
inundated 46% of the time; whereas, the vegetated intertidal habitat was inundated only 25% of
the time. The decreased inundation time and phytoplankton import of the intertidal habitats,
relative to the subtidal habitats, did not translate to decreased TPOC import into the intertidal
habitats (Table 5 and Fig 6B). The vegetated subtidal habitat imported the greatest TPOC
annually (-1.7 x 10s gC); the other three habitats were similar relative to TPOC import (Table 5
and Fig. 6B). All four habitats imported dissolved inorganic nitrogen and the annual TDIN
imported was correlated to the annual phytoplankton mass imported as phytoplankton remove
nitrogen from the water column (Fig. 6A and 6D).
Flux: The exchange of dissolved and particulate materials between coastal marshes and
adjacent environments are governed by their historical geomorphology, basin configuration, and
hydroperiod (Childers et al. 1993; Rozas 1995). The hydroperiod of an individual marsh is
unique and has a significant influence on the horizontal and vertical material exchanges
(Vorosmarty and Loder 1994). Despite the differences in geomorphology and hydroperiod
among marshes, it is useful to compare and contrast the flux characteristics among marshes as a
way to identify spatial or temporal patterns (Childers 1992). The material exchange estimates
generated for the littoral zone of the Goodwin Islands NERR were compared to annual flux
estimates derived from other studies (Table 8).
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 41
-------
0)
M
3)
r-l
42
-------
TableS. Summary of marsh ecosystem flux estimates assembled from published literature. Negative flux
(-) denotes net import to marsh, positive flux (+) denotes net export from marsh.
Location
Carter Creek, Virginia
Carter Creek, Virginia
Size
(ha)
10
10
POCflux
gCyr-1
1.17xl07
na
DOC flux
gCyr1
-2.5 x 106
na
DIN flux
mmoleN yr"1
-2.6 xlO7
-2.0 x 107
Literature
Source
Axelradetal. 1976
Wolavere/a/. 1983
Rhode River, Maryland
Low marsh
13
na
na
-8.6 x 106
Conelletal. 1992
Ely Creek, South Carolina
Marsh 12
-2.1 xlO7 1.8 xlO7 -1.2 xlO7
Dameetal. 1991
Goodwin Islands Marsh
Model
85
-1.4xl07 -l.OxlO7 -1.5xl07
This Study
The Goodwin Islands marsh model predicted a net POC import to the marsh annually
(Table 8). The difference between the POC export at Carter Creek and the POC import predicted
using the Goodwin Islands model was attributed to the absence of an upland connection at the
Goodwin Islands NERR and the proximity of the terrestrial boundary at Carter Creek (Axelrad et
al. 1976). Like the Carter Creek marsh on the York River, Virginia, the Rhode River, Maryland
marsh has an upland connection (Correll et al. 1992).
Conclusions
The Chesapeake Bay Littoral Zone and Wetlands Models employed a series of simulation
models to calculate annual carbon production, nitrogen demand, and water column dynamics in
the littoral zone of the lower Chesapeake Bay. These models also investigated potential change in
habitat and ecosystem properties. Environmental problems that can be explored include: potential
effects of changing water quality on productivity in the eelgrass community; the possible effects
that significant changes in the distribution and abundance of eelgrass might have on primary
production and nitrogen uptake in the subtidal habitats; and the potential effects of changes in
mean sea level on intertidal productivity and material exchange properties.
The current models were designed to be coupled to coarser-scale models of water quality
in the Chesapeake Bay watershed (Cerco 1993). Of course, many important physical and
biogeochemical processes are currently not present in the models. The development of the
sediment state variables and processes and linkages to the overlying water column must be
included to better investigate production and material cycling in shallow and irregularly flooded
littoral zone habitats. Phosphorus dynamics, the contribution of living and dead plants to DOC
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 43
-------
production and exchange, and the nitrogen relationships of sediment microalgae would greatly
improve the biogeochemical portions of the models. Secondary productivity within the different
littoral zone habitats should be included as a vehicle to transfer energy and nutrients between the
autotrophs and higher trophic levels. This would provide additional mechanisms to link the
habitats in time and space, as well as better address other living resources issues facing the
management community (e.g. higher trophic level dynamics and controls).
Data Needs
The dynamics of 37 different state variables can be represented by these four littoral zone
habitat models (Figure 1). The output of only a few of these state variables was validated in this
summary. Although one of the objectives of this modeling project was to organize data relevant
to Chesapeake Bay littoral zone ecology, another was to identify information that was lacking.
Additional data are required in several areas, including: the annual variation in the productivity
and biomass of sediment microalgae in all habitats; the relationships between sediment microalgal
production and the effects of plant canopy shading; the spatial and temporal dynamics of dissolved
organic carbon in shoal waters of lower Chesapeake Bay; the processes of gross and net
photosynthesis, nitrogen uptake, and internal carbon and nitrogen translocation in Spartina
alterniflora; and the horizontal exchange of dissolved and particulate materials between the
vegetated intertidal marsh and the surrounding habitats.
Management Implications and Future Directions
The Lower Chesapeake Bay Littoral Zone and Fringing Wetlands Ecosystem Model links
water quality and living resource dynamics in the littoral zone habitats that are essential to the
survival of important juvenile and adult fishes and invertebrates. This model can investigate the
potential effects of habitat alteration, water quality conditions, relative sea level variation, or
watershed practices on the ecology of the littoral zone. Specific questions that can be addressed
within this modeling framework are:.
What is the relationship between habitat area (i.e., changes in area covered by SAV or marsh)
and water quality or plant production?
What is the impact of short-term, high frequency or seasonal pulses (of a water column
constituent such as total suspended solids, chlorophyll, dissolved inorganic nitrogen, or of
seasonally immigrating grazers and predators) on plant production within SAV and marsh
habitats? Is there a difference in modeled outcome based on the frequency of pulsed events? Is
there a critical time of year for water or habitat quality?
What is the relationship between plant production and area of habitat to the potential trophic
transfer to, or production of, invertebrates and fishes?
How do predictions for littoral zone habitats from the large-scale Chesapeake Bay Program
Water Quality model compare with a small-scale model of SAV production and littoral zone
water quality? If they differ, why and what are the implications of this?
44 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
References
Axelrad, D. M, Moore, K. A., and Bender, M. E. 1976. Nitrogen, phosphorus, and carbon flux in
Chesapeake Bay marshes (OWRT Project B-027-VA Bulletin 79). Blacksburg, Virginia:
Virginia Institute of Marine Science.
Bach, H. K. 1993. A dynamic model describing the seasonal variations ingrowth and the
distribution of eelgrass (Zostera marina L.) I. Model theory. Ecological Modeling 65:31-
50.
Baird, D. and Ulanowicz, R. E. 1989. The seasonal dynamics of the Chesapeake Bay ecosystem.
Ecological Monographs 59:329-364.
Batuik, R. A., Orth, R.J., Moore, K.A., Dennison, W.C., Stevenson, 1C., Staver, L.W., Carter,
V., Rybicki, N.B., Hickman, R.E., Kollar, S., Bieber, S., and Heasly, P. 1992.
Chesapeake Bay Submerged Aquatic Vegetation Habitat Requirements and Restoration
Targets: A Technical Synthesis: U.S. EPA, Chesapeake Bay Program.
Blum, U., Seneca, E.D. and Stroud, L.M. 1978. Photosynthesis and respiration ofSpartina and
Juncus salt marshes in North Carolina: Some models. Estuaries 1:228-238.
Blum, L. K. 1993. Spartina alterniflora root dynamics in a Virginia marsh. Marine Ecology
Progress Series 102:169-178.
Boon, J. D. and Byrne, R. J. 1981. On basin hypsometry and the morphodynamic response of
coastal inlet systems. Marine Geology 40:27-48.
Boumans, R. M. J. and Sklar, F. H. 1990. A polygon-based spatial (PBS) model for simulating
landscape change. Landscape Ecology 4:83-90.
Borum, J., Murray, L., and Kemp, W.M. 1989. Aspects of nitrogen acquisition and conservation
in eelgrass plants. Aquatic Botany 35:289-300.
Boyer, J. N., Christian, R. R., and Stanley, D. W. 1993. Patterns of phytoplankton primary
productivity in the Neuse River estuary, North Carolina, USA. Marine Ecology Progress
Series 97:287-297.
Buzzelli, C. P. 1991. Sediment Inorganic Nitrogen Stocks and Root-Rhizome Ammonium Uptake
by Eelgrass (Zostera marina L.) in the lower Chesapeake Bay: M. A. Thesis, School of
Marine Science, College of William and Mary, Williamsburg, VA.
Buzzelli, C.P. 1996. Integrative Analysis of Ecosystem Processes and Habitat Patterns in the
Chesapeake Bay Littoral Zone: A Modeling Study of the Goodwin Islands National
Estuarine Research Reserve. PhD. Dissertation, School of Marine Science, College of
William and Mary, Williamsburg, VA.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 45
-------
Buzzelli, C.P. and R.L. Wetzel. 1996. Modeling the Lower Chesapeake Bay Littoral Zones and
Fringing Wetlands: Ecosystem Processes and Habitat Linkages. II. Model Sensitivity
Analysis, Validation, and Estimates of Ecosystem Processes. Special Report No. 335 in
Applied Marine Science and Ocean Engineering, Virginia Institute of Marine Science,
Gloucester Point, VA.
Buzzelli, C.P., Wetzel, R.L., and Meyers, M.B. 1995. Modeling the lower Chesapeake Bay
littoral zone and fringing wetlands: Ecosystem processes and habitat linkages. I.
Simulation Model Development and Description. Special Report No. 334 in Applied
Marine Science and Ocean Engineering. School of Marine Science-Virginia Institute of
Marine Science, College of William and Mary.
Cerco, C.F. 1993. Three-dimensional eutrophication model of Chesapeake Bay. Journal of
Environmental Engineering 119(6): 1006-1025.
Cerco, C.F. and T. Cole. 1994. Three-dimensional eutrophication model of Chesapeake Bay:
Volume 1, main report. Technical Report EL-94-4, United States Army Engineer
Waterways Experiment Station, Vicksburg, Mississippi.
Childers, D. L. 1992. Fifteen years of marsh flumes-A review of marsh-water column interactions
in Southeastern U.S. estuaries. INTECOL Fourth International Wetlands Conference,
Elsevier Press.
Childers, D. L., H.N. McKellar, R.F. Dame, F.H. Sklar, and E.R. Blood. 1993. A dynamic
nutrient budget of subsystem interactions in a salt marsh estuary. Estuarine Coastal and
Shelf Science 36:105-131.
Colijn, F., and V.N. deJonge. 1984. Primary production of microphytobenthos in the Ems-Dollard
estuary. Marine Ecology Progress Series 14:185-196.
Correll, D. L., T.E. Jordan and D.E. Weller. 1992. Nutrient flux in a landscape: Effects of
coastal land use and terrestrial community mosaic on nutrient transport to coastal waters.
Estuaries 15:431-442.
Costanza, R., F.H. Sklar, and M.L. White. 1990. Modeling coastal landscape dynamics.
BioScience 40:91-107.
Daehnick, A. E., Sullivan, M. J., and Moncreiff, C. A. 1992. Primary production of the sand
microflora in seagrass beds of Mississippi Sound. Botanica Marina 35:131-139.
Dai, T., and Wiegert, R. G. (in press). Estimation of the primary productivity ofSpartina
alterniflora using a canopy model. Ecography.
46 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Dame, R. F., J.D. Spurrier, T.M. Williams, B. Kjerfve, R.G. Zingmark, T.G. Wolaver, T.H.
Chrzanowski, H.N. McKellar, and FJ. Vernberg. 1991. Annual material processing by a
salt marsh-estuarine basin in South Carolina, USA. Marine Ecology Progress. Series
72:153-166.
Dame, R. F. and P.O. Kenny. 1986. Variability ofSpartina alterniflora primary production in the
euhaline North Inlet estuary. Marine Ecology Progress Series 32:71-80.
de Jonge, V. N. and F. Colijn. 1994. Dynamics of microphytobenthos biomass in the Ems estuary.
Marine Ecology Progress Series 104:185-196.
Dennison, W. C., R.J. Orth, K.A. Moore, J.C. Stevenson, V.C. Carter, S. Kollar, P.W.
Bergstrom, and R.A. Batuik. 1993. Assessing water quality with submersed aquatic
vegetation. Bioscience 43:86-94.
Drake, B. G. and M. Read. 1981. Carbon dioxide assimilation, photosynthetic efficiency, and
respiration of a Chesapeake Bay salt marsh. Journal of Ecology 69:405-423.
Eiser, W. C. and B. Kjerve. 1986. Marsh topography and hypsometric characteristics of a South
Carolina salt marsh basin. Estuarine Coastal and Shelf Science 23:595-605.
Friedrichs, C. T. and D. G. Aubrey. 1994. Uniform bottom shear stress and equilibrium
hypsometry of intertidal flats. In Mixing Processes in Estuaries and Coastal Seas (ed.
C. Pattiaratchi). Washington, D.C.: American Geophysical Union.
Giurgevich, J.R. and E.L. Dunn. 1979. Seasonal patterns of CO2 and vapor exchange of the tall
and short height forms ofSpartina alterniflora Loisel. in a Georgia salt marsh.
Oecologia 43:139-156.
Gould, D. M. and E.D. Gallagher. 1990. Field measurements of specific growth rate, biomass,
and primary production of benthic diatoms of Savin Hill Cove, Boston. Limnology and
Oceanography 35:1757-1770.
Gross, M. F., M.A. Hardisky, P.L. Wolf, and V. Klemas. 1991. Relationship between
aboveground and belowground biomass ofSpartina alterniflora (Smooth Cordgrass).
Estuaries 14:180-191.
Heck, K. L. and T.A. Thoman. 1984. The nursery role of seagrass meadows in the upper and
lower reaches of Chesapeake Bay. Estuaries 7:531-540.
Hopkinson, C. S. and J.P. Schubauer. 1984. Static and dynamic aspects of nitrogen cycling in the
salt marsh graminoid Spartina alterniflora. Ecology 65:961-969.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 47
-------
lizumi, H. and A. Hattori. 1982. Growth and organic production of eelgrass (Zostera marina L.)
in temperate waters of the Pacific coast of Japan. III. The kinetics of nitrogen uptake.
Aquatic Botany 12:245-256.
Joint, I. R. 1978. Microbial production of an estuarine mudflat: Estuarine Coastal and Shelf
Science 7:185-195.
Keller, A. A. 1988. An empirical model of primary productivity (14C) using mesocosm data along
a nutrient gradient. Journal of Plankton Research 10(4):813-834.
Keller, A. A. 1989. Modeling the effects of temperature, light, and nutrients on primary
productivity: An empirical and a mechanistic approach compared. Limnology and
Oceanography 34:82-95.
Kenworthy, W. J. and G.W. Thayer. 1984. Production and decomposition of the roots and
rhizome of seagrasses, Zostera marina and Thalassia testudinum, in temperate and
subtropical marine ecosystems. Bulletin of Marine Science 35(3):364-379.
Kneib, R. T. and S.L. Wagner. 1994. Nekton use of vegetated marsh habitats at different stages
of tidal inundation. Marine Ecology Progress Series 106:227-238.
Kuo, A. Y. and K. Park. 1995. A framework for coupling shoals and shallow embayments with
main channels in numerical modeling of coastal plain estuaries. Estuaries 18:341-350.
Mallin, M. A. 1994. Phytoplankton ecology of North Carolina estuaries. Estuaries 17:561-574.
Malone, T. C., W.M. Kemp, H.W. Ducklow, W.R. Boynton, J.H. Turtle, and R.B. Jonas. 1986.
Lateral variation in the production and fate of phytoplankton in a partially stratified
estuary. Marine Ecology Progress Series 32:149-160.
Mendelssohn, I. A. 1973. Angiosperm production of three Virginia marshes in various salinity and
soil nutrient regimes. M.A. Thesis, School of Marine Science, College of William and
Mary, Williamsburg, Virginia.
Moncreiff, C. A., M.J. Sullivan, and A.E. Daehnick. 1992. Primary production dynamics in
seagrass beds of Mississippi Sound: the contributions of seagrass, epiphytic algae, sand
microflora, and phytoplankton. Marine Ecology Progress Series 87:61-171.
Moore, K. A., J.L. Goodman, J.C. Stevenson., L. Murray, L. and K. Sundberg. 1994. Chesapeake
Bay Nutrients, Light, and SAV: Relations Between Variable Water Quality and SAV in
Field and Mesocosm Studies: U.S. EPA Chesapeake Bay Program.
Morris, J. T. 1982. A model of growth responses by Spartina alterniflora to nitrogen limitation.
Journal of Ecology 70:25-42.
Morris, J. T., R.A. Houghton and D.B. Botkin. 1984. Theoretical limits of belowground
48 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
production by Spartina alterniflora: An analysis through modeling. Ecological Modelling
26:155-175.
Morris, J. T. and P. Bradley. 1990. Influence of oxygen and sulfide concentration on nitrogen
uptake kinetics in Spartina alterniflora. Ecology 71:282-287.
Murray, L. 1983. Metabolic and Structural Studies of Several Temperate Seagrass Communities
with Emphasis on Microalgal Components. Ph.D. Dissertation, School of Marine Science,
College of William & Mary, Williamsburg, VA.
Murray, L. and R.L. Wetzel. 1987. Oxygen production and consumption associated with the
major autotrophic components in two temperate seagrass communities. Marine Ecology
Progress Series 38, 231-238.
Neckles, H. A. 1990. Relative Effects of Nutrient Enrichment and Grazing on Ephiphyton-
Macrophyte (Zostera marina L.) Dynamics. Ph.D. Dissertation, School of Marine
Science, College of William & Mary, Williamsburg, VA.
Neikirk, B.E.B. 1996. Exchanges of dissolved inorganic nitrogen and dissolved organic carbon
between salt marsh sediments and overlying tidal water. M. A. Thesis, School of Marine
Science, College of William and Mary, Williamsburg, VA.
Ornes, W. H. and D.I. Kaplan. 1989. Macronutrient status of tall and short forms of Spartina
alterniflora in a South Carolina salt marsh. Marine Ecology Progress Series 55:63-72.
Orth, R. J. and K.A. Moore. 1986. Seasonal and year-to-year variations in the growth of Zostera
marina L. (eelgrass) in the lower Chesapeake Bay. Aquatic Botany 24:335-341.
Pezeshki, S. R., R.D. DeLaune and W.H. Patrick. 1987. Gas exchange characteristics of Gulf of
Mexico coastal marsh macrophytes. Journal of Experimental Marine Biology and Ecology
111:243-253.
Pinckney, J. and R. Zingmark. 1993 a. Modeling the annual production of intertidal benthic
microalgae in estuarine ecosystems. Journal of Phycology 29:396-407.
Pinckney, J. and R. Zingmark. 1993b. Biomass and production of benthic microalgal communities
in estuarine sediments. Estuaries 16:887-897.
Rizzo, W. M., G.J. Lackey and R.R. Christian. 1992. Significance of euphotic, subtidal sediments
to oxygen and nutrient cycling in a temperate estuary. Marine Ecology Progress Series
86:51-61.
Roman, C. T., and K.W. Able. 1988. Production ecology of eelgrass (Zostera marina L.) in a
Cape Cod salt marsh-estuarine system, Massachusetts. Aquatic Botany 32:353-363.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 49
-------
Roman, C. T., K.W. Able, M.A. Lazzari and K.L. Heck. 1990. Primary productivity of
angiosperm and macroalgae dominated habitats in a New England salt marsh: A
comparative analysis. Estuarine Coastal and Shelf Science 30:35-46.
Rozas, L. 1995. Hydroperiod and its influence on nekton use of the salt marsh: A pulsing
ecosystem. Estuaries 18(4):579-590.
Sand-Jensen, K. and J. Borum. 1991. Interactions among phytoplankton, periphyton, and
macrophytes in temperate freshwaters and estuaries. Aquatic Botany 41:137-175.
Short, F. T. and C.P. McRoy. 1984. Nitrogen uptake by leaves and roots of the seagrass Zostera
marinaL. BotanicaMarina, 27:547-555.
Smith, K. K., R.E. Good and N.F. Good. 1979. Production dynamics for above and belowground
components of a New Jersey Spartina alterniflora tidal marsh. Estuarine. and Coastal
Marine Science 9:189-201.
Spinner, G. P. 1969. Serial Atlas of the Marine Environment. In: The wildlife wetlands and
shellfish areas of the Atlantic coastal zone, vol. Folio 18. New York City American
Geophysical Union.
Stevenson, J. C., L. G. Ward, and M. S. Kearney. 1988. Sediment transport and trapping in marsh
systems: Implications of tidal flux studies. Marine Geology 80:37-59.
Strahler, A. N. 1952. Hypsometric (area-altitude) analysis of erosional topography. Bulletin of
Geological Society of America 63:1117-1142.
Sullivan, M. J. and Moncreiff, C. A. 1988. Primary production of edaphic algal communities in a
Mississippi salt marsh. Journal of Phycology 24:49-58.
Sundback, K. and Graneli, W. 1988. Influence of microphytobenthos on the nutrient flux between
sediment and water: a laboratory study. Marine Ecology Progress Series 43:63-69.
van Lent, F. and J.M. Verschuure. 1994. Intraspecific variability of Zostera marina L. (eelgrass)
in the estuaries and lagoons of the southwestern Netherlands. I. Population dynamics.
Aquatic Botany 48:31-58.
Vorosmarty, C. J. and T.C. Loder. 1994. Spring-Neap tidal contrasts and nutrient dynamics in a
marsh dominated estuary. Estuaries 17:537-551.
Wetzel, R. L. and C.S. Hopkinson. 1990. Coastal ecosystem models and the Chesapeake Bay
Program: philosophy, background and status, pp7-23. IN: M. Haire and E.G. Krome (es.)
Perspectives on the Chesapeake Bay 1990: Advances in Estuarine Sciences.
CBP/TRS41/90, CRC, Inc., Gloucester Point, VA.
50 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Wetzel, R.L. and M.B Meyers. 1994. Ecosystem process modeling of submerged aquatic
vegetation in the lower Chesapeake Bay. Final report (1993) to USEPA Region III
Chesapeake Bay Program Office.
Wetzel, R. L. and H.S. Neckles. 1986. A model ofZostera marina L. photosynthesis and growth:
Simulated effects of selected physical-chemical variables and biological interactions.
Aquatic Botany 26:307-323.
Wolaver, T. G., J.C. Zieman, R.L. Wetzel and K.L. Webb. 1983. Tidal exchange of nitrogen and
phosphorus between a mesohaline vegetated marsh and the surrounding estuary in the
lower Chesapeake Bay. Estuarine Coastal and Shelf Science 16:321-332.
Modeling the Lower Chesapeake Bay Littoral Zone and Fringing Wetlands 51
-------
Trophic Control of SAV Responses to Nutrient Enrichment:
Simulation Modeling of Mesocosm Experiments
R. Bartleson, L. Murray, W. M. Kemp
Horn Point Environmental Laboratory
University of Maryland
Cambridge, MD
Background
Nutrient enrichment of coastal waters promotes growth of planktonic and epiphytic
algae, which tend to inhibit production and survival of estuarine submerged aquatic vegetation
(SAV) and Bay grasses (Kemp et al. 1983; Twilley et al. 1985; Coleman and Burkholder 1995;
Duarte 1995; Short etal. 1995). In many shallow eutrophic habitats, epiphytic algal growth
appears to be of primary importance for regulating abundance of SAV (Borum 1985; Duarte
1995). In some environments, invertebrate grazing can effectively reduce the effects of epiphytes
on SAV. However, the relative effectiveness of herbivorous grazing control on epiphytes tends
to vary with season, region, and feeding mode of the grazer populations (Howard 1982;
Bronmark 1985; Neckles et al. 1993). In addition, significant changes in the mortality of these
grazers, which may result from predation or altered environmental conditions (e.g., Lubbers et al.
1990), can regulate the ability of grazers to control epiphyte growth. Ultimately, differences in
trophic structure of the community associated with submersed plants can radically alter the plant
responses to changes in nutrient levels.
Relating SAV growth and survival to both in situ environmental conditions and
habitat requirements (e.g., nutrient concentrations) and to the effects of nutrient load reduction
strategies (e.g., 40 % reduction goals; Dennison et al. 1993) are of particular interest to
environmental decision-makers. Although a variety of mesocosm experiments (e.g., Twilley et
al. 1985, Neudnorfer and Kemp 1990) have demonstrated the potential importance nutrient
enrichment on Bay grass production, it is difficult to relate experimental nutrient loading rates
and concentrations to conditions observed in the field. Size limitations for most mesocosm
experiments preclude inclusion of full trophic structure as it occurs in nature. Current
experiments in the Environmental Protection Agency (EPA) Multiscale Experimental Ecosystem
Research Center (MEERC) program at the University of Maryland Center for Environmental and
Estuarine Studies (UMCEES) are examining the effects of variations in trophic structure on
susceptibility of SAV populations to stress from nutrient enrichment.
One way to relate results from mesocoms experiments to field conditions is with the
aid of numerical ecosystem process models. Several models of SAV production have examined
responses to nutrient enrichment and grazing regulation of epiphyte growth (Wetzel and Neckles
1986; Kemp et al. 1995; Madden and Kemp 1996). These models have not simulated dynamics
of the grazers themselves and feedback interactions among grazers, epiphytes and SAV. In
addition, no experiments or models have explicitly considered how herbivorous grazing may be
controlled by resident predator populations. This model explores some of the consequences of
interactions between fish and grazers at different nutrient concentrations. Results should
improve understanding of how nutrient concentrations affect the survival of SAV.
52 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Model Development
The current model has several modifications from the original SAV model (e.g.,
Kemp et al. 1995; Madden and Kemp 1996). The model equations were developed from
empirical relationships based on laboratory and field data from mesohaline Chesapeake Bay,
other estuaries, and theoretical functions. The model structure is shown in Figure 1. Because
internal nutrient storage may affect competition between SAV and algae, plant nitrogen (N),
including nonstructural nitrogen, was modeled separately from carbon (C). Nonstructural
carbohydrate pools of SAV biomass were also included in the model to allow realistic
C-allocation. The model was designed to track N, phosphorus (P), oxygen (0), and C through
the system compartments. It has 19 state variables: SAV leaf C and N; SAV root C and N; SAV
nonstructural C and N; epiphytic algae; benthic algae; phytoplankton; amphipods; fish; infauna;
sediment labile organic carbon; sediment refractory organic carbon; water column dissolved
inorganic N and P; sediment porewater dissolved inorganic N and P; and dissolved oxygen. Data
from MEERC mesocosm experiments (e.g., Sturgis and Murray 1997) were used to calibrate the
model.
The model used a time step of 15 minutes to capture diel effects of light on nutrients,
production and respiration. Integration was performed using fourth order Runga-Kutta
numerical integration. The model was calibrated with field data and data collected from 1 m2
mesocosms in a 1995 experiment (Figure 2). The model experiments examined the interacting
effects of trophic
complexity and
nutrient enrichment
(Figure 3). The model
was run at three
nutrient addition rates
(1, 15, and405Ml'1
d"1) for each of three
trophic complexity
conditions: no grazers;
with grazers; and with
grazers and fish.
Grazers (1.5 gC)
were added at week 2,
and fish (2.5 g C)
were added at week 5
at densities in the
mid-range
(1-2 g C m'2) reported
from shallow-water
systems (e.g. Vimstein
etal. 1983, Kemp et
al. 1983).
SAV Mesocosm Modei
From Chopt
Pump
Figure 1.
Gereral structure of SAV model is shown. Some material
flows are left out for simplicity. N is modeled in stoichicmetiy with
Cejocept in SAV where structural and nonstructural N and Care
modeled separately.
Trophic Control of SAV Responses to Nutrient Enrichment
53
-------
1995 Mesocosm Experiment and Model
2°i ; 1 Results and Discussion
Baseline Experiment
The model output
matched the oxygen and
dissolved inorganic nitrogen
data from the 1995
mesocosm experiments quite
well (Fig. 2). Simulations of
community production and
respiration, as well as growth
rate and SAV biomass
measurements, were
comparable with mesocosm
experimental data, but
biomass of epiphytes and
grazers diverged from model results. The model predicted that grazers would control epiphyte
biomass in the high-grazer treatment and that epiphyte biomass would be higher in the low-grazer
treatment.
Simulations of Interactions Between Nutrient Addition
Rate and Number of Trophic Groups
No Grazers Grazers Added Grazers and Fish Added
Low Nutrients
12,
10
Week
Figure 2
6.
SEDORQN
SAV .- - "
. ' SEDORGN
- tiiN' 'Epjphyifc' Xlcw.'' 'Srar»ra
SAV .
SED one N
G'rj airs' .....
Moderate Nurients
12r-
SED OflG N .
SAV,..
DIN eiazeo
SEDORQN-
..'Epiphytic Algw SAV ..
Grucrs
High Nutients
Figure 3
54
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
The measured epiphyte biomass, however, was not significantly different among treatments.
This discrepancy could have resulted from an increased regeneration rate of P in the high-grazer
treatment, which led to enhanced growth of N-fixing epiphytic blue-green algae. It appearred
that these blue-greens were not effectively controlled by amphipod grazing. Phosphorus was not
explicitly included in this model, but expeimental data revealed that P was quite low in the
treatments and control, so regeneration rates would be important in controlling growth rates of
the autotrophs. In another study, blue-green algae dominated the epiphytic algal community
when amphipods and isopods were present in a nutrient enriched eelgrass microcosm (Neckles et
al. 1994).
Trophic complexity experiment
Grazers can significantly increase or decrease nitrogen retention (Fig. 3, panels d and
e) in an SAV mesocosm. This is due to the stimulating effect of grazers on primary production,
through nutrient regeneration, and the negative effect of reducing epiphytic algal biomass to low
levels. When producer biomass is sufficient, increased nitrogen cycling and production should
reduce the rate of outflow. Adding fish, however, increased (fig. 3, panels e and f) nitrogen
cycling with moderate nutrient additions and decreased (panels h and i) nitrogen cycling with
high nutrient additions. Results were dependent on the coefficients picked for the model. For
example, if a lower grazing rate was used in the moderatre nutrient, grazer added scenario (panel
e), epiphyte biomass may have remained high, giving a model result more like the fish addition
scenario (panel f). Sensitivity analysis will show the coefficients that have the most influence on
model behavior.
Using the mesocosms to improve the model
The epiphytic algae did not decrease as expected in the amphipod treatments, because
the filamentous blue-green algae were not noticeably grazed. In the next generation of the
model, blue-green algae will be added as an explicit state variable. Since water-column P was
found in more limiting concentrations than N, it was necessary to add P to the model. This will
improve the model's ability to accurately predict rates of primary production. Fish will be added
to the model to help determine the number that will reduce, but not eliminate the amphipods. In
this way, the model may be used to improve the experiment.
Trophic Control of SAV Responses to Nutrient Enrichment 55
-------
References
Borum, J. 1985. Development of epiphytic communities on eelgrass, Zostera marina,
along a nutrient gradient in a Danish estuary. Mar. Biol. 87: 233-241.
Bronmark, C. 1985. Interactions between macrophytes, epiphytes and herbivores: an
experimental approach. Oikos 45: 26-30.
Coleman, V., and J. Burkholder. 1995. Response of microalgal epiphyte communities to
nitrate enrichment in an eelgrass (Zostera marina) meadow. 31: 36-43
Dennison W.C., Orth R.J., Moore K.A., Stevenson J.C., Carter V., Kollar S., Bergstrom P.W.,
Batiuk R.A. 1993. Assessing water quality with submersed aquatic vegetation: habitat
requirements as barometers of Chesapeake Bay health. Bioscience 43:86-94.
Duarte C.M. 1995. Submerged aquatic vegetation in relation to different nutrient regimes.
Ophelia 41:87-112
Howard, R. K. 1982. Impact of feeding activities of epibenthic amphipods on
surface-fouling of eelgrass leaves. Aquat. Bot. 14: 91-97
Kemp, W. M., W.R. Boynton, R.R. Twilley, J.C. Stevenson, J. C. Means. 1983. The
decline of submerged vascular plants in upper Chesapeake Bay: Summary of results
concerning possible causes. Mar. Technol. Soc. J. 17: 78-89
Kemp, W. M., W. R. Boynton, A. J. Hermann. 1995. Simulation models of an esituarine
macrophyte ecosystem, pp. 262-278, In: B. Patten and S. E. J0rgensen (eds.)
Complex ecology. Prentice Hall, Englewood Cliffs, NJ.
Lubbers, L., W. R. Boynton and W. M. Kemp. 1990. Variations in structure of estuarine fish
communities in relation to abundance of submersed vascular plants. Mar. Ecol. Progr. Ser.
65: 1-14.
Madden, C. J. and W. M. Kemp. 1996. Ecosystem model of an estuarine seagrass
community: Calibration and simulation of eutrophication responses. Estuaries. (In
press).
Neckles, H. A., R.L. Wetzel, RJ. Orth. 1993. Relative effects of nutrient enrichment and
grazing on epiphyte-macrophyte (Zostera marina L.) dynamics. Oecologia 93:
285-295
Neckles, H. A., E.T. Koepfler, L.W. Haas, R.L. Wetzel, R.J. Orth. 1994. Dynamics of
epiphytic photoautotrophs and heterotrophs in Zostera marina (Eelgrass)
microcosms: responses to nutrient enrichment and grazing. Estuaries 17: 597-605
56 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Patten, S. E. Jorgensen and S. Auerbach (eds.) Complex ecology. Prentice Hall,
Englewood Cliffs, NJ.
Short, F., D. Burdick and J. Kaldy. 1995. Mesocosm expeiments quantify the effects of
eutrophication on eelgrass, Zostera marina. Limnol. Oceanogr. 40: 730-749.
Sturgis and Murray 1997. Scaling of nutrinet input to submersed macrophyte
mesocosms. Submitted to Mar. Ecol. Progr. Ser.
Twilley, R. R, W.M. Kemp, K.W. Staver, J.C. Stevenson, 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
Wetzel, R. L. and H. Neckles. 1986. A model of Zostera marina L. photosynthesis and
growth: Simulated effects of selected physical-chemical variables and biological
interactions. Aquat. Bot. 26: 307-323.
Trophic Control of SAV Responses to Nutrient Enrichment 57
-------
A Simulation Model of Submersed Aquatic Vegetation (SAV) and Water
Quality in Chesapeake Bay
Christopher J. Madden and W. Michael Kemp
Horn Point Environmental Laboratory
University of Maryland
Cambridge, MD
Background
A simulation model of the shallow littoral zone (<3 m) of the Patuxent River was created to
investigate responses of submersed aquatic vegetation (SAV) to increased nutrient and sediment
inputs. The objective of the modeling effort was to explore mechanisms responsible for the recent
decline of SAV in Chesapeake Bay and to provide input to the management process about
conditions required for habitat restoration. The model is currently calibrated to represent habitat
dominated by Potamogeton perfoliatus, a well-studied grass formerly prevalent in the mesohaline
estuary and found in abundance along the flanks and in many tributaries of the mainstem Bay.
During the 1970s and 1980s, the distribution of Potamogeton and other submersed plants was
severely reduced by a combination of environmental perturbations (Batiuk et al. 1992, Madden
and Kemp 1996). Current management plans are directed at restoring SAV and increasing
abundance and distribution of Potamogeton to former levels, using among other things,
environmental benchmarks for several water quality characteristics which will hopefully lead to
habitat improvement (Batiuk et al. 1992). The model described here was developed with
assistance from the US Environmental Protection Agency and Chesapeake Bay Program Office,
Annapolis, Maryland, to aid managers with refining restoration guidelines and facilitate the
incorporation of littoral processes into other water quality models (Cerco and Cole 1993). This
model is also being used by ecosystem researchers to test hypotheses and predict ecosystem
responses to natural, anthropogenic, and management changes.
The model focuses on the relationship of light, nutrients and submersed vegetation. State
variables for above-ground (leaves and stems) and below-ground (roots and rhizomes) plant
biomass, phytoplankton, epiphytes, and paniculate sediment carbon are included. The model has
a time-step of 0.25 day and a simulation length of one to five years. The underwater light
environment as influenced by sediment, phytoplankton, and epiphytes is strongly emphasized in
model equation formulation, and four levels of PAR (photosynthetically active radiation) are
calculated for each time-step: PAR at the surface, at 0.25% of total depth, at the average canopy
depth, and at the SAV leaf surface. The model is used to examine how availability of nutrients
and light control growth and productivity of SAV plants, phytoplankton and epiphytes in
mesohaline Chesapeake Bay. Other environmental variables tracked by the model include
suspended sediment concentration, total water column turbidity, phytoplankton chlorophyll
concentration, and dissolved inorganic nitrogen concentration.
In previous model reports (Kemp et al. 1995, Madden and Kemp 1996), model results focused
on the relationship between dissolved inorganic nitrogen (DIN) and SAV biomass. These results
indicated that increases in DIN of 40% over baseline (1976) levels were sufficient to eliminate
58 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
SAV beds after three years. Increased DIN concentrations enhanced growth of algae, especially
epiphytes attached to SAV leaves, which reduced light at the SAV leaf surface below the
threshold required for net productivity in rooted aquatics. Reductions in water column nutrient
concentrations (in the model) did not result in rapid restoration of SAV biomass during the first or
even second year of simulated restored conditions. This lack of immediate response primarily was
attributed to the lag time required to regenerate substantial below-ground root biomass. That is,
the model showed that although water column light conditions may be favorable for SAV growth,
if the below-ground structure required to regenerate above-ground shoots was not sufficient,
appreciable growth rates would not be attained.
Model analysis also revealed a hysteresis in the relationship between ambient nutrient levels
and the amount of SAV biomass present. This result indicated that despite a strong correlation
between the amount of SAV biomass and ambient nutrient concentrations, levels of SAV
productivity and biomass depended on whether recent nutrient concentrations were trending
upward or downward. Previously higher nutrient concentrations stressed the SAV system and
significant damage to habitat structure may limit growth of SAV, even under declining or
acceptably "low" nutrient conditions. P. perfoliatus, with significant root/rhizome biomass,
tended to integrate the history of plant health over a longer period in the below-ground structure
than the above-ground structure. After periods of eutrophy resulting from high nutrient
concentrations, the amount of below-ground biomass available was reduced compared with plants
exposed to less eutrophic conditions. The model should quantify such interactions, predicting
system responses based on current and previous conditions. It should also predict expected SAV
recovery time under a variety of conditions.
New model experiments involved sediment and nitrogen feedbacks. Model experiments were
performed to elucidate the nature of these feedbacks between SAV beds and the water column.
Through such hypothesis testing, the efficacy of habitat criteria in restoring or protecting desirable
submersed aquatic plants can be explored. Understanding these relationships is important to
environmental resource managers because many environmental variables in question are used as
habitat criteria that relate water quality improvements to SAV restoration.
Model Development
The SAV model was comprised of systems of simultaneous differential (finite-difference)
equations, solved using a second-order Runge Kutta numerical integration scheme on a
microcomputer platform. A time-step (dt) of 6 hours permitted resolution of diel light,
photosynthesis and plant respiration patterns, and in the future, will permit the addition of tidal
and organism movement signals. Rate equations for physiological processes were derived using
relevant literature and empirical relationships from Patuxent River, Choptank River, or other
Chesapeake Bay tributaries. The baseline condition of the model was established at 1968, when
SAV was in significant abundance in the lower Patuxent River. Simulation experiments were
performed to project light and nutrient conditions in the bed under a variety of nutrient loading
rates. The model tracked biomass in units of organic carbon; phytoplankton, SAV, epiphyte, and
detritus stocks were calculated in mg C nr2. TheRedfield atom ratio of 106:16:1 for C:N:P, and
106:16 for C:Si for diatoms, was assumed for plant tissue. This ratio was used to index nutrient
uptake to carbon flow. N and P in the biota were also reported in mg nv2. Initial model
configuration simulated annual patterns for Potamogeton perfoliatus; the model is currently under
Simulation Model of SAV and Water Quality in Chesapeake Bay 59
-------
expansion to simulate additional species by varying rate coefficients, initial conditions and model
equations.
The model was optimized for performing specific experiments on the subsurface light regime
and SAV production. It describes a 1-3 m unstratified water column overlaying a benthic
community with which the water column community interacted through sedimentation, diffusive
fluxes, and nutrient translocation. Target variables included above-ground and below-ground
Potamogeton biomass, epiphytic growth on plant leaves, two groups of phytoplankton, and
organic material in sediments. Exchange between the grassbed in the model and the adjacent deep
river channel occurred via lateral advective flow through the model boundary. Flows were
regulated by a 10% mass water exchange flow and concentration gradients resulting from
processes in the model. These gradients were primarily driven by photosynthetic production,
metabolism, and the removal of epiphytes and phytoplankton by herbivorous grazers. Model
forcing functions included temperature, PAR, inorganic nitrogen loading, and suspended
particulates.
Growth in autotrophic compartments was controlled by maximum temperature-dependent
photosynthetic rates and was reduced by turbidity from phytoplankton, suspended particulates,
and epiphytes. Twilley et al. (1986) demonstrated that the attenuation of PAR at the leaf surface
followed an exponential decline with increasing epiphyte biomass. The relationship between
inorganic nutrient concentrations and photosynthesis followed Michaelis-Menton kinetics.
Zooplankton and herbivorous fish controlled grazing rate of phytoplankton and a community
of composite herbivorous grazers consumed epiphytes. Predator grazing on epiphytes and two
phytoplankton groups were described by separate Ivlev functions for each predator, which
included a preference function when there was a choice of prey. There is no direct grazing of
plant leaves in the model, as this loss for Potamogeton was generally deemed insignificant in
natural systems (Stevenson 1988).
Light in the model was calculated in five stages, starting from daily incident PAR, supplied by
forcing function. Number of daylight hours per day was determined using algorithms from Kirk
(1983). Incident PAR (variable name=PAR,) was partitioned into means for time intervals 0600-
lOOOh, 1000-1400h, and 1400-1800h. PAR, was reduced 10% at the water surface by a mean
scattering and reflection coefficient (PARj). Attenuation of light with depth was modeled with a
Beers-Lambert function (Kirk 1983) resulting in exponential decay of PAR intensity with depth.
The equation:
PAR^PAR^e'VZ) (1)
yielded light intensity in uEin m'V1 at any depth (Z), derived from PARj at the surface. Down-
welling attenuation (KD) was the sum of component attenuation from phytoplankton (Kc),
suspended particulates (Ks), and water itself (Kw) (Kirk, 1983). Equations of the form: KC= 0.054
* Chl°«7+ 0.0088 * Chi (Kroner and Nixon 1978); Ks= 0.0396*SPM+0.39 (Twilley etal. 1985);
and Kw= 0.03 (Kirk 1983) determined components of attenuation in units of m'1 for K, ug L'! for
Chi, and mg L'1 for suspended particulate material (SPM). Active chlorophyll concentrations
were based on chl:c ratios for diatoms of 1:50 and for other phytoplankton, 1:90 (Vant 1991,
McBrideetal. 1993).
PAR intensity at 30 cm (PAR,) represents integrated mean quantum irradiance for the 1 m
water column and for calculations of phytoplankton photosynthesis rates. Measurements by
Goldsborough and Kemp (1988), indicate mean leaf length for P. perfoliatus of 0.25 m, setting
60 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
the depth of the plant canopy (ZP) in the model at a constant 0.75 m for calculating irradiance at
canopy depth (PAR4) and photosynthesis rates for epiphytes. Light reaching submersed plant
leaves (PAR,) is used to calculate photosynthetic rates for submersed plants, was attenuated
through the epiphyte layer, based on epiphyte density on leaf surfaces as in Twilley etal. (1985):
PAR, = P AR, * e(C 32"c"-42'2 5AE» (2)
The quantity 2.5*AE = the density of epiphytes per unit plant surface area (mg cm*2), based on a
carbon to dry weight ratio of 1:2.5.
Difference equations for autotrophic compartments (Table 1; end of this article) incorporated
inputs of photosynthetic production; advective import; and in the case of rooted plants, growth of
regenerative shoots. Loss terms included mortality; advective export; respiration; exudation; and
in the case of rooted plants, translocation of carbon to roots and rhizomes. Biomass-specific
photosynthesis (Ps) for autotrophs was a function of maximum specific growth rate P,^, related
to temperature by an Arrhenius function (Odum 1983), adjusted by nutrient (N) and light (L)
limitation function:
Ps= P.. * f(N, L) (3)
where Ps and u are in units of mg C (g C d)"1; and f(N, L) ranges from 0-1, determined by nutrient
and light availability.
For submersed plant production, an exponential growth function with Q10 of 2.15 was used:
PAQ=DAO * 0.25 "(e'00633*T-Gmt)-1) (4)
where T = water temperature; G^,, the critical temperature coefficient, was set at 12° C,
calibrated to match the rate and timing of growth initiation of P. perfoliatus in Chesapeake Bay
tributaries (Stevenson et al. 1993). Such a temperature coefficient has been used in various plant
models (cf. Verhagen and Nienhuis 1983) to calibrate seasonal growth patterns. The self-
limitation term DAG is a density-dependent function determined by calibration calculated to apply
the effect of crowding and self-shading in the canopy:
D^HAG/cL^)2 (5)
AG is above-ground submersed plant biomass in g m2 and ct,^ a constant = 200 g C nr2. A
similar density-limited growth expression was used in the epiphyte production function:
DAE = (1^*0.17 *e(o°2*T) (6)
where the density dependent self-limiting effect for algae is proportional to carbon density on
plant leaves (mg C cm"2 of plant) as follows:
D^l-p^/d^)2 (7)
DESA = density of epiphytes on host leaves and d^^ is a constant = 20 mg C cm"2. Maximum
densities for both submersed plant and epiphyte were derived by inspection of values attained
Simulation Model of SAV and Water Quality in Chesapeake Bay 61
-------
under nutrient and light saturating conditions in field and laboratory studies (Staver 1984,
Twilley et al. 1985).
The phytoplankton maximum growth function is defined by a temperature-dependent
formulation similar to that employed by Kremer and Nixon (1978) in a model of Narragansett
Bay phytoplankton:
GP = 0.59e(0-0633*(T-Gcrit)) (8)
where Gp is maximum photosynthetic rate in mg C mg C1 d"1, and T is Celsius temperature;
growth is modified by addition of a specific temperature coefficient, where G^ =: 1 ° C for
diatoms and 8° C for other (summer) phytoplankton. The higher temperature coefficient for
summer plankton delays growth and formation of peak biomass until early June; diatoms, with
lower temperature optima, peak in March. Self-shading in phytoplankton was modeled
implicitly via the plankton chlorophyll contribution (Kc) to total water column attenuation
(KD) (Huisman and Weissing 1994).
For all autotrophs, maximum temperature-based photosynthetic rates were downward-
adjusted by limitation terms describing instantaneous nutrient and light-based carbon uptake.
Carbon fixation rates for epiphytes and phytoplankton were modeled as functions of water
column nutrient concentrations following Michaelis-Menton kinetics. Photosynthesis related
to irradiance by rectangular hyperbola (e.g. Dennison and Alberte 1986) using different
photosynthetic parameters (P^ and alpha) for each producer group to generate characteristic
P vs I relationships. The minimum of all potential rates (nutrient and light) acted as the
limiting effect on growth.
F(N,L) = MIN(N , P , Si , PAR; ) (9)
N+K. P+K,, Si+K,; PARs+Ik
This minimum value formulation selected a single limiting rate (either nutrient or light) by
comparison of availability versus requirement for all four resources, (Kremer and Nixon
1978), and assumed no multiplicative limitation through nutrient interaction (cf. DiToro et al.
1971).
For nutrient uptake kinetics by submersed plants, the standard Monod (1942) formulation
was modified to account for dual sources of nutrients: rooted plants, with access to pools in
sediment pore waters as well as water column, incorporate both root and leaf uptake
pathways (Kemp et al. 1995) as follows:
U= Nw + k'K (10)
K + (Nw + k'Ns)
where U (relative nutrient uptake rate) represents total uptake potential from both the
sediment and water column pools. Nw and Ns represent water column and pore-water nutrient
concentrations, respectively; K is the half-saturation constant for uptake of Nw; and (k'/K) is
the reciprocal of the half-saturation for Ns (Kemp et al. 1995). It can be seen that when the
water column concentration (Nw) is zero, the equation reduces to the simple Michaelis-
Menton equation for uptake through roots. When the sediment pool concentration is equal to
zero, the equation becomes which describes uptake kinetics through leaves alone.
Intermediate conditions between these extremes produced an aggregate uptake rate based on
62 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
the relative proportions of half-saturation concentrations represented by pool concentrations.
This formulation, therefore, simultaneously calculated uptake via both pathways, based on the
relative proportion of demand satisfied by each pool, and the model inputs that combined
proportion (0-1) to individual equations for N, P, and Si. The minimum function (eq 9)
selected the term limiting to photosynthesis, as for the phytoplankton nutrient kinetics
described above. Kinetic coefficients were chosen from literature values from controlled
experiments by Izumi and Hattori (1982) and Thursby and Harlin (1982). Model runs
produced patterns of DIN uptake by leaves and roots similar to experimental observations of
Thursby and Harlin (1982) for Potamogeton sp.
Results and Discussion
Calibration
Above- and below-ground SAV biomass patterns generated by model baseline runs
calibrated well with available data from Chesapeake Bay and were representative of SAV
dynamics in mesohaline estuaries such as the Patuxent and Choptank Rivers (Kemp et al.
1984, Twilley et al. 1985, Lubbers et al. 1990). The model exhibited stable behavior under a
variety of initial and forcing conditions, producing simulations closely replicating plant growth
patterns measured in laboratory studies using Chesapeake Bay plants (Staver 1984,
Neundorfer and Kemp 1993), as well as in field studies from sites around the world
(Cambridge et al. 1986, Izumi and Hattori 1982, Huisman and Weissing, 1994, Burkholder et
al. 1994). The model was robust over long simulation times of up to five years, although it
was strongly responsive to small changes in key variables such as DIN concentrations and
epiphyte growth rate. Monte Carlo analysis is being conducted to assess model sensitivity to
each variable in a systematic manner. Because the model operates in a spatially averaged
mode, and does not have a spatial component, the expansion or contraction of SAV
distribution zones cannot be explicitly produced by the model in its current form. However,
spatial information can be inferred from predictions of success or failure of submersed
vegetation in the 1 m2 spatial unit.
Experiments
The model was used to assess ecosystem relationships that would be relevant to
management activities, including the development of habitat criteria for increasing baygrass
distributions and densities. The response of model state variables to elevated nutrient
concentrations simulated the effects of eutrophication, including higher phytoplankton and
epiphyte growth and biomass, and the decline in SAV productivity as observed in Chesapeake
Bay. When comparing nutrient concentrations required to elicit reductions in growth
response or elimination of SAV in the model with those in the real system, model output was
found to parallel observations in the Patuxent and Choptank Rivers, as well in various
mesocosm and pond experiments (Kemp et al. 1984, Lubbers et al. 1990, Neundorfer and
Kemp 1995, Twilley et al. 1986). Elevated DIN concentrations of 40% above baseline led to
extremely low P. perfoliatus density and eventual loss of SAV.
Simulation experiments were conducted to examine critical feedback relationships in the
SAV habitat. Two are reported here, including the effect of SAV on water column nutrient
levels and the effect of SAV on water column turbidity. The nutrient experiment was
Simulation Model of SAV and Water Quality in Chesapeake Bay 63
-------
Influence of SAV on Nitrogen Dynamics
Inamato
DINIOBlng
% reduction in
UN by SAV
40%
200*
29%
performed by subjecting the
baseline SAV model system
to increasing nutrient inputs
and comparing the resulting
SAV levels and dissolved
nutrient concentrations with
those in a similarly treated,
but unvegetated site.
Except for the absence of
above-ground SAV leaves,
the modeled "bare" site was
identical in all respects to the
vegetated site, including
phytoplankton
concentrations., sediment
nutrient levels, potential for
epiphyte colonization, etc.
Initially, SAV plants were
set at densities of 0 g nr2 leaf
biomass and 45 g nr2
root/rhizome biomass in
both "vegetated" and
"unvegetated" sites, but
unvegetated sites were not
allowed to initiate above
ground growth. Both
habitats were exposed to nutrient loads of 0%, 20%, 40%, 100%, and 200% above baseline.
The effect of SAV on water column nutrient levels was evaluated by determining mean
nutrient concentrations predicted by the model in vegetated sites during the growing season
(May-November) and comparing them with nutrient levels in unvegetated sites during the
same period. The experimental simulation was carried out for a period of two years, ensuring
model stability and reducing treatment sensitivity to initial conditions (Figure 1).
Results showed a number of interesting behaviors related to the way in which nutrients
were processed by the SAV community. Vegetation significantly affected the nutrient levels
in the beds. When present, SAV significantly reduced DIN concentrations during summer,
even under the highest nutrient loading treatments. During Year 1, at all nutrient levels, DIN
concentrations were about 50% lower in the waters overlying SAV beds, compared with
unvegetated sites, indicating that after a certain lag time, most of the available new nutrient
inputs were being taken up by plants. During Year 1, SAV began to decline during late
summer in higher nutrient treatments and in Year 2, SAV biomass continued a progressive
downward trend until the bed was essentially eliminated in the highest nutrient treatments.
With reductions of SAV in Year 2, DIN concentrations progressively increased 'with
additional DIN loading and the difference in DIN concentration over bare and vegetated sites
diminished. In the +200% treatment, the difference between sites was insignificant.
Because the SAV root/rhizome component was impacted on a longer time scale than the
above-ground component, loss of below-ground material in some cases was not evident until
Figure 1
64 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
the second year. In Year 1, high initial SAV biomass growth potential skewed the comparison
of vegetated and unvegetated treatments because calculations of average nutrient
concentrations included the spring growth phase, during which epiphytes did not have a
chance to affect biomass production. In Year 1, due to healthy below-ground biomass, SAV
grew strongly in spring even under high nutrient treatments and was able to take up significant
water column nitrogen. In Year 2, the long-term effects of DIN loading were more apparent,
and spring growth of SAV was affected by the previous year's poor growth and low below
ground reserves. Running the model through Year 2 was thus essential to reduce the
influence of initial conditions (i.e. the high potential for SAV growth) from results.
In all treatments with SAV present, there was a significant regeneration of DIN in fall and
winter due to senescence. During fall and winter, DIN concentrations in water overlying
vegetated sites exceeded that in unvegetated sites, a pattern which was most marked in the
highest loading treatments. This suggested that below-ground decomposition and recycling
contribute to nutrient pools, while at the same time, SAV was not present to intercept the
sediment nutrients as they were released to the water column.
A second experiment involved the effect of SAV on suspended sediment concentrations in
the bed. Ward etal. (1985) showed that SAV leaves could significantly reduce water current
speeds within SAV beds, permitting sediments to fall out of suspension. Roots and rhizomes
reduced the tendency for water motion to resuspend bottom sediments. Using empirical data
on concentrations of suspended paniculate material (SPM) in Choptank River Potamogeton
beds of different leaf densities, a relationship was developed which represented the reduction
of SPM in waters overlying the grassbed:
SPM = SPMa * 1-(0. 1 *e
-------
Interestingly, applying proposed habitat criteria for nitrogen in the model automatically
generated levels for other variables that approximated habitat criteria levels: the model
predicted that chlorophyll a and KD would reduce to 15 mg nr3 and 1.5m'1 respectively if
water column nitrogen concentration were reduced to to 10 uM, as called for in Tier 1
restoration criteria. Based on model predictions, Tier I habitat criteria should be appropriate
for achieving goals established for SAV in the mesohaline Bay and may enable regrowth of P.
perfoliatus in distributions observed in the lower Patuxent and Choptank Rivers during the
1960-1970s. As demonstrated, however, complexities in several feedback relationships, such
as loss of sediment trapping ability of SAV, may make synergistic impacts more intractable
than anticipated. Ongoing model studies will contribute to further understanding of these
relationships and assist in refining critical thresholds for habitat quality.
66 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
_c
en
V
1
CO
0>
**
eS
tn
t_i
<2
en
_O
"J3
T3
O
0
.Is
-a
c
cS
en
O
°+3
CO
§3
O
g s
0) >->
c c
j
2 S
f2"§
r* en
4-»
o
*
^s
o
00
1
Q
o,
o
<
Q
a,
1
Q
of
c/
F^
1
Q
a.
v<
tS
1
Q
a,
5"
1
S
£
O
1
Q
a.
V3
t3
*
^2
o
00
1
o
4"
1
2
i
o
CL
S
I
2
X
w
i
g
s2
i
o
0.
S
t
o.
<]
^*-4
+
C/3
fe
.+
H^^
T3
i
<->
O
OH
II
S
o
OH
^-*
T3
*
/-»v
^
JD
00
i
5?
*|
i
U
«?
0<
!3
O
U
<
^
+
9
J?
en .5°
en OH
^3 N|*^'
£ +
.2 ^
PQ -S
U i
£>^
t^
W II
If
T3
3
2
O
o
_o
<
^^^
en
en
C
£
o
: Leaf/Stem Bi
CQ
S
4-J
"O
*
1
|
1
CA
H
1
o
H
i
D
««
(Ji
1
3
H
+
§
<
W3
OH
v»^
+
-*-»
"O
i
*->
O
<
ll
^-/
O
<:
^
*n
^j
c
o
O
_o
*-»
a.
fll
.>
*-*
o
V
T^
W
CQ
II
'«5
c"
_O
*-»
O
3
a
o
UN
a.
_o
4-1
CD
.C
^_»
>>
en
O
41
O
JH
cx
II
en
OH
ooplankton, menhaden, or epiphytic grazers);
utput; Sd=sedimentation; and Sb= leaf substrate loss.
N 0
^ >
en .£,
° ts
M ^
.£ T3
N «J
S II
oo o
i<
O <=>
^ S
*.> 3
O t!
, 0 g
\i^
o -^
i t^~i
^
§ S
"§£»
Sf
J3 O,
cx pq
feT3
o o
en en
en en
CS ctS
|l
S CQ
NO
O
CQ
S
;* Q o 5
C« OH OH
-------
s s s
u u o
bO (50 GO
ESS
o o
o o
o o
o o
T es
4
O T3
CCo>
3 3 ~
2 2 -S
t_ tm CC
tso 60-r
O ^ ^l
a^ s
^2 ^2 o
^ -a
o o
tC
o
en
en
.^
pa
06
CO
a
os
U
o
K>
"«
T3
O
s
£
VI
O
o o
< CQ
oe
-------
References
Batiuk, R. A, R. J. Orth, K. A. Moore, W. C. Dennison, J. C. Stevenson, L. W. Staver, V. Carter,
N. B. Rybicki, R. E. Hickman, S. Kollar, S. Bieber, and P. Heasly. 1992. Submerged aquatic
vegetation habitat requirements and restoration targets: A technical synthesis. CBP/TRS
83/92. Annapolis, MD 186 pp.
Burkholder, J. M., H.B Glasgow and J. E. Cooke. 1994. Comparative effects of water column
nitrate enrichment on eelgrass Zostera Marina, shoalgrass Halodule wrightii, and
widgeongrass Ruppia maritima. Marine Ecology Progress Series 105:121-138.
Cambridge, J. L., A. W. Chiffings, C. Brittan, L. Moore, and A. J. McComb. 1986. The loss
of seagrass in Cockburn Sound, Western Australia II: Possible causes of seagrass decline.
Aquatic Botany 24:269-285.
Cerco, C., and Cole, T. 1993. Three-dimensional eutrophication model of Chesapeake Bay.
Journal of Environmental Engineering 119(6): 1006-1025
Dennison, W. C. and R. S. Alberte. 1986. Photoadaptation and growth of Zostera marina L.
(eelgrass) transplants along a depth gradient. Journal of Experimental Marine Biology and
Ecology 98:265-282.
DiToro, D. M., D. J. O'Connor, R. V. Thomann. 1971. A dynamic model of phytoplankton
populations in the Sacramento-San Joaquin Delta. Advances in Chemistry Series 106. p
131-180.
Goldsborough, W. J., and W. M. Kemp. 1988. Light responses of a submersed macrophyte:
Implications for survival in turbid waters. Ecology 69:1775-1786.
Huisman, J. and F. J. Weissing. 1994. Light-limited growth and competition for light in well-
mixed aquatic environments: An elementary model. Ecology 75(2):507-520.
Izumi, H. and A. Hattori. 1982. Growth and organic production of eelgrass (Zostera marina)
in temperate waters of the Pacific coast of Japan. III. The kinetics of nitrogen uptake.
Aquatic Botany 12:245-256.
Kemp, W.M., W.R. Boynton, J.C. Stevenson, and L.G. Ward. 1984. Influence of submersed
vascular plants on ecological processes in up Chesapeake Bay. In: V.S. Kennedy (ed.)
The Estuary as a Filter. Academic Press. New York. Pp. 367-393.
Kemp, W. M., W. R. Boynton, and A. J. Hermann. 1995. Simulation models of an estuarine
macrophyte ecosystem. In: B. C. Patten (ed.) Complex Ecology. Prentice Hall,
Englewood Cliffs, NJ. p 262-277.
Simulation Model of SAV and Water Quality in Chesapeake Bay 69
-------
Kirk, J. T. 0. 1983. Light and photosynthesis in aquatic ecosystems. Cambridge Univ. Press.
Cambridge, UK. 401 pp.
Kremer, J. and S. W. Nixon. 1978. A Coastal Marine Ecosystem: Simulation and. Analysis.
Springer-Verlag, New York. 217 pp.
Lubbers, L., W. R. Boynton, and W. M. Kemp. 1990. Variations in structure of estuarine fish
communities in relation to abundance of submersed vascular plants. Marine Ecology
Progress Series 65: 1-14.
Madden, C. J. and W. M. Kemp. 1996. Ecosystem model of an estuarine plant community:
Calibration and simulation of eutrophications responses. Estuaries 19(2B):457-474.
McBride, G. B, W. N. Vant, J. E. Cloern, J. B. Liley. 1993. Development of a model of
phytoplankton blooms in Manukau Harbor. NIWA Ecosystems Publication No. 3.
Hamilton, New Zealand.
Monod, J. 1942. Recherches sur la croissance des cultures bacteriennes. Paris: Herman et Cie.
Neundorfer, J. V. and W. M. Kemp 1993. Nitrogen versus phosphrus enrichment of
brackish waters: responses of the submersed plant Potamogeton perfoliatus and its
associated algal community. Marine Ecology Progress Series 94:71-82.
Odum, H. T. 1983. Systems Ecology: An Introduction. Wiley and Sons. New York. 644 pp.
Staver, K. 1984. Responses of epiphytic algae to nitrogen and phosphorus enrichment and
effects on productivity of the host plant, Potamogeton perfoliatus L. in estuarine waters.
MS thesis, Univ. MD, College Park, MD.
Stevenson, J. C. 1988. Comparative ecology of submersed grass beds in freshwater, estuarine,
and marine environments. Limnology and Oceanography 33:867-893.
Stevenson, J. C., L. W. Staver, K. W. Staver. 1993. Water quality associated with survival of
submersed aquatic vegetation along an estuarine gradient. Estuaries. 16(2): 346-361.
Thursby, J. E. and M. M. Harlin 1982. Leaf-root interaction in the uptake of ammonium by
Zostera marina. Marine Biology 72:109-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. Marine Ecology Progress
Series 23:179-191.
Twilley, R. R., G. Ejdung, P. Romare, and W. M. Kemp. 1986. A comparative study of
decomposition, oxygen consumption and nutrient release for selected aquatic plants
occurring in an estuarine environment. Oikos 47:190-198.
70 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Vant, W. N. 1991. Underwater light in the Northern Manukau Harbor, New Zealand.
Estuarine, Coastal and Shelf Science. 33:291-307.
Verhagen, J. H. G, and P. H. Nienhuis. 1983. A simulation model of production, seasonal
changes in biomass distribution of eelgrass (Zostera marina) in Lake Grevelingen. Marine
Ecology Progress Series 10:187-195.
Ward, L. G., W. M. Kemp, and W. R. Boynton. 1985. The influence of waves and seagrass
communities on suspended particulates in an estuarine environment. Mar. Geol. 59:85-
103.
Simulation Model of SAV and Water Quality in Chesapeake Bay 71
-------
Simulation Model of Biogeochemical Processes in Ma rsh
Mesocosms
Christopher J. Madden, Jennifer Zelenke, and J. Court Stevenson
Horn Point Environmental Laboratory
University of Maryland
Cambridge, MD
Background
The high degree of complexity and variability in coastal marshes makes the use of models an
important means of synthesizing data and understanding the biogeochemistry of marsh systems.
In mesocosms established at the Multiscale Experimental Exosystem Research Center (MEERC)
facility at Horn Point Environmental Laboratory (HPEL), University of Maryland, a series of 12
artificial marshes containing sediments and several types of vegetation will be analyzed both
experimentally and with modeling techniques. Research on these mesocosm communities will
determine the effects of time and community complexity on system functioning. Experiments will
involve the manipulation of groundwater nitrogen inputs, plant community diversity, hydroperiod
and water depth to determine their effects on nutrient processing, marsh productivity and plant
succession. Results will be extrapolated using numerical models to help understand natural marsh
system functioning. The experiments will be conducted over a period of several years beginning
in 1996.
The primary purpose of the modeling effort will be to synthesize, interpret and extend the
results of the mesocosm experiments. The plan is to develop the model concurrently or in advance
of experiments, giving researchers tools for testing hypotheses, predicting outcomes of
experiments, and even aiding in designing treatments and methodologies. Because model and
experimentation will be interactive, feedback from each should be of assistance with refining the
other.
Hydrological, biological and chemical processes will be examined using the model and
relationships of these processes will be interpreted across spatial scales. Results of mesocosm
experiments will be continuously assimilated into the model and results of modeling experiments
will be extrapolated to the Chesapeake Bay ecosystem. The model is in the early phase of
development and will be described here in general terms, along with a presentation of preliminary
results.
Model Development
A working six-compartment model of the MEERC marsh mesocosms has been developed and
initial diagnostics are being run. Initial efforts focus on modeling a single plant species, creating
hydrological and geochemical sectors around the biological structure of the well-known marsh
cordgrass species, Spartina alterniflora. Subsequently, sectors describing additional plant species
and, eventually, consumers will be incorporated. This scheme gives the advantage of simplifying
early stages of modeling as the model framework (forcing functions, physical dimensions) and
72 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
sectors (hydrology, sediment chemistry) that will be common to all biological sectors are
formalized.
Tidal and groundwater forcing functions in the mesocosms will have significant impact on
both the actual mesocosm and the model processes, and special emphasis has been placed on
accurately capturing hydrological parameters of the mesocosm. The marsh mesocosms are raised
at one end to create a 5 ° slope, allowing partial inundation by an artificial tide entering at the
lower end. We determined, in consultation with the experimental team, that a two-cell model
(high marsh and low marsh) should be established to represent the non-flooded and flooded
portions of the mesocosm. The model will use a floating boundary between cells, its position
dependent on the physical extent of flooding in the mesocosms, which is determined by
experimental treatment. Because flooding is a variable and dominant factor in sediment
chemistry, plant productivity and animal response, it is of paramount importance to capture the
spatial extent of the water coverage on each tidal cycle in the model.
The biological submodel consists ofSpartina alterniflora above-ground and below-ground
carbon, as well as associated leaf and root litter. Biological sectors will be expanded to include
Eleocharis sp., Spartinapatens, Hibiscus sp., and Scirpus olneyi plant populations. The
difference equations for each of the plant species will be similar to the equations for S.
alterniflora:
(1) SHOOT(t) =
SHOOT(t - dt) + (Sht_Prod + Upjrans - Sht_Graze -
Sht_Mort - Sht_Rsp - H - Dnjrans) * dt
(2) ROOT(t) =
ROOT(t - dt) + (Dnjrans - Rt_Mort - Rt_Graze - Rt_Resp - Upjrans) * dt
where above-ground biomass SHOOT (g C nr2) is the sum of shoot production by photosynthesis
(Sht_Prod), upward translocation from roots (Upjrans), minus losses from direct grazing
(Sht_Graze), mortality (Sht_Mort), respiration (Sht_Rsp), harvesting (H) for samples, and
downward translocation to roots (Dnjrans). Root biomass, ROOT (g C nr2), is the sum of
downward transport of photosynthesis (Dnjrans), minus mortality (Rt_Mort), direct loss to
grazing (Rt_Graze), root respiration (Rt_Resp) and translocation of stored carbon to shoots
(Upjrans).
The geochemical submodel includes nitrogen, phosphorus, sulfur and oxygen dynamics.
Hydrologic input via groundwater is forced using experimental design specifications. Sediment
porosity, bulk density and diffusivity of mesocosm soils are being analyzed and will be
incorporated into the model as data become available. Presently, groundwater functions are based
on empirical measures of water transit time and groundwater export rates in selected MEERC
mesocosms. Model parameters include a time unit of 1 day, a time-step (dt) of 0.25, a simulation
length of 1 yr. The spatial unit of the model is m"2, and forcing functions include air temperature,
PAR (photosynthetically active radiation), groundwater input rates, and groundwater nitrate
concentrations. Boundary conditions included opaque, impermeable walls containing two
horizontal cells for high and low marsh. Two vertical layers consisted of sediments and overlying
water. Mesocosms were open and exposed to natural sunlight, air and precipitation.
Simulation Model of Biogeochemical Processes in Marsh Mesocosms 73
-------
Results
Calibration
The model has initially been developed for a monospecific stand ofSparfina alterniflora,
calibrated to produce growth and biomass patterns characteristic of temperate US East coast
marshes. Two datasets were used to calibrate the Spartina model (Figure 1): a South Carolina
salt marsh (Morris and Haskin 1990) and an Eastern Shore, Maryland salt marsh in Monie Bay
(Stribling and Cornwell 1995). The model generates reasonable values for above-ground
biomass, showing growth initiation in early spring, reaching a peak in late summer at about 175 g
C nr2 before dieoff occurs. Both model and field data indicate that a significant amount of
standing stock (approx. 1 10 g C nr2) remains in place at the end of the growing season and
persists into winter. This material is assumed to be mostly highly refractory standing dead
material. To reduce standing material to winter levels in the field (about 40 g C nr2), a special
function will be required in the model.
As yet, the forcing (e.g. wind, rains, floods) which actually removes this material in the
field, and its schedule, is not determined. This does however, illuminate an interesting point
observed when comparing model results to the mesocosm biomass patterns: at the end of the first
growing season in the actual mesocosm marsh community, a large percentage of standing dead
biomass persisted throughout the
winter into spring. In other words,
the mesocosm behaved more like the
model than like the natural marsh,
which generally has low standing
biomass at the end of winter.
Apparently, a physical forcing is
absent in the mesocosms which is
responsible for removing standing
dead material in the field each year.
Model-mesocosm analysis may
provide information concerning this
E
O
I
01
i
I
200.00
100.00
\
\
MODEL OUTPUT
./
\
S. Carolina
Morris & Haskin 1990
Monie Bay, MD
StribKng 1995
o oo
F M A M J J
O
observation.
Figure 1 . Calibration of plant biomass to literature values.
Biomass sampling
The marsh model has been used to aid design of mesocosm experiments a priori to predict
effects of experimental treatments on marsh productivity. Initially, it was used to estimate the
expected range of standing biomass in the area dominated by S. alterniflora in the mesocosms and
to quantitatively assess appropriate scales and schedules for biomass sampling. The procedure by
which productivity will be determined in the mesocosms is destructive biomass harvesting. This
procedure introduces some error into the experiments. One objective of the modeling exercise
was to determine the maximum rate of destructive sampling that S. alterniflora plantings could
support, while sustaining minimal impact on plant growth. The model was run with different
levels of biomass removal on a variety of schedules to observe the projected impact on standing
stock. Experimenters intended to use the information to adjust sampling rates to optimize for
plant biomass, replicating the nonharvested condition. Results indicated that monthly sampling of
74
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
approximately 5% of the above-ground biomass would not cause large deviations in standing
stock from the unimpacted condition.
An important point needs to be raised concerning this experiment. Successful duplication
of unimpacted patterns of biomass density does not necessarily mean that marsh processes are
equivalent in both harvested and unharvested conditions. As revealed in the model, following
each harvesting, a brief increase in plant productivity occurred, probably in response to increased
light, reduced nutrient depletion, and lowered density-related inhibition. This resulted in a quick
"filling in" of marsh standing stock, but the acceleration of processes (in the model) associated
with the additional production caused changes in nutrient utilization, oxygenation of the root
zone, and decomposition rates. Thus, although standing stocks in the mesocosm may be
minimally impacted by a particular harvest schedule, other processes may in fact be strongly
impacted. The model is currently being configured to analyze harvest impacts on all relevant
geochemical processes and a suite of responses will be generated.
Nitrogen Transformations
Model experiments were performed to determine how processing of nitrogen may be
affected by changes in the rate of gross primary production of S. alterniflora (Figure 2).
Denitrification represents a permanent removal of nitrogen from the system via conversion of
nitrate to the diatomic gaseous form (NO3 to N2(g)) and is an important term in the nutrient
budgets of the mesocosms. Denitrification rates are dependent on sediment oxygen demand,
temperature and nitrate concentrations. A series of experimental treatments using a range of
Spartina biomass levels and productivity values was performed using the model. Simulated
denitrification rates increased in the modeled system with increasing plant biomass and
productivity. Because loss rates of nitrogen via denitrification often exceeded the supply of
nitrate from groundwater inputs, other sources of nitrate (in addition to groundwater inputs)
probably accounted for modeled denitrification rates. Much of the denitrification loss was
apparently supported by decaying plant tissue which produced an increase in remineralized
ammonium and was subsequently oxidized to N03 via nitrification (NH4 to NO3). Increases in
primary productivity strengthened the coupling of the denitrification and nitrification processes by
increasing ammonium supply to the oxidized upper sediments. The importance of this mechanism
is currently quantitatively being investigated using the model. Both processes by which
denitrification might increase (increased ammonium supply and reduced redox potential in the
lower root zone) are dependent upon rates of vascular plant production. The interaction of
biology and geochemistry has far-reaching management importance in marsh environments and
the proper interpretation of the mesocosm experiments using modeling techniques will hopefully
enable the extrapolation of experimental results to the natural system.
Discussion
The mesocosm model demonstrated reasonable behavior for U.S. temperate East Coast
salt marshes and calibration of the model for the mesocosm system is currently underway. The
model is useful as an interactive tool for designing experiments. Its utility with assisting
experimental design was demonstrated; levels of manipulations and measurement techniques and
expected impacts on the mesocosm system were calculated. In other experiments, results
suggested that permanent nitrogen loss through denitrification was dependent on productivity of
the marsh plants, but the importance of nitrogen removal versus nitrogen storage in plant biomass
Simulation Model of Biogeochemical Processes in Marsh Mesocosms 75
-------
has yet to be determined. The model should prove useful in examining these questions.
Development will continue with linkage of the hydrologic submodel to biogeochemical
components. As empirical data
become available, incorporation into a
the model will enable better e.oo
convergence of the two (experimental s
and model) systems and scaling of "j
mesocosm processes to the natural
system will follow.
g 3.00
9
0.00
FMAMJJ ASOND
I
|
^j
5.00
2.50-]
a
0.00
4.0
0.00 2.50
Nitrification mg/m2/day
5.00
Figure 2. Modeled responses of denitrifcation to changes in S. altemitlom productivity. Gross
production rate multiplier is indicated and was used to mimic relative ecosystem health, (a) ch
denitrification in response to primary productivity (b) shifting of nilrification-denitrification in resp
to primary production.
76
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
References
Kaplan, W, I. Valiela, and JM Teal. 1979. Denitrification in a saltmarsh ecosystem. L&O 24(4):
726-734.
Morris, JT, and B Haskin. 1990. A 5-year record of aerial primary production and stand
characteristics of Spartina alterniflora. Ecology 71(6): 2209-2217.
Thompson, SP, HW Paerl, and MC Go. 1995. Seasonal patterns of nitrification and denitrification
in a natural and a restored salt marsh. Estuaries 18(2): 399-408.
White, DS and BL Howes. 1994. Translocation, remineralization, and turnover of nitrogen in the
roots of Spartina alterniflora (Gramineae). American Journal of Botany: 81(10):
1225-1234.
Simulation Model of Biogeochemical Processes in Marsh Mesocosms 77
-------
Zooplankton Process Model for Mesohaline Chesapeake Bay
Christopher J. Madden and W. Michael Kemp
Horn Point Environmental Laboratory
University of Maryland
Cambridge, MD
Background
Analysis of the role of zooplankton processes in trophic and nutrient cycling is an important
component of the ecosystem process modeling strategy for Chesapeake Bay. Acartia tonsa and
Eurytemora affinis are the two most abundant zooplankton species in the mesohaline reaches of
the Bay mainstem, as well as in the lower parts of many Maryland tributaries, particularly the
mesohaline Patuxent River. The two species differ in their responses to seasonal changes,
resulting in winter-spring dominance by Eurytemora and summer dominance by Acartia.
Zooplankton represent a significant and variable source of predation to lower trophic levels,
including herbivory of phytoplankton, consumption of microzooplankton and protozooplankton,
and cannibalism of juvenile mesozooplankton. Resolution of grazer population cycles in
Chesapeake Bay is important to understanding overall system functioning and is imperative for
properly modeling the ecology of the system. Therefore, a stand-alone model of zooplankton
interactions was developed which can be used to test the role of zooplankton in carbon and
nutrient cycling in Chesapeake Bay.
This simple nutrient-phytoplankton-zooplankton (N-P-Z) model was conceptualized with
zooplankton compartments aggregated into three state variables: Eurytemora sp., Acartia sp., and
microzooplankton (Figure 1). These compartments dynamically interact with lower trophic
levels. Predation from higher trophic levels (notably Anchovy) occurs in this model via seasonally
variable forcing functions. The model will be used to examine phytoplankton-zooplankton
interactions, simulate trophic carbon and nutrient transfers, interpret empirical observations of
abundance and distribution, quantify effects of predation on zooplankton, and predict the response
of zooplankton to water quality changes via bottom-up effects. Simulations will examine the
functional role of zooplankton in nutrient regeneration, carbon transfer, anoxia and phytoplankton
dynamics to assess the importance of zooplankton in the environment from a management
perspective. This model will be linked to a series of similar models developed to examine aspects
of zooplankton ecology and function in the mainstem Chesapeake Bay and reaches of the other
tributaries.
A second model (See Madden etal., A Stage-Structured Model..., this report) was created
to examine factors influencing zooplankton community development. Environmental factors affect
generation time, abundance-biomass relationships, organism condition, and size-dependent
predator-prey relationships. Such a model requires stage-specific modeling of the zooplankton
community and will be used to evaluate phytoplankton size relative to zooplankton size as a
means of determining feeding behavior and top-down effects on primary production. This
approach also yields a size-frequency spectrum of the zooplankton community, which is of
relevance in their role as food items for planktivorous fish. (See Madden et al., A Stage-
78 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Aggregated Chesapeake Bay Zooplankton Dynamics (PNZ) Model:
Growth and Abundance Controls
Fish* (1-EXfr 0.01 -(C6-8))
EXpO.ISTTemplC-IOh
% recruit= EN/ (p- EXpO-15-(TemplC-1Cy
Structured Model..., this
report.)
Together, the two
modeling approaches
synthesize Chesapeake Bay
monitoring and modeling
activities to assimilate and
interpret new data. The
models will help direct
monitoring activities which
will improve data collection
efficiency and productivity.
Model analysis can help
determine sampling
locations, as well as temporal
and spatial scales for data
collection. Ecosystem
processes models can take
information gained from
small-scale sampling and
extrapolate to large-scale
processes.
Model Development
A simple numerical model
of nutrient-phytoplankton-
zooplankton interactions is
being developed using the
STELLA II simulation
environment on a RISC-based microcomputer. The strategy for developing process models of
Chesapeake Bay and tributaries is to calibrate separate models for representative segments of the
habitats and then combine the models as confidence is gained in the calibration of each
component. Model development is currently underway for tidal fresh, transitional, and
mesohaline habitats in the Patuxent River. These models will be linked via hydrologic exchanges.
The transitional and mesohaline segments will be represented by identical model structure, with
variations in parameters, time constants, and hydrology, and the input functions specific to each
area (e.g. nutrient loading). The tidal fresh version of the Patuxent River model includes
additional freshwater species (e.g. Bosmina sp)
The mesohaline Patuxent River model represents a spatially averaged 1 m2 area of water
column and benthic sediment in the mesohaline open water region of the Patuxent River. The
zooplankton process model described here divides zooplankton into two species groups (Acartia
and Eurytemord), and a "microzooplankton" state variable which aggregates nauplii and
Feeding= Ax =*CAR 1 - EXP" °-00017" (° 4'fo
-------
copepodid stages from both species, plus other small zooplankton and protozoa. The water
column is treated as a single homogeneous layer in the littoral model. In the deep channel
version, where there is two-layered structure, each state variable has upper and lower water-
column components. Two versions of the mesohaline model were developed: one for the shallow,
single-layered water column in the littoral zone and another for the deeper, channel region where
stratification in the river results in formation of a two-layer water column structure. In the
shallow model, the system is treated as vertically averaged over the 1-3 m water column. In the
two-layer version, the upper layer is fixed at 6 m and {he bottom at 10m. For other reaches of
the Patuxent River (transitional, tidal fresh) and other tributaries, the vertical dimensions of the
model were adjusted to the mean bottom depth.
In the shallow littoral model, the water column is treated as a homogeneous parcel which
exchanges nutrients and carbon with the benthos and sediments to a depth of 15 cm. In the two-
layer model, the upper and lower water column exchange phytoplankton and particles by sinking;
zooplankton via bidirectional migration; and nutrients via mixing across the pycnocline. As in the
single layer model, the lower water column exchanges nutrients and carbon with the benthos and
sediments to a depth of 15 cm.
Carbon transfer is tracked in units of mg C, with conversion of C to N using Redfield
stoichiometry of 7.2:1. A time-step of 0.25 day, a base simulation length of 1 year, and a second-
order Runge-Kutta integration scheme were used in the model. Simulation runs were conducted
for single and multi-year periods.
The model consists of 12 state variables, each represented by a difference equation calculating
total carbon biomass from summed inputs and outputs over each time-step. State variables
include diatoms, flagellates and chlorophytes, Acartia tonsa and Eurytemora qffinis,
microzooplankton, benthic filter feeders, dissolved inorganic nitrogen (DIN), and oxygen.
Modeling of the dominant zooplankton functional groups allows a distinction based on differences
in food supply and temperature response resulting in a cold weather dominance of Eurytemora
and a warm weather dominance of Acartia. Phytoplankton were also modeled as two separate
state variables, with functional distinction based on differences in nutrient and light response.
Temperature, nitrogen concentration and underwater light intensity regulate specific
photosynthesis rates. Phytoplankton growth rate is determined by a photosynthesis versus
irradiance (P vs I) relationship and Michaelis-Menton nutrient kinetics, with N representing the
limiting agent. Maximum growth rate for phytoplankton P^ was defined by a temperature-
dependent formulation (Kremer and Nixon 1978):
P^ 0.59 e«M)633*
-------
where PAR= light intensity or photon flux density at any depth (Z) is derived from PAR,, the
incident light entering the water column. Downwelling attenuation (KD) is summed from
component attenuations of phytoplankton (Kc), suspended particulates (Ks), and water (Kw).
Equations use the form: KC= 0.054 * Chl0667+ 0.0088 * Chi, Ks= 0.0396*SPM+0.39, and Kw=
0.03. K is in units of m"1, Chi in micro-g L"1, and suspended particulate material (SPM), mg L"1.
Self-shading effects on phytoplankton photosynthesis rates are accounted implicitly by the
phytoplankton contribution (K^) to total water column attenuation (KD). Active chlorophyll
concentrations are based on chl:c ratios for diatoms of 1 : 50 and for other phytoplankton, 1:90.
Maximum temperature-based photosynthetic rates are downward-adjusted by limitation terms
describing instantaneous nutrient and light-based carbon uptake. Carbon fixation rates for
phytoplankton are modeled as functions of water-column nutrient concentrations following
Michaelis-Menton kinetics for nitrogen (N), phosphate (P), and silica (Si), using half-saturation
constants of 1.5, 0.9, and 1.5 M, respectively (Monod 1942). Photosynthesis by phytoplankton is
related to irradiance (PARs, mean irradiance in the water column) by rectangular hyperbola P-I,
with a light saturation onset parameter (1^ of 90 Bin m"2 s"1. The minimum of all nutrient and
light - related rates was taken as the overall limiting rate for phytoplankton growth:
F(N, L) = MIN (N_ , P , Si , PAR; ) (3)
N+K. P+Kp Si+IQ
Growth in the zooplankton compartments is the sum of grazing inputs, modified by an
assimilation coefficient minus losses to predation, respiration, sinking and mortality. Conversion
of zooplankton biomass to abundance of adult individuals per m3 used a conversion constant of 5
micro-g carbon per adult (Roman, pers. comm.). Grazing is characterized by an Ivlev function
based on predator and prey abundances. A separate Ivlev formulation was used for each
zooplankton species for grazing on phytoplankton, incorporating temperature and oxygen effects
(example for Acartia):
Gap=ACAR*(l - (e(-0-<°'«^ + 0<*°PH°)) * Tf * (02/ (1 + O2)) (4)
where Gap is the ingestion rate by Acartia of phytoplankton, ACAR = the biomass of adult
Acartia; 0.001 is the Ivlev coefficient; DIA = the biomass of diatoms and OP = the biomass of
other phytoplankton (mg m"3), k is a preference function based on the seasonal availability of each
phytoplankton group, where DIA is 50% more likely to be consumed during spring-early summer
(DAY<190) and 25% more likely to be consumed in late summer though winter; Tf = an
Arrhenius temperature function with Q10 of approximately 2 and of the form:
Tf=SQRT(e(Temp*ao69)*0.5) (5)
and the final term is an oxygen relationship which inhibits feeding at low oxygen concentrations,
having a 1 mg L"1 half-saturation coefficient. A similar grazing function was used to represent
zooplankton predation on microzooplankton (including juveniles):
Gam=ACAR*0.075*(l - (e^117*^-5)))* Tf * (027 (1+02)) (6)
Zooplankton Process Model for Mesohaline Chesapeake Bay 81
-------
where Gam is the ingestion rate byAcartia of microzooplankton, MZ = the biomass of small
zooplankton in mg m"3. Mesozooplankton grazed on detrital material according to the following
relationship:
Gad=ACAR*0.005*(SPOC/(1000+SPOC)) * Tf * (O2/ (1+O2)) (7)
where Gad is the ingestion rate byAcartia of detritus, SPOC is suspended particulate organic
carbon (mg m'3).
The feeding relationships for Eurytemora were similar, using an Ivlev coefficient (0.007) and
the temperature coefficient adjusted to suppress feeding by 50% after mid-year. This formulation
has the effect of enhancing Eurytemora feeding and growth during the early part of the year and
restraining it during mid-year whenAcartia dominates.
Losses from the zooplankton compartments include respiratory expenditures, natural
mortality, and grazing. Respiration includes a basal metabolic rate that is temperature-driven, plus
an active respiration rate related to feeding rate, as 15% of the carbon ingested, of the form:
Rz=ACAR*0.055*Tf+0.15*I (8)
where Tf is the above Arrhenius function and 1= the ingestion rate (sum of grazing rates from eqs.
4, 6, 7). Losses of zooplankton carbon due to grazing by planktivorous fish (anchovy) are related
to zooplankton abundance from the model calculations and fish abundance provided by forcing
function following the familiar Ivlev formulation:
Gfz=F*R3*0.11*(l
where Gfz is grazing rate ofAcartia by anchovy, F = anchovy abundance; R3 = a randomizing
term to impart a 25% variability around the forcing function; -0.01 is the literature value Ivlev
coefficient, empirically determined and ACAR = prey abundance (mg m"3).
Zooplankton egg production was related to feeding rate and temperature following the
relation:
Egg=I*e-°-01'a)*eT*ao84 (10)
where Egg is the number of eggs produced per mg of zooplankton biomass, I is the ingestion rate
of the zooplankton. Recruitment to naupliar, then copepodid and adult stages, are functions of a
characteristic development time (maximum=24 d), modified by temperature and losses to natural
mortality and grazing.
Results and Discussion
The model is currently in the mid-development phase and has undergone successful
preliminary calibration. Initial results showed stable model behavior, reflecting appropriate
seasonal dynamics and feeding relationships. Calibration datasets for zooplankton were taken
from mesohaline and Patuxent River monitoring data from 1964 and 1985 (Flemer 1970, Heinle
1974, MDE 1988). Calibration data for phytoplankton were from several sources reporting
82 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
measurements in the Patuxent and Choptank Rivers in 1970 and the 1980s (Flemer 1970,
Boynton etal. 1982, Flemer etal. 1983, MDE 1988, Magnien etal. 1991, Kemp etal. 1981,
1995).
Seasonal phytoplankton and zooplankton abundances in the single layer water column model
reproduced the spring period of dominance by diatoms and Eurytemora, yielding to summer
dominance by flagellates (other phytoplankton) and Acartia. Model output for Eurytemora and
Acartia adults followed calibration data from 1964 reasonably well. A peak in Acartia abundance
in April-May in three calibration datasets was reproduced by the model. A peak in August
predicted by the model corresponded well with a small summer peak in the 1964 dataset for the
mesohaline Patuxent. The biomass peak also coincided in timing, but not in magnitude, to the
1985 mesohaline and transitional Patuxent dataset.
These results showed the model to be stable and within the range for Chesapeake Bay
zooplankton populations. Further model refinement will likely improve matching of the recent
trend of higher zooplankton abundance in summer and fall. The annual pattern of Eurytemora
abundance was captured by the model, with a strong peak in March, declining to about the level
of Acartia abundance during summer and fall. Eurytemora does not appear to exhibit a secondary
peak, as does Acartia during late summer.
Seasonal phytoplankton and zooplankton abundances show that during an annual cycle
baseline run, phytoplankton components of the aggregated model in both the upper and lower
water column exhibited similar seasonal patterns of diatom and flagellate dominance as in the
shallow littoral model. Some significant differences in vertical distribution were noted. In both
phytoplankton groups, upper water layer compartments were generally 10-25% higher in mean
algal biomass.
In the stratified model, Eurytemora and Acartia biomass in upper and lower water columns
also reflected similar patterns as in the shallow model. However it was noteworthy that the
second Acartia peak, which developed in late summer, occurred only in the upper water column.
Additional analysis is required, but preliminary assessment ruled out food limitation, because in
late summer, the upper and lower phytoplankton compartments were almost equally abundant.
Numerical abundance of Acartia naupliar and copepodid stages in the upper layer followed adult
patterns of abundance, indicating that during times of increased food availability and adult growth,
there was increased egg production. The sum of all naupliar stages was about an order of
magnitude greater than the adult density, and the sum of copepodid stages about 40% of the
nauplii. The order of magnitude of copepodid abundance agreed well with the 1964 calibration
data from Patuxent for most of the year, although the summer peak exceeded calibration data by
about a factor of four and the spring peak in the model precedes the data by about one month.
Model-based analysis enables quantification of cause-effect relationships and the examination
of emergent properties in the data. The aggregated two-layer model was used to analyze top-
down control of phytoplankton biomass by varying the amount of predation pressure from
Anchovy on the adult zooplankton components. Anchovy predation (a forcing function) was
varied by changing predator abundance in successive runs to equal 10, 40, 80, 100, 120% of the
baseline (1988) case (Figure 2), while maintaining constant feeding pressure offish on
zooplankton. Increased predation pressure on zooplankton reduced grazing on phytoplankton.
Approximately 10 mg C m"3 increase in phytoplankton biomass was realized per 10% increase in
anchovy predation due to top-down control. A doubling of Anchovy biomass resulted in an
Zooplankton Process Model for Mesohaline Chesapeake Bay 83
-------
increase in mean spring phytoplankton biomass of 80%. Model analysis focusing on the
relationships between turbidity and submerged aquatic vegetation (SAV) evaluated effects of
increased phytoplankton abundance on the underwater light regime. Model results suggested that
as phytoplankton chlorophyll a doubled, water column attenuation (KD) increased from 1.1 m"1 to
1.5 m"1. This is relevant to management of living resources, as subsurface light is a critical factor
for SAV survival. According to the SAV model, by doubling water column chlorophyll, PAR at
the plant canopy declined by about 50% (204 to 97 Ein m"2 s'1). A halving of PAH would severely
inhibit SAV growth. Thus, albiet using forced top predator data, the model demonstrates a strong
trophic cascade effect from top-down control. Use of the model to explore the magnitude of
these effects, and the relative dependence of phytoplankton concentration on top-down versus
bottom-up control will yield information useful to management of the mesohaline Bay.
Mean spring phytoplankton biomass vs. anchovy abundance
(preliminary model results)
600
£ 550
O
o»
e SOD
450
I 400
jt.
c
p. 350
.c
°- 0
0 10 20 30 40 50 60 70 80 90 100 120
Anchovy abundance (% 1988 level)
Madden & Kemp 1995
Figure 2
84
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
References
Boynton, W. R., W. M. Kemp, C. W. Keefe. 1982. A comparative analysis of nutrients and other
factors influencing estuarine phytoplankton production. In: V. S. Kennedy (ed). Estuarine
Comparisons, p 69-89.
Flemer, D. A. 1970. Primary productivity in the Chesapeake Bay. Chesapeake Science 11:117-
129.
Flemer, D. A., MacKiernan, G., W. Nehlsen, and W. K. Tippie. 1983. Chesapeake Bay: A profile
of environmental change. USEPA. 200 pp.
Heinle, D. R. 1974. An alternative grazing hypothesis for the Patuxent estuary. Chesapeake
Science 15:146-150.
Kemp, W. M., A. Hermann, and W. R. Boynton. 1981. Resource dynamics and ecology of
submerged aquatic vegetation in Chesapeake Bay: A modeling approach to demonstrate
resource management concepts. Univ. MD Center Environ. Estuar. Stud. Red No. HPEL-81-
214, Cambridge, MD.
Kemp, W. M., W. R. Boynton, and A. J. Hermann. 1995. Simulation models of an estuarine
macrophyte ecosystem. In: B. Patten, (ed.) Complex Ecology. Prentice Hall.
Kirk, J. T. O. 1983. Light and photosynthesis in aquatic ecosystems. Cambridge Univ. Press.
Cambridge, UK. 401 pp.
Kremer, J. and S. W. Nixon. 1978. A Coastal Marine Ecosystem: Simulation and Analysis.
Springer-Verlag, New York. 217 pp.
Magnien, R. E, R. Eskin, R. Hoffman, and T. Parham. 1991. Water quality characterization
report for the 1991 re-evaluation of the Chesapeake Bay nutrient reduction strategy. 107 pp.
MDE-Maryland Department of the Environment. 1988. Chesapeake Bay Water Quality
Monitoring Program Mesozooplankton Component. Versar, Inc.
Monod, J. 1942. Recherches sur la croissance des cultures bacteriennes. Paris: Herman et Cie.
Roman, M. Personal communication. Horn Point Environmental Laboratory, Cambridge, MD.
Zooplankton Process Model for Mesohaline Chesapeake Bay 85
-------
A Stage-Structured Model of Zooplankton Dynamics in Chesapeake Bay
Christopher J. Madden, W. Michael Kemp, and Flemming M0hlenberg
Horn Point Environmental Laboratory
University of Maryland
Cambridge, MD
Background
A model of zooplankton-phytoplankton interactions was developed to investigate aspects of
carbon flow and trophic transfer relationships in Chesapeake Bay. In estuarine systems, the
zooplankton play an extremely important role in determining nutrient recycling rates,
phytoplankton size, species composition, productivity, and rate of transfer of carbon to higher
trophic levels. However, the zooplankton community is highly dynamic and the range of impacts
can be highly variable, both seasonally and inter-annually. This zooplankton community dynamics
model was developed as part of the Chesapeake Bay Ecosystem Process Model initiative as a
means of studying and quantifying the impact of the zooplankton on lower trophic levels, nutrient
cycling, and size-frequency distributions of their prey. The important feature of this model is the
highly disaggregated zooplankton community, which is monitored by developmental stage.
Generally, estuarine simulation models follow the zooplankton using one or two aggregated state
variables, which follow juvenile and adult components. Madden et al., Zooplankton Process
Model..., this report, describes an aggregated, three-component model of the Chesapeake Bay
zooplankton community which includes state variables for Eurytemora affinis, Acartia tonsa, and
generalized microzooplankton.
This second zooplankton modeling effort was initiated to address specific questions about
how zooplankton community development affects upper and lower trophic levels, and to provide
guidance for monitoring efforts in Chesapeake Bay. It is expected that this model will also
interface with bioenergetic models of Chesapeake Bay fish communities (Figure 1), allowing
estimates of food availability and preferability for desirable fish species under a variety of
environmental conditions. Ultimately this sort of model linkage could permit the prediction of
potential fish production based on primary production rates. The current high-resolution
zooplankton model will continue to be developed in parallel with the three-compartment model
and both will be used in tandem to address different questions related to research and
management.
The Chesapeake Bay zooplankton community is highly diverse in both function and size.
Functionally, zooplankton community feeding rates vary depending on abundance and species
composition of the zooplankton, water temperature, prey type and abundance. Size of the
zooplankton community ranges from sub-mm for nauplia to several mm for adults of some
species. The importance of predator-prey size relationships is emphasized in this model, as
selective feeding will cause shifts in size and in species composition in prey. Tracking of
zooplankton size distribution in the model is also important for zooplankton-fish trophic
interactions. Fish gape size relative to prey size is critically important to year-ckss success
(Gushing 1958). The match-mismatch hypothesis of Gushing has been borne out in numerous
studies indicating the sensitivity of cohort success to coincidence of gape size and prey complex at
86 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
CONCEPTUAL MODEL OFZOOPLANKTON TROPHIC PROCESSES
COPEPOD
Interlace w/
Rsh Bnenergetic
^TeggpfQd
%= % maturation d1
e= excretion
m= mortality
r= respi ration
c= cannibalistic predation
Madden and Kemp 1995
Figure 1
the appropriate point in larval and fingerling development.
The Chesapeake Bay zooplankton database presents an excellent resource for development
and calibration of the model. Ongoing monitoring programs provide biweekly or monthly
information on species abundances at several sites which will be used for model calibration and
validation. At present, the model is uncalibrated. However, it serves an important role as a
heuristic tool for testing hypotheses even in its current state.
The Chesapeake Bay ecological model featuring zooplankton trophic dynamics is directed
toward answering the following questions:
1) What is the quantitative importance of zooplankton as food for higher trophic levels?
2) What are the critical dimensional considerations in predator-prey relationships in
Chesapeake Bay? What size resolution is needed for effective modeling and sampling?
3) What level of model aggregation is needed to achieve reliable model results and
calibration?
4) What is the quantitative impact of zooplankton on Chesapeake Bay ecosystem trophic
dynamics? On nutrient cycling? On oxygen dynamics?
Stage-Structured Model of Zooplankton Dynamics in Chesapeake Bay
87
-------
Several model simulations are planned using zooplankton biomass at different size-frequency
distributions to demonstrate the importance of feeding impacts on phytoplankton. Similarly, the
species (and size) composition of the phytoplankton will be varied in simulations, while holding
total biomass constant. This will allow observation of the impact of size structure on zooplankton
productivity.
The resolution of stage-dependent zooplankton development and number of zooplankton life
stages required to realistically model total zooplankton biomass/abundance will be investigated
using this model. It will explicitly and dynamically quantify relationships between zooplankton
and higher trophic levels, such as planktivorous fish and filter feeders. This model uses size-
frequency distributions and stage-specific abundances to estimate food availability to
planktivorous fish. The size structure of the zooplankton community determined by model
calculations will, thus, yield the means to estimate potential fish production based on fish
bioenergetic relationships. The high degree of disaggregation of the zooplankton compartment
will allow complete analysis of size-regulated transfer efficiencies to piscivorous fish. This aspect
of the model represents a first step toward linking the ecosystem process models with fish
bioenergetic models and linking nutrient loading to high level living resources.
Model Development
A numerical model was developed using the STELLA II simulation environment on a RISC-
based microcomputer. The model consists of a set of 20 difference equations representing state
variables for small and larger (>12 micrometer) planktonic algae; ciliates; rotifers; the
zooplankter, Acartia tonsa; a benthic filter feeder; dissolved inorganic nitrogen (DIN); and
sediment paniculate organic carbon (POC). Phytoplankton are modeled as two separate state
variables, allowing a distinction between functional groups based on differences in nutrient and
light response. This strategy allows the model to reflect seasonal succession of spring dominance
by diatoms and summer dominance by a mixed flagellate and chlorophyte community. A time-
step of 0.25 day, a simulation length of 1-3 yr, and a Runge-Kutta 2 integration scheme are used
for model simulations.
The model is spatially averaged, simulating aim2 unit in the mesohaline littoral region of
Patuxent River. Simulations can be run with the water column varying from 1-3 m in depth. The
water column in each case is modeled as a single homogeneous layer mixed to the bottom. The
shallow depth of the system indicates that benthic filtering of zooplankton could represent a large
mortality fraction, and a benthic community is included in the model. Benthic nutrient and carbon
compartments are calculated based on a uniform sediment depth of 15 cm. Carbon is tracked in
units of mg C nr3 and conversion of C to N follows Redfield stoichiometry of 7.2:1 for
phytoplankton.
Increases in the phytoplankton compartment include advective import and growth. Growth is
modeled as a function of nutrient and light availability. Temperature regulates maximum
photosynthetic rate and nutrient uptake is expressed using Monod functions for each nutrient.
Light regulates phytoplankton photosynthesis in the model based on P-I (photosynthesis -
irradiance) relationships specific for each phytoplankton group. Losses from the phytoplankton
community are due to advective export, sinking, natural mortality, and grazing by zooplankton.
Additional details of how phytoplankton processes are modeled are discussed in Madden and
Kemp, Simulation Model of SAV and Water Quality..., this report.
88 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
The zooplankton comprise 13 state variables, representing eggs, six naupliar, and six
copepodid developmental stages. Stage-specific feeding and mortality rates are partitioned among
each size-class, permitting refined accurate tracking of stage population and size-classes. This
allows accurate accounting of size-dependent trophic transfers, as bivalve and fish predation rates
can be highly dependent on the size-structure of zooplankton prey. Zooplankton predation on
phytoplankton is partially dependent on predator size, and the phytoplankton group that is
accessible to the zooplankton predation is determined by the relative size of the zooplankton
stage to the phytoplankton equivalent spherical diameter.
Zooplankton growth is a function of recruitment, temperature, food concentration,
assimilation efficiency, and diapause egg hatch rate. Recruitment through each of the 12 naupliar
and copepodid stages is regulated by temperature and food abundance and converted to biomass
using an average mass conversion factor for each stage; conversion factors range from 0.05-4.2
micro-g individual'1 as follows:
Egg=0.05 Cl=0.33
Nl=0.05 C2=0.46
N2=0.05 C3=0.75
N3=0.06 C4=1.8
N4=0.08 C5=2.5
N5=0.12 C6=4.2
N6=0.2
Egg production is a function of temperature and food available to the female C6 population:
Egg=0.5*AC6*5*ln(AP-40.5)*Tf (1)
Egg= #eggs produced m"3 d"1; AC6 is the total adult Acartia population; AP is the total prey
biomass available to adults (mg m"3); and Tf = an asymptotic sigmoidal temperature function
ranging from 0-1 for temperatures ranging from 5-25° C.
Recruitment (maturation of a developmental stage) by zooplankton of each size class to the
next higher class occurs based on a stage-specific development rate. This rate is translated to a
maturation time period (MT in days) during which the entire cohort occupies the specific class.
Each day, a proportion of the cohort moves to the next class based on the reciprocal of the
maturation period (ie. if MT-5 days, 1/5 of the cohort population advances daily). The rate of
development of each zooplankton stage is regulated by different factors, depending on age. For
Acartia eggs (AE) and first naupliar stage (AN1), development is dependent on temperature only:
Dl=0.85*(X)*Tf (2)
where Dl=proportion of cohort (Egg or Nl) advancing to next stage (d"1); X=number of
individuals in cohort (# m"3); and Tf is the temperature function described above.
For higher stages, development is dependent on temperature and food concentration:
D2=XDf (3)
Stage-Structured Model of Zooplankton Dynamics in Chesapeake Bay 89
-------
where D2=proportion of cohort (N2 through N6) advancing to next stage (d"1); Df is a
development fimction based on temperature and total prey availability, defined as follows:
Df=(TP0.45*(APn-40)/(APn+100)) (4)
where Df = the fraction of cohort maturation (d"1); Tf is defined above and AP=Total prey
(mg m"3) available to nauplia, assumed to be the smaller phytoplankton fraction plus ciliates. In
the larger zooplankton (C1-C6), APn is replaced by APc, the total prey complex of the
zooplankton, which includes APn plus larger phytoplankton forms and rotifers. Thus it is possible
(and the norm) to have each developmental stage maturing at different rates in the model. We
expect this represents what is actually occurring in nature. The above expression yields a
Michaelis-Menton type function for development rate ranging from 0-0.45 d"1 with a threshold at
40 mg m"3, half-saturation at 100, and 95% maximum at about 2500. This is the percentage of
the cohort state variable (numerical abundance) that passes to the next stage each day. The
maximum rate of 0.45 corresponds to a 2.2 day development time for the stage, realizing a 45%
recruitment rate for all individuals in the stage per day. Lower prey concentrations result in
slower growth, extending recruitment time to infinity. At maximum rates, the entire life cycle is
completed in about 24 d from egg to adult.
Zooplankton losses occur from natural mortality, benthic grazing, zooplankton
cannibalism, ciliate grazing, respiration, sinking, carbon exudation, and advective export. Losses
from each stage combine natural mortality, mortality from filter feeders, and planktivorous fish
grazing terms. Because the empirical relationship between food concentration and Acartia
development rate and recruitment was developed using the number of individuals successfully
recruited to the next stage, the formulation yields an extra accumulation of production.
Respiration is calculated from these results for purposes of balancing nutrient and carbon flows in
other parts of the model, but is not used in calculating carbon loss from the Acartia stage, since
that is already incorporated in the empirical calculation of cohort recruitment.
Grazing loss (Mzgr) from the zooplankton compartment via planktivorous fish (anchovy)
predation relate predator and prey abundance following an Ivlev formulation:
Mzgr=Fa*R3*0.11*(l-e(-°01*ACAR)) (5)
where Fa represents predator abundance (g C m"3); R3 is a randomizing term to impart a 25%
variability around the anchovy forcing function; and ACAR = total adult zooplankton abundance.
Water column physics are incorporated in the model using relationships of current velocity
and wind speed to calculate turbulent energy dissipation rate in a shallow water column.
Turbulence in the model has two functions: it modulates predator-prey contact rates via an
empirical relationship and it determines the depth of the mixed layer. Turbulent dissipation (TD)
is calculated as a quadratic of wind speed, modified by a turbidity index and the air-water thermal
gradient as follows:
TD= ((0.134*Wind2-0.000347*Wind+0.00045) (6)
Depth (m)
where Wind= wind velocity in m s"1.
90 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
An additional term, [ PAR'7*ln(Kd+l)+ (TW-TA)*2E'S ], is subtracted from the above
expression to account for the effect of increasing water stability with increasing paniculate
concentration. Kd is the downwelling attenuation coefficient (measuring turbidity) in nr1; PAR=
mean integrated PAR in Ein d"1; TW-TA= the gradient between water and air Celsius temperatures.
An increase in thermal mass of the upper water column contributes to water column stability.
When the model is run using 1-3 m water column depths, the stability term is not relevant because
the mixed layer always reaches the bottom, however, future expansion of the model into the
deeper tributary channels and mainstem Bay requires this term (in conjunction with a sigma term
to account for salinity driven density differences) to adjust the depth of the mixed layer and the
degree of contact of the benthic community with the water column.
Results and Discussion
This model is in the process of being calibrated for Chesapeake Bay. Initial pre-calibration
runs indicated that carbon flows balanced properly and that the model was stable and within
normal ranges for a temperate estuarine system. Forcing functions in the model for preliminary
diagnostic execution runs were established using a generalized nitrogen loading cycle similar to
that for Chesapeake Bay and initial conditions for state variables based on data from Flemer
(1970). The model was run for one, two and three years until stable. Output agreed well with the
range expected for temperate estuarine systems.
In standard diagnostic run of the model, a strong pulse of DIN (by forcing function)
persisting from January through March months entered the system via runoff from winter
precipitation and spring runoff. Water column DIN concentration increased in response to
freshwater inputs, peaking in March at 20 uM. Chlorophyll a initially remained at low initial
levels, but bloomed in mid-April as light and temperature increased to thresholds required for
phytoplankton (particularly diatom) growth; chlorophyll attained values of more than 60 //g I/1
during the spring bloom phase. After the bloom phase, phytoplankton were severely reduced due
by top-down control. Analysis of the spring phytoplankton decline indicated that initially, losses
in March were due to grazing by pelagic herbivores (Acartia, ciliate and rotifer populations), and
benthic filter feeders. In early summer, however, DIN often fell below half-saturation
concentrations for both spring and summer phytoplankton functional groups. Because the
diatoms dominated during spring bloom and smaller phytoplankton were negligible, the adult
zooplankton were responsible for all of the pelagic grazing losses through the zooplankton
pathway. Juveniles were limited to consuming the smaller phytoplankton fractions and detrital
paniculate carbon (POC) due to size constraints. Summer phytoplankton group losses were
distributed across all zooplankton groups because the smaller average cell size permitted grazing
access by all zooplankton (by model parameterization).
During the spring phase of increasing DIN concentration and phytoplankton abundance, a
period of rapid growth in the benthic macrofauna compartment occurred, reflected in an increase
in bivalve biomass by 50% within a month. After this initial growth period, bivalve biomass
steadily declined throughout the summer and fall, despite regularly responding to increases in
phytoplankton biomass throughout the summer with short bursts of filtration activity. These short
bursts of filtration effectively reducing phytoplankton concentrations and stimulating brief periods
of bivalve growth. These episodes occurred periodically, at 5-10 d intervals, as pulses of primary
production responded to releases of regenerated nitrogen and intervals of grazing pressure by
Stage-Structured Model of Zooplankton Dynamics in Chesapeake Bay 91
-------
bivalves decreased as particle concentrations fell below feeding thresholds. Riverine loading of
water and nutrients (via forcing function) was reduced during the summer months and dissolved
nutrient concentrations remained low, often limiting phytoplankton growth through November.
Under conditions of low summer nitrate input and efficient benthic grazing, model calculations of
water column turbidity declined during the growing season, with chlorophyll concentrations
ranging to below 10 ug I/1. In preliminary runs without filter feeders, the model showed high
levels of phytoplankton productivity, which were ultimately became by nitrogen and biomass in
excess of 65 \ig I/1, indicating a large degree of top-down control exerted by benthic organisms
when present.
An experiment to determine the effect of prey type on zooplankton size distribution was
performed by comparing abundances of all stages of the zooplankton (Acartia sp.) population
under conditions of 50% reduced diatoms, versus 50% reduced chlorophytes. The term
"chlorophytes" was used to represent the smaller summer dominant functional group of
phytoplankton in the model. Comparison of results of a model run having reduced spring diatom
abundance with the reduced chlorophyte run revealed that spring-dominant diatoms seemed less
important in determining the ultimate success of the zooplankton community in summer.
Reasons for a stronger reduction in one model population over the other are complex, but
a simple explanation may present itself in comparison of output of the juvenile stages in both runs.
The initial spring increase in zooplankton abundance was nearly identical in both graphs,
indicating that during spring, juvenile zooplankton were probably not limited by food
concentrations and may be controlled by grazing, temperature, or initial conditions of adult
population and egg abundance. Adult zooplankton seemed slightly reduced in the spring period
under the reduced diatom condition. As the summer bloom developed, the two treatments
markedly diverged in the adult stages: reduced chlorophytes led to reduced abundance in all
zooplankton stages in summer compared with the reduced diatom run. Obviously, because
diatoms were concentrated in spring, the forced reduction resulted in lower food availability
during spring only; whereas, forced reduction of chlorophytes reduced summer food availability.
This was the period when zooplankton were most active and respiratory expenditures were
exponentially higher than in spring. This simple model result may indicate an important feature of
food availability controlling zooplankton populations via bioenergetics. This represents an area of
ongoing study using the model.
Examination of model output of zooplankton abundances by life stage provides insight to
rates of zooplankton maturation and development, and may address the above concept of
energetics. One important diagnostic of appropriate model functioning is the ability of changes in
abundance of younger stages to propagate through the zooplankton stage state variables. Model
output indicated that changes in earlier stages were expressed in later stages with appropriate time
delays. Moreover, the lags between individual peaks in abundance shorten from spring to
summer, then began to lengthen toward fall as temperature and food availability altered the rates
of development of each stage. This is simply illustrated with lines drawn which connect cohort
peaks in response to population events at various times through the year. The first spring event
noted propagated through the population slowly, requiring more than 10 d to reach C-6, while the
summer events require 1-3 d at most. These rates reflected total development times (egg to adult)
of well over 20 d during the cooler months to less than 10 d in mid-summer. The rates of
turnover corresponding to these development times were indicative of the energy required from
food to maintain the population and may partially explain the sensitivity of the zooplankton
population to food concentration in summer in the experiment described above. Of course, the
92 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
population curves were shaped by other forces in addition to development time, such as
differential grazing based on size, but the likelihood was that such top-down forces will adjust the
size of a peak, but the bottom-up forces control the existence of the peak. The model appeared
useful for addressing such questions of trophic control in detail.
Management Implications
Initial results indicated that the highly resolved zooplankton model presented here will be
an important complement to the aggregated zooplankton model. Each will allow the investigation
of different questions. Eventually, this model should provide a means of linking size-frequency
distributions of phytoplankton to feeding efficiencies of zooplankton. The ability to track
individual zooplankton stages gives the opportunity to dynamically calculate the size-frequency
distribution of the zooplankton community. This will allow the determination of the prey and size
classes available to subsequently higher trophic levels such as planktivorous fish.
Stage-Structured Model of Zooplankton Dynamics in Chesapeake Bay 93
-------
References
Gushing, D. H. 1958. The effects of grazing in reducing the primary production: A review.
Cons. Explor. Mer 144:73-75.
Flemer, D. A. 1970. Primary productivity in the Chesapeake Bay. Chesapeake Science 11:117-
129.
Nomenclature
EGG= AcartiaEggs
AN1-AN6= Acartia Nauplia, stages 1 through 6
AC1-AC6= Acartia Copepodids, stages 1 through 6, AC6= Adult Acartia=Total
Acartia biomass
AP= Total prey biomass
Tf= Temperature function
Df= Development function
TD= turbulent energy dissipation rate
AEclr= Acartia egg clearance rate by filter feeders
AEmort= Acartia egg natural mortality rate
Xn= Numerical abundance of Stage n
n= Stage number, for N1-C6 n=l-12
Aclr= Acartia mortality clearance rate by Filter feeders
Amort= Acartia natural mortality rate
Astrv= Acartia starvation mortality
Adens= Acartia density dependent mortality
Afg= Acartia grazing mortality by fish
K= Stage specific modification coefficient for reduction of grazing mortality rate due to
increased avoidance of Filter feeder clearance
Stage mortality modify= -0.0377*x, where x=
Nl -0.0377
N2 -0.0754
N3 -0.1131
N4 -0.1508
N5 -0.1885
N6 -0.2262
Cl -0.2639
C2 -0.3016
C3 -0.3393
C4 -0.3770
C5 -0.4147
C6 -0.4524
94 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay:
Model Simulations of Responses to External Perturbations
W. Michael Kemp and Richard D. Bartleson
Horn Point Environmental Laboratory
University of Maryland
Cambridge, MD
Background
Many of the important biological and chemical processes occurring in estuarine ecosystems
are controlled by exchanges between planktonic and benthic subsystems (Kemp and Boynton 1990).
In most estuaries, autotrophic processes, which tend to be limited to the upper portion of the water
column, depend in part on nutrients recycled from benthic metabolism (Kemp and Boynton 1984;
Boynton and Kemp 1985). On the other hand, heterotrophic activities, which are concentrated at the
sediment surface, depend on the deposition of paniculate organic matter (POM) from the overlying
water (Hargrave 1973). Quantitatively, nutrient recycling from the benthos is typically sufficient to
support more than half of the primary production in the overlying water (Nixon 1981); whereas,
growth of benthic animals and metabolism of sediment bacteria in deeper portions of estuaries is
largely dependent on inputs of organic matter from plankton above (Parsons et al. 1984). The
mechanistic interactions that control exchanges of nutrients and POM in estuaries involve a complex
assemblage of geochemical, biological and physical processes. These interactions are not readily
understood using the information typically generated in field and laboratory experiments.
In recent years, understanding of key processes controlling pelagic-benthic interactions has
increased. Dynamic behavior of the benthic subsystem is controlled by numerous factors including
inputs of POM from above, sediment redox conditions, as well as macrofaunal grazing and
bioturbation (Kemp and Boynton 1990). Rates of deposition of POM are regulated by rates of
phytoplankton production and, equally important, by rates of consumption by planktonic bacteria,
protozoa and copepods that live in the overlying water (Hargrave 1973). Rates of nutrient recycling
are regulated by bacterial decomposition, availability of oxygen, and the activities of benthic
macrofauna (Nixon and Pilson 1983). In temperate estuaries such as Chesapeake Bay, seasonal
patterns of temperature (and particularly rates of vernal warming) play a major role in determining the
proportion of the algal production consumed in the water versus deposited to the sediments.
Seasonal temperature patterns also control rates at which POM deposited in winter-spring is
metabolized in spring-summer (Westrich and Berner 1984).
Chesapeake Bay is a partially mixed estuary with seasonally varying stratification and reduced
exchange between surface and bottom waters. Nutrient loading to the estuary, from the atmosphere
and watershed, results in elevated production and consumption of organic matter, which can lead to
summertime depletion of oxygen from bottom waters (Officer et al. 1984). Nutrient enrichment can
also affect the exchange of POM and nutrients between planktonic and benthic subsystems.
Increased phytoplankton production associated with higher nutrient loading tends to decrease the
proportion of that production deposited to the benthos (Oviatt et al. 1986). This may be because,
with higher production rates, planktonic consumers become more effective at removing
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay 95
-------
phytoplankton before they sink to the bottom (Margrave 1973). Under eutrophic conditions, the
structure of the planktonic trophic web may change markedly (Landry 1977), possibly leading to
shifts in the timing, quantity and quality of POM deposition (Smetacek 1984).
Eutrophication and associated bottom water hypoxia can lead to dramatic shifts in redox
conditions and in microbial metabolism (Kemp et al. 1992), as well as seasonal mortality of benthic
animals (Holland et al. 1977) and associated bioturbation of sediments. These hypoxic conditions
also inhibit the coupled benthic process of nitrification-denitrification, which provides a natural
mechanism for gaseous removal of excess nitrogen from the estuary (Kemp et al. 1990). With the
reduction of denitrification under eutrophic/hypoxic conditions, nitrogen is rapidly recycled back to
the overlying water where it can support more phytoplankton growth. Given the diversity of
interacting processes, it is challenging to predict how changes in nutrient loading to estuarine
ecosystems affect the nature of exchanges between water column and sediment subsystems.
Numerical simulation models have successfully dealt with the complexity of interactions
characterizing coastal marine ecosystems (e.g., Kremer and Nixon 1978). These models often
combine populations of estuarine organisms into aggregated variables, as well as emphasize
mechanistic detail concerning factors controlling ecosystem processes. Such models have been used
to examine various aspects of aquatic ecology, including: planktonic trophic interactions (Pace etal.
1984); planktonic nutrient cycling controls on production (Scavia 1979); dynamics of organic carbon
and dissolved oxygen (Summers 1985); factors regulating POM deposition from planktonic systems
(Andersen and Nival 1988); sediment nitrogen cycling (Billen and Lancelot 1988); and benthic
suspension-feeding effects on phytoplanktonic biomass (Gerritsen et al. 1994). None of these models
focused on pelagic-benthic coupling by combining processes associated with simultaneous flows of
organic carbon, inorganic nutrients and dissolved oxygen.
An ecosystem simulation model was used to examine the role of planktoriic-benthic
interactions and system responses to bottom-up and top-down external perturbations in the
mesohaline region of Chesapeake Bay. The primary objectives were to investigate the system
responses to changes in: 1) POM deposition; 2) bottom-water oxygen concentrations; and 3)
bottom-up and top-down external perturbation.
Model Development
Overall Model Structure
The model was developed using finite difference equations which describe the temporal
rates of change of the system state variables. The simulations were conducted on an Apple
Macintosh computer using the STELLA simulation software (High Performance Systems, Inc.).
Simulation time for an annual cycle was approximately 7 minutes when integrations were
performed by Euler rectangular approximation using a time-step of 3 hours.
This model simulated ecological processes for an average square meter area of the mesohaline
Bay. This square meter was separated into three layers in the vertical dimension: upper photic layer
of the water column; lower aphotic layer of water; and sediments. The total volume modeled was
12.19 x 109 m3. Average depth of the upper layer was 10 m in deep water and 4.4 m where depth
was less than 10 m. The lower layer depth averaged 6.08 m. The surface area was 527 x 10s m2.
Mass movement of materials into and out of the system, and through the pycnocline, was determined
by mass balance calculations and adjacent upstream and downstream salt concentrations.
96 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
The planktonic-benthic interactions model included 40 state variables detailing planktonic and
benthic processes in the mesohaline portion of the Chesapeake Bay. The water column was divided
by the pycnocline into two sections. The ecosystem was defined in terms of its state variables,
external forcing functions and interactions among variables. The equations used to define this model
employed well-established relationships from the scientific literature. For example, ammonium effects
on nitrification were described using a hyperbolic Michaelis-Menten expression (Henriksen and Kemp
1988). The use of such widely accepted relationships in structuring the model insured a high degree
of generality, making the model applicable for a broad range of conditions and Bay regions.
Pelagic state variables were represented in both upper and lower layers. Phytoplankton (PP)
were divided into two groups, large diatoms (DIA) and other algae dominated by nannoplankton
(OP); these two groups had distinctly different seasonalities and physiological characteristics.
Phytoplankton assimilated dissolved inorganic nitrogen (DIN) and were grazed by copepods (ZP) and
protozooplankton (MP). Protozooplankton fed mostly on smaller phytoplankton (OP) (Heinbokel
1978) but also grazed diatoms, bacteria (Sherr and Sherr 1987) and suspended particulate organic
carbon (POC). Copepods also consumed protozoa and POC. Bacterioplankton (BP) assimilated both
dissolved organic matter (DOM) generated from phytoplankton exudation, as well as ZP and MP
excretion and sloppy feeding. BP also assimilated DIN in varying stoichiometries, depending on the
DOM sources. Ctenophores and fish fed on copepods and phytoplankton. Respiration in the water
column consumed oxygen.
The benthic subsystem was averaged over the upper 10 cm of sediments, where most of
the biological activity was concentrated. Particulate organic matter (POM) was separated into
labile (LPOM) and refractory (RPOM) fractions, based on the ability of the organic matter to be
decomposed (Billen et al. 1988). POM was deposited to the sediments from the overlying water
in proportion to the POC and phytoplankton pools. The POM was consumed by deposit-feeding
macrofauna and by both aerobic and anaerobic bacterial respirations. Rate coefficients for
biological utilization of labile and refractory fractions were based on geochemical experiments
(e.g. Westrich and Berner 1984). Suspension-feeding macrobenthos fed on POC in the bottom
layer of the water column and in the portion of the surface layer that was in direct contact with
the benthos. Both macrofaunal groups excreted feces to POM and nitrogenous wastes to
dissolved inorganic nitrogen (DIN), as well as consumed oxygen in respiration. These metabolic
processes were described by hyperbolic functions of food consumption and biomass. Iron and iron
sulfide complexes were included because they played a role in sediment oxygen flux. Benthic
algae played a role in oxygen production and nutrient uptake at the surface of shallow euphotic
sediments.
Forcing functions included photosyntheticly active radiation (PAR) (Pers. comm., Tom
Fisher, Horn Point Environmental Lab, University of Maryland), particulate carbon, Si
concentrations, water temperature and salinity. For the most part, these data were generated
from 1986-1989, for the same region of the mainstem Bay (between the mouths of the Choptank
and Patuxent Estuaries). Studies were supported by the Maryland Department of the
Environment (MDE) Bay Monitoring Program, National Oceanic and Atmospheric
Administration (NOAA)/Sea Grant, National Science Foundation, and the Environmental
Protection Agency (EPA)/MDE Sediment Data Collection Program.
Process Description
Phytoplankton carbon fixation was modeled as a function of temperature, light, and nutrients.
The effect of temperature was assumed to be exponential with a QIO of 1.7 for diatoms and 2 for PO.
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay 97
-------
This formulation was used by Kremer and Nixon (1978) and supported by several published values
(Fasham et al 1983, Bannister 1974). Enzyme inhibition occurs at high temperatures among single
species, but the assumption was made that it didn't occur among the whole assemblage of
phytoplankton at normal summer temperatures. Because phytoplankton size varied inversely with
temperature (Malone et al. 1991), an allometric formulation was used for this relationship. The
phytoplankton community is capable of adapting to the seasonal change in temperature and growth
may be independent of temperature (Sheridan and Ulik 1976, Davison 1987). Maximum growth rates
and temperature coefficients for PD and for PO are calibrated within the range of reported values
(Tailing 1957, Eppley 1972, Ojala 1993, etc.). Production of picoplankton (< 1-3 //M) is low (1 ^g C
//g Chi'1 h'1) below 20° C but exceeds 10 //g C //g Chi'1 h'1 at higher temperatures (Malone et al. 1990).
The effect of light on maximum photosynthesis was calculated by using a formulation for
depth-integrated photosynthesis modified for the effect of self-shading. This expression did not
incorporate photoinhibition, which was assumed to be unimportant in the relatively turbid waters of
the mid-bay. An attenuation coefficient for self-shading of 0.002 per mg phytoplankton C was used
(Steemann-Nielsen 1962). Light availability at the base of the pycnocline was assumed to be negligible
and if mixing rate was small, production also would be negligible.
The effect of dissolved nitrogen (DIN) on photosynthesis was modeled as a hyperbolic
function with a half-saturation coefficient (K,) of 10 mg nv3 for PO (Goldman and Glibert 1983
obtained a K, value of 7 //g I'1) and 30 mg nr3 for PD (Scavia 1980). Diatom growth additionally was
controlled by silica concentrations (forcing function) with a K, of 0.2 //M (Conley and Malone 1992).
Inputs of phytoplankton were determined by assuming a ratio of 70 mg C mg Chr1 (Malone 1982).
Phytoplankton losses occurred due to respiration; grazing by fish, copepods and protozoa
(Capriulo and Carpenter 1983, Burkhill etal. 1987, Gallegos 1989); sinking; exudation; natural
mortality; and downstream export. Respiration was modeled as a function of temperature (Scavia et
al 1976), biomass and production. Reported specific respiration rates range from 0.02 to 1.2 d'1
(Geider 1992), or about 10% of maximum photosyntetic rate (Parsons et al. 1984). Respiration rates
for flagellated species may be quite high, relative to the diatoms due to their active nature (Geider and
Osborne 1989). Grazers of phytoplankton included copepods, protozoa, and menhaden. Copepods
preferentially graze diatoms and, although protozoans preferred other algal taxons (Burkhill et al.
1987), they still ingested significant amounts of diatoms. The model assumed that copepods preferred
diatoms 5 to 1 over PO. Copepod grazing loss included a percentage due to inefficient feeding.
Menhaden biomass was a forcing function that peaked at 50 mg C nv3 in November 1989 (Brandt
1990) and consumed from 5 to 40% of their weight per day, depending on temperature.
To establish initial values for model variables and equation coefficients, data were compiled
from published and unpublished reports on research in the mid-Bay. For the most part, these data
were generated from 1986-1989, for the mainstem Bay between the mouths of the Choptank and
Patuxent Estuaries, in the same studies cited above for the forcing functions.
98 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Results and Discussion
Seasonality of POM deposition and bottom O2 balance
Phytoplankton dominated the sinking participate matter in spring; whereas, detrital matter
dominated in summer through fall (Figure 1). This was partly due to the fact that the spring diatoms
were large and sank faster, and partly due to the deficiency of grazing in the spring. Two scenarios
examined the effect of timing of POM Carbon Sedimentation
deposition: (1) a 50% increase in inputs
of POM from upstream in winter-spring;
and (2) the same increase in summer-fall.
If sedimentation rate was increased,
production decreased, so the effect seen
wasn't exclusively the result of POM
deposition change. With increased POM
total sedimentation decreased, so the
increased winter sedimentation rate
resulted in lower sedimentation; this
resulted in slightly less ammonium flux
and higher lower-layer dissolved oxygen
concentrations. The
nitrification-denitrification equation is
parameterized to fit some seasonal data,
but effects of oxygen on this process are
not well described. Thus, related
simulations may not be realistic.
The seasonal O2 balance showed that inputs from mixing become important in the summer
when upstream concentrations were reduced (Figure 2). Pelagic respiration dominated the losses in
the summer.
Nutrient input/output balance
The seasonal nitrogen-balance from model output (Figure 3) showed that most of the
nitrogen came from river flc in the spring and from regeneration in the summer, with most of the
M
M j
Month
Figure 1
Lower Layer O2 Inputs 3utputs
Figure 2
Inputs
E3 production
E3 upstream
[5] upper layer
Outputs
gj downstream
13 upper layer
S Pelage resp
ffl Berthic uptake
net change
regeneration occurring in the water
column. As nutrient inflow increased,
outflow increased to a lesser extent.
During the summer, the effect of
increased nutrients wasn t seen in the
outflow, due to the high demand.
Responses to variations in nutrient
inputs
Increasing landward or seaward
DIN inputs by 100 mg m"3 d~l removed N
limitation. Production was controlled by
temperature, dissolved silica
concentrations and light availability.
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay
99
-------
M
Regardless of the source of additional DIN, phytoplankton biomass, sedimentation and lower-layer
dissolved oxygen all changed by the same amount. A difference may have been detectable at a lower
addition rate.
Effect ofbenthic bioturbation
The physical activities ofbenthic deposit feeders influence the oxygen concentration in the
sediments in proportion to their feeding rate. When we increased the effect of bioturbation on
sediment oxygen the nitrification-denitrification rate increased. Previous experimental and modeling
_..,... _. . studies have reported similar effects of
Dissolved Nitrogen Budget ^^ bun-owing on nitrification
(Henriksen and Kemp 1988). NH4+ and
sulfide fluxes were not affected by
bioturbation.
Effects offish grazing
Increased menhaden ingestion
reduced summer phytoplankton slightly in
summer and substantially during October
(Figure 4, left panels). Dissolved nitrogen
concentration in the upper layer was
increased slightly in early summer and
October. Sedimentation changed only
slightly before October, when it decreased
with phytoplankton biomass. Increased
anchovy feeding rates resulted in reduced copepod populations and increased spring phytoplankton
biomass and sedimentation. However, increased anchovy feeding had little effect on model behavior
the rest of the year. If fish were hypothesized to control copepod biomass and predation is decreased,
there was a decrease in spring phytoplankton and sedimentation, but an increase in fall phytoplankton
and sedimentation. Although grazing increased, the increased nutrient release stimulated production.
Effect ofbenthic suspension feeders
If oxygen concentrations remained above 2 mg I"1 the standing stock ofbenthic suspension-feeders
(oysters, clams, tunicates and polychaetes) could be much higher than the levels in the base run. The
increased activities of suspension-feeding would result in increased benthic metabolism and would
impact the dissolved oxygen concentration. We ran two scenarios to examine this impact. In one
scenario, biomass and ingestion rates of both shallow and deep suspension feeders were increased. In
the second, just the shallow suspension feeder biomass and ingestion rate were increased. Both
scenarios reduced phytoplankton biomass and sedimentation (Figure 4, right panels). When both
suspension-feeder groups were increased, the summer lower-layer oxygen concentrations were
basically the same as the base run. Thus, although sedimentation decreased, total input to the deep
sediments was not significantly changed because the suspension-feeders contributed feces to the
sediments. Lower layer oxygen increased slightly in the second scenario since inputs to the deep
sediments were increased.
M J
Month
Figure
100
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Conclusions and Future Directions
There are many interactions to consider when modeling a complex natural system, such as
Chesapeake Bay. By keeping track of carbon and providing some feedback effects, the results of this
model were more useful than "back of the envelope" computations. In the future, flows will be
determined by hydrodynamic
models and the spatial resolution Increased Menhaden Ingestion Increased FF Biomass, Growth
will be increased.
The model shows that
various processes can affect
sedimentation of organic carbon
and dissolved oxygen in the
deeper waters of the Bay.
Increasing nitrogen did not
substantially increase hypoxia,
even though sedimentation
increased 30%. This may be due
to reduced sediment respiration
rate at low oxygen
concentrations, or the lack of
feedback from the southern
portion of the Bay. Decreased
nitrogen inflows also resulted in
only a small increase in lower
layer dissolved oxygen. A larger
change occurred if particulates
were also reduced.
Fish, whether feeding on
benthos, copepods or plankton,
had little effect on product' :in or
sedimentation; whereas, be thic suspension feeders potentially had a significant effect. Because beds
of submerged aquatic vege; '.on (SAV) can also remove paniculate material from the water column,
it is possible that increased V biomass will have a similar influence. Freshwater inflow can affect
both POM deposition and n; Jng rates, which have a significant impact on lower-layer dissolved
oxygen concentrations.
Figure
Pianktonic-Benthic Interactions in Mesohaline Chesapeake Bay
101
-------
References
Andersen, V., and P. Nival. 1988. A pelagic ecosystem model simulating production and
sedimentation of biogenic particles: Role of salps and copepods. Mar. Ecol. Prog. Ser. 44:
37-50.
Bannister, T. T. 1974a. Production equations in terms of chlorophyll concentration, quantum
yield, and upper limit to production. Limnol. Oceanogr. 19: 1-12.
Billen, G., S. Dessery, C. Lancelot, and M. Maybeck. 1989. Seasonal and interannual variations of
nitrogen diagenesis in the sediments of a recently impounded basin. Biogeochem. 8: 73-100.
Boynton, W. R. and W.M. Kemp. 1985. Nutrient regeneration and oxygen consumption by
sediments along an estuarine salinity gradient. Mar. Ecol. Progr. Ser. 23: 45-55
Brandt, S.B. 1990. Acoustic quantification offish abundance in the Chesapeake Bay: applications
to power plant siting and evaluation. Progress Report to: Maryland Department of Natural
Resources, Power Plant Research Program, Tawes State Office Building, B-3, Annapolis,
Maryland 21401 {UMCEES} CBL 90-047
Burkill, P. H., R.F. Mantoura, C. Llewellyn, and C.A. Owens. 1987. Microzooplankton grazing and
selectivity of phytoplankton in coastal waters. Mar. Biol. 93: 581-590
Capriulo, G. M., and E.J. Carpenter. 1983. Abundance, species composition and feeding impact of
tintinnid microzooplankton in central Long Island Sound. Mar. Ecol. Prog. Ser. 10: 277-288
Davison, I. R. 1987. Adaptation of photosynthesis inLatninaria saccharina (Phaeophyta) to
changes in growth temperature. J. Phycol. 23: 273-283.
Eppley, R. W. 1972. Temperature and phytoplankton growth in the sea. Fishery Bull. 1063-1085.
Fasham, M. J. R., P. M. Holligan, and P. R. Pugh. 1983. The spatial and temporal development of
the spring phytoplankton bloom in the Celtic Sea, April 1979. Prog. Oceanogr. 12: 87-145.
Gallegos, C. L. 1989. Microzoplankton grazing on phytoplankton in the Rhode River, Maryland:
nonlinear feeding kinetics. Mar. Ecol. Prog. Ser. 57: 23-33
Geider, R. J. 1992. Respiration: Taxation without representation. In: Falkowski, P. G. and
Woodhead, A. D. (eds.) Primary productivity and biogeochemical cycles in the sea. Plenum
Press, New York, p. 333-360
Geider, R. J., and B.A. Osborne. 1989. Respiration and microalgal growth: A review of the
quantitative relationship between dark respiration and growth. New Phytol. 112: 327-340
Gerritsen, J., A. F. Holland, and D. E. Irvine. 1994. Suspension-feeding bivalves and the fate of
primary production: An estuarine model applied to Chesapeake Bay. Estuaries 17: 403-416.
102 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Goldman, J. C., and P.M. Gilbert. 1983. Kinetics of inorganic uptake by phytoplankton. In:
Carpenter, E. J. Capone, D. G. (eds.) Nitrogen in the marine environment. Academic Press,
New York, p.
Hargrave, B. T. 1973. Coupling carbon flow through some pelagic and benthic communities. J.
Fish. Res. Board Can. 30: 1317-1326
Heinbokel, J. F. 1978. Studies on the functional role of tintinnids in the Southern California Bight.
I. Grazing and growth rates in laboratory cultures. Mar. Biol. 47: 177-189.
Henriksen, K., and W. M. Kemp. 1988. Nitrification in estuarine and coastal marine sediments:
Methods, patterns and regulating factors. In 1. H. Blackburn and J. S0rensen [eds.], Nitrogen
cycling in coastal marine environments. John Wiley
Holland, A. F., N.K. Mountford, and J.A. Mihursky. 1977. Temporal variations in upper bay
mesohaline benthic communities: I. The 9-m mud habitat. Chesapeake Sci. 18: 370-378
Kemp, W. M., and W.R. Boynton, W. R. 1984. Spatial and temporal coupling of nutrient inputs to
primary production: The role of particulate transport and decomposition. Bull. Mar. Sci. 35:
522-535
Kemp, W. M., and W.R. Boynton. 1990. Benthic-pelagic interactions: Coupling of nutrient input
and recycling processes to plankton production and oxygen depletion. In: Leffler, M. Smith, D.
(eds.) Dissolved oxygen dynamics in Chesapeake Bay. Maryland Sea Grant College Publ.,
College Park, MD, p. 1-40
Kemp, W.M., P.A. Sampou, J.C. Cafrrey, M. Mayer, K. Henriksen, and W.R. Boynton. 1990.
Ammonium recycling versus denitrification in Chesapeake Bay sediments. Limnol. Oceanogr. 35:
1545-1563
Kemp, W.M., P.A. Sampou, J. Garber, J. Turtle, and W.R. Boynton. 1992. Seasonal depletion of
oxygen from bottom waters of Chesapeake Bay: roles of benthic and planktonic respiration and
physical exchange proct es. Mar. Ecol. Prog. Ser. 85: 137-152
Kremer, J. N., and S.W. NL n. 1978. A coastal marine ecosystem simulation and analysis. 1 ed.
Ecological Studies. Springer-Verlag, New York
Landry, M.R. 1977. A reviev of important concepts in the trophic organization of pelagic
ecosystems. Helgolander viss. Meeres. 10: 8-17
Malone, T.C. 1982. Phytoplankton photosynthesis and carbon-specific growth: light-saturated rates
in a nutrient-rich environment. Limnol. Oceanogr. 27: 226-235
Malone, T.C., and H.W. Ducklow. 1990. Microbial biomass in the coastal plume of Chesapeake
Bay: Phytoplankton-bacterioplankton relationships. Limnol. Oceanogr. 35: 296-312
Malone, T. C., H.W. Ducklow, E.R. Peele, and S.E. Pike. 1991. Picoplankton carbon flux in
Chesapeake Bay. Mar. Ecol. Prog. Ser. 78: 11-22
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay 103
-------
Nixon, S.W. 1981. Remineralization and nutrient cycling in coastal marine ecosystems. In:
Estuaries and nutrients. B. J. Nielson. and. L. E. Cronin (eds.) Humana, p. 111-138
Nixon, S.W., and M.E.Q. Pilson. 1983. Nitrogen in estuarine and coastal marine ecosystems. In:
Carpenter, E. J. Capone, D. G. (eds.) Nitrogen in the marine environment. Academic Press,
New York, p. 565-649
Odum, H. T. 1971. Environment, power and society, John Wiley Publ.
Officer, C.B., R.B. Biggs, J.L. Taft, L.E. Cronin, M.A. Tyler, and W.R. Boynton. 1984.
Chesapeake Bay Anoxia: Origin, development and significance. Science 223: 22-27
Ojala, A. 1993. Effects of temperature and irradiance on the growth of two freshwater
photosynthetic cryptophytes. J. Phycol. 29: 278-284.
Oviatt, C.A., A. A. Keller, P. A. Sampou, L.L. Beatty. 1986. Patterns of productivity during
eutrophication: Amesocosm experiment. Mar. Ecol. Progr. Ser. 28: 69-80
Pace, M. L., J.E. Glasser, L.R. Pomeroy. 1984. A simulation analysis of continental shelf food
webs. Mar. Biol. 82: 47-63
Parsons, T.R., M. Takahashi, andB. Hargrave. 1984. Biological Oceanographic Processes. Third
Edition ed. Pergamon Press, New York
Scavia, D., B. J. Eadie, and A. Robertson. 1976. An ecological model for the Great Lakes. In W.
R. Ott [eds.]. Environmental modeling and simulation. EPA 600/9-76-016 (U. S. Environmental
Protection Agency).
Scavia, D. 1979. Examination of phosphorus cycling and control of phytoplankton dynamics in
Lake Ontario with an ecological model. J. Fish. Res. Board Can. 36: 1336-1346
Scavia, D. 1980. An ecological model of lake Ontario. Ecol. Modelling 8: 49-78
Sheridan, R.P., and T. Ulik. 1976. Adaptive photosynthesis responses to temperature extremes by
the thermophilic cyanophyte Synechococcus lividus. J. Phycol. 12: 255-261
Sherr, E. B., and B. F. Sherr. 1987. High rates of consumption of bacteria by pelagic ciliates.
Nature 325: 710-711
Smetacek, V. 1984. The supply of food to the benthos. In: Fasham, M. J. (eds.) Flows of energy
and materials in marine ecosystems. Theory and practice. Plenum Press, New York, p. 517-547
Steemann-Nielsen, E. 1962. On the maximum quantity of plankton chlorophyll per surface unit of a
lake or the sea. Int. Rev. ges. Hydrobiol. 47: 333-338
Summers, J. K. 1985. A simulation model of carbon and oxygen dynamics in a reservoir. Ecol.
Model. 28: 279-309
104 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Tailing, J. F. 1957. Photosynthetic characteristics of some freshwater plankton diatoms in relation
to underwater radiation. New Phytol. 56: 29-50
Westrich, J.T., and R.A. Berner. 1984. The role of sedimentary organic matter in bacterial sulfate
reduction: The G model tested. Limnol. Oceanogr. 29: 236-249
Planktonic-Benthic Interactions in Mesohaline Chesapeake Bay 105
-------
Fish Bioenergetics:
Relating Nutrient Loading to Production of Selected Fish Populations
Jiangang Luo, Stephen B. Brandt, Kyle J. Hartman,
Melinda A. Gerken, and Michael T. Weimer
Great Lakes Center, Buffalo State College
State University of New York
Buffalo, NY
Background
Fish production depends on two factors: (1) changes in the number of individuals in the
population over time, which are caused by mortality, migration, and recruitment; and (2)
growth rates of individuals within the population. The growth rate of an individual fish is a
highly pliant, species- and size-specific response to habitat conditions and food availability.
Because fish growth is influenced by both natural variability and anthropogenic causes of
ecosystem change, it has been used as an indicator of aquatic ecosystem health. The ability of
a fish to achieve its maximum growth potential is a relative measure of fish well-being and its
ultimate survival. Growth rate is also closely linked to reproductive potential, because larger
females often produce more and larger eggs, which can enhance larval survival (Zastrow et al.
1991; Monteleone and Houde 1990).
Fish growth rates are very sensitive to water temperature and food supply (Bartell et al.
1986). Thus, bioenergetics models are useful for evaluating the effects of changes in
temperature and prey abundances on growth, consumption, and ultimately, trophic interactions
(Brandt and Hartman 1993). Bioenergetics models are simply energy budgets that can be used
to estimate how much food was consumed, based on observed growth rates, thermal history,
and diets of fish. They simulate growth rates for fish "grown" under different conditions of
diet, prey availability, or water temperature. Bioenergetics models were used to evaluate the
growth rate potential of a fish under given environmental conditions.
In the years ahead, Chesapeake Bay Program managers will require an integrated
ecological framework that can guide decisions and priorities. Information from modeling and
analysis of ecosystem processes enables managers to direct efforts toward an ecosystem-based
management program. Management scenarios produced from the ecosystem models will be
applied over the next two years to answer specific scientific and resource management
questions. Questions concerning physical habitat alterations, living resource exploitation, and
nutrient input reduction will be addressed. For example, such models could be used to
consider how implementation of the Baywide 40% nutrient reduction strategy will change the
production of recreationally important fishes in Chesapeake Bay (Kemp et al. 1992). The fish
bioenergetics models provide detailed analyses of how environmental conditions (e.g.,
temperature, food, dissolved oxygen) affect fish production and consumption. Bioenergetics
modeling techniques were recently applied to fish species inhabiting estuaries such as the
106 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Hudson River and Chesapeake Bay (Brandt and Kirsch 1993; Luo and Brandt 1993; Hartman
and Brandt 1993).
A key aspect of the ecosystem modeling program is effective integration of the
modeling approaches. The ecosystem process models use nutrient loading and other inputs to
predict the quality and quantity of food and habitat available to fish populations. The fish
bioenergetics models use these outputs as inputs to predict the potential production of striped
bass, bluefish, weakfish, bay anchovy, menhaden, spot, and white perch. The models will be
combined to directly incorporate ecological feedbacks associated with top-down control of fish
predators on their prey and on other components of the ecosystem, such as SAV and water
quality. The simultaneous evolution of the two types of mechanistic models (fish bioenergetics
models and ecosystem process models) for Chesapeake Bay provides an unprecedented
opportunity to merge and link the two essentially different modeling frameworks into a unified
management tool which encompasses ecosystem processes from nutrients and habitat quality to
the highest level living resources.
Atlantic menhaden
Model Development
The Atlantic menhaden (Brevoortia tymnnus) is an abundant and commercially
important fish in the Chesapeake Bay estuary and nearshore habitats of the eastern United
States (U.S. Nat. Mar. Fish. Serv. 1978). Menhaden account for nearly half the total east
coast fishery harvest (Peters and Schaaf 1991). The menhaden is also consumed by many
commercially and recreationally important species such as bluefish, striped bass, and weakfish
(Hartman and Brandt 1995a). On the other hand, the menhaden consume mostly
phytoplankton; it converts primary production (phytoplankton and plant detritus) directly into
fish production (growth and reproduction) (Peter and Schaaf 1991, Rippetoe 1993). The
Chesapeake Bay is the major nursery ground of Atlantic menhaden. Menhaden spawn in
coastal ocean water during later fall and winter. The larvae enter Chesapeake Bay in winter
and spring and use Chesapeake Bay as an important nursery ground in summer and fall.
The main objective of this study was to link an Atlantic menhaden bioenergetics model
with the 3-D Water Quality Model in Chesapeake Bay. The model can be used to evaluate
growth rate potential of youi;_-of-year (YOY) Atlantic menhaden in the Bay, to quantify the
habitat quality of YOY Allan, sc menhaden in the Bay, and to estimate the carrying capacity of
YOY Atlantic menhaden in the Bay.
Menhaden Foraging Model
Menhaden of 50 mm or larger are mostly filter feeders on phytoplankton (Rippetoe
1993). The consumption rate of a filter feeding organism can be estimated as:
con = (phy * vol) * eff (1)
vol = gap * u (2)
Fish Bioenergetics 107
-------
where, con = the consumption rate (g/s); phy = the phytoplankton density (g/m3); vol = the
filtration rate (mVs); eff = the phytoplankton retention efficiency; gap = mouth open area
(m2); and u = swimming speed (m/s).
The swimming speed of menhaden was modeled as function of body length and water
temperature (Figure 1). The filtration retention efficiency (Figure 2) was modeled as function
of body length. Atlantic menhaden of different sizes were sampled in Chesapeake Bay and
brought back to the laboratory. The height, the width of the mouth and the total length of
menhaden were measured with a caliber to nearest 0.1 mm. The shape of the mouth was a
close approximate of an ellipse. Therefore, the height and width was used to estimate the
mouth open area as an ellipse, then we modeled the mouth open area as an exponential
function of total body length (Figure 3).
2.5
2
s:
J» 2.0
1.5
1.0
0.5
3
0.0
Swimming speed
5 10 15 20 25
Temperature
30
Figure 1
0.50
0.40
s-
§ 0.30
.a
"5
S 0.20
ss
£L
0.10
0.00
Gill-raker retention efficiency
so eo so 120
Menhaden total length {mm}
150
Figure 2
Atlantic Menhaden Gapes
600
0 50 )00 150 200 250 300
Total Length (mm)
Figure 3
Bioenergetics Model
The bioenergetics model for Atlantic
menhaden was developed with parameters
derived by Rippetoe (1993). In the model, the
YOY menhaden were started on Julian day 152
(June 1) at 50 mm total length. The growth rate
was modeled as function of temperature, its
physiological constraint, and phytoplankton
biomass density. Temperature and fish
physiology determined the maximum potential
consumption of phytoplankton by menhaden, but
phytoplankton density determined the amount of
food that a menhaden could retain by filtering the
water through its gill rakers. The maximum
108
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
consumption for a 1 g menhaden at the optimal temperature for consumption was 1.294 g/g/d.
The exponent for the weight dependence of the consumption was -0.312. The temperature-
dependence of consumption was defined by the Thornton-Lessem algorithm (Thornton and
Lessem, 1978), which included two sigmoid curves: one fit the increasing portion of the water
temperature dependence curve; the other fit the decreasing portion. The optimal temperature
for maximum consumption was about 28 °C and the maximum limit of temperature for
consumption was about 33 °C. Respiration rate was modeled as an allometric function of body
weight, water temperature, fish activity level, and specific dynamic action. The rate for
standard respiration of a Ig menhaden at the optimal respiration temperature was 0.003301 g
O2 /g/d. The exponent for the weight dependence of respiration was -0.2246. The
temperature-dependence of respiration was defined by a non-linear function in Kitchell et al.
(1977). The slope for temperature dependence of respiration was 2.07. The optimum
temperature for respiration was 33 °C and the maximum limit of temperature for respiration
was 36 °C. The activity multiplier was modeled as a function water temperature. Specific
dynamic action coefficient, defined as the metabolic cost of digestion, absorption and
deposition of consumed energy, was
equal to 10% of assimilated energy.
Egestion was modeled as a constant
20% of consumed food, and
excretion was modeled as a constant
10% of assimilated food.
Frequency distribution
Cumulative distribution
100
§ 80
Spatially Explicit Model of Growth
Rate Potential
Traditional models of
population production are normally
based on average conditions over
large areas and assume homogeneity
and constancy of the environment.
As fragmentation and decrease of
natural habitat occurred, scientists
and managers realized that traditional
modeling techniques were n; longer
reliable tools to assess and nanage
any complex ecosystem. Recent
work suggested that spatial
patchiness of density-dependent and
density-independent processes could
significantly affect population
processes, including predator-prey
interactions, trophic efficiency, and
system-level production (Kareiva and
Anderson 1988; Possingham and
Roughgarden 1990).
60
~5 40
g 20
June 1
-004
100
SO
60
40
20
-0 02
0 00 0 02 0 04
0.06 0 08
July 1
-004 -0.02
100T
go - October i
60-
40
20
OOC 002 004
0.06 0 08
-0.04
-0.02
0 00 0,02 0 04
Growlh rate potential (g/g/d)
0 06 0 OS
Figure 4
Fish Bioenergetics
109
-------
A three-dimensional Chesapeake Bay geographic information system (GIS) program
(developed at State University of New York (SUNY) as part of the EPA project) was used to
link the ecological processes and fish bioenergetics models to the Chesapeake Bay Water
Quality Model. This spatially explicit model demonstrated how the spatial patterning of the
environment affected species-specific consumption and growth rates. In spatially explicit
modeling, the aquatic habitat was modeled as an explicit feature of environment by
subdividing the pelagic habitat into small homogeneous units that defined a cube for a
geographical coordinate system and water depth. Each cubic cell represented a small volume
of water that was characterized by a specific set of attributes, including prey density, prey
size, water temperature and dissolved oxygen (DO) that were measured in the field or were
simulated (e.g. Chesapeake Bay Water Quality Model).
The conceptual framework of the three-dimensional fish growth and production model
was similar to the two-dimensional, spatially explicit models (Brandt et al. 1992, Brandt and
Kirsch 1993, Luo and Brandt 1993, Mason and Patrick 1994). Fish production was
determined by the relationship of the supply of prey resources to the amount of prey that the
fish required for growth (predator demand). This relationship depended on the particular
temperature-dependent physiological needs and growth rate potential of the predator and on the
ability of the predator to make use of the prey supply. Process-oriented simulation models of
the same model structure were run in each cell, but were parameterized differently, according
to the habitat conditions of each cell and the specific size and species of fish being modeled.
A foraging submodel estimated consumption rate by converting measured prey densities and
sizes into prey availability for the predator. The bioenergetics submodel estimated the growth
rate of the predator from consumption rate and was based on the physiology of the predator
and habitat conditions in the cell. The result was a 3-D representation of Atlantic menhaden
growth rate potential for each day of the year. The growth rate potential integrated the
physiological response and the needs of the predator and, thus, could be interpreted in the
context of habitat quality. Volumetric analyses of the growth rates over the entire Bay defined
the proportion of the Bay volume that could support various levels of fish growth. The
frequencies below zero growth were cumulated and plotted with the dates to demonstrate the
seasonal change of menhaden growth habitat.
Spatially Explicit Carrying Capacity Estimation for Menhaden
Carrying capacity is the maximum number of individuals or inhabitants that a given
environment can support without detrimental effects. Traditionally, carrying capacity was
estimated as a singular value of a given environment. For example, we would say "the
carrying capacity of herring in the North Sea", or "the carrying capacity of striped bass in
Chesapeake Bay". This type of estimation usually could be derived from quantitative study of
energy flow in a whole ecosystem. It required information on standing stocks of the living and
nonliving ecosystem components, diets of the feeding species, and rates at which ingested
materials were used, transferred among various entities in the food web. M^ost early studies of
ecosystem dynamics assumed processes occurred over a homogeneous environment both
spatially and temporally (e.g. Baird and Ulanowicz 1989). Recent studies in spatial ecology
(Brandt et al. 1992, Brandt and Kirsch 1993, Luo and Brandt 1993, Mason and Patrick 1994)
indicated that spatial heterogeneity could have significant effects on overall results, because
most ecological processes occurred within a short time-period and at small spatial scales. For
110 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
example, a school of planktivore fish will not interact with a patch of zooplankton 10 km
away.
In this study, the carrying capacity of the Atlantic menhaden was estimated as
spatially and temporally explicit features of the environment, based on the spatially explicit
model of growth rate potential. The carrying capacity (CC,., , g/m3 ) was estimated for each
cell (/, 1 - 4073) and for each day (/', 1 - 365) as:
where, fp = the fraction of phytoplankton production (10%) that can be consumed by the
menhaden; pfy = daily phytoplankton production and biomass ratio (Nixon 1981) on day/;
phyfj = hytoplankton biomass density (g/m3 ) for cell i and on day/; cons,y = weight specific
consumption (g/g/d) for cell / and on day/; ffj(g) = a growth rate dependent sealer for cell i
and on day/ (Figure 4); and f^(do) = a dissolved oxygen dependent sealer for cell i and on
day/.
Results of Menhaden Modeling
Growth Rate Potential
100
80 _
60 r
40
"
Growth rate sealer
1.0
0.8
0.6
0.4
0.2
0.0
-O.C
/ ^
- / -
/
/
K15 0.000 0.005 0.010
Uan 1Mar 1May iJul iSep iNov iJan
Time (daytfnonifi)
Growth rate potential (cft/d)
Figure 6
Atlantic menhaden growth rate
potential varied over space and time. Growth
rate potential ranged from -0.02 to 0.08 g/g/d
over the entire bay in early summer (Figure 5,
top panel), and more than 80 % of the Bay
had positive growth rate potential. On July 1
(Figure 5, middle panel), maximum growth
rate potential was less than 0.04 g/g/d and
Figure 5
Fish Bioenergetics
111
-------
50% of the Bay had negative growth rate
potential (i.e. menhaden would loose weight if
they remained in these areas). On October 1
(Figure 5, bottom panel), growth rate potential
ranged from -0.01 to 0.04 g/g/d and over 70% of
the Bay had positive growth rate potential.
If it was assumed that menhaden would
only stay in cells with growth rate greater than a
given value, the percentage of those cells would
define the quantity of menhaden habitat. At
menhaden growth rates greater than 0.0 g/g/d,
which represents a positive growth rate (Figure
6), the percentage of menhaden habitat in
Chesapeake Bay ranged from 20 to 90 % of the
Bay and the average growth rate ranged from
0.01 to 0.02 g/g/d during summer and fall. A 50
mm (about 0.5 g ) menhaden growing under these
conditions would only reach 10 g (98 mm) by the
end of the year. At growth rate greater than
0060
a
S 0 050 -
£
0040 -
i 0030 -
I 0020 -
a
g 0.010 -
3 0000
-0010
50
5 40
«
g 30
c
o 20
I 10
iMay Uul iSep
Time (day/nonih)
iMay Uul iSep
Time (day/month;
Figure 8
Figure 7
0.005 g/g/d (Figure 7), the percentage of
habitat ranged from 10 to 75 % of the Bay and
the average growth rate ranged from 0.015 to
0.025 g/g/d during summer and fall. The same
fish would reach 20 g (120 mm) by the end of
the year. At growth rates greater than 0.01
g/g/d (Figure 8), the percentage of habitat was
reduced to 5 to 50 % of the Bay and the
average growth rate increased to 0.02-0.035
g/g/d. A 50 mm (about 0.5 g ) menhaden
growing under this condition would reach about
50 g (157 mm) by the end of the year. Trawl
survey data from Virginia Institute of Marine
Sciences (VIMS), (Bonzek el al. 1992) showed
that the sizes of YOY menhaden ranged from
120 mm to 160 mm in November and
December. Therefore, the menhaden habitat
could be defined as the cells with growth rate
potential equal or greater than 0.005 g/g/d.
112
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Figure 9
Carrying Capacity
Results indicated a great degree of
patchiness in carrying capacity, with
carrying capacity varying over space and
time. High carrying capacity occurred in
early summer, which sharply declined in
mid-June. A second peak occurred in mid-
July and gradually declined in later summer
and fall. In early summer (Figure 9, top),
carrying capacity ransedfrom 0 to 3.5 g/m3
over the entire Bay; 75 % of the bay haa
carrying capacity greater than 2.0 g/m3. On
July 1 (Figure 9, middle), 60% of the bay
had no carrying capacity; whereas, 25 % of
the Bay had carrying capacity greater than
3.0 g/m3. In fall, the capacity differences
between cells decreased; 30% of the bay
had no carrying capacity and 99% of the
bay had carrying capacity less then 2.0 g/m3
(Figure 9, bottom). Another measurement
of menhaden habitat quantity was the
percentage of the Bay with greater than zero
carrying capacity. During summer and fall,
the minimum habitat volume occurred in mid-June and the maximum habitat volume occurred
in mid-September (Figure 10, top).
Integrating the carrying capacity of each cell over the entire Bay, for each day,
produced the total carrying capacity (on weight base) of Chesapeake Bay for each day (Figure
10, middle). The weight-based carrying capacity was not a good measurement in an
ecological sense, however, because mortality and migration were operating on an individual
organism. The individual-based carrying capacity can be derived by dividing the weight-based
carrying capacity with the average weight of fish on each day (Figure 10, bottom).
Discussion
Results indicated large spatial and temporal variations in growth rate potential and
carrying capacity for the Atlantic menhaden in Chesapeake Bay. Compared with the
traditional estimate of carrying capacity (a singular value for a given environment), the
spatially explicit carrying capacity had many advantages. Spatial patchiness of growth rate
potential and carrying capacity characterized the quality and quantity of the Atlantic menhaden
habitat in Chesapeake Bay. Temporal changes in growth rate potential and carrying capacity
portrayed the dynamics of the Chesapeake Bay ecosystem.
Comparing predicted growth with the growth of Atlantic menhaden observed from
field, suggested that menhaden should choose water bodies where growth rates were no less
then 0.005 or 0.01 g/g/d. Observed differences in menhaden sizes (20 g to 50g, Figure 11,
bottom) at the end of year could be contributed to small differences in growth rates (0.003-
0.004 g/d/d, Figure 11, middle).
Fish Bioenergetics
113
-------
Cumulative distributor
Comparing predicted carrying capacity in
Chesapeake Bay from this study with estimated
recruitment of YOY (at age 0.5) menhaden
calculated for the entire East coast (Ahrenholz et al.
1987), suggested that Chesapeake Bay could nurse
the recruits of the entire Atlantic menhaden stock.
Ahrenholz et al. (1987) estimated the number of
recruits at age 0.5 for the entire stock from 1955 to
1979 from a virtual population analysis method.
The maximum number of recruits was 18.6 billion in
1958, and the minimum number of recruits was 1.5
billion in 1967. During the late 1970's, the number
of recruits ranged from 6 to 10 billion. The current
Atlantic menhaden landings were very close to
those of late 1970's. Therefore, current number of
recruits can be assumed at a similar level, about 10
billion. Spawning peak of Atlantic menhaden
occurs in December, therefore, age 0.5 occurs in
June. Results from this study showed a maximum
carrying capacity of over 150 billion in early June
and a minimum of 20 billion in mid-June (Figure
10, bottom), which was greater than the maximum
recruits (18.6 billion) of the entire Atlantic
menhaden stock. This minimum dip of carrying
capacity is the bottle-neck of menhaden recruitment level. No matter how large the carrying
capacity was in earlier and later months, this bottle-neck will limit the recruitment level.
The gradual decline of carrying capacity in later summer and fall can be described by
the natural succession of the ecosystem. As the menhaden grow larger in the fall, each
individual fish will consume more food and require more space than the fish in early season.
Therefore, the system would hold fewer individuals later in the season. As the season
progressed, many fish would be eaten by large predators such as striped bass, bluefish and
weakfish, or would migrate out of the Bay. The decline of capacity can be quantified by an
exponential function: jVt = N0 e"8t (Figure 10 bottom).
1
"o
5
o
c
1
o
e
1
5
o
I
1
5
s
s
1
I
°"
100
80
60
40
20
0
0
100
80
60
40
20
0
0
100
60
60
40
20
0
June 1
-
_
... A
\ . . . y s- ....
0 05 10 15 20 25 30 3.5 40 45 5.
July 1
-
r" '
-
rl /
I , J, V
0 05 10 15 20 25 30 35 40 45 5
October 1
_
- . ' '
\ A
L . / \
00 0.5 10 1.5
20 25
30
Carrying Capacity (gfm*J
Figure 10
35 40 4.5 50
114
Ecosystem Models of Chesapeake Bay: 1994-1996
-------
growth rate >0.01 growth rate >0.005
Currently proposed nutrient reductions in Chesapeake Bay may not decrease the
recruitment level of the entire Atlantic menhaden population and might increase the
recruitment level. If overloaded nutrients (mostly, nitrogen and phosphate) produced excess
primary production that sank to the bottom in spring, two negative effects would occur in
summer. First, excess primary production would deplete the limiting nutrients (e.g. silicon)
and cause the sharp decline in primary production in summer. Second, the decaying of the
over-produced primary production on the bottom could deplete the water column oxygen
concentration in summer, thus, degrading the habitat quality and decreasing the quantity of
habitat available to fish. Therefore, reducing nutrient loadings in spring may increase the
primary production and dissolved oxygen level in Chesapeake Bay in summer. This would
reduce the bottle-neck effect on the menhaden recruitment level and might increase the
recruitment level.
The modeling approach described above is still at its frontier development stages with
many simplified assumptions and no time dynamics. Information from this type modeling and
analysis of ecosystem processes may enable managers to direct their efforts toward an
ecosystem-based management program, considering model predictions of the effects of habitat
quality changes and resource management options on different trophic levels in the ecosystem.
In the future, the bioenergetics model and the spatially explicit model need to dynamically
linked with the Water Quality Model to
better relate nutrient levels to habitat
volumes and prey production.
Bay Anchovy Model
A major revision on the STELLA
bay anchovy model (Anchoa mitchilli) was
made (eg. Luo and Brandt 1993). Most
progress was made at the population level of
the model (immigration and emigration) and
simulation time. In the previous version of
the model, the model could only simulate a
one-year period, because only the growth of
individual growhorts (fish hatched within a
time-period; 15 days in this case) was
modeled for a year. Also, recruitment
processes were not included. In the new
version, the growth of each growhort was
followed for up to two years (Figure 2).
Spawning, immigration and emigration were
added to the model. The model can be now
run for a number of consecutive years. The
recruitment for each growhort was assumed
to occur at the same time each year. At
73 0.060
§ 0.050
I 0.040
| 0.030
1 0.020
°> 0.010
g 0.000
2 -0.010
Uan 1Mar
1May Uul ISep
Time (day/monlh)
Figure 11
1Nov Uan
Fish Bioenergetics
115
-------
present, the level of recruitment was set at a constant value. In the next step, the recruitment
level will be modeled as a biomass density-dependent function.
Spatial Planktivore Forage Model
Planktivorous fish can be very effective at altering the size distribution and abundances
of zooplankton in aquatic systems (Brooks and Dodson 1965, Hutchinson 1971, Sprules 1972,
Lynch 1979, O'Brien 1979, Vanni 1986, MacDonald et al 1990). Changes in size structure
of the zooplankton community can permeate throughout the food web and change production
dynamics and rates of nutrient cycling (Stein et al. 1988). Predicting the size-selectivity of a
planktivore becomes critical to evaluating the effect of planktivory on the zooplankton
community as well as the relationship of zooplankton resources to planktivore production.
Many studies have sought to determine prey size selection or preference by comparing
the size composition of zooplankton in the diet to the size composition of ambient prey
measured in the field (Baird and Hopkins 1981, Scott and Murdoch 1983, Mikheev 1984,
Magnhagen 1985, Main 1985, Collie 1986, Grover and Olla 1986, Khadka and Ramakrishna
1986, Schmitt 1986, Confer et al 1990, Forrester et al. 1994). Other studies have suggested
that apparent size and encounter frequency were the determinant of prey size selection (Werner
and Hall 1974; Confer and Blades 1975; O'Brien et al. 1976, 1985; Eggers 1977, 1982;
Pastorok 1981). A typical index of feeding electivity was usually calculated in the form of an
"odds ratio" or "forage ratio" (Fleiss 1973, Jacobs 1974, Chesson 1978, Gabriel 1979, Johnson
et al. 1990) that compares the percent of the size class in the plankton to that in the diet.
Such indices, however, are not really an electivity index from the fish's perspective; they are
more appropriately called size-difference indices. The actual prey size distribution perceived
by a fish may be very different from that measured.
In this study, a spatial simulation model was developed from the fish's perspective of
the size distribution of prey. Field data from the planktivorous bay anchovy, (Anchoa mitchilli)
was used to test the model. The bay anchovy is the most abundant planktivorous visual feeder
in Chesapeake Bay and the nearby coastal areas ( Johnson et al. 1990, Klebasko 1991, Luo
1991). Previous studies showed that patches of high consumption by the bay anchovy could
occur even though the average consumption was relatively low (Klebasko 1991, Luo and
Brandt 1993). First, the model performance was tested with hypothetical prey size
frequencies. Then, model-predicted size distributions of the prey were compared with the diet
of the bay anchovy.
Results and Discussion
Bay anchovy predation was simulated with a 3-D spatial model that also provided an
animated visualization during simulation runs. The input size-frequency distributions of prey
was compared with the size-frequency distributions of prey encountered, selected, and
captured under different assumptions. The size-frequencies of prey encountered did not differ
with number of time-steps (Figure 6A) and differed only slightly with prey densities (Figure
6B), light attenuation coefficient (Figure 6C), capture efficiencies (Figure 6D), predator's
speeds, and reactive distances (Figure 7). The distance distribution of prey (Figure 8) did
differ and most prey were located somewhere in the middle distance of its visual field (Figure
8 A) depending on types of prey size distribution (Figure 8D).
116 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Results showed that reactive distance and predator speed affected the size-frequency of
prey encountered (Figure 7). An increase of reactive distance and predator speed will increase
the encounter rate of prey (Gerritsen and Strickler 1977). However, it has not been
demonstrated how changes in reactive distance and predator speed will affect the size
frequency of prey encountered. If the reactive distance was held as a constant and predator's
speed increased, the predator would encounter more, larger prey (Figure 7A-C). But, if
predator's speed was held as a constant and the reactive distance increased, the predator would
encounter more, smaller prey (Figure 7D). An encounter was determined by both visual
acuity (function of reactive distance) and predator swimming speed.
The size-frequencies of prey eaten were compared for assumptions of "Random choice"
(RC) and "apparent size choice" (AC) for different prey densities, reactive distances and
predator's speeds were compared (Figure 9). There was no difference between the two
choices at different prey densities, if the reactive distance and predator's speed were assumed
to be equal (Figure 9A-B). However, there were significant differences between the two
choices when the predator's speed was below the reactive distance (Figure 9C-D). The prey
distance distribution between the two assumptions of prey choices was also compared (Figure
8). At low prey density (N=250, Figure 8A), there was no obvious differences of prey
distance distribution between the two assumptions. At higher prey densities (Figure 8B, C), if
the predator selected prey according to the "AC" assumption, more prey would be selected at a
closer distance than if the predator selected prey by random choice. At low prey density, the
predator rarely saw two prey at the same time therefore no choice was needed; the predator
attacked what it saw. In contrast, at higher prey density, the predator might see two or more
prey at the same time and must choose; the predator attacked the preferred prey. From an
energenetics point of view, a predator should use the "AC" model to select prey. Selecting
prey at closer distance not only saved energy in reaching the prey but also could increase the
capture efficiency. This phenomenon may help explain why larval fish need higher prey
densities to survive (Lasker 1975, 1978; Houde 1978) even though the expected encounter
rates are high enough for survival at lower prey densities.
Prey size selection by the bay anchovy was evaluated by running the model with the
zooplankton data collected in mid-Chesapeake Bay during April, May, August, and October
1990 with different combinations of swimming speed, reactive distance, and capture
efficiency. Then, the predicted size distribution was compared with the dietary size
distribution of the bay anchovy collected at the same site and time. The Kolmogorov-Smirnov
goodness of fit test (Table 1) was used. The dietary prey size distributions were from bay
anchovies ranging in size from 40 to 60 mm total length (Figure 10). The model predicted
diet size frequencies quite well for most of the observed patterns of prey size selection by the
bay anchovy even though errors could have occurred in zooplankton sampling (missing small
size by net, missing large size by pump), the patchiness of zooplankton, or in measuring
preserved zooplankton, in measuring partially digested prey from fish stomachs.
This 3-D spatial foraging model (SFM) has good potential for the study of planktivore
and plankton interactions in aquatic environment. This model will provide a link between
anchovy predation and size-dependent mortality of zooplankton. Future model developments
include changing the assumption of randomness regarding prey and predator movements,
allowing predation to alter the prey density and size distribution through time, incorporating
other environmental parameters (e.g. temperature which may affect swimming speeds) and
Fish Bioenergetics 117
-------
running the model in a spatially heterogeneous environment (Brandt and Kirsch 1993, Mason
and Patrick 1993) .
Table 1. Kolmogorov-Smimov goodness of fit test for comparing model predicted prey
size frequencies with prey size frequencies found in bay anchovy stomachs in April, May,
August, and October, 1990 in mid-Chesapeake Bay. Ed is the reactive distance (mm) to a 1-
mm prey; y is bay anchovy speed (mm/s); CE is the capture efficiency; x is prey size (mm);
dmax is the test statistics (Zar 1984); values in parenthesis are P values. A P value less than
0.05 indicates that predicted frequency is statistically different from the observed frequency.
Rd
100
100
100
100
200
200
200
200
₯
100
100
50
50
100
100
150
150
CJE
1
l-0.5x
1
l-0.5x
1
l-0.5x
1
l-0.5x
April May
dmax rimax
8.66(0.21) 11.98(0.051)
13.51(0.015) 7.07 (0.40)
18.52 (0.001) 4.24 (0.70)
24.15(0.001) 8.67(0.25)
9.88(0.11) 8.70(0.25)
14.71 (0.004) 3.97 (0.75)
11.68(0.06) 8.50(0.26)
18.02 (0.001) 4.45 (0.70)
August October
dmax dmax_
20.05(0.001) 11.02(0.075)
15.72 (0.004) 7.24 (0.40)
10.00(0.10) 7.03(0.40)
7.68(0.35) 14.17(0.01)
16.06(0.003) 11.10(0.075)
11.51 (0.065)7.04(0.40)
16.34(0.003)9.34(0.16)
10.31 (0.08) 7.93 (0.30)
Incorporation of Fish Bioenergetics Model into the Ecosystem Process Models
The Atlantic menhaden and bay anchovy models were linked with the ecosystem
process models developed at the University of Maryland (see SAV, plankton-benthos and
zooplankton models in this document) . Because the ecosystem process models used different
mass balance units (carbon) than the fish bioenergetics model (calories), unit conversions had
to be made before the models could be coupled. For example, to convert zooplankton to fish,
the biomass of zooplankton in mg carbon/m3 was first converted into gram wet weight /m3 by
40% carbon from dry weight and 15% dry weight from wet weight. Then, gram wet weight
/m3 was converted into calories/m3 by the caloric density of zooplankton (average for
copepoda was 650 cal/g). In reverse, wet fish weight was converted to dry weight at 20%,
then the dry fish weight was converted into carbon at 40 %. These conversion coefficients
were set as constant for convenience, but in reality, they varied both seasonally and
ontogenetically.
In this combined planktivore-ecosystem model, bay anchovy production was linked
with copepod production by relating copepod abundance to the proportion of maximum
consumption that a bay anchovy could achieve in a day. Menhaden production was linked
with phytoplankton production by relating phytoplankton abundance to the proportion of
maximum consumption that a menhaden can obtain in a day. Fish grazing had negative effects
118 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
on both phytoplankton and zooplankton production by taking a proportion of daily production,
and positive effects by returning the egestion and excretion into the particular organic matter
(POM) and dissolved organic matter (DOM) pools.
In preliminary simulations, the fish-ecosystem model ran for a four-year time period.
The results showed that all variables followed consistent annual cycles, indicating a steady
state of the model (Fig. 3, Fig. 4, Fig. 5). The dissolved inorganic nitrogen concentration
(DIN) peaked in early spring, followed by diatoms (DIA), which peaked in mid spring. Other
phytoplankton (OP) then peaked in late spring and late summer (Fig. 3). The zooplankton
biomass density peaked in spring and fall (Fig. 4), which matched quite well with the observed
data in 1984 and 1985. The bay anchovy population had two peaks (Fig. 2), a small one in
spring as a result of immigration; and a larger one in fall as a result of production. The
Atlantic menhaden population biomass sharply increased in summer and peaked in fall as a
result of rapid growth of individual fish. It then quickly decreased in late fall and winter as a
result of emigration (Fig. 1). At the bottom of the water column, the deposit feeder biomass
peaked in summer and the benthic filter feeder peaked in fall (Fig. 5).
In the next steps of the project, the piscivore bioenergetics models for striped bass and
bluefish will be coded into the STELLA models. The planktivore-ecosystem model will also
be refined. The piscivore models will be incorporated into the planktivore-ecosystem model to
form the piscivore-planktivore-ecosystem model. When this model is completed,
managegement simulations will be able to examine such features as bottom-up controlling
process by reducing the nutrients input and top-down controlling processes by increasing the
piscivore abundance (or reducing the fishing mortality).
Striped Bass Model
Effect ofHypoxia on Striped Bass Growth and Consumption
Laboratory experiments were conducted to test the effect of hypoxia on striped bass
(Morone saxatilis) growth and consumption (as part of a MD Sea Grant-sponsored project).
Low (dissolved oxygen) DO can be a limiting factor in many aquatic ecosystems (Coutant
1985; 1990). In particular, fish distributions, growth, consumption, and metabolic rates can
be significantly affected by the amount of DO available. In systems with a long history of
cultural eutrophication, such as Chesapeake Bay, low DO can limit striped bass production
during the summer when stratification occurs. Striped bass consumption and growth were
measured over a two-week time period, at three oxygen levels and four temperatures. Both
consumption and growth decreased significantly as DO levels decreased. Some interactions
between temperature and oxygen concentrations were present. Laboratory results were used to
develop equations to modify existing bioenergetics models of striped bass in Chesapeake Bay.
The bioenergetics models were then used to examine trade-offs between temperature, prey
availability, and oxygen levels across Chesapeake Bay. These results will be applied in the
spatially explicit model to test how the intensity of hypoxia in Chesapeake Bay will affect the
striped bass growth and habitat volume.
Fish Bioenergetics 119
-------
Field Verification of a Bioenergetics Model for Striped Bass
Verification of bioenergetics modeling with independent field estimates of consumption
or growth produced mixed results, often differing by 50-200% between field and model
estimates. One outcome of these mixed verification results was a recommendation that model
results for untested models be accepted with caution pending field validations. In earlier work,
a bioenergetics model for striped bass in Chesapeake Bay was developed and laboratory-
validated. Independent laboratory experiments indicated the striped bass bioenergetics model
provided good estimates of growth and consumption. During October 1993, age-0 striped bass
were collected with bottom trawls every 2.5 hours over two 24-h periods from near Pooles
Island in Chesapeake Bay. Fish were frozen in water for later analysis. In the laboratory,
stomach content weights were estimated using standard methods. Daily ration estimates were
calculated with a gastric evacuation model for age-0 striped bass derived from laboratory
experiments. Field measures of consumption and growth were similar (within 20%) to those
estimated with the bioenergetics model. The results of this field verification suggested that the
Hartman and Brandt (1995) age-0 striped bass model was a robust model that provided realistic
estimates of growth and consumption in field applications.
The successful testing of this bioenergetics model for striped bass with field estimated
consumption and growth will allow this model to be used with confidence by managers and
researchers. This model can be used by managers and researchers to predict the feeding
requirements of striped bass with commonly obtained information such as the growth of these
fish. The bioenergetics models will allow managers to determine if sufficient prey fishes are
present to support stocking of striped bass into reservoirs. Such information will also allow a
better understanding of the importance of striped bass in aquatic food webs.
Changes in Habitat Suitability for Bluefish and Striped Bass in Mid-Chesapeake Bay
Fish growth rate potential was previously defined as a measure of the habitat suitability
for fish (Brandt et al 1992; Brandtand Kirsch 1993; Goyke and Brandt 1993). Growth rate
potential is the expected growth rate of a predator if placed in a particular volume of water
with known physical and biological characteristics. Fish growth rate potential combines
foraging and bioenergetics models with spatially explicit field measures of temperature, prey
density and prey size to arrive at the potential growth rate a predator may experience if
constrained within the spatial cell. Thus, growth rate potential (GRP) provides a measure of
habitat quality from the fish's perspective. However, previous work on spatial modeling of
fish GRP examined only single sizes (adults) of predators (Brandt 1992; Brandt and Kirsch
1993; Goyke and Brandt 1993; Mason et al. 1995).
As fish age, the type and size of prey typically changes (Mathur and Robbins 1971;
Knight et al. 1984; Setzler-Hamilton and Hall 1991; Hartman and Brandt 1995a). Many fish
species begin life as planktivores and switch to a piscivorous feeding mode with aging.
Maximum prey size and size ranges of prey consumed also usually increase with age in
piscivores (Knight et al. 1984; Hartman and Brandt 1995a). Due to these changes in prey
sizes with age, the available prey field may differ among different sized predators of the same
species. Thus, depending upon the size distributions of prey, food availability may differ
dramatically among different sized predators.
120 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Striped bass, Morone saxatilis, and bluefish, Pomatomus saltatrix, represent an
interesting ecological contrast. Both also support valued sport and commercial fisheries
(Setzler et al. 1980; Mid-Atlantic Fishery Management Council 1989). Striped bass and
bluefish are the dominant predators in Chesapeake Bay and represent two families
(Percichthyidae and Pomatomidae, respectively) from the order Perciformes). Striped bass are
resident in Chesapeake Bay until about age-6 (Kohlenstein 1981); whereas, bluefish are
seasonal residents, using the bay as a nursery and growing area during summer (Norcross et
al. 1974; Kendall and Walford 1979; McBride and Conover 1992). Feeding strategies also
differ among the two species bluefish are cruising predators with teeth that permit searing
detention (Baird 1873; Bigelow and Schroeder 1953); striped bass feed both as sit-and-wait
and cruising predators (Setzler et al. 1980). Both predators commonly feed on the same
species such as bay anchovy (Anchoa mitchitti), Atlantic menhaden (Brevoortia tyrannus) and
spot (Leiostomus xanthurus)(Grant 1962; Homer and Boynton 1978; Friedland et al. 1988;
Hartman and Brandt 1995a,c).
Growth rates of fish are mediated by body size, temperature, and prey availability.
Because water temperatures and available prey sizes and abundance (prey spectrum) are
heterogeneously distributed, the conditions for growth of a predator will also be patchily
distributed. Similarly, temperatures and prey spectrum changes seasonally, so conditions for
growth of a predator may differ among seasons. Due to the complex interaction of spatial,
seasonal, biotic, abiotic, and ontogenetic changes in growth conditions for fish, an adequate
evaluation of habitat quality for any fish species would require a seasonal and ontogenetic
analysis of fish growth rate potential. In this study, spatially explicit modeling was used to
combine previously developed bioenergetics models of bluefish and striped bass (Hartman and
Brandt 1995b) with quantitative acoustic measurements of prey spectra to evaluate habitat
suitability for the two piscivores in the mid-Chesapeake Bay.
Results and discussion
The perception of habitat quality, or growth rate potential, changes with fish species,
fish size, and season. Within a season and species, spatial maps of fish growth rate potential
substantially differed across fish sizes (Fig. 6). This was largely due to different prey size
ranges that resulted in different prey spectrums and, in part, to differing thermo-physiological
constraints on consumption and metabolic rates (Hartman and Brandt 1995b). Small (10 g)
predators usually had the lowest percentages of cells that supported positive growth, probably
due to the low availability of the smaller prey sizes these predators could consume. The large
predators consumed a wider range of prey sizes which led to a wider prey field. Thus, larger
predators tended to have higher percentages of cells that supported positive growth rate
potential (Fig. 7, 8). The exception to this generalization was during winter, when 10 g striped
bass had slightly higher possible growth rates and a greater percentage of cells with positive
growth than 2000 g striped bass.
Water column averages of growth rate potential suggested that fall and spring provided
the best, and summer and winter the worst overall habitat quality for striped bass and bluefish
(Table 2). For each predator, highest averages of water column GRP occurred in fall. Striped
bass average GRP was +0.0039 (10 g) to +0.0209 g g1 d~' (100 g) in fall. Bluefish average
GRP was -0.0099 g g1 d"1 for 10 g, but increased to +0.0187 g g'1 d'1 for 2000 g fish (Table
2). Poorest average growth was in summer for 10 g fish (-0.0146 to -0.0327 g g"1 d"1), but all
Fish Bioenergetics 121
-------
predators had poorest average growth of any season in summer. In winter and spring, 2000 g
predators had higher average GRP than smaller conspecifics, with average spring GRPs of
+0.0034 to +0.0096 g g"1 d"1 estimated for bluefish and striped bass, respectively (Table 2).
Table 2. Seasonal whole water-column averages of growth rate potential (x 100) for 10,
100, and 2000 g striped bass and bluefish in Chesapeake Bay. Standard deviations of mean
values (x 100) are given parenthetically.
Predator Winter Spring Summer Fall
Striped Bass
10 g
100 g
2000 g
Bluefish
10 g
100 g
2000 g
-0.08
(0.33)
-0.09
(0.25)
-0.06
(0.12)
-0.70
(0.14)
-0.38
(0.07)
-0.17
(0.03)
-0.03
(0.65)
+0.05
(0.60)
+0.96
(1.11)
-0.94
(0.54)
-0.48
(0.34)
+0.34
(0.52)
-1.46
(0.67)
-1.06
(0.47)
-0.26
(0.86)
-3.27
(0.73)
-1.69
(0.59)
-0.29
(1.20)
+0.39
(1.84)
+2.09
(2.28)
+ 1.54
(1.12)
-0.99
(1.79)
+0.25
(1.76)
+ 1.87
(1.23)
Measures of habitat quality (spatial growth potential) supported the concept that fish
use of the Bay may be driven by both thermal physiology and prey availability. Striped bass
were much better suited to the thermal regimes and prey fields available during winter and
spring than were bluefish. Growth potential for both species during winter was limited by low
water temperatures. Bluefish of any size would lose weight in all cells during winter;
whereas, some small growth was possible for striped bass. Growth rate was also suppressed
for bluefish during spring, with a small volume of water supporting growth for 10 or 100 g
bluefish (Fig. 7). Positive GRP was available in 25 % of cells for 10 and 100 g striped bass
and 2000 g fish had positive GRP in over 65 % of all cells (Fig. 8).
Given the differing growth maps for striped bass and bluefish during winter and spring,
it was not surprising that bluefish are not found in Chesapeake Bay during those seasons (when
spring data were collected; Norcross et al. 1974; Kendall and Watford 1979) and that striped
bass are annual residents of the Chesapeake Bay (Kohlenstein 1981). Bluefish overwinter in
the coastal waters off North Carolina, from November to May, where temperatures, and in all
122 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
likelihood prey availability, exceed those found in Chesapeake Bay during that time (Norcross
et al. 1974; Kendall and Watford 1979).
Spatial maps of growth rate potential of small (10 g) fish may be less representative of
growth conditions for these predators (than for larger predators) in Chesapeake Bay. First, we
have considered only pelagic fish as possible prey for small striped bass and bluefish. Actual
diets of these small predators often include a large variety of abundant prey not accessible to
our hydroacoustic instruments such as benthic invertebrates, shrimps, and some benthic fishes
(Boynton et al. 1981; Friedland et al. 1988; Hartman and Brandt 1995a). Second, the
frequency of the acoustic transducer (120 kHz) did not effectively sample very small
organisms, thus any prey less than 12 mm were not included in the prey field. Therefore,
prey fields of the small fish in the model were probably less accurate than for larger fishes
which fed primarily on larger pelagic fishes (Grant 1962; Schaefer 1970; Texas Instruments
1976; Homer and Boynton 1978; Gardinier and Hoff 1982; Rulifson and McKenna 1987;
Hartman and Brandt 1995a).
Water temperatures and prey densities combined to limit the maximum growth rate and
ranges of growth rates for striped bass and bluefish. Within a season, water temperatures fell
within a restricted range (±. 4°C). Thus, GRP was moderated more by the absolute density of
prey than by temperatures within a season. Across seasons, maximum prey densities within a
cell were similar (although spatial patterning and seasonal prey densities differed greatly, prey
densities of over 5 g nT3 were found in all seasons) and the resulting GRP was limited by
temperature. This suggested that across seasons, temperature limited the maximum growth
rate potential of predators; but, within a season, the distribution of prey density limited the
proportion of cells providing favorable growth conditions.
Expansion of the concept of fish growth rate potential to include different sized
predators with careful consideration of each predator's size-specific prey field improved
understanding of ecosystem function and the importance of spatial patterning in the physical
and biological habitat upon available habitat for a species throughout it's ontogenetic
development. These models allowed evaluation of habitat quality in a comparative nature,
across species and sizes of fish. This modeling has shown that habitat quality and the spatial
patterning of habitat in the environment (as assessed by fish growth rate potential) can be
quantified from the fish's perspective and that it changed seasonally and with age for striped
bass and bluefish.
Weakfish
Spatial Analyses of Predator Growth Potential: Biotic and Abiotic Influences
Chesapeake Bay is a nursery and growing area for weakfish. In summer, hypoxic
mesohaline waters may affect distributions, migrations, and growth of these fish. Growth is
influenced by biotic variables (e.g. prey availability and size) and abiotic variables
(temperature, D.O., salinity). These biotic and abiotic factors vary in time and space. A
spatial model was developed that used field-derived measures of abiotic variables and
hydroacoustically derived prey size, density, and distribution arrays. These arrays were
combined with weakfish bioenergetics and foraging models to predict growth potential. Fish
growth potential represented the growth rate of a fish if constrained within a spatial cell (about
Fish Bioenergetics 123
-------
0.5 m depth by 30 m cross-sections) with prey and physical attributes inherent to that cell.
Modeling results showed DO may play a key role in summer growth patterns of weakfish,
particularly age-0 fish. Weakfish growth potential (all ages) was generally highest near the
surface where anchovy densities were greatest. Anchovies are a primary food source for
weakfish (Hartman and Brandt 1995a). A sequential analysis of the physical data showed that
DO and temperature were more important than salinity in influencing growth potential. These
results suggested that DO levels might limit weakfish production; therefore weakfish habitat
could be improved by continued nutrient control programs in Chesapeake Bay.
The spatial modeling in this study permited a better understanding of the importance of
temperature, salinity, DO, and food availability upon habitat quality for weakfish in
Chesapeake Bay. The spatial maps of fish growth rate potentials illustrated how habitat quality
compared across different sections of the water column and the importance DO levels, relative
to available temperatures if you were a weakfish! These types of color maps help to
illustrate the importance of reducing nutrients for future management of coastal waters.
White perch
Trophic Dynamics: Diet, Growth, and Bioenergetics
A bioenergetics model for the white perch (Morone americand) is currently being
developed. The white perch is an abundant benthic feeder in the freshwater and mesohaline
portions of the Chesapeake Bay estuary and may have significant effects on the macrobenthic
community of the bay. White perch are fished commercially and recreationally in many
coastal waters. They provide a prey resource to larger piscivores, such as striped bass and
bluefish. To define the role of the white perch in the Bay's food web, diet, growth, and
bioenergetics of this species were examined. White perch were sampled at four different
seasons of the year (April, June, October, and December 1995) over a 24-hour period with
bottom trawls. Fish lengths, wet weights, and ages were assessed. The main diets of white
perch consisted of amphipods (Gammarid spp.), copepods, and cladocerans. Consumption
estimates will be derived from a bioenergetics model and compared to direct measures of diel
and seasonal predation rates on various prey types. Overall, white perch appear to play a
dominant role in the benthic food web interactions in Chesapeake Bay.
Future Objectives
During the next year, the principle objective will be to emphasize the development of
biologically complex and spatially explicit models. An important task is to couple the fish
bioenergetics models with the ecosystem process models and with the water quality model. The
connections of ecosystem process and fish bioenergetics models accomplished during previous
years will be defined and used to assess the effects of patchiness and spatial resolution on
model results and predictions for water quality and habitat improvement. A range of policy
scenarios will be simulated using the individual and integrated models to provide resource
managers with a unique tool for relating fish production to proposed and planned nutrient
124 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
reductions. An assessment will be made to determine the water quality improvements needed
to achieve targeted habitat and species restoration goals and to identify the associated sensitive
biological components. The combination of ecosystem process models and fish bioenergetics
models will provide the capability to predict potential fish production from nutrient loading
and to simulate management scenarios.
Fish Bioenergetics 125
-------
References
Ahrenholz, D.W., W.R. Nelson, and S.P. Epperly. 1987. Population and fishery
characteristics of Atlantic menhaden, Brevoortia tyrannus. Fish. Bull., U.S. 85:569-
600.
Aksnes, D. L., and J. Giske. 1993. A theoretical model of aquatic visual feeding. Ecol.
Model. 67:233-250
Baird, R. C. and T. L. Hopkins. 1981. Trophodynamics of the fish Valenciennellus
tripunctulatus. 2. Selectivity, grazing rates and resource utilization. Mar. Ecol. Prog.
Ser. 5:11-19.
Baird, R.M. and R.E. Ulanowicz. 1989. The seasonal dynamics of the Chesapeake Bay
ecosystem. Ecol. Monogr. 59:329-364.
Baird, S. F. (1873). Natural history of some of the more important food fishes of the south
shore of New England. EL The bluefish, (Pomatomus saltatrix) (Linn.). Rep. U.S.
Fish. Comm. for 1871-1872: 235-252.
Bartell, S.M., I.E. Breck, R.H. Gardner and A.L. Brenkert. 1986. Individual parameter
perturbation and error analysis of fish bioenergetics models. Can. J. Fish. Aq. Sci.
43:160-168.
Beamish, F. W. H. 1978. Swimming capacity. Pages 101-189 in W. S. Hoar, and D. J.
Randall, editors, Fish Physiology. 7. Locomotion. Academic Press, New York, USA.
Bigelow, H. B., and W.C. Schroeder. 1953. Fishes of the Gulf of Maine. U.S. Fish Wild.
Ser., Fish. Bull. 53.
Blaxter, J. H. S., and R. S. Batty. 1985. Herring behaviour in the dark: responses to
stationary and continuously vibrating obstacles. J. Mar. Biol. Assoc. U.K. 66:1031-
1049.
Boynton, W. R., Polgar, T. T., and H.H. Zion. 1981. Importance of juvenile striped bass
food habits in the Potomac Estuary. Trans. Am. Fish. Soc. 110:56-63.
Bonzek, C.F., P.J. Geer, J.A. Colvocoresses and R.E. Harris. 1992. Juvenile finfish and
blue crab stock assessment program. Virginia Institute of Marine Science. Special
science report No. 124. Vol. 1992.
Brandt, S. B. (1992). Acoustic quantification of fish abundance in the Chesapeake Bay.
Final Report to Maryland Department of Natural Resources, Power Plant Topical
Research Program, [UMCEES]CBL 92-054, Annapolis.
126 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Brandt, S.B., D.M. Mason and E.V. Patrick. 1992. Spatially explicit models offish growth
rate. Estuaries 17(2): 23-3 5.
Brandt, S.B. and K.J. Hartman. 1993. Innovative approaches using bioenergetics models: Future
applications to fish ecology and management. Trans. Amer. Fish. Soc. 122:731-735.
Brandt, S.B. and J. Kirsch. 1993. Spatially explicit models of striped bass growth in the
mid-Chesapeake Bay. Trans. Amer. Fish. Soc. 122:845-869.
Brooks, J. L. and S. I. Dodson. 1965. Predation, body size, and composition plankton.
Science 150:28-35.
Chesson, J. 1978. Measuring preference in selective predation. Ecology. 59:211-215.
Collie, J. S. 1987. Food selection by yellowtail flounder (Limandaferrugined) on Georges
Bank. Can. J. Fish. Aquat. Sci. 44:357-367.
Confer, J. L. and P. I. Blades. 1975. Omnivorous zooplankton and planktivorous fish.
Limnol. Oceanogr. 20:571-579.
Confer, J. L, E. L. Mills, and L. O'Bryan. 1990. Influence of prey abundance on species and
size selection by young yellow perch (Percaflavescens). Can. J. Fish. Aquat. Sci.
47:882-887.
Coutant, C.C. 1990. Temperature-oxygen habitat for freshwater and coastal striped bass in a
changing climate. Trans. Am. Fish. Soc. 119:240-253.
Coutant, C.C. 1985. Striped bass, temperature, and dissolved oxygen: a speculative
hypothesis for environmental risk. Trans. Am. Fish. Soc. 114:31-61.
Day, J. W., Jr., Hall, C. A. S., Kemp, W. M., Yanez-Araniciba, A. (1989). Estuarine
Ecology. John Wiley and Sons, Inc. New York.
Eggers, D. M. 1982. Planktivore preference by prey size. Ecology. 63:381-390.
Eggers, D. M. 1977. The nature of prey selection by planktivorous fish. Ecology. 58:46-
59.
Fleiss, J. L. 1973. Statistical methods for rates and proportions. John Wiley and Sons, New
York.
Forrester, G. E., J. G. Chace and W. McCarthy. 1994. Diel an density-related changes in
food consumption and prey selection by brook charr in a New Hampshire steam. Env.
Bio. Fish. 39:301-311.
Fish Bioenergetics 127
-------
Friedland, K. D., Garman, G. C., Bedja, A. J., Studholme, A. F., OUa, B. (1988).
Interannual variation in diet and condition in juvenile bluefish during estuarine
residency. Trans. Am. Fish. Soc. 117: 474-479.
Gabriel, W. L. 1979. Statistics of selection. In: Lipovsky, S. J., Simenstad, C. A. (eds).
Gutshop 1978: fish food habits studies. Proc. 2nd Pacific Northwest technical
workshop. Washington Sea Grant, Seattle, p. 62-66.
Gardinier, M. N., Hoff, T. B. 1982. Diet of striped bass in the Hudson River Estuary.
N. Y. Game Fish J. 29: 152-165.
Gerritsen, J. and J. R. Strickler. 1977. Encounter probabilities and community structure in
zooplankton: a mathematical model. J. Fish. Res. Board Can. 34:73-82.
Goyke, A. P., and S.B. Brandt. 1993. Spatial models of salmonine growth rates in Lake
Ontario. Trans. Am. Fish. Soc. 122: 870-883.
Grant, G. C. 1962. Predation of bluefish on young Atlantic menhaden in Indian River,
Delaware. Ches. Sci. 3:45-47.
Grover, J. J. and B. L. Olla. 1986. Morphological evidence for starvation and prey size
selection of see-caught larval sablefish, Anoplopomafimbria. Fish. Bull. 84:484-489.
Harding, L W. Jr. 1994. Long-term trends in the distribution of phytoplankton in Chesapeake
Bay: roles of light, nutrients and streamflow. Mar. Ecol. Prog. Ser. 104:267-291.
Hartman, K. J. 1994. Striped bass, bluefish, and weakfish in the Chesapeake Bay: energetics,
trophic linkages, and bioenergetics model applications. PhD thesis. The University of
Maryland, College Park, MD, USA. 336 pp.
Hartman, K. J., and S.B. Brandt. 1995a. Trophic resource partitioning, diets and growth of
sympatric estuarine predators. Trans. Am. Fish. Soc. 124: 520-537.
Hartman, K. J., and S.B. Brandt. 1995b. Comparative energetics and the development of
bioenergetics models for sympatric estuarine piscivores. Can J. Fish. Aq. Sci. 52:
1647-1666.
Hartman, K. J., and S.B. Brandt. 1995c. Predatory demand and impact of striped bass,
bluefish, and weakfish in the Chesapeake Bay: applications of bioenergetics models.
Can J. Fish. Aq. Sci. 52: 1667-1687.
Hartman, KJ. and S.B. Brandt. 1993. A bioenergetics model for age 0 striped bass:
evaluation of parameter accuracy and precision. Trans. Amer. Fish. Soc. 122:912-926.
128 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Hofmann, E. E. and J. W. Ambler. 1988. Plankton dynamics on the outer southeastern U.S.
continental shelf. Part n: A time-dependent biological model. J. Mar. Res. 46:883-
917
Homer, M., and W.R. Boynton. 1978. Stomach analysis of fish collected in the Calvert
Cliffs region, Chesapeake Bay - 1977. Final Rep. to MDNR, PPS Progress,
Annapolis, MD. Ref. No. UMCEES 78-154.
Houde, E. D. 1978. Critical food concentrations for larvae of three species of subtropical
marine fishes. Bull. Mar. Sci. 28:395-411.
Hutchinson, B. P. 1971. The effect of fish predation on the zooplankton of ten Adirondack
lakes, with particular reference to the alewife, Alosapseudoharengus. Trans. Am.
Fish. Soc. 100:325-335.
Jacobs, J. 1974. Quantitative measurement of food selection: a modification of the forage
ratio and Ivlev's electivity index. Oecologia (Berl.) 14:413-417.
Jobling, M. 1982. Food and growth relations of the cod, Gadus morhua L., with special
reference to Balsfjorden, north Norway. J. Fish Biol. 21: 357-363.
Johnson, W. S., D. M. Allen, M. V. Ogburn and S. E. Stancyk. 1990. Short-term predation
responses of adult bay anchovies Anchoa mitchilli to estuarine zooplankton availability.
Mar. Ecol. Prog. Ser. 64:55-68.
Kareiva, P. and M. Anderson. 1988. Spatial aspects of species interactions: The wedding of
models and experiments, pp. 38-54, In: Hastings, A. (ed.), Community Ecology,
Springer-Verlag, New York.
Kemp, W.M., S.B. Brandt, W.R. Boynton, R. Bartleson, J. Hagy, J. Luo and C.J. Madden.
1992. Ecosystem models of the Patuxent River Estuary relating nutrient loading to
production of selected fish populations. Final Rept, Year 1. Maryland Dept. Nat.
Res. Tidewater Admin., Chesapeake Bay Research and Monitoring Division, Annapolis
MD 21401.
Kendall, A. W., and L.W. Walford. (1979). Sources and distribution of bluefish Pomatomus
saltatrix, larvae and juveniles off the east coast of the United States. Fish. Bull. 77:
213-227.
Khadka, R. B. and R. T. Ramakrishna. 1986. Prey size selection by common carp (Cyprinus
carpio var. communis) larvae in relation to age and prey density. Aquaculture. 54:89-
96.
Fish Bioenergetics 129
-------
Kitchell, J.F., DJ. Stewart and D. Weininger. 1977. Applications of a bioenergetics model
to yellow perch (Perca flavescens) and walleye (Stizostedion vitreum). J. Fish. Res.
Bd. Can. 34:1922-1935.
Kitchell, J. F., DJ. Stewart, and D. Weininger 1977. Applications of a bioenergetics
model to yellow perch (Perca flavescens) and Walleye (Stizostedion vitreum vitreum).
J. Fish. Res. Board. Can. 34:1922-1935.
Klebasko, M. J. 1991. Feeding Ecology an daily ration of bay anchovy (Anchoa mitchilli) in
the mid-Chesapeake Bay. MS. thesis, Univ. of Maryland, College Park. 103 p.
Knight, R. L., F.J. Margraf, and R.F. Carline. 1984. Piscivory by walleyes and yellow
perch in western Lake Erie. Trans. Am. Fish. Soc. 113:677-693.
Kohlenstein, L. C. (1981). On the proportion of the Chesapeake Bay stock of striped bass
that migrates into the coastal fishery. Trans. Am. Fish. Soc. 110:168-179.
Lasker, R. 1975. Field criteria for survival of anchovy larvae: the relation between inshore
chlorophyll maximum layers and successful first feeding. U. S. Fish. Bull. 73:453-
462.
Lasker, R. 1978. The relation between oceanographic conditions and larval anchovy food in
the California Current: Identification of factors contributing to recruitment failure.
Rapports et Proces-Verbaux des Reunions du Conseil International pour 1'Exploration
de la Mer 173:212-230.
Lippson, A. J. 1985. The Chesapeake Bay in Maryland - an atlas of natural resources.
Johns Hopkins University Press, Baltimore, Maryland.
Luo, J. and S.B. Brandt. 1993. Bay anchovy, Anchoa mitchilli, production and consumption
in mid-Chesapeake Bay based on a bioenergetics model and acoustical measures of fish
abundance. Mar. Ecol. Prog. Ser. 98:223-236.
Luo, J. 1991. The life history of the bay anchovy, Anchoa mitchilli, in Chesapeake Bay. Ph. D.
dissertation, College of William and Mary, Williamsburg, VA. 108 p.
Lynch, M. 1979. Predation, competition, an zooplankton community structure: an experimental
study. Limnol. Oceanogr. 24:253-274.
MacDonald, M. E., L. B. Crowder and S. B. Brandt. 1990. Changes in Mysis and Pontoporeia
populations in southeastern Lake Michigan: A response to shifts in fish community.
Magnhagen, C. 1985. Random prey capture or active choice? An experimental study on prey
size selection in three marine fish species. Oikos. 45:206-216.
130 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Main, K. L. 1985. The influence of prey identity and size on selection of prey by two marine
fishes. J. Exp. Mar. Biol. Ecol. 88:145-152.
Mason, D. M. and E. V. Patrick. 1993. A model for the Space-time dependence of feeding
for pelagic fish populations. Trans. Am. Fish. Soc. 122:884-901.
Mason, D. M., Goyke, A., and S.B. Brandt. 1995. A spatially-explicit bioenergetics
measure of habitat quality for adult salmonines: Comparison between Lakes Michigan
and Ontario. Can. J. Fish. Aq. Sci. 52: 1572-1583.
Mathur, D., and R.W. Robbins. 1971. Food habits and feeding chronology of young white
crappie, Pomoxis annularis, Rafinesque, hi Conowingo Reservoir. Trans. Am. Fish.
Soc. 100:307-317.
McBride, R. S., andD.O. Conover. 1992. Recruitment of young-of-the-year bluefish
Pomatomus saltatrix to the New York Bight: variation in abundance and growth of
spring- and summer-spawned cohorts. Mar. Ecol. Prog. Ser. 78: 205-216.
McHugh, J. L. 1967. Estuarine nekton. In: Lauff, G.. H. (ed.), Estuaries. Am. Assoc.
Advanc. Sci., Publ. No. 83, Washington, DC, pp 581-620.
Mid Atlantic Fishery Management Council, Atlantic States Marine Fisheries Commission,
National Marine Fisheries Service, New England Fishery Management Council, South
Atlantic Fishery Management Council. 1989. Fishery management plan for the
bluefish fishery. 76 pp.
Mikheev, V. N. 1984. Prey size and food selectivity in young fishes. J. Ichthyol. 24:66-76.
Miller, T. J, L. B. Crowder, and J. A. Rice. 1993. Ontogenetic changes in behavioural and
histological measures of visual acuity in three species of fish. Environmental Biology
of Fishes 37:1-8.
Monteleone, D.M. and E.D. Houde. 1990. Influence of maternal size on survival and growth
of striped bass Morone saxatilis (Walbaum) eggs and larvae. J. Exp. Mar. Biol. Ecol.
140:1-11.
Nixon, S. W. 1981. Freshwater inputs and estuarine productivity. In R. D. Cross and D. L.
Williams (Eds), Proceedings of the National Symposium on Freshwater Inflow to
Estuaries. U.S. Fish and Wildlife Service, Office of Biological Service, FWS/OBS-
81/04, Washington, D. C., Vol. 1, pp 31-57.
Norcross, J. J., S.L. Richardson., W.H. Massmann, E.B. Joseph. 1974. Development
of young bluefish (Pomatomus saltatrix) and distribution of eggs and young in
Virginia coastal waters. Trans. Am. Fish. Soc. 103:477-497.
Fish Bioenergetics 131
-------
O'Brien, W. J. 1979. The predator-prey interaction of planktivorous fish and zooplankton.
Am. Sci. 67:572-581.
O'Brien, W. J., B. Evans and C. Luecke. 1985. Apparent size choice of zooplankton by
planktivorous sunfish: exceptions to the rule. Env. Bio. Fish. 13:225-233.
O'Brien, W. J., N. A. Slade, and G. L. Vinyard. 1976. Apparent size as the determinant of
prey selection by bluegill sunfish (Lepomis Macrochirus). Ecology. 57:1304-1310.
Officer, C.B., R.B. Briggs, J.L Taft, L.E. Cronin, M.A. Tyler, and W.R. Boynton. 1984.
Chesapeake Bay anoxia: origin, development and significance. Science 223:22-27.
Park, S. K. and K. W. Miller. 1988. Random Number Generators: Good ones are hard to
find. Communication of the ACM, 31:1192-1205.
Pastorok, R. A. 1981. Prey vulnerability and size selection by Chaoborus larvae. Ecology.
62: 1311-1424.
Peters, D. S. and W. E. Schaaf. 1981. Empirical Modelof the Trophic Basis for Fishery
Yield in Coast Waters of the Eastern USA.
Possingham, H.P. and J. Roughgarden. 1990. Spatial population dynamics of a marine
organism with a complex life cycle. Ecology 71:973-985.
Rippetoe, T.H. 1993. Production and energetics of Atlantic menhaden in Chesapeake Bay.
Master's thesis, Univ. of Maryland, College Park. 113 pp.
Rulifson, R. A., McKenna, S. A. 1987. Food of striped bass in the upper Bay of Fundy,
Canada. Trans. Am. Fish. Soc. 116:119-122.
Schaefer, R. H. 1970. Feeding habits of striped bass from the surf waters of Long Island.
N.Y. Game Fish J. 17:1-17.
Schmitt, P. D. 1986. Prey size selectivity and feeding rate of larvae of the northern anchovy,
Engraulis mordax Girard. Rep. CCOFI. 27:153-161.
Scott, M. A. and W. W. Murdoch. 1983. Selective predation by the backswimmer,
Notonecta. Limnol. Oceanogr. 28:352-366.
Setzler, E. M., plus eight others. 1980. Synopsis of biological data on striped bass,
Morone saxatilis (Walbaum). NOAA Tech. Rep., NMFS Circular 433. FAO
Synopsis No. 121. 69pp.
132 Ecosystem Models of Chesapeake Bay: 1994-1996
-------
Setzler-Hamilton, E. M., Hall, L., Jr. 1991. pp 13-1 to 13-25, In: Habitat Requirements
for Chesapeake Bay Living Resources. S. L. Funderburk, J. A. Mihursky, S. J.
Jordan, and D. Riley (eds.), Chesapeake Research Consortium, Solomons, Maryland.
Sprules, W. G. 1972. Effects of size-selective predation and food competition on high altitude
zooplankton communities. Ecology 53:375-386.
Stein, R. A., S. T. Threllkeld, C.D. Sandgren, W. G. Sprules, L. Persson, E. E. Werner,
W. E. Neill, and S. I. Dodson. 1988. Size-structured interactions in lake
communities. In: Complex interactions in lake communities, Carpenter, S .R. (ed).
161-179. Springer-Verlag New York.
Texas Instruments. 1976. Predation by bluefish in the lower Hudson River. Final Report
to Consolidated Edison Company of New York.
U.S. National Marine Fisheries Service. 1978. Fishery statistics of the United States, 1977.
U.S. National Marine Fisheries Service Current Fishery Statistics, No. 7500. 112 pp.
Vanni, M. J. 1986. Fish predation and zooplankton demography: indirect effects. Ecology
67:337-354.
Werner, E. E. and D. J. Hall. 1974. Optimal foraging and the size selection of prey by the
bluegill sunfish (Lepomis macrochirus). Ecology 55: 1042-1052.
Zastrow, C.E., and E.D. Houde, and L.G. Morin. 1991. Spawning, fecundity, hatch-date
frequency and young-of-the-year growth of bay anchovy Anchoa mitchilli in mid-
Chesapeake Bay. Mar. Ecol. Prog. Ser. 73:161-171.
Fish Bioenergetics 133
-------
-------
------- |