EPA
United States
Environmental Protection
Agency
Industrial Environmental Research
Laboratory
Research Triangle Park NC 27711
EPA-600/7 80 O68
March 1980
TVA
Tennessee Valley
Authority
Division of Energy
Demonstrations and Technology
Chattanoog^T^ 3 7401
EOT 1O1
ORNL
Oak Ridge
National Laboratory
Environmental
Sciences Division
Oak Ridge, TN 37830
A Partial Differential
Equation Model of Fish
Population Dynamics and
Its Application in
Impingement Impact
Analysis
Interagency
Energy/Environment
R&D Program Report
-------
RESEARCH REPORTING SERIES
Research reports of the Office of Research and Development, U.S. Environmental
Protection Agency, have been grouped into nine series. These nine broad cate-
gories were established to facilitate further development and application of en-
vironmental technology. Elimination of traditional grouping was consciously
planned to foster technology transfer and a maximum interface in related fields.
The nine series are:
1. Environmental Health Effects Research
2. Environmental Protection Technology
3. Ecological Research
4. Environmental Monitoring
5. Socioecbnomic Environmental Studies
6. Scientific and Technical Assessment Reports (STAR)
7. Interagency Energy-Environment Research and Development
8. "Special" Reports
9. Miscellaneous Reports
This report has been assigned to the INTERAGENCY ENERGY-ENVIRONMENT
RESEARCH AND DEVELOPMENT series. Reports in this series result from the
effort funded under the 17-agency Federal Energy/Environment Research and
Development Program. These studies relate to EPA's mission to protect the public
health and welfare from adverse effects of pollutants associated with energy sys-
tems. The goal of the Program is to assure the rapid development of domestic
energy supplies in an environmentally-compatible manner by providing the nec-
essary environmental data and control technology. Investigations include analy-
ses of the transport of energy-related pollutants and their health and ecological
effects; assessments of, and development of, control technologies for energy
systems; and integrated assessments of a wide range of energy-related environ-
mental issues.
EPA REVIEW NOTICE
This report has been reviewed by the participating Federal Agencies, and approved
for publication. Approval does not signify that the contents necessarily reflect
the views and policies of the Government, nor does mention of trade names or
commercial products constitute endorsement or recommendation for use.
This document is available to the public through the National Technical Informa-
tion Service, Springfield, Virginia 22161.
-------
EPA-600/7-80-068
TVA EDT-101
March 1980
A Partial Differential Equation
Model of Fish Population
Dynamics and Its Application
in Impingement Impact Analysis
by
P.A. Hackney and T.A. McDonough and D.L DeAngelis and M.E. Cochran
TVA, Office of National Resources ORNL, Environmental Sciences Division
Norris, Tennessee 37828 Oak Ridge, Tennessee 37830
EPA Interagency Agreement No. IAG-D8-E721-BE
Program Element No. INE624A
EPA Project Officer: Theodore G. Brna
TVA Project Director: Hollis B. Flora II.
Industrial Environmental Research Laboratory
Office of Environmental Engineering and Technology
Research Triangle Park, NC 27711
Prepared for
U.S. ENVIRONMENTAL PROTECTION AGENCY and TENNESSEE VALLEY AUTHORITY
Office of Research and Development Division of Energy Demonstrations and Technology
Washington, DC 20460 Chattanooga, TN 37401
-------
ABSTRACT
This study was undertaken to (1) develop a model describing fish
populations as a function of life process dynamics and facilities which
impose additional mortality on fish and (2) improve objective impingement
impact prediction. The mathematical model developed accounts for hatching,
growth, and mortality as functions of time and permits computer simulation
of impingement impact. It also accounts for the genetic and environmental
heterogeneity effects on the growth of a cohort of fish. Gizzard shad data
collected by TVA were used to corroborate the model.
Simulated impingement impacts for the steam-electric generating plant
and reservoir studied were much less than could be measured in field
studies. For a 10-fold increase over observed impingement losses, the
model predicted that gizzard shad stock levels would fall by less than
10 percent for any age group. Similarly, the model with a 100-fold
increase over the observed losses, predicted that age IV gizzard shad
stock levels were reduced about 65 percent from baseline values, with
younger age groups showing less response. Model simulations revealed
that current levels of intake-induced mortality reduced the total numbers
of gizzard shad in each age class by less than 1 percent. These findings
show little effect to a species having high natural mortality, but cannot
be generalized to other species having significantly different natural
mortality patterns.
This report was submitted in fulfillment of Task 4 Subagreement 21
of the interagency agreement between TVA and EPA (TV-41967A, EPA-IAG-D5-
0721) under the sponsorship of the U.S. Environmental Protection Agency.
This report covers the period October 1, 1978, to December 17, 1979, with
work completed as of March 15, 1980.
111
-------
DISCLAIMER
This report was prepared by the Tennessee Valley Authority and has
been reviewed by the Office of Environmental Engineering and Technology,
U.S. Environmental Protection Agency and approved for publication.
Approval does not signify that the contents necessarily reflect the views
and policies of the Tennessee Valley Authority or the United States
Environmental Protection Agency nor does mention of trade names or
commercial products constitute endorsement or recommendation for use.
IV
-------
CONTENTS
Abstract iii
Disclaimer iv
Figures vi
Tables viii
Abbreviations and symbols x
1. Introduction 1
2. Conclusions 3
3. Recommendations 4
4. Materials and methods 5
The model 5
Data source 16
5. Experimental procedures 42
Baseline conditions 42
Zero plant mortality 42
Ten-fold mortality increase 43
One-hundred fold mortality increase 43
6. Results and discussion 44
Bibliography 51
Appendices
A. Analysis of model 53
B. Use of computer program 75
1. General information 76
2. Description of input data 77
3. Example application of the program 81
C. Computer program 88
v
-------
FIGURES
Number Page
1 The hatching of a cohort of fish through time
(the reproduction season) 6
2 A cohort of fish at a given instant in time 7
3 The same cohort of fish in Figure 2 shown at four
instances in time 9
4 Another representation of the cohort of fish in Figure 3
to better show growth 10
5 Mortality in a cohort of fish 11
6 Relationship of hatching, growth, and mortality in a
cohort of fish 12
7 Size distributions at three successive time intervals ... 14
8 Location of Cumberland Steam-Electric Plant on
Barkley Reservoir 18
9 Intake and associated structures, Cumberland Steam-
Electric Plant 19
10 Distribution of gizzard shad in Barkley Reservoir,
1974-1978 21
11 Total length - scale radius relationship of gizzard
shad from Barkley Reservoir, 1974-1978 22
12 Comparison of back calculated growth rates of gizzard
shad for the 1969 through 1976 year classes 24
13 Estimated length frequency distributions of gizzard
shad in 1977 and 1978 cove rotenone samples 25
14 Mean length and mean length increment of gizzard
shad in 1977 and 1978 cove rotenone samples 28
15 Yearly length increment of gizzard shad in relation
to total length 29
16 The relationship between back calculated first and
second year length increment and gizzard shad
population density 30
VI
-------
Number Pagt
17 Three-dimensional graph of the logarithm of gizzard
shad numbers impinged through time (May 1975 to
April 1976) by 25 mm groups 32
18 Mean length of the 1974, 1975 and 1976 year classes
of gizzard shad impinged on the intake screens at
Cumberland Steam-Electric Plant by month 33
19 Log length - log weight relationship plotted against
length for gizzard shad collected in cove rotenone
samples in Barkley Reservoir, 1974-1978 34
20 Density of gizzard shad in Barkley Reservoir cove
rotenone samples (average of 1977 and 1978 samples)
by age class 35
21 Relationship between survival rate of gizzard shad
and population density 37
22 Increase in mortality rate with increasing length 38
23 Length frequency distribution in cove samples compared
with open water areas in the Crooked Creek
Study Area 40
24 Comparison of age structure in cove rotenone samples
with results of simulations run under four test cases
(simulated populations under baseline and zero plant
mortality situations are represented by the same
line) 47
25 Length frequencies simulated using baseline conditions
compared with cove rotenone estimates 48
26 Length frequencies of the simulation in which the
cohort is divided into nine subcohorts with
different growth rates, using baseline conditions,
compared with rotenone estimates 49
A.I The integral surfaces defined by Eqs. (A.17) and (A.18).
The intersection of these surfaces is the
characteristic curve 57
A.2 The surface N(s,t) as defined by Eq. (A.21). It is
composed of characteristic curves 58
A.3 A hypothetical recruitment rate, B(t), between the times
t=0 and t=T 61
s
A.4 Plots of B.{t - (l/g0)ln(s/s_} as functions of time, t,
for several values of size, s, in units of millimeters.
The reproduction function, B.(t) is a truncated normal
(Figure A.3) 62
vii
-------
Number
A.5
A.6
A.7
BQ{t - (l/g0)ln(s/s0}/gQs) as function of size, s, for
three values of time, t. The reproduction function,
Bn(t), is a truncated normal (Fig. A.3)
N(s,t) from Eq. (A.32) as a function of size, s, for
several values of time, t, and for arbitrarily chosen
parameter values
The mean size in the cohort (black dots) and the standard
deviation in size (white dots) as functions of time
from Eq. (A.32). The parameter values have been chosen
arbitrarily
63
66
67
Vlll
-------
TABLES
Number Page
1 Mean estimated total length (mm) at each annulus for
gizzard shad collected in Barkley Reservoir
(1974-1978) 23
2 Estimated age distribution of gizzard shad (Dorosoma
cepedianum) in cove rotenone samples for the year
1977 26
3 Age distribution of gizzard shad (Dorosoma cepedianum)
in cove rotenone samples for the year 1978 27
4 Analysis of variance for effects of the logarithm of
length and population density on the logarithm of
weight 31
5 Percent of survival and mortality of each age class
between 1977 and 1978 36
6 Cove rotenone population density estimates, yearly
impingement rates, (August to July), estimates of
total, impingement and natural mortality rates,
1974-1977 41
7 Length frequency of gizzard shad in cove rotenone
samples by age class compared with the results of
simulations run under four test cases 45
B.I The input data cards relevant to the example. The
meaning of the individual cards is given in
section 2 83
B.2 The form in which the input data is printed out by
the program 84
B.3 Predicted numbers of larvae per 1000 cubic meters in
the first eight 0.5 millimeter length classes for
twenty time periods 85
B.4 Predicted numbers of larvae per 1000 cubic meters in
the first eight 5.0 millimeter length classes for
twenty time periods 86
B.5 Predicted total population numbers per 1000 cubic meters
for twenty time periods 87
IX
-------
LIST OF ABBREVIATIONS AND SYMBOLS
ABBREVIATIONS
cm
ha
kg
m
MWe
mg
mm
-- centimeters
-- hectares
— kilograms
— meters
-- megawatt electrical
-- milligrams
-- millimeters
SYMBOLS
B(s,t)
B0(t)
As
At
G(s,t)
G(s)
KSO,S)
J(sQ,s)
N(s,t)
N
N
t+1
t+1
s
S0
s
max
t
TL
T
s
Z(s,t)
Z(s)
fecundity going into size class s
fecundity, all going into size class s»
increment in size
increment in time
growth rate as a function of size and time
growth rate, assumed dependent only on size
growth rate constant
integral, defined in Eq. (A.23)
integral, defined in Eq. (A.24)
size distribution function
number of fish at start of year
number of fish at end of year
number of young-of-the-year fish at end of year
scale radius
size
minimum size of fish in the cohort
maximum size of fish in the cohort
time
total length
length of spawning period
mortality rate coefficinet as a function of size and time
mortality rate, assumed dependent only on size
-------
ZQ — constant mortality rate
Z.f(s,t) ~~ fishing mortality rate coefficient
Z (s,t) -- natural mortality rate coefficient
XI
-------
SECTION 1
INTRODUCTION
Section 3l6(b) of the Clean Water Act has focused much attention on
the impingement of fish on water intake screens at power plants. Although
well documented, impingement is a poorly understood phenomenon. In many
instances, the species and sizes of fish affected are known to have swim-
ming speeds sufficiently greater than needed to avoid entrapment. While
it is frequently possible to relate unusually heavy impingement of some
species to cold shock, heavy losses of cold tolerant species (or cold
sensitive species during periods of warm water temperatures) continue to
puzzle investigators.
It is interesting to note that trawls fished for commercial fish
stocks are generally towed at speeds well below either burst or sustain-
able swimming velocities of the target species. Therefore, as a first step
in attempting to understand the impingement phenomenon, it may be produc-
tive to think of water intake screens as "stationary trawls." Water is
pulled through this "stationary trawl" rather than the trawl being pulled
through the water. Indeed, intake screen "catches" have been compared to
conventional trawl catches in order to gain some perspective of the
magnitude of impingement losses.
Assessments required by Section 3l6(b) of environmental impacts owing
to impingement have taken many forms. In some cases, although no attempt
was made to relate impingement mortality to the affected population, it
was the opinion of the investigator(s) that estimated impingement losses
were not an adverse impact to the fish community (Duke Power Company,
Undated). Other investigations have related impingement losses to estimated
standing stocks in the affected water body, and while noting that some
fraction of the stock was killed, similarly opined that such losses did
not constitute an adverse impact (Commonwealth Edison Company, 1977).
Another approach (Tennessee Valley Authority, 1977a) translated
impingement losses into reservoir area of lost production (i.e., for an
estimated standing stock of 10,000 fish/ha coupled with impingement of
2,000,000 fish annually, a loss of 200 ha of production is predicted.
For a 20,000-ha reservoir, this is 1 percent of the total area). A fourth
approach compared impingement losses with commercial catch for a given
species (Moseley, et al., 1975). Since numbers of impinged fish are always
much less than the commercial catch, impact is presumed to be small or
negligible. Although most of these "impact assessment" techniques put
impingement losses into at least a reasonable perspective, none truly
addresses impact; i.e., the actual or predicted modification of fish
communities by an extraneous mortality source.
In the case of new or proposed facilities which take in considerable
volumes of water, preoperational studies of the near- (and/or far-) field
fish community.can be compared to similar studies during the operational
phase of the facility. Although this approach would appear ideal for
determining the impact of such a facility to the resident fish community,
-------
it is fought with uncertainty. Documented changes in the far-field fish
community during an extended study period may not necessarily be due to
plant operation since changes through time are frequently observed in
waters where facilities do not exist (Tennessee Valley Authority, 1978).
Additionally, sampling techniques may not be sufficiently sensitive to
reliably determine subtle far-field impacts.
While determination of impacts in the above situation seems tenuous
at best, in the case of existing facilities for which preoperational data
are unavailable, criteria for assessing operational impacts to the far-
field fish community are even more obscure. The use of documented changes
during a protracted study period is obviously subject to the earlier given
criticisms for such an approach. It is clear that the existing fish com-
munity "is what it is" under the operating regime of the water intake
facility. What is not clear is what the far-field fish community in
question would be in absence of the plant.
If observable adverse modifications of fish communities cannot
reasonably be coupled with cause-effect relationships owing to plant
operation, or if the fish community was not studied prior to operation of
the facility, the investigator is reduced to two approaches: (1) offering
an "expert" opinion as to probable impact, or (2) modeling the fish commu-
nity (or population) and the effect of plant operation on it.
It is in this latter vein that this project was undertaken. A model
which describes fish populations in terms of both life process dynamics and
their interaction with facilities which impose additional mortality is
much needed. The purpose of this task was 2-fold: (1) to develop such a
model, and (2) to remove as much of the subjective process of impingement
impact prediction as possible.
-------
SECTION 2
CONCLUSIONS
The parameters needed by the model were estimated from data of a
gizzard shad population in Barkley Reservoir. The model was then used
to simulate the total population and length distribution under (1) baseline
or current conditions, (2) the use of assumed zero plant mortality, (3) a
10-fold increase in plant mortality, and (4) a 100-fold increase in plant
mortality.
The model proved to be effective in simulating the observed total
population numbers, though its predictions of the length distributions
were narrower than those of the observed population. Model simulations
indicated that present plant-induced mortality lowers total numbers in
each age class by less than 1 percent. A 10-fold increase in plant mor-
tality is predicted to cause at worst a 10-percent decrease of age class
IV. A 100-fold increase in plant mortality would have significant effects
according to the model, lowering the number of fish in age class IV by
about 65 percent.
Apparently, the gizzard shad population in Barkley Reservoir is
virtually unaffected by impingement at current levels. However, these
results cannot be generalized to other species which do not possess great
compensatory powers.
-------
SECTION 3
RECOMMENDATIONS
The unexpectedly low impingement impacts predicted for gizzard shad
in Barkley Reservoir were at least in part due to the compensatory powers
possessed by this species. Relatively short lifespan, density dependent
growth, and high natural mortality all contributed to greatly reduce the
effects of impingement losses for this species in Barkley Reservoir.
Therefore, these results probably cannot be generalized to other species
and the model should be tested with a relatively long-lived species which
does not typically possess great compensatory powers.
Model predictions for impingement impact levels approximating those
in both cooling ponds and rivers where virtually the total flow is entrained
should be compared with actual data from such situations. This would pro-
vide a valid test of the model's predictive abilities. Should relatively
great modifications be required to achieve accurate predictions, they could
well lead to new approaches in modeling the field of fish population
dynamics.
The model is more flexible and has greater capabilities than were used
in this study. It could be used to predict the total impact of intake-
related mortalities (entrainment as well as impingement) since the model
permits inclusion of all sizes and ages of fish if data are available.
Therefore, it is recommended that the model be tested using all sizes of
a particular fish species affected by intakes.
-------
SECTION 4
MATERIALS AND METHODS
THE MODEL
All of the individuals of a particular fish species in a given body
of water collectively form a unit termed a population. These individuals
undergo three distinct processes: birth, growth, and death. Although
these processes are properties of individuals, taken collectively they form
the basis for what is known as dynamics of the population.
Birth
Virtually every current fish population model assumes that all indi-
viduals of a cohort are born (or hatch) at the same time. This obviously
simplifies the determination of age, growth, and mortality rates through
time. However, hatching typically occurs over a period of time as depicted
in Figure 1. The essential elements of this figure are the time when hatch-
ing begins, builds to a peak, declines, and finally ceases. This period
of hatching may be from a few days to several months in length, depending
upon the species in question. In the latter instance, the simplifying
assumption that all hatch on the same day is not only unwarranted, it is
in serious error.
Returning to Figure 1, note that all members of a cohort are shown
in this graph, and as such the collective property (birth or hatching)
depicted becomes a population process. This graph is actually a 3-variate
or 3-dimensional figure of numbers through time for which the third variable,
size, is fixed. In this case, the "fixed size" is length at hatching. It
could just as easily have been "fixed" for some other length of "significance";
i.e., the length(s) at which a cohort becomes vulnerable to impingement.
However, the choice of length at hatching has two particularly important
attributes. The area beneath the curve in Figure 1 is the total number
hatched. This is determined by integrating an equation describing this
curve and the slope of this integral form is the hatching rate.
The concepts of three dimensionality and fixing of one variable are
important as they will form the basis for mathematical description of
the model developed in this report.
Growth
As mentioned previously, growth is a property of individuals. When
growth is considered collectively for a cohort, it becomes a cohort or
population process.
Figure 2 is a commonly used graphic representation of cohort length-
frequency. Numbers and length are seen to vary for a "fixed" or particular
instant in time. Although the members of the cohort are all in the same
year of life, they are not of exactly the same age as noted earlier.
While the scatter of observed lengths is in part due to differences in age
-------
HATCHING
5
3
Z
I
INITIATION PEAK CESSATION
TIME
Figure 1. The hatching of a cohort of fish through time
(the reproduction season).
-------
V)
D
Z
LENGTH
Figure 2. A cohort of fish at a given instant in time,
-------
(thus a longer period of growth), differential growth of individuals is
also a factor. The positive skew is a typical pattern although others,
such as biomodality, occur.
If a cohort is followed through time, the pattern shown in Figure 3
will typically emerge. Here time is effectively "fixed," at several points
to study growth. A more common graphical representation of this phenomenon
is shown in Figure 4. In this instance numbers are "fixed." The scatter
of lengths in Figure 4 is the same as that for the respective groups in
Figure 3. Growth is simply determined in Figure 4 as the average length
for the various age groups.
Mortality
Although death or mortality of an organism also is clearly a property
of individuals, it too can be modeled as a population process. Determina-
tion of mortality rate is readily accomplished from the idealized data of
Figure 3. The total number of cohort members extant in time periods I
through IV is determined simply as the area under the respective curves
(length-frequency distributions). This can be plotted as in Figure 5,
the resulting curve being a mortality rate which describes decline in
cohort numbers through time.
The Model
The length-frequency distributions of Figure 3 are erected on Figure 4,
matching their common length ranges and age groups, a 3-dimensional figure
will be produced. This is shown in Figure 6. Figure 1 has also been imposed
on this graph and is the far left distribution. Note that its "plane of
reference" differs from the other distributions. The variables along the
three axes are now free to vary simultaneously.
This then is the conceptual model. Hatching, growth, and mortality
are linked through time in a framework that combines conventional concepts
and data in a model which is completely new to the field of fish population
dynamics. It is a model which, while data dependent, is extremely flexible,
readily understandable in concept, and does not depart from traditional
fish population dynamic concepts.
To be useful, the model must be quantitative as opposed to the
dimensionless values used to develop it up to this point. This will
require the mathematical description developed in subsequent sections.
MATHEMATICAL DESCRIPTION OF THE CONCEPTUAL MODEL
Formulation of a partial differential equation model
Partial differential equations were probably first applied to popula-
tion dynamics by von Foerster (1958). Such equations are useful in describ-
ing a population through time, not only in terms of numbers of individuals,
but also age or size distributions of the population. Description of
length distributions is especially important since most fisheries are con-
cerned with fish numbers and sizes. Also, the inclusion of age or size
-------
LENGTH
Figure 3. The same cohort of fish in Figure 2
shown at four instances in time.
-------
AGE (OR TIME)
Figure 4. Graphic representation of the cohort of fish in
Figure 3 to better show growth.
10
-------
CO
j L
AGE (OR TIME)
Figure 5. Mortality in a cohort of fish.
11
-------
LENGTH
~
Figure 6. Relationship of hatching, growth, and mortality
in a cohort of fish.
-------
structure implies a time-lag effect in the population, essential in properly
describing the dynamics. Partial differential equation models have been
used to describe the age and size structures of populations by Trucco
(1965a,b), Oldfield (1966), Weiss (1968), Sinko and Streifer (1967, 1971),
Langhaar (1972), Oster and Takahashi (1974), Rubinow (1973), Rotenberg
(1975), Griffel (1976), and Van Sickle (1977), among others.
Before developing the model, two population characteristics of
interest are described, mortality and growth of the individual fish.
The rate of growth of an organism can often be approximated by the
differential equation
= G(s(t),t) , (1)
where s(t) is some measure of organism size (e.g, length) as a function of
time, and G(s,t) is a function of both size and time. When G(s,t) takes
some simple functional form, independent of t, for example,
G(s) = g()(l - ^— )s , (2)
max
where g- and s are constants, then Eq. (1) yields the solution,
u rndx
*o(t - V
s-s e
0 max
_ _
- tt ft - t 1
(s - SA) + sAe80U V
max 0 0
where tQ is some initial time and s,. is the initial size.
The mortality of fish in a single cohort is typically modeled as,
= -Z{s(t),t}N(t) , (4)
where N(t) is the number of organisms and Z{s(t),t} is a mortality func-
tion, in general dependent on both time and the average size of organisms
in the cohort. For simple cases, Eq. (4) is easily solved. For example,
if Z{s(t),t} is merely a constant, Zn, then
-Z (t - t )
N(t) = NQe ° ° (5)
where NO is the number in the cohort at time to.
Size distribution, N(s,t), means the number of organisms per unit
size class at a given time. An equation will be derived to describe the
change in N(s,t) through time. Consider, at time t, the number of organ-
isms in a size class of width As centered around the size s, and denote
this number by N(s,t) (see Fig. 7). Assume a very short interval of time,
At, passes. What is the new distribution at t + At? Conceptually, the
number of fish in size class s and time t + At can be symbolized as
13
-------
ORNL-DWG 79-9243
ND(L+AL,t+At)
At
Figure 7. The three curves represent size distributions at three successive
time intervals. The distributions N(s,t) is a continuous distribution,
while N-(L,t) is a discrete approximation of this distribution.
-------
N(s,t + At) = N(s,t)
- {number of fish lost to mortality}
- {number of fish lost because of growth to the next
larger size class}
+ {number of fish entering because of growth from
the next smaller size class}
+ {newly reproduced fish of size s}
+ {immigration from other populations}
- {emigration to other populations}. (6)
Ignoring immigration and emigration, the conceptual model will be
expressed in mathematical terms. The following assumptions and equations
form the basis for this model.
Number of Fish Lost to Mortality--
Assume this is proportional to the number of fish in size class s,
N(s,t), the mortality coefficient, Z(s,t), and the time interval, At.
Number of Fish Lost Because of Growth to the Next Larger Size Class--
Assume this is proportional to the number of fish in size Class s,
N(s,t), and the fraction of that size which grows to the next larger size
class, G(s,t) (At/As). To help understand this term, note that G(s,t) At
represents the amount of a particular size class population that moves
from one size class to the next during interval At, so G(s,t) (At/As) is
the fraction of a given size class population that moves.
Number of Fish Gained (Recruited) Because of Growth from Next Smaller
Size Class--
Analogous to the above, this is G(s - As,t)(At/As)N(s-As,t).
Number of Fish Gained in Size Class s from Reproduction—
This is a product of the number of reproductives, N (t), the
fecundity, M(t), and the fraction, T(s), of the fecundity going into
size class s; i.e., B(s,t) = Nre (t)M(t)T(s).
Combining the above terms yields
N(s,t + At) = N(s,t) - AtZ(s,t)N(s,t) - (At/As)G(s,t)N(s,t)
+ (At/As)G(s - As,t)N(s - As,t) + B(s,t)At , (7)
for the numbers of size s at time t + At (Fig. 7). To obtain a govern-
ing equation from Eq. (7), rearrange (7) as
N(s.t + At) - N(s,t) + G(s,t)N(s,t) _ G(s - As.t)N(s - As.t)
At As As
= -Z(s,t)N(s,t) + B(s,t) (8)
15
-------
and take the limits
lim N(s,t + At) - N(s,t) _ 3N(s,t)
At->0 At at
(G(s.t)N(s.t) - G(s - As,t)N(s - As,t) = &
At-»0 As ~ 9
to obtain
9N^>t) + |^ (G(s,t)N(s,t)} = -Z(s,t)N(s,t) + B(s,t). (11)
This equation must be supplemented with initial conditions on N(s,t)
at the initial time, t_. In the remainder of this report, we shall let
t0 = 0 and assume that this instant precedes production of any of the
cohort, so that
N(s,0) = 0 (for all s) . (12)
Equation (11) is the general cohort model used in this report.
Appendix A shows how this equation can be solved analytically, at least
to the point where only certain integrals are left to be done numerically.
This appendix also presents certain limiting cases in which the solutions
are completely analytic. The remaining appendices describe the computer
programs that numerically solve the general model in which both the growth
rate, G(s,t), and the mortality rate, Z(s,t), are size and time dependent
and in which the reproduction (or recruitment) rate, B(s,t), varies over
a spawning period.
Appendix A also describes two further generalizations of the model.
First, rather than being necessary to prescribe a reproductive rate, B(s,t),
it is possible simply to prescribe an initial size distribution, N_(s)
at time t,.. Second, rather than each fish in the cohort being given the
same growth rate, G(s,t), there is a provision for dividing the cohort
into subcohorts, each subcohort with its own growth rate, G.(s,t).
DATA SOURCE
A large volume of computer-based fisheries data is available for
Barkley Reservoir and TVA's Cumberland Steam-Electric Plant. As a con-
sequence, this system (plant and reservoir) was chosen as a test case for
the model. The gizzard shad (Dorosoma cepedianum) was chosen as an
example species because of its common occurrence in impingement samples
and widespread distribution. Also, it is not as subject to winter kill
as the threadfin shad (D. petenense) and this, coupled with greater
longevity, results in more nearly stable populations.
The plant, reservoir, and attendant gizzard shad population are
described in greater detail below.
16
-------
BARKLEY RESERVOIR
Barkley Reservoir is a mainstream impoundment on the lower Cumberland
River. Barkley Dam, located at Cumberland River Mile (CuRM) 30.6 near
Grand Rivers, Kentucky, was closed in 1964. The reservoir extends approxi-
mately 103 km to Cheatham Dam (CuRM 148.7). At normal full pool (108 m
above MSL) Barkley Reservoir has a surface area of 23,440 ha and a shore-
line length of 1,615 km. Barkley Reservoir is connected with Kentucky
Reservoir (a Tennessee River impoundment) by Barkley Canal, allowing
integrated operation of the two reservoirs for navigation, flood control
and hydroelectric power production.
CUMBERLAND STEAM-ELECTRIC PLANT
Physical Features
Cumberland Steam-Electric Plant is located on the south bank of
Barkley Reservoir in Stewart County, Tennessee, at Cumberland River Mile
(CuRM) 103 (Figure 8). This facility is the Tennessee Valley Authority's
largest fossil fuel plant, having two units rated at 1,300 MWe each. The
first unit began commercial operation in March 1973 and the second unit
began operation in November 1973.
Condenser cooling water is provided for this plant thorugh a 200-m
long intake channel located perpendicular to the shoreline (Figure 9).
A 339-m long skimmer wall located at the mouth of the intake channel
allows selective withdrawal of cooling water from the lower reaches of
the reservoir. Maximum total condenser and auxiliary flow is 120 m3/sec.
The intake pumping station at Cumberland Steam-Electric Plant con-
sists of eight condenser cooling water pumps which draw water through
sixteen intake chambers. At the mouth of each intake chamber are trash
racks which keep out large debris. Behind the trash racks are the vertical
traveling screens. These are a series of 2.3 m by 0.6 m interlocking
screen panels attached to endless chains which rotate on sprockets
located at the top and bottom of this intake structure. The screen panels
consist of 12-gauge galvanized wire having 9.5 mm square openings.
Impingement
Impingement of fish on the screens at cooling water intakes consti-
tute a potential adverse ecological impact. In response to the passage
of Section 316(b) of the Clean Water Act, the Tennessee Valley Authority
began a weekly fish impingement monitoring program at Cumberland Steam-
Electric Plant in August 1974. Each Tuesday, all vertical traveling
screens in operation were rotated and washed clean of fish and debris.
Twenty-four hours later, the screens were again cleaned and the impinged
fish flushed into a sluice pipe leading to a catch basin. All fish were
separated by species into 25 mm length groups and total numbers and
weights were recorded by species and 25-mm group.
17
-------
ORNL-DWG 80-9341
-BARKLEY DAM
BARKLEY
CANAL
KENTUCKY
RESERVOIR
kilometer
KENTUCK^
TENNESSEE
CUMBERLAND
CITY
BARKLEY
RESERVOIR
CUMBERLAND RIVER
BARKLEY RESERVOIR
0 5 10
i . i . . i i
WELLS CREEK
CUMBERLAND STEAM PLANT-
CLARKSVILLE
CHEATHAM
DAM
Figure 8. Location of Cumberland Steam-Electric Plant on Barkley Reservoir.
-------
ORNL-DWG 80-9340
LEGEND:
+ RIVER MILE
DIRECTION OF FLOW
feet
0 1000 2000 3000
SKIMMER
WALL
INTAKE
STRUCTURE
CUMBERLAND
STEAM PLANT
Figure 9. Intake and associated structures at Cumberland Steam-Electric Plant.
-------
Results of this monitoring program have been presented previously
(Tennessee Valley Authority, 1977a). This study estimated the total
numbers impinged yearly and, using cove rotenone population estimates,
computed plant exploitation rates for each species. Two additional
studies, McDonough and Hackney (1978) and McDonough and Hackney (1979)
related impingement rates to physical factors (temperature and water
elevation), plant intake design, and life history of the species.
GIZZARD SHAD
Estimates of gizzard shad population density and size structure in
Barkley Reservoir were made each summer from 1974 through 1978.using cove
rotenone samples (Figure 10). Field procedures for treatment and collec-
tion of data followed standard cove rotenone techniques.(Tennessee Valley
Authority 1977b). Coves were blocked with 1.3 cm bar mesh nets and 5 per-
cent emulsifiable rotenone was applied at a concentration of 1.0 mg/liter.
All fish collected during a 2-day period were identified to species and
grouped by 25 mm length classes. Numbers and weights were recorded for
each length class. Ten coves were sampled each year.
In 1974, 14,390 gizzard shad weighing 692.4 kg/hectare were collected
in cove rotenone samples in Barkley Reservoir (Figure 14). The population
was dominated by large fish with low numbers of young-of-year fish. Both
numbers and biomass decreased steadily to a low of 4,873 fish weighing
191.7 kg/hectare in 1976. Numbers then increased greatly (15,369 fish/
hectare) with only a slight increase in biomass (229.6 kg/hectare),
apparently the result of a strong year class in 1977. In 1978 numbers
per hectare decreased to 12,906.4 fish per hectare while biomass increased
to 365.7 kg/hectare, apparently due to growth of the abundant 1977 year
class.
Estimation of Age, Growth and Mortality
Gizzard shad collected from gill nets, hoop nets, electrofishing and
cove rotenone samples were used in age and growth studies. Scales were
taken from an area near the front of the dorsal fin above the lateral line.
Scale impressions were produced on cellulose acetate strips using a roller
press and then projected using an Ederback model 2700 microprojector (40 X).
Between September 1974 and August 1978, scale samples were taken from 679
gizzard^shad for the determination of age and growth characteristics.
Back calculations of length at annulus formation provided information on
growth rates of the 1969 through 1977 year classes.
A plot of total length to scale length (40 X) indicated that the
relationship was linear (Figure 11). The regression formula: TL = 85.86 +
0.83 SR where TL = total length and SR = scale radius was fitted using least
square procedures. The distance from the focus to each annulus was
substituted for SR in the above equation in order to calculate the length
of the fish at the time of annulus formation.
Calculated mean length of fish at annulus formation was determined
for each age group collected for the 1969 through 1977 year classes (Table 1)
Growth was most rapid during the first year and then decreased. Weighted
-------
15,000"
10,000-
5,000-
1 0,000.
5,000.
UJ
^ 10,000.
a
I 5,000.
cc
UJ
o.
UJ 15,000.
CO
3 10,000.
5,000.
15,000-
1 0,000-
5,000-
1974
-Tl-t-r
1
-L
1975
••••^•a
1976
i n ii i i , lii
1977
~~|^-, -J— h-i_
1978
JTT-Th-^
100 *
i CMT
i, I. £ }- fc £ >
§ | § § § * *
PER HECTARE
Lograms per Hectare) distribution
i Barkley Reservoir, 1974-1978.
•i-l -H
23
p,-« w ^ S W
2 B TJ
QC 0 M
O-H CO
h,W Q * S
^^ i-H -H
mmt CO 00
^ O 4H
.500 *
o
F— 1
cu
•1,500 j?
.1 ,000
.500
00 300
2 -
»TLI "s* T:^
oo>
21
-------
350
300
250
O
Z
LU
200
o
150
100
50
100
150
200
250
SCALE RADIUS
Figure 11. Total length - scale radius relationship of gizzard
shad from Barkley Reservoir, 1974-1978. Dots represent
means of 10 mm length groups; vertical lines represent
the mean plus and minus two standard errors.
22
-------
mean increments for years I-IV were 144.2, 41.1, 15.7, 12.5, and 28.5. The
increment for age V is of doubtful accuracy, because this estimate is based
upon a single fish. Lee's phenomenon of apparent change in growth rates
was evident.
TABLE 1. MEAN ESTIMATED TOTAL LENGTH (mm) AT EACH ANNULUS FOR
GIZZARD SHAD COLLECTED IN BARKLEY RESERVOIR (1974-1978)
Age Group
I
II
III
IV
V
Number
of Fish
237
266
122
19
1
I
149.47
143.55
136.97
133.79
144.08
II
190.13
177.44
169.29
191.48
III
202.20
192.93
206.46
IV
212.54
232.24
V
243.05
Numbers
Grand mean
645
408
142
20
1
Length
Increments
144.19
144.19
185.32
41.13
200.99
15.67
213.53
12.54
242.05
28.52
Growth rates tended to increase from 1969 through 1976 (Figure 12), per-
haps in response to decreasing stock density. Barkley Reservoir was impounded
in 1964. Since fish population density has been shown to decrease following
impoundment, the increasing growth rates may be a result of decreasing stock
density.
Scales were taken from representative individuals in each 25 mm length
group (> 100 mm) collected in 1977 and 1978 cove rotenone samples for
determination of fish age. Using the proportion of individuals of each age
in the subsample, age composition for each 25 mm size group was determined
(Tables 2 and 3). Fish smaller than 100 mm were assumed to be young-of-
year. Figure 13 shows length frequencies by age groups and for the entire
population.
23
-------
300
250
200
UJ
150-
100
AGE CLASS
Figure 12. Comparison of back calculated growth
rates of gizzard shad for the 1969
through 1976 year classes.
24
-------
LU
OC
o
lu
OC
UJ
Q.
10,000-
1,000-
100
10
1
0.1
10,000
'1,000
100
10
1
o.H
1977
50
100
150
200
250
1978
300
100
150
200
250
300
TOTAL LENGTH,
Figure 13. Estimated length frequency distribution of gizzard
shad in 1977 and 1978 cove rotenone samples.
25
-------
TABLE 2. ESTIMATED AGE DISTRIBUTION OF GIZZARD SHAD (DOROSOMA CEPEDIANUM)
IN COVE ROTENONE SAMPLES FOR THE YEAR 1977
Length (mm)
Minimum Maximum
Age Groups in Subsample
0 I II III Fish/ha
Calculated Age Representation
0 I II III
26 50
51 75
76 100
101 125
126 150
151
176
201
226
251
276
175
200
225
250
275
300
TOTAL
9
3 1
18
24
24
3
3
11
20
2
1
24
9,506
2,377
866
28
293
1,329
836
1 103
3 3
15,369
.8 24.8
.6 9,506.6
.5 2,377.5
.0 866.0
.1 21.1
.3
.1
.9
.0
.9
.2
.4 12,796
7
293
1,181
573
12
2,068
.0
.3
.4
.9
.9
.5
147
263
85
1
498
.7
.0
.8
.6
.2
.3
4.
2.
6.
3
3
6
These age composition and length frequency distributions were used to
estimate growth and mortality rates. The length distributions and age sub-
samples from 1977 and 1978 cove rotenone samples were averaged and used to
determine growth characteristics. A pattern similar to that found from
back calculation of lengths at each annulus was evident (Figure 14). Growth
was rapid during the first year, then declined during the second and third
year. The fourth year estimated growth was unexpectedly high; however, this
estimate was based upon a single fish.
Yearly length increment was plotted against length (Figure 15). Young-
of-year gizzard shad averaged 77.91 mm and were approximately three months
old when the samples were taken. Therefore, young-of-year fish were growing
at an estimated 312 mm per year during this period. The experimental
equation: Length Increment = 504.32 length ' was found to fit this
relationship very well (excluding growth during the fourth year for reasons
mentioned previously), having a correlation coefficient of 0.990.
Population density appeared to have an effect on growth rate. Both
length at annulus formation, and the change in young-of-year length distri-
bution through time for shad impinged at Cumberland Steam-Electric Plant
were found to decrease as population density increased. The length-weight
relationship, an indirect indicator of growth differences was also found
to decrease with population density.
Back calculated length at first annulus formation and length increment
between the first and second annulus both tended to decrease with increas-
ing population density (Figure 16). The 1974 year class showed a high
growth rate despite high population density. The estimated growth rate
for 1974, however, was based upon only five individuals while estimates
of growth of the 1979, 1976, and 1977 year classes were based on 49, 144,
and 156 individuals, respectively.
26
-------
TABLE 3. AGE DISTRIBUTION OF GIZZARD SHAD (DOROSOMA CEPEDIANUM) IN COVE ROTENONE SAMPLES
FOR THE YEAR 1978 AS DETERMINED BY SCALE READINGS OF SUBSAMPLES FROM EACH
SIZE GROUP.
Length (mm)
Maximum Minimum
26
51
76
101
126
151
176
201
226
251
276
301
TOTAL
50
75
100
125
150
175
200
225
250
275
300
325
Age Group in
0 I
1 1
13 22
7 39
42
50
8
II
1
12
32
27
4
Subsample
III
1
3
6
3
1
IV Fish/ha
2
2
3
1
1
1
12
474.3
,527.5
,353.5
684.2
,386.1
,648.9
,055.2
711.7
56.8
6.7
1.2
0.3
,906.4
Calculated Age Representative
0
474.3
2,527.5
2,353.5
342.1
1,257.7
245.6
7,200.7
I
342.1
2,128.4
1,368.2
820.7
428.7
12.0
5,100.1
II
35.1
234.5
274.4
40.4
2.7
587.1
III
8.5
4.5
4.0
1.2
.1
18.3
IV
.1
.1
-------
300
TOTAL LENGTH
250
O
LU
200
LU
150
100
50
/ LENGTH INCREMENT
I II
AGE CLASS
Ml
IV
Figure 14.
Mean length and mean length increment of gizzard
shad in 1977 and 1978 cove rotenone samples.
28
-------
300
z
LU
2
LU
fiC
O
200
O
UJ
b
CC
LU
100
50
100
150
200
TOTAL LENGTH, MM
Figure 15. Yearly length increment versus total length.
29
-------
155-
2 150-
Z
111
5
Hi 145'
DC
£ ^
?E 55
O
UJ
50-
45-
AGE 0 »1974
• 1976
• 1975
• 1977
~7
/
• 1975
•1974
AGE I
• 1973
• 1976
~t.
5,000 10,000
NUMBER PER HECTARE
15,000
Figure 16. The relationship between back calculated first and
second year length increment and gizzard shad
population density.
30
-------
The effect of population density on young-of-year gizzard shad growth
rate was also reflected in the length-frequency distributions of fish
impinged at Cumberland Steam-Electric Plant. Figure 17 is a 3-dimensional
graph of the length distributions of gizzard shad impinged from May 1975
through April 1976. Numbers are expressed as logarithms to permit greater
detail for all sizes. Young-of-year fish first became large enough to be
impinged in late July or early August. Growth after this period resulted
in a gradual increase in length distributions through the remainder of the
year. Although this same pattern was observed each year, the rate of
growth differed each year, apparently due to differences in stock density.
The average length of young-of-year fish was computed for each month for
the 1974, 1975, and 1976 year classes using length frequency analyses to
separate young-of-year fish from older individuals. The growth rate
increased each year (Figure 18), while population densities decreased.
In addition to its usefulness in converting length to weight, the
length-weight relationship is also useful as an indicator of population
growth differences. The relationship was computed for gizzard shad col-
lected during the first day of cove rotenone samples. Fish were grouped
into 25-mm groups, and numbers and weight of individuals in each size
group were recorded. Mean weight per individual in each size group was
computed by dividing the total weight by the number in the size class.
The relationship: Log (weight) = -4.65 + 2.828 x log (length) describes
this relationship having an R2 of 0.97. The deviations from this regres-
sion for each year were plotted against population density (Figure 19).
In years when numbers were low, the weight of fish at a given size tended
to be higher than in years when densities were high. Sequential F-test
(Draper and Smith, 1966) of the relationship between mean weight and popu-
lation density (with length already in the model) indicates that this
relationship is highly significant (Table 4).
TABLE 4. ANALYSIS OF VARIANCE FOR EFFECTS OF THE LOGARITHM OF LENGTH
AND POPULATION DENSITY ON THE LOGARITHM OF WEIGHT
Source of
Variation
Degrees of
Freedom
Sum of
Squares
Mean
Square
Prob F
Corrected Total 471
Regression 2
due to log ..(length) 1
due to density given 1
Iog1() (length)
Residual 469
471.0480
216.1010
215.9516
0.1494
6.9470
108.0504
215.9516
0.1494
0.0148
7,294.60 0.0001
14,579.11 0.0001
10.08 0.0016
The average of the age distributions in the 1977 and 1978 cove rotenone
samples was used to estimate annual total mortality rates. Numbers of
gizzard shad (on log scale) for each age is plotted against age. The
slope of the descending right limb describes the instantaneous mortality
rate (Figure 20). Mortality rates were found to increase each year. Total
31
-------
Figure 17. Three-dimensional graph of the logarithm of gizzard
shad numbers impinged through time (May 1975 to
April 1976) by 25 mm groups (N = numbers, T = time
in weeks, S = size).
32
-------
140
120
10°
UJ
-I
iy
80
1975
40 J
N
MONTH
M
•1974
Figure 18, Mean length of the 1974, 1975 and 1976 year classes
of gizzard shad impinged on the intake screens at
Cumberland Steam Electric Plant by month.
M
33
-------
0.06-
0.04-
0.02-
j
o
o-
LU
DC
-0.02-
-0.04-
2
i
• 1976
i
• 1975
1
H978
'
i
»1974
•1977
5000
10,000
15,000
NUMBER PER HECTARE
Figure 19. Residuals of the log length - log weight
relationship as a function of population
density.
34
-------
10,000
1,000
UJ
cc
u
*
cc
S
CC
UJ
10-
1-
0.1
II
III
IV
AGE CLASS
Figure 20. Density (number per hectare) of gizzard shad in
Barkley Reservoir cove rotenone samples (average
of 1977 and 1978 samples) by age class.
35
-------
annual mortality rates were also computed by comparing the change in num-
bers of fish in each year class collected in the 1977 cove rotenone samples
to abundance of the same year class in the 1978 samples (Table 5). Percent
mortality rose as fish age increased: 60 percent for the young-of-year
fish in 1977, 72 percent for age I, 96 percent for age II, and 98 percent
for age III.
TABLE 5. PERCENT OF SURVIVAL AND MORTALITY OF EACH AGE CLASS
BETWEEN 1977 AND 1978
Year
Class
1977
1976
1975
1974
Number
in 1977
12,796.0
2,068.5
498.3
6.6
Number
in 1978
5,100.1
587.1
18.3
0.1
Survival
Rate
0.3985
0.2839
0.0367
0.0152
Mortality
Rate
0.6015
0.7162
0.9633
0.9848
Instantaneous
Mortality Rate
0.920
1.259
3.304
4.189
The percent of fish surviving from one year to the next was estimated
for each year using the formula:
N
- N
s =
t+1
where S is percent survival, N is the total number of fish at some
time, N is the total number of fish one year later, and N_ is the
t+1 ut + 1
number of young-of-year fish one year later. The age 0 of fish were separated
from older fish using subsamples of each size group as described previously
for the 1977 and 1978 samples. Length frequency analyses were used to iden-
tify young-of-year fish in other years. The survival rate was found to
decrease with increasing population density (Figure 21). Survival from 1977
to 1978 was high with high population density; however, the population in
1977 was dominated by young-of-year fish, which have a higher survival rate
than older fish.
Mortality rates were found to increase in proportion to increasing fish
length (Figure 22). The equation: Z = -1.7926 + 0.0207 length was found to
be the best linear description of this relationship, explaining 81 percent
of the variation. Second-year mortality, however, was over estimated while
third-year mortality was under estimated. The increase in mortality rate
with length (or age) is unexpected, mortality would be expected to be high-
est during the first year, due to heavy predation. As shad become larger,
the mortality rate would be expected to decrease because they are no longer
available as prey to most predatory fish (Bodola 1965).
An apparent increase in mortality with increasing age would result from
migration of older individuals to open water areas. Bodola (1964) found
that in western Lake Erie young-of-year and age I fish were close to shore
in mid-summer; usually in shallow water quite close to shore. As gizzard
shad became older, they were found to inhabit deeper water.
36
-------
60-
•1976
50-
40-
1977
30-
o
in
1975
• 1974
20-
10-
5,000 10,000
NUMBER PER HECTARE
.15,000
Figure 21. Relationship between survival rate of gizzard
shad and population density.
37
-------
5-
4-
UJ
3-
2-
cc
O
1-
100
200
300
TOTAL LENGTH, MM
Figure 22. Increase in mortality rate with increasing length.
38
-------
In order to examine the question of how well cove samples represent
open water areas of a reservoir, the Reservoir committee of the Southern
Division of the American Fisheries Society sampled Crooked Creek Bay (85
hectares) in Barkley Reservoir. This embayment was divided into cove and open
water areas and treated with rotenone. The mouth of the bay was blocked with
nets and the open water was divided by nets into three areas. All coves in
the embayment were also blocked and the larger coves were subdivided into
compartments.
Length distributions of gizzard shad collected from the small coves
and ends of the large coves were compared with the open water areas (Figure
23). No evidence of migration of older fish from coves to open water was
found. Most length groups between 63.5 mm and 317.5 mm were overestimated
by the cove samples except for the 91.4 mm to 114.3 mm and the 218.4 mm to
241.3 mm groups. The open water samples tended to have more individuals
smaller than 63.5 mm and larger than 317.5 mm, probably because the large
sample area increased the likelihood of encountering individuals of these
sizes.
The convex shape of the catch curve may be due to high population
density. Michaelson (1970) and Anderson (1973) examined the dynamics of
balanced and unbalanced bluegill (Lepomis macrochirus) populations in ponds.
Balanced populations were defined as those with a satisfacotry relationship
between the fish population and its food supply (Anderson 1973). Balanced
populations were characterized by high and relatively stable recruitment of
young-of-year, while unbalanced populations were characterized by wide
variations in year class strength. They also found high mortality of young
fish in balanced populations. Unbalanced populations had low natural mor-
tality in young fish but high mortality for age II and older. Anderson (1973)
discussed the similarities between characteristics of pond bluegill popula-
tions and gizzard shad populations where wide fluctuations in year-class
strength may also be indicative of imbalance for reservoir populations of
this species.
Estimates of annual impingement were determined by integrating the area
between successive pairs of sample dates over a period from August 1 to
July 31 for each year. The beginning of this period approximates both the
first occurrence of a year class of gizzard shad on the intake screens and
the cove rotenone sample period.
Between August 1974 and July 1975, an estimated 165,050 (7.04 per
hectare) shad were impinged (Table 6). Impingement estimates decreased to
95,103 (4.06 per hectare) between August 1975 to July 1976. From August 1,
1976, to March 23, 1977, an estimated 190,267 fish were impinged. Expanding
this figure to a full year yielded an estimated 296,828 (12.66 per hectare)
gizzard shad impinged. This is probably an overestimate since fall and
winter are usually the periods of highest impingement.
39
-------
10,000
OPEN WATER
COVES
1,000-
UJ
CC
o
UJ
CC
Q.
CC
UJ
CO
100-
10
1-
0.1-
100 200
TOTAL LENGTH, MM
300
Figure 23. Length frequency distribution in cove samples
compared with open water areas in the Crooked
Creek embayment of Barkley Reservoir.
40
-------
TABLE 6. COVE ROTENONE POPULATION DENSITY ESTIMATES, YEARLY
IMPINGEMENT RATES,* ESTIMATES OF TOTAL, IMPINGEMENT
AND NATURAL MORTALITY RATES, 1974-1977
Numbers Per Hectare
Yearly Mortality Rate
Year
Total Initial
Population
Age 1+ at
End of Year
Numbers Impinged
Per Year
Total
Plant
Natural
1974
1975
1976
MEAN
14,390.09
10,484.93
4,872.59
9,915.87
3,582.92
3,380.84
2,573.40
3,179.05
7.04
4.06
12.66
7.92
1.3904
1.1318
0.6384
1.1375
0.0020
0.0012
0.0049
0.0025
1.3884
1.1306
0.6335
1.1351
* From August to July.
Using abundance estimates obtained from cove rotenone samples, it was
estimated that in the 1974-1975 season, 0.049 percent of the gizzard shad in
Barkley Reservoir were impinged. This figure decreased to 0.039 percent in
1975-1976, and increased to 0.206 percent in 1976-1977. Using averages for
this period: Z (the total mortality rate) was 1.1375, Z (the natural mor-
tality rate) was 1.1351, and Zf (the plant fishing mortality rate) was 0.0025.
41
-------
SECTION 5
EXPERIMENTAL PROCEDURES
As mentioned earlier, impingement impacts owing to fish mortality at
existing facilities are not readily assessed using conventional hypothesis
testing techniques. The existing fish community "is what it is" under the
operating regime of the facility. What it was or "would be" in the absence
of the plant is not an observation normally available to the experimentalist.
The impact of impingement mortalities on fish populations can be mani-
fested in three general ways or plateaus of impact. There may be: (1) no
impact to the fish population if growth and survival mechanisms are suffi-
ciently great to compensate impingement losses, (2) depression to some lower
level of abundance where compensatory mechanisms resist further decreases,
or (3) eventual extirpation of the fish population at some time. A fourth
possibility exists at the community level. Fish species for which individuals
are not normally impinged in numbers sufficiently high to adversely affect
their own populations may interact with species which are adversely impacted
by impingement. Community structure (e.g., predator-prey interactions) could
thus be modified as large impingement losses for one species are translated
through the fish community. This is a problem which cannot be addressed with
a single species model such as developed here except in the trivial case where
no impact to the population under study is predicted (where interacting
populations are presumed also to be unaffected).
With these factors in mind, a model was constructed which utilized an
estimated fishing mortality rate (i.e., the impingement rate) for Cumberland
Steam-Electric Plant in conjunction with the natural mortality rate estimated
for the gizzard shad population of Barkley Reservoir. This procedure, as a
result of deliberately varying the fishing rate coefficient, was used to
estimate the effects of the plant on the gizzard shad population as described
below. However, it is important to note that these efforts do not represent
experiments in the classical sense of hypothesis testing. Rather, they are
hypothetical experiments which are beyond the ability of the investigators
to perform and for which a model has been developed to predict outcomes.
BASELINE CONDITIONS
The term "baseline conditions" as used here refers to the population
as it exists with operation of the plant. This was accomplished by first
running the model with the mortality coefficients applicable for the gizzard
shad population as it exists with the operating plant. This provided a
reference with which to make assessments considering other operational
regimes for the platn.
ZERO PLANT MORTALITY
In this scheme, the model was run with the plant fishing mortality
coefficient set equal to zero. Thus, the gizzard shad population predicted
by the model is that expected for Barkley Reservoir if Cumberland Steam-
Electric Plant was not present or if impingement mitigation was undertaken
42
-------
and 100 percent effective. The difference between the zero mortality and
baseline condition estimates is the predicted impact of the plant on the
gizzard shad population.
TEN-FOLD MORTALITY INCREASE
For this experiment, the plant was assumed to impinge 10 times as
many gizzard as are presently estimated to be killed, with other baseline
parameters unchanged. This is roughly equivalent to siting a similar
facility on a lake one-tenth the size of Barkley Reservoir or, increasing
plant size by a factor of 10. This translates into about 0.9 MWe/ha, a
normal ratio for cooling ponds but considerably lower than ratios for
plants located on multipurpose reservoirs. Plant intake volume would be
about 10 percent of the reservoir per day. The prediction obtained in this
case perhaps reflects the impact of impingement alone to gizzard shad
populations in cooling ponds.
ONE-HUNDRED-FOLD MORTALITY INCREASE
This experimental run was similar to the previous test except that the
plant was assumed to impinge 100 times as many gizzard shad as are currently
estimated to be killed. Plant-to-reservoir ratios in this case were 9 MWe/ha
with a daily intake volume equivalent to 97 percent of the reservoir volume.
These ratios are not approached for plants located on reservoirs and are
rarely found in flowing streams, and then usually for only a small portion
of the year.
43
-------
SECTION 6
RESULTS AND DISCUSSION
The total numbers of fish in the cohort and cohort length distributions
were first simulated for baseline or current conditions at yearly intervals
(Table 7). The total numbers of fish simulated by the model for the first
year of life were generally very similar to cove rotenone sample estimates
(Figure 24). Simulated second- and third-year numbers were somewhat lower
than the actual measured numbers, while fourth-year model predictions were
very close to the numbers observed.
The model also predicted changes in the length distribution through
time. The predicted length distributions followed the same pattern as the
observed size distribution; that is, they become narrower through time.
This phenomenon is expected if the growth rate is decreasing as a function
of length (see DeAngelis and Mattice 1979, Bodola 1965, Ricker 1975), and is
known as compensatory growth. The model length distributions tended to be
narrower than the observed distribution (Figure 25). Despite these differences
in the variances of the length distributions between the model and observa-
tions, the mean length at the end of each year predicted by the model was
close to that observed except in the case of age IV. This is easily explained
by the fact that there was only a single individual of age IV sampled (see
Materials and Methods section).
The differences in the variances of the real and simulated size distri-
butions are probably due to the fact that there is a great deal of hetero-
geneity within the actual population (i.e., differences in growth rates
among individuals due to genetic and environmental variability) that was
not incorporated in the simulation shown in Figure 25. However, the model
is structured to permit inclusion of these differences within the cohort
in an approximate way by allowing division of the cohort into a number (up
to 20) of subcohorts, each with its own growth rate G.(s,t). Because of
this poor fit, a simulation was made in which the cohort was divided into
nine subcohorts, distributed roughly normally about the mean growth rate,
ranging from 0.6 to 1.4 of the mean value. The results of this simula-
tion are shown in Figure 26, where the numbers are the sums over the sub-
cohorts, these sums being the discrete approximation to Eq. (A.49) in the
appendix. Note that the spread in growth rates has resulted in an increase
in variance of the simulated size distribution so that it more closely
approximates the observed distribution. Since there is little available
information on the true spread in growth rates within a cohort, these
results can only be viewed as indicative.
The model was next run with plant mortality set equal to zero.
Estimated total numbers differed by less than 1.0 percent from the baseline
case for every age class and length distributions were identical in appear-
ance for the two cases. These results indicate that the effect of the plant
on a gizzard shad cohort is negligible.
44
-------
TABLE 7. LENGTH FREQUENCY OF GIZZARD SHAD IN COVE ROTENONE SAMPLES BY AGE CLASS
COMPARED WITH THE RESULTS OF SIMULATIONS RUN UNDER FOUR TEST CASES
Length Group
25-49 50-74 75-99 100-124
COVE ROTENONE SAMPLES
Age 0 249.550 6017.150 2365.500 704.600
Age I 70.50
Age II
Age III
Age IV
BASELINE CONDITIONS
INITIAL 249.550 6017.150 2365.500 704.600
Age I
Age II
Age III
Age IV
ZERO PLANT MORTALITY
INITIAL 249.550 6017.150 2365.500 704.600
Age I
Age II
Age III
Age IV
TEN-FOLD INCREASE
INITIAL 249.55 6017.15 2365.50 704.60
Age I
Age II
Age III
Age IV
ONE -HUNDRED -FOLD INCREASE
INITIAL 249.550 6017.150 2365.500 704.600
Age I
Age II
Age III
Age IV
125-149 150-174
700.350 140.580
1006.750 851.500
14.910
700.350 140.580
3014.121
700.350 104.580
3021.291
704.35 104.58
2953.860
704.350 104.580
2271.787
-------
TABLE 7. (Continued)
Length Group
175-199 200-224 225-249 250-274 275-299 300-324
COVE ROTENONE SAMPLES
Age 0
Age I 971.380 485.580 14.180
Age II 220.770 282.160 60.570 2.120 0.180
Age III 6.560 5.150 3.180 0.520 0.100
Age IV 0.100
BASELINE CONDITIONS
INITIAL
Age I 590.572 11.300
Age II 198.835 1.538
Age III 4.203 0.010
Age IV 0.095
ZERO PLANT MORTALITY
INITIAL
Age I 591.977 11.327
Age II 199.782 1.546
Age III 4.233 0.010
Age IV 0.096
TEN-FOLD INCREASE
INITIAL
Age I 578.765 11.074
Age II 190.964 1.477
Age III 3.956 0.010
Age IV 0.087
ONE-HUNDRED-FOLD INCREASE
INITIAL
Age I 445.122 8.517
Age II 112.955 0.874
Age III 1.800 0.004
Age IV 0.031
-------
10,000-
1,000-
100
LLJ
CC
!^ 1O-
OC
LLJ
O- 1
CC
LLJ
CO
0.1-
0.01-
-ACTUAL
-BASELINE AND ZERO
PLANT MORTALITY
-•TEN FOLD INCREASE
-ONE HUNDERED FOLD INCREASE
II
III
IV
AGE CLASS
Figure 24. Comparison of age structure in cove rotenone samples
with results of simulations run under four test cases
(simulated populations under baseline and zero plant
mortality situations are represented by the same line)
47
-------
ioc
1
AGE IV
100
1
AGE III
AGE II
ACTUAL
SIMULATED
INITIAL POPULATION
100
150
200
250
300
TOTAL LENGTH, MM
Figure 25. Length frequencies simulated using baseline conditions
compared with cove rotenone estimates.
48
-------
LJ
tr
o
UJ
I
a:
LJ
a_
a:
UJ
m
100
1
100
1
10,000
100
1
ORNL-DWG 80-9342
I I
AGE IY H
AGE III _
' ACTUAL
_ SIMULATED
10,000
100
10,000
100
1
AGE I
INITIAL
POPULATION
I ~
0 50 100 150 200 250
TOTAL LENGTH (mm)
300 350
Figure 26. Length frequencies of the simulation in which the cohort
is divided into nine subcohorts with different growth
rates using baseline conditions, compared with cove
rotenone estimates.
-------
Next, 10- and 100-fold plant mortality rates increases were assumed
in the model. The 10-fold increase reduced total cohort numbers by less
than 10 percent in every age class. The 100-fold plant mortality increase
resulted in significant reductions in total numbers, approximately 65
percent in age class IV, which was most affected.
50
-------
BIBLIOGRAPHY
Anderson, R. 0. 1973. Application of Theory and Research to Management of
Warm Water Fish Populations. Trans. Amer. Fish. Soc. 102:164-170.
Bodola, A. 1965. Life History of the Gizzard Shad, Dorosoma cepedianum
(Le Sueur), in Western Lake Erie. Fish Bull., U.S. Fish and Wldlf. Svc.
Washington, DC. ''
Commonwealth Edison Company. 1977. Impingement and Entrainment Investigation.
Section 6. In: Annual report for fiscal year 1975, Lake Sangchris
Project. IL. Nat. Hist. Surv., Urbana, IL.
Courant, R. and D. Hilbert. 1966. Methods of Mathematical Physics Vol II
Wiley, NY.
DeAngelis, D. L., and J. S. Mattice. 1979. Implications of a Partial
Differential Equation. Math. Biosc. 47:271-285.
Draper, N. R. and H. Smith. 1966. Applied Regression Analysis. John Wiley
and Sons, Inc., NY. 407 pp.
Duke Power Company. Undated. Fish Impingement Study, Allen Steam Station.
Duke Power Company.
Griffel, D. H. 1976. Age-dependent Population Growth, J. Inst. Math Appl
17:141-152.
Hildebrand, F. G. 1962. Advanced Calculus for Applications, Prentice-Hall
Englewood Cliffs, NJ.
Langhaar, H. L. 1972. General Population Theory in the Age-time Continuum,
J. Franklin Inst. 293:199-214.
McDonough, T. A. and P. A. Hackney. 1978. An Analysis of Factors Affecting
Impingement at Cumberland Steam Plant. Proc. Int. Symp. in Env. Effects
of Hydraulic Eng. Works. Knoxville, TN. pp. 49-57.
McDonough, T. A. and P. A. Hackney. 1979. Relationship of Threadfin Shad
Density and Size Structure to Impingement at a Steam-electric Plant.
Proc. Annu. Conf. Southeast. Assoc. Game and Fish Comm. In Press.
Michaelson, S. M. 1970. Dynamics of Balanced and Unbalanced Bass-bluegill
Populations in Ponds in Boone County, Missouri. M.A. Thesis, University
of Missouri. 67 pp.
Moseley, F. N., B. J. Copeland, L. S. Murray, T. S. Jinnette, and P. T. Price.
1975. Further Studies on the Effects of Power Plant Operation on Cox Bay,
Texas. Central Power and Light Company, Corpus Christi, Texas. 153 pp.
51
-------
Oldfield, D. G. 1966. A Continuity Equation for Cell Populations, Bull.,
Math. Biophys. 28:545-554.
Oster, G. and Y. Takahashi. 1974. Models for Age-specific Interactions in a
Periodic Environment, Ecol. Monogr. 44:483-501.
Ricker, W. E. 1975. Computation and Interpretation of Biological Statistics
of Fish Populations. Bull Fish. Res. Board. Can. 119. 300 pp.
Rotenberg, M. 1975. Equilibrium and Stability in Populations Whose
Interactions Are Age-specific, J. Theoret. Biol. 54:207-224.
Rubinow, S. I. 1973. Mathematical Problems in the Biological Sciences, SIAM,
Philadelphia.
Sinko, J. W. and W. Streifer. 1967. A New Model for Age-size Structure for a
Population, Ecology 48:910-918.
Sinko, J. W. and W. Streifer. 1971. A Model for Populations Reproducing by
Fission, Ecology. 52:331-335.
Tennessee Valley Authority. 1977a. 316 (a) and (b) Demonstration - Cumberland
Steam Plant. Vol. 5: Effects of the Cumberland Steam Plant Cooling Water
Intake on the Fish Populations of Barkley Reservoir. Tennessee Valley
Authority. Morris, Tennessee. 84 p.
Tennessee Valley Authority. 1977b. 316 (a) and (b) Demonstration Cumberland
Steam Plant. Vol. 4: Effects of Thermal Discharges from Cumberland
Steam Plant on the Fish Populations of Barkley Reservoir. 234 pp.
Tennessee Valley Authority. 1978. Preoperational Fisheries Report for the
Sequoyah Nuclear Plant. Tennessee Valley Authority. Norris, TN. 179 pp.
Trucco, E. 1965. Mathematical Models for Cellular Systems. The Von Foerster
Equation. Part I, Bull. Math. Biophys. 27:285-304.
Trucco, E. 1965. Mathematical Models for Cellular Systems. The Von Foerster
Equation. Part II, Bull. Math. Biophys. 27:385-304.
Van Sickle, J. 1977. Analysis of a Distributed-Parameter Population Model
Based on Physiological Age, J. Theoret. Biol. 64:571-586.
Von Foerster, H. 1959. Some Remarks on Changing Populations. In: The
Kinetics of Cellular Proliferation, F. Stohlman, Jr., Ed., Grune and
Stratton, NY. pp. 382-407.
Weiss, G. H. 1968. Equations for the Age Structure of Growing Populations,
Bull. Math. Biophys. 30:427-435.
52
-------
APPENDIX A: ANALYSIS OF MODEL
53
-------
1. Solution of the partial differential equation for time-invariant
growth and mortality rates
Solution of Eq. (11) is obtainable by the method of characteristics.
See, for example, Hildebrand (1962), Courant and Hilbert (1966), or Van
Sickle (1977). To apply this method, it is easiest to first make the further
assumptions that G(s,t) = G(s) and Z(s,t) = Z(s) and that all recently
hatched fish have the same size s = s». Then B(s,t) can be written
B(s,t) = B0(t)6(s-s0) , (A.I)
where 6(s-s,.) is the Dirac delta-function, defined by the properties
6(s-S()) = 0 for sfsQ (A.2)
(b
6(s-s0)f(s)ds = f(sQ) for (sa < S() < sb) , (A.3)
'Sa
where f(s) is an arbitrary function. The Dirac delta-function may be visu-
alized as a very sharp spike at s=sn, having unit area.
The replacement of B(s,t) by B0(t)6(s-sQ) reduces Eq. (11) to
^ + =| {G(s)N(s,t)} = - Z(s)N(s,t) + B.(t)6(s-sn) (A.4)
w OS U U
or
3N(s.t) . „, x 9N(s,t) ,„, , ^ aG(s),iT, ..>_..,, ^Mf s /-. CN
9t' + G(s) gs? = ~{Z(s) + -g^-}N(s,t) + B0(t)6(s-sQ) (A.5)
To determine the boundary condition, integrate both sides of Eq (A.4)
over a very small region surrounding s=s0;
{G(s)N(s,t)} + Z(s)N(s,t)
e
S0 - 2
- B0(t)6(s-s0)J ds = 0 . (A.6)
Thus, taking each term in (A.6), the following results are obtained:
54
-------
£
S0 - 2
s
S0 " 2
3N
3N(s t)
ds £ £ jrr—^ = 0
at 6 o
(G(s)N(s,t)} ds = G(s)N(s,t)
* G(s )N(s ,
(A.7)
(A.8)
S0 " 2
Z(s)N(s,t) ds = £Z(sQ)N(s0,5) -> 0
t* \J
(A.9)
£
S0 " 2
B0(t)6(s-s0) ds = BQ(t)
(A.10)
The result (A.8) follows since it is assumed that N(s,t) = 0 for all s < s~.
See Eq. (A.3) for the derivation of (A.10). Therefore, using these results
in Eq. (A.6), the boundary conditions on N(s,t) are found to be
N(s0,t) = B0(t)/G(S()). (A.11)
The theory of partial differential equations (e.g., Hildebrand 1962)
implies that the general solution of Eq. (A.5) is of the form F(u.,u2) = 0,
where F(u-,u2) is any function relating u^ and u , and where u (N,s,t) = C.
and u0(N,s,t) = C% are solutions of any two differential equations that imply
dt _ ds
1 " G(s) ' JZ^T
-dN
dG(sTTN '
ds
(A.12)
where C. and C» are constants of integration.
In the present case, in which boundary conditions are posed on the
s = SQ plane, it is convenient to choose to integrate the equations formed
by terms 1 and 2 and by terms 2 and 3. Performing the integrations, there
result
55
-------
u.CN.s.t,) = t - I(sn,s) = C.
(A.13)
and
where
u2(N,s,t) = ln(N) + In {G(s)/G(sQ)} + J(sQ,s) = ln(C2):
(A.14)
KSO,S) =
ds'
G(s')
(A.15)
and
J(sQ,s) =
Z(s') ds'
G(s')
bo
Equations (A.13) and (A.14) can be expressed as
t = I(s0,s) + Cj
and
G(sn) -J(sn,s)
\r — ft r " •> "
(A.16)
(A.17)
(A.18)
These two equations each define a family of integral surfaces in (N,s,t)-
space (Fig. A.I). The intersection of a given pair of surfaces from these
families defines a curve called the characteristic curve.
What is now necessary is to specify the particular solution in the
relation F(u ,u2) = 0, or, equivalently, F(C ,€„) = 0. This relation must
be such that the boundary value condition, Eq. (A.11), is satisfied. Choose
as F(C1,C-) the equation
C2 = lyCjJ/GCs^ , (A. 19)
or, from Eq. (A.17),
C2 = BQ{t - I(s0,s)}/G(s0) . (A.20)
It can be seen from Eq. (A.18) that, using this assumption, the particular
solution is
N(s,t) =
Bn{t - l(sn,s)}e
-J(sQ,s)
Q 0
(A.21)
56
-------
ORNL-DWG 80-9042
\
Figure A.I. The integral surfaces defined by Eqs. (1-17) and (A.18).
The intersection of these surfaces is the characteristic
curve.
57
-------
ORNL-OWG 80-9046
N
Figure A.2. The surface N(s,t) as defined by Eq. (A.21). It is composed
of characteristic curves.
58
-------
(Fig. A.2). The surface defined by Eq. (A.21) is composed of characteristic
curves satisfying Eq. (A.19). It can be seen that for s = sn, N(s,t) =
B0(t)/G(sQ). °
The simplicity of the above analytic solution hinges on the time
invariance of G(s) and Z(s) and on the approximation (A.I) for B(s,t).
Cases where these assumptions need not hold will be considered later.
2. Particular solutions of the equation
Any attempt to model empirical data requires that specific forms of
G(s), Z(s), and Bfl(t) be used. Fortunately, G(s) can often be approximated
by very simple functions. Two cases are given below.
Case 1
The early growth of many fish is approximately exponential. Hence,
G(s) = 8()s , (A.22)
where gn is a constant. If the mortality rate, Z(s), is assumed to be
constant, Z», then I(s_,s) and J(s0,s) become
s
ds
;rr = — ln(s/S())
(A.23)
J(sQ,s) =
ds'Z
0
80S
'0
Equation (A. 21) then becomes
4 in(s/so}
(A.24)
N(s,t) =
Bn{t -
)ln(s/s
B{t -
-(Z0/g0)
(A.25)
Some discussion of the meaning of this expression may be helpful at
this point. First, to interpret BQ{t - (l/g0)ln(s/sQ)}, recall that BQ(t)
59
-------
is reproduction as a function of time. For most fish species of interest,
spawning is confined to a relatively short period of time during the year,
perhaps a few weeks. Denote by TO = 0 the beginning of the spawning period
and the end by T . Over this period, B(t) will have some distribution, which
s
can perhaps be very crudely approximated by a truncated normal distribution
(Fig. A. 3).
The expression BQ{t - (l/g0)ln(s/sQ)} in Eq. (A. 25), however, is a
function of t - (l/g0)ln(s/s0) , and hence of both t and s. Consider this
expression from two points of view. First, consider B~{t - (l/g,.)ln(s/s-) }
as a function of time for various values, s., of s (s. ^ s,,). Several such
curves are plotted in Fig. A. 4. Note that as s. increases, the curve moves
to the right. If t. signifies the time when cohort members first reach
size s., then note that as s. increases, the time interval between successive
i' i '
values of times, t. and t. , decreases since growth is accelerating.
One can also consider the expression Bft{t - (l/g..)ln(s/s0)} from the
standpoint of variation with respect to s at fixed times, t , t~ , ..., t
(Fig. A. 5). Each curve is the size distribution for the particular asso-
ciated time, t.. Note that the size distribution spreads through time
since the ratio of the upper to lower limits, s /s
UT)
, ,
.LOW
is
'low
(A.26)
In Fig. A.5, Bn{t - (1/g )ln(s/sn)} has been divided by gns, as in Eq. (A.25),
0
0
(T
so that in the absence of mortality (Z~ = 0), population number is conserved;
i.e., it can be shown by a change of variable that
up
Bn{t - (l/gn)ln(s/sn)}
V
ds =
BQ(x) dx
(A.27)
'low
0
for all values of t. The factor (s/s,J
-(Z0/g0) _
in Eq. (A.25) is a quantity
between 0 and 1.0, representing the fraction of surviving fish.
Case 2
Growth characteristics for individual organisms of many species have
the appearance of logistic curves when examined over the total life spans;
60
-------
ORNL-DWG 80-9045
t
B(t)
0
>s
t, TIME
Figure A.3. A hypothetical recruitment rate, B(t), between the times
t=0 and t=T .
s
61
-------
ORNL-DWG 80-9043
1200
1000
800
£ 600
400
200
S= 10
\
S = 30 S = 50
T
= 70
1
10 20 30 40
t, TIME
50
60
Figure A.4.
Plots of B_{t - (l/g-)ln(s/s-} as functions of time, t, for several values of size,
s, in units of millimeters. The reproduction function, B/>(t) is a truncated normal
(Figure A.3).
-------
ORNL-DWG 80-9044
OJ
(A
1200
1000
800
600
400
200
t= 20
t = 35
I I
I I
t = 50
I I
20 40
60 80
s, SIZE
100 120 140
Figure A.5. BQ{t " (l/gQ)ln(s/s())}/8Qs) as function of size, s, for three values of time, t.
Tfie reproductive function, BQ(t), is a truncated normal (Fig. A.3).
-------
some examples are given by DeAngelis and Mattice (1979). Early growth in
the life span is approximately exponential, but later slows down and even-
tually plateaus. Growth rate in these cases can be approximated by the
function
G(s) =
max
(A.28)
where s is a constant representing the upper limit of size. Again, if
IT13X
Z(s) = Z_, then I(s_,s) and J(sn,s) are easily computed using the method
of partial fractions. Note that
1 1 .1
g
ov
max
)g
- s)
(A.29)
so that
. s
ds'
s - s s
max
(A.30)
= { - r- ln(s_ - s') + i- ln(s')}
max
s(s
= — In
~
max 0
80 " S0(smax - s)
Then
s(s - sn)
0
(A.31)
N(s,t) =
s
max
gn(s - s)s
sn(s
0 max
s(s - i
max
s)
5o}_
B{t - I(s0s)} (A.32)
The size distribution, N(s,t) exhibits interesting behavior. Assume
that B(t) has a truncated normal distribution,
64
-------
B(t) =
0 t < 0
b exp -(t - 0.5T )2/b^ 0 < t < T (A.33)
u si s
0 T < t
s
where b- and b are constants. In Fig. A.6 the size distribution of the
cohort is computed for several values of t and for a set of hypothetical
parameter values. In this example, Z_ is assumed to be zero for all values
of s, so N(s,t) is given by Eq. (A.32) with Z- = 0. Time is measured in
years from initial production of the cohort. Note that unlike the preceding
example, the initial broadening of the size distribution is followed by a
narrowing through time (Fig. A.7). The mean size follows a logistic curve,
as expected. The standard deviation increases at first, then asymptotes and
finally decreases as the sizes of individuals in the cohort approach their
upper limit.
The variations in standard deviation are somewhat unexpected but can
easily be interpreted. In the initial phase of growth, when size is
2 2
accelerating (d s/dt > 0), the largest fish (those produced earliest)
always grow at faster rates than smaller individuals, leading to a broaden-
2 2
ing of the size distribution. Later, when d s/dt < 0, the smallest fish
grow faster than the larger ones and tend to catch up in size, causing the
size distribution to become narrower.
A rough mathematical interpretation for the occurrence of the peak in
the standard deviation curves close to the point of the maximum growth rate
can be made. The size distribution should be broadest in the region in which
the argument of B, t - I(s.,s), changes most slowly with respect to s. For
an increment in I(s ,s), AI(sQ,s) = As/G(s), the argument is least sensitive
when G(s) attains its maximum value.
3. Generalization to time-varying coefficients
At best, Eqs. (A.22) and (A.28) may be reasonable approximations of the
growth rate under special conditions. They will certainly not he good
approximations when environmental parameters such as temperature and food
availability change significantly over time scales of interest. Seasonal
changes in temperature greatly influence growth rates. Hence, we can hardly
expect the assumptions that G(s,t) = G(s) and Z(s,t) = Z(s) to be valid in
real cases, so solutions developed in the preceding two sections are not
sufficient.
Despite the insufficiency of Eq. (A.21) for describing cohort dynamics
when mortality and growth rate are functions of time, it can be extended to
approximate these cases. If Z(s,t) and G(s,t) change relatively slowly with
respect to time, say on a scale of months, then Eq. (A.21) for N(s,t) may be
valid over short periods of time. Solutions for these short time periods can
then be pieced together.
65
-------
600
500
t 400
ORNL-DWG 78-22347R
LU
Q
CC
UJ
CD
300
200
100
0
0
= 4.0
= 5.0
= 8.0
A
40
80 120
LENGTH (mm)
160
Figure A.6. N(s,t) from Eq. (A.32) as a function of size, s, for several values of time, t, and
for arbitrarily chosen parameter values.
-------
ORNL-DWG 78-22354
E
E
e>
Ld
160
120
80
40
4 6
TIME (years)
8
10
8 E
Ld
Q
a
or
|
I
CO
Figure A. 7.
The mean size in the cohort (black dots) and the standard deviation in size (white
dots) as functions of time from Eq. (A. 32). The parameter values have been chosen
arbitrarily.
-------
Assume that the spawning period has ended and that N(s,t.) has been
computed for time t. from Eq. (A.21). Assume that time t. represents the end
of one month and that during the next month G(s,t) and Z(s,t) are different
from their original values. The equation for N(s,t) during the next month is
}- + Z(s)} N(s,t)
+ N0(s)6(t-t.) (A.34)
where N_(s) is the size distribution at the end of the first month.
One can write down the solution for N(s,t) using the method of char-
acteristics. However, it is useful to go through an explicit, detailed
solution for N(s,t), since this will help keep better track of the
mathematical steps.
Introduce the variable t, defined by
T = t + I(s) (A.35)
where
I(s) = J ds|
This substitution permits the left hand side of Eq. (A.34) to be written as
a total derivative with respect to t since
dN _ 3t 3N(s,t) 3£ aN(s,t) _ 3N(s,t) , , 3N(s,t) ,
dr ~ at at at as at ^*J as • ^
Therefore, Eq. (A.34) can now be written as
dN = _{dG(sl + z(g)} N(g>t) + N (s)6(t_t ). (/
Q L Q S v X
Equation (A.36) can be solved to give
N(s,t) = e"R(T) I eR(Tl)N0(s)6(t - t±) dl' (A.38)
where
} df . (A. 39)
R(t) = (Z(s) +
68
-------
In Eqs. (A.38) and (A.39), s and t must be expressed in terms of t using
Eq. (A.35). In general, an explicit analytic expression for s in terms of
T is not possible since I(s) will not always be integrable analytically.
Therefore, it may be necessary to obtain s in terms of t numerically from
Eq. (A.35), which we call s = F(t - t). Equations (A.38) and (A.39) can
be reexpressed as
N(s,t) =
NQF(T' - t)6{T' - I(s) -
dt' (A.40)
and
R(t) =
t
dt1 .
(A.41)
Equation (A. 40) integrates to
N(s,t) =
x
G{F(t
- t)} -Q{t.
, t
- t)}
(A.42)
where
t + I(s)
JZ{F(T' - t)}ldt'
(A.43)
t. + Ks)
Special case
Consider the special case where
G(s) = v(l - ) s and
max
Z(s) = Z
0
In this case, I(s) has an analytic form,
v s - s
max
as does F{I(s) + t. - t} ,
69
-------
- t}
e -.-> + Vtj
s smax . (A.44)
max
R{t±+ I(s), t + I(s)} = - (Z0+ v)(t - t ) + 2 In { 1
-v(t - t.)
s - s
max
max
- 2
= - (v + ZQ)(t - t.) - 2 In {- - i - } . (A. 45)
max
-v(t - t.)
(s - s) + s e
Finally, N(s,t) is
-Up + v)(t - t.) f Snax "I
N(S>t) = C I (s -.)*..-""'I']
v max ' '
XN° f~77
Is + (
\ VI t - t . J
s - s) e v i
max
From Equation (A.42), N(s,t) can be calculated for all times in month i,
up to time t = t. .. Then the size distribution, N(s,t.+1) can be used
as the new initial value function, N. -(s), from which N(s,t) can be
calculated for the succeeding month by the process outlined above.
4. Size-dependent mortality
The general formula for N(s,t) given by Eq. (A.21) and its extension
to situations of time-varying mortality and growth rate in the preceding
section should enable consideration of realistic populations. Of special
interest in this report will be the influence of various types of mortality
on cohort population number and length distribution.
It is useful to distinguish among the various types of mortality
expected to occur; for example, between fishing mortality, Zf(s,t), and
natural mortality, Z (s,t). Natural mortality can itself be divided into
different classes. Suppose the effects of a particular predator on the
species being modeled are of interest. It will be useful to separate the
70
-------
mortality, Z ..(s,t), caused by this predator from the remaining natural
mortality, Z
What makes Zf(s,t) and Z -(s,t) of interest is that these mortality
sources are likely to act on restricted size class ranges. In particular,
Z j(s,t) will depend on the size ratios of predator to prey.
Assume, for example, that Z .(s,t) has the form
Z ,(s,t) = C ,(s) (A.47)
nl pred
where
Smax,p
s
c
where ^ fi,(s ,t) is the number of predators of size s , s is the maximum
size at which a predator can devour a prey of length s, s the maximum
max,p
size of predators in the population, and C(s /s) the likelihood of a
predator of size s devouring a prey of size s in an encounter.
5. Distribution of parameter values
Another source of potential inaccuracy in the model is the assumption
of certain types of uniformity in the cohort populations. For example, in
the particular case for which N(s,t) is given by Eq. (A.32), the parameters
sr\> ^m<.«» an<* 8n are treated as constants having the same values for every
U IDciX U
member of the population. This is an approximation, for certainly these
parameters will differ for individual organisms due to both genetic and
environmental variation. This variation can be taken into account in our
model at the cost of some complication. Suppose for example that newly
produced organisms do not all have the same growth rate coefficient, g-,
but that the sizes have some general distribution, F(gQ). The size
distribution, N(s,t), through time is given by
N(s,t) =
o"
80
F(g())N(s,t;g0) dg() , (A.49)
where g' and g" are the lower and upper bounds, respectively, on the
distribution of gfl, and N(s,t;gQ) represents N(s,t) for a particular value
of g . If Eq. (A.32) is the solution for N(s,t), and distribution functions
are available for s and SA as well as for gn, then one can generalize
max U u
71
-------
Eq. (A.49) to a triple integral over the distribution functions for SA,
s , and grt.
max' 60
In general, one must use numerical methods to evaluate the integral
(A.49). However, in some special cases the integral can be performed
analytically, so that the effects of variation in a parameter value within
the cohort can be investigated in more detail. For example, suppose all
reproduction takes place close to a single instant, t=0, so that BQ(t) can
be approximated by Bft_6(t). Assuming N(s,t) is given by Eq. (A.32) and
also assuming that the parameters SA, s , and gA do not (at this point)
U IftHX U
vary among the members of the cohort, then from
N(s,t) =
max
V max
- s)s
SA(S - s)
0 max
s(s - SA)
max 0
V«0
B0()6{t - I(s0,s)} , (A.50)
with I(s ,s) given by Eq. (A.28), the sizes of all surviving fish in the
cohort will the the same at all times and be described by a logistic function.
The total number of fish, N(t), at any time t is
N(t) =
N(s,t) ds = BnAe
(A.51)
where s' and s" are the upper and lower limits on sizes in the cohort at
time t. This can be proven by substituting N(s,t) given by Eq. (A.50) into
the integral of Eq. (A.51) and integrating, after making the change of
variable from s to u, where
u = t - I(s0,s) = t - (l/gQ)ln
Eliminating ds in favor of du by means of
s(s - SA)
max 0
s,.(s - s)
Ov max
ds =
g0(u-t)
grtsrts (s - s_)e j
60 0 max max 0 du ,
' S0) + S0e
g0(u-t)
(A.52)
(A.53)
and then using
s =
s s_e
max 0
80(u-t)
(A.54)
72
-------
s - s =
max
s (s SA)
max max 0
(A.55)
(s - s.) + s..e
max 0 0
g (u-t)
and
S0(smax " S)
S(smax ' S0)
=
(A.56)
the integral (A.51) becomes
u'
ZQ(u-t) -Z t
due0 B006(u)=B00e
u'
(A.57)
Now let one of the parameters, say g-, have a normal distribution,
F(gQ), about the mean, gfl;
F(g0) = - l-ur exp { -(g - g~0)2/2p2} .
° 1/2 ° °
(A. 58)
The size distribution, N(s,t), can be calculated from
N(s,t) =
00
H.(s,t;g0)F(g0)
(A.59)
Substituting N(s,t) from (A.50) and F(gQ) from (A.58) into (A.59), and
making a change of variables from g» to u using
u = t - — ln(D) (A.60)
0
dgfl
du = ;f ln(D)
«0
(A.61)
where
D =
s(s - s_)
max 0
sn(s - s)
0 max
(A.62)
it follows that
73
-------
N(s,t) = - B°°Smax
(27t)1/2p(s_ - s)s
max
-Z (t - u)/ln(D)
ln(D) - 8()(t - u)}2
n " i-'-L^~"--
- u)'
D\J 4.1J \ L. U; ft \ .
e r 5(u) du
u - t
-00
(A.63)
or, using
-Z t/ln(D) -Znt
D ° = e ° , (A.64)
-{ln(D) - g()t}2/(2p2t2) -ZQt
B00Smaxe _ e . (A. 65)
In a similar manner, it is possible to compute the size distribution,
N(s,t), given parameter distributions in either s or sn.
ffldX U
6. Density-dependent growth rate
It is reasonable to expect that in sufficiently crowded populations
the growth rate of individuals may be reduced because of food resource
limitations. There are very few empirical data available on which to base
models of density-dependent growth. Nonetheless, a simple relationship
between the actual growth rate, G'(s,t), the growth rate under uncrowded
conditions, G(s,t), and the total population size, N-(t) can be postulated;
G'(s,t) = G(s,t)/{l.O + pNQ(t)} , (A. 66)
where p is a constant coefficient.
74
-------
APPENDIX B: USE OF THE COMPUTER PROGRAM
75
-------
1. General information on the program
The purpose of the computer program is to solve the basic partial
differential equation, Eq. (11), for the size distribution, N(s,t), through
time. The growth rate, G(s,t), and the mortality rate, Z(s,t), depend on
the size (e.g., length), s, and on the time of the year, t. The reproduc-
tion rate, B(s,t) = B(t), can be allowed to vary during the spawning
period. The numerical formulas, Eqs. (A.21) and (A.42), are used to compute
N(s,t). Note that GA(J), ZA(J), and BA(J) in the computer program correspond
to discrete arrays of the mathematical variables G(s,t), Z(s,t), and B(t),
respectively, in the text. To obtain G(s,t), Z(s,t), and B(t) from GA(J),
ZA(J), and BA(J), interpolation is used; e.g.,
G(s,t) = GA(J) + (GA(J+1) + GA(J))*((s - SIZE(J))/(SIZE(J+1) - SIZE(J))),
where SIZE(J) < s < SIZE(J+1), and the array, GA(J), is assumed to represent
the growth rate at size SIZE(J) at time t. (The array GA(J) can be reset
at specified times as described below.)
Some simplifying assumptions are built into the program in its present
form to keep it from becoming too complex. One of the main simplifications
is that the temporal variations in GA(J) and ZA(J) are assumed to be dis-
continuous rather than smooth. For example, GA(J) for a particular month
(or whatever time period is chosen) will be constant in time during the
whole month (through varying with s, or the computer variable J), but will
change for the next month. Secondly, all new recruits to the cohort are
assumed to have the same length, SQ. These assumptions can be relaxed,
but at the cost of appreciable complication in the computer program.
Below, the present section outlines the general operation of the
computer program, section 2 describes the setup of the input data cards,
section 3 is an example of a specific application of the program, and
Appendix C is a printout of the program.
The computer program is divided into a MAIN PROGRAM and a subroutine,
SUBROUTINE TRAP. The purposes and operation of each of these is discussed
in some detail below.
Main Program
The first task of the MAIN PROGRAM is to read in all initial input
data and print these out. Only data on later temporal changes in the
growth and mortality rates are read in later in the program.
The principal DO-loop in the MAIN PROGRAM is DO-loop 900. In general
it is assumed that the cohort can be divided into several subcohorts, each
with a different growth rate, G(s,t). The DO-loop 900 sums over all of
the subcohorts.
The simulation of each subcohort size distribution through time is
performed by DO-loop 700 in the MAIN PROGRAM. Within this DO-loop are
two main subloops. The first of these, DO-loop 100, uses Eq. (A.21) to
calculate the part of N(s,t) resulting from reproduction during the
76
-------
preceding time interval. Within DO-loop 100, a call is made to SUBROUTINE
TRAP to calculate I(s0,s) and J(sQ,s). The second loop, DO-loop 500,
calculates N(s,t) resulting from an initial size distribution at the
beginning of the preceding time period using Eq. (A.42). Within DO-loop 500
is a main subloop, DO-loop 450, which calculates R{t. + I(s), t + I(s)}
(see Eq. A.41).
Values of the size distribution, N(s,t), through time are stored in
the array SDIST(I,J), and then printed out near the end of the MAIN PROGRAM.
Subroutine Trap
SUBROUTINE TRAP is called from DO-loop 100 of the MAIN PROGRAM. It
calculates I(sft,s) and J(sft,s) for every value of s. The trapezoidal
method is used in evaluating the integrals. These integrals are used in
DO-loop 100 to evaluate Eq. (A-21). The values of I(s~,s) from s0 to s
are stored and then later inverted numerically in DO-loop 500 to obtain s
in terms of values of I(s ,s). This yields the function s = F{l(sQ,s)}.
Later in DO-loop 500, the values t. - t + I(s~,s) are used as the argument
of F, to give ?{t± - t + I(SG,S)} (see Eq. A.42). Then a linear interpola-
tion technique is used to compute N.(F{t. - t + KS.S,)}), which is necessary
in evaluating Eq. (A.42).
2. Input data
Card A
Input parameters: NSIZES, NSIZEC, NBIRTH, NCHNGE, NRUNTM, Nil, NIII, NGROW
Format: 815
NSIZES = Number of size classes into which the cohort is divided.
NSIZEC = Number of sizes at which input data on size-dependent mortality
and growth rates are given.
NBIRTH = Number of time intervals during the spawning period at which
the instantaneous numbers spawned per month are given. If there
is no reproduction, set NBIRTH = 0.
NCHNGE = Number of times during the projected course of the simulation at
which the mortality and growth rates change. Since these
quantities are usually assumed to change month-to-month, NCHNGE
is a measure of the number of months in the projected simulation.
NRUNTM = Number of time steps in the simulation.
Nil, = These integers control the printing out of detailed information
NIII on the computations, which may be useful for diagnostic purposes.
For example, if we set Nil = 6, certain computational details
77
-------
will be printed out every time the index I in DO-loop 700 is a
multiple of six. For other values of I, the WRITE statements
are skipped. The other integer, NIII, controls other WRITE
statements. To understand and utilize these computational
printouts, the program users will have to familiarize themselves
with the details of the program.
NGROW = Number of subcohorts having different growth rates.
Card B
Input parameter: TO
Format: 2E10.0
TO = Time of the beginning of the simulation
Card C
Input parameters: DELSIZ, DELCLA, DELSZA
Format: 4E10.0
DELSIZ = Length of size classes into which the cohort is divided. It
seems best to make DELSIZ as small as practical for more
accuracy.
DELCLA = Length of size intervals between which mortality and growth
rate data as functions of size are given.
DELSZA = Length of desired size-class printout; obtained by summing
over groups of size classes of width DELSIZ. For example,
DELSIZ may be equal to 0.5 millimeters, but if the user wants
to print out results of N(s,t) for 5.0 millimeter, then the
user must set DELSZA =5.0.
Card D
Input parameter: SIZE(l)
Format: E10.0
SIZE(l) = Size at time of reproduction of cohort members.
Card E
Input parameters: ZNA(I), 1=1, NSIZEC
Format: 7E10.0
ZNA(I) = Size-specific natural mortality for size class I.
78
-------
Cards F
Input parameters: ZIA(I), I=1,NSIZEC
Format: 7E10.0
ZIA(I) = Size-specific impingement mortality for size class I.
Cards G
Input parameters: GA(I), I=1,NSIZEC
Format: 7E10.0
GA(1) = Size-specific growth rate for size class I.
Cards H
Input parameters: BA(I), I=1,NBIRTH
Format: 7E10.0
BA(I) = Instantaneous reproduction rate (numbers per month, for example)
through time during the spawning period. If there is no
reproduction, leave one blank card here.
Cards I
Input parameters: TBRT(I), I=1,NBIRTH
Format: 7E10.0
TBERT(I) = Times during the spawning period at which numbers spawned
per unit time (e.g., month) are given. If there is no
reproduction, leave one blank card here.
Cards J
Input parameters: TIMEA(I), I=2,NCHNGE
Format: 7E10.0
TIMEA(I) = Times during the projected simulation run at which the
growth and mortality rates are changed. TIMEA(l) is
automatically set to 0.0.
Cards K
Input parameters: TIMEPL(I), I=1,NRUNTM
Format: 7E10.0
79
-------
TIMEPL(I) = Times at which the size distribution is computed and printed
out. It is necessary to at least have one value of TIMEPL(I)
paired to each value of TIMEA(I). Actually it is best to
set TIMEPL(I) to a value very slightly smaller (say a fraction
of a day) than the paired value of TIMEA(I).
Card L
Input parameter: DENDPA
Format: E10.0
DENDPA = Coefficient of effects of total population number on the growth
rate (p in Eq. 75).
Cards M
Input parameters: FRACA(I), I=1,NGROW
Format: 7E10.0
FRACA(I) = Constant setting growth rate of subcohort I. GA(J) = FRACA(I)*
GAA(J) is the growth rate of subcohort I at size J.
Cards N
Input parameters: FRACB(I), I=1,NGROW
Format: 7E10.0
FRACB(I) = Fraction of the total cohort in subcohort I.
Card 0
Input parameters: NINIT, NUINT
Format: 215
NINIT: If NINIT = 1, there is an initial size distribution at time t = TO;
otherwise all the size classes initially have zero population.
NUINT: The first NUINT size classes are assigned initial values greater
than or equal to zero.
Cards P
Input parameters: SNUMIN(I), 1=1,NUINT
Format: 7E10.0
SNUMIN = Initial population (t = TO) at size class I
80
-------
Card Sets Q, R, S
A new sequence of these three sets is read in every time the current time,
T, exceeds the next value of TIMEA(I).
Cards Q
Input parameters: ZNA(I), I=1,NSIZEC
Format: 7E10.0
ZNA(I) = Size-specific natural natural mortality for size class I.
Cards R
Input parameters: ZIA(I), I=1,NSIZEC
Format: 7E10.0
ZIA(I) = Size-dependent impingement mortality for size class I.
Cards S
Input parameters: GA(I), I=1,NSIZEC
Format: 7E10.0
GA(I) = Size-specific growth rate for size class I.
3. Example application of the program
As an example showing how to set up input data for the computer program,
let us consider measurements of larval crappie in an arm of Pickwick Reservoir
on the Tennessee River in 1976 (TVA 1976, Hackney and Webb 1977, DeAngelis
et al. 1979).
The basic information read as input data is shown in Table 1. It
includes the numbers of new recruits (per 1000 cubic meters) into the 5.0
millimeter length class as a function of time in months (Cards H and I).
Data on growth rates and natural mortality rates as functions of length
are also assumed known (Cards E and G), and are assumed to be constant
through time, so NCHNGE = 2 and TIMEA(2) = a very large number. No initial
length distribution is assumed, so NININT = 0 (Card M) and no number
density-dependence of growth rate is assumed, so DENDPA = 0.0 (Card L).
The size of the length classes, DELSIZ, is set at 0.5 millimeters (Card C).
The length distribution simulation is to be printed out 20 times at the
specified values of TIMEPL(I) (Cards K). No density-dependence is assumed
in the growth rate (Card L), and only one subcohort (the whole cohort) is
assumed (Card M and Card N).
The input data are printed out by the program as shown in Table 2.
Selected results are shown in Tables 3, 4, and 5. Table 3 shows the
numbers of larvae per 1000 cubic meters in the first 10 0.5-millimeter
81
-------
length classes (out of the 200 length classes actually printed out by the
program) for all 21 times. Table 4 shows the numbers in the first 10
5.0 millimeter length classes (of the 20 actually printed out by the
program) at all times. Finally, Table 5 shows the total population at all
times.
82
-------
00
200 10 20 3 20 20 20 1 A
B
C
D
5
i?
•.j
F
P
G
G
H
H
H
I
I
I
J
K
K
K
L
0.0
0.5
5.0
2.1
2.1
0.0
0.0
5.45
10.9
0.0
0.0
0.0
0.0
1.^5
3.50
100.
0.0
1.75
3.50
0.0
1.0
1.0
0 0
40.
2.1
2.1
0.0
0.0
54.0
10.9
2.5
0.0
0.0
0.25
2.0
3.75
0.25
2.0
3.75
5.0
2.1
2.1
0.0
0.0
95.0
10.9
3.75
0.0
0.0
0.50
2.25
4.00
0.50
2.25
4.00
2.1
0.0
10.9
10.9
0.0
0.0
0.0
0.75
2.5
4.25
0.75
2.5
4.25
2.1
0.0
10.9
10.9
50.
0.0
0.0
1.0
2.75
4.50
1.0
2.75
4.50
2.1
0.0
10.9
250.
0.0
0.0
1.25
3.0
4.75
1.25
3.0
4.75
2.1
0.0
10.9
187.5
0.0
1.5
3.25
5.00
1.5
3.25
Table B.I. The input data cards relevant to the example. The meaning
of the individual cards is given in section 2.
-------
PARTIAL DIFFERENTIAL EQUATIOH HODEL OF A FISH COHORT SIZE DISTRIBOTIOH
HSIZES= 200 HSI£BC= 10 HBIRTH= 20 HCHHGE= 3
SIZE CUSS
HOSTILITY RATE
GROWTH
HROHTH= 20 SII= 20 HIII= 20
HIHIT= 0
DELSIZ= 0.50000000E 00 DEI.CLA= 0.40000000E 02
DELSZA=
5.0000
0.50000000E 01
0. 45000000E 02
0,85000000E 02
0, 12500000E 03
0. 16500000E 03
0. 20500000E 03
0.24500000E 03
0.28500000E 03
0.32500000E 03
0, 36500000E 03
0.20999994E 01
0. 2099999UE 01
0.2099999HB 01
0.2099999UE 01
0.2099999HE 01
0.20999994E 01
0.2099999UE 01
0.20999994E 01
0.20999994E 01
0.2099999HE 01
0.51HI99998E 01
0.5KOOOOOOE 02
0.95000000E 02
0.10900000E 02
0. 10900000E 02
0. 10900000B 02
0.10900000E 02
0. 10900000E 02
0.10900000E 02
0. 10900000E 02
PiRiHETEB CHAHGES OCC08 AT
TIHEA= 0.0
TISEA= 0.10000000S 03
TIHEA= 0.0
c»
TIHES FOR RHICH SIZE DISTRIBOTIOK IS CALCULATED
TIHE
BIRTH RATE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
0.0
0.25000000E 00
0.50000000E 00
0.75000000F 00
0.100000001! 01
0.12500000E 01
0.15000000E 01
0.17500000F 01
0.20000000B 01
0.22500000? 01
0.25000000? 01
0.27500000^! 01
0.30000000K 01
0.32500000* 01
0.3500000Gk 01
0.3750000CB 01
O.UOOOOOOOB 01
O.U250000CE 01
0.4500000CE 01
O.H7500000E 01
0.0
0. 25000000E 00
0.50000000E 00
0.75000000E 00
0. 10000000E 01
0.12500000E 01
0. 15000000B 01
0.17500000E 01
0.20000000E 01
0. 22500000B 01
0.25000000E 01
0.27500000E 01
0.30000000E 01
0.32500000E 01
0.35000000E 01
0.37500000E 01
0.40000000E 01
0.42500000E 01
O.H5000000E 01
0.47500000E 01
DEUDPA= 0.0
0,0
0.25000000E 01
0.37500000B 01
0.0
0.50000000E 02
0.25000000E 03
0.18750000E 03
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
Table B.2. The format in which the input data is printed out by the
program.
-------
ESTIMATED SIZE DISTBIBOTIOHS
SIZES
Ln
TIHE
0.0
0.2500
0.5000
0.7500
1.0000
1.2500
1.5000
1.7500
2.0000
2.2500
2.5000
2.7500
3.0000
3.2500
3.5000
3.7500
4.0000
4.2500
4.5000
4.7500
5.0000
0.0
0.4587
0.6881
0.0
9.1743
45.8716
34.4037
0.0
0.0
0.0
0.0
0.0
0.0
0*0
0.0
0.0
0.0
0.0
0.0
0.0
5.5000
0.0
0. 2395
0.4981
0.2083
4, 7899
26.7261
31.8458
10.4128
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
6.0000
0.0
0.0867
0.3296
0.2993
1.7336
12.6592
26.4571
14.9671
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
6.5000
0.0
0.0006
0.2225
0,3325
0.0120
4.4928
22.2095
16.6234
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
7.0000
0.0
0.0
0.1318
0.2437
Q.0689
2.6369
14.1025
14.4785
3.4426
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
7.5000
0.0
0.0
0.0738
0.1855
0.1122
1.4755
8.8742
13.0159
5.6120
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
8.0000
0.0
0.0
0.0303
0..1362
0.1362
0.6053
4.8426
11.3513
6.8112
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
3.5000
0.0
0.0
0.0022
0.1012
0.1468
0.0443
2.1794
9.9555
7.3420
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
9.0000
0.0
0.0
0.0001
0.0708
0.1202
0.0228
1.4225
7.3607
6.7179
1.0605
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
9.5000
0.0
0.0
0.0
0.0473
0.0969
0.0389
0.9467
5.2521
6.1421
1.9439
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
Table B.3. Predicted numbers of larvae per 1000 cubic meters in the first
eight 0.5 millimeter length classes for twenty time periods.
-------
BSTIH&TBD SIZE DISTRIBUTIONS
SIZES
00
TIHB
0.0
0.2500
0.5000
0.7500
1.0000
1.2500
1.5000
1.7500
2.0000
2.2500
2.5000
2.7500
3.0000
3.2500
3.5000
3.7500
4.0000
4.2500
4.5000
4.7500
5.0000
0.0
0.7855
1.9765
1.6218
16.3911
94.5733
147.2839
103.4173
36.0678
3.0044
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
10.0000
0.0
0.0
0.0
0.0508
0. 3203
0.4518
1.1442
9.9611
29.9345
26. 8748
6. 4241
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
15.0000
0.0
0.0
0.0
0.0
0.0014
0.0867
0.1787
0.1066
1.8343
9.5982
11.5279
3.8952
0.0037
0.0
0.0
0.0
0.0
0.0
0.0
0.0
20.0000
0.0
0.0
0.0
0.0
0.0
0.0009
0.0480
0.0889
0.0463
1.0238
5.1141
5.3876
1.4140
0.0
0.0
0.0
0.0
0.0
0.0
0.0
25.0000
0.0
0.0
0.0
0.0
0.0
0.0
0.0037
0.0363
0.0491
0.0782
0.9855
3.4136
2.5960
0.2150
0.0
0.0
0.0
0.0
0.0
0.0
30.0000
0.0
0.0
0.0
0.0
0.0
0.0
0.0000
0.0106
0.0288
0.0194
0.2114
1.3156
2.0872
0.9713
0.0002
0.0
0.0
0.0
0.0
0.0
35.0000
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0008
0.0145
0.0215
0.0177
0.3429
1.4176
1.1481
0.1131
0.0
0.0
0.0
0.0
0.0
40.0000
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0000
0.0059
0.0141
0.0078
0.1186
0.6963
0.9630
0.3890
0.0
0.0
0.0
0.0
0.0
45.0000
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0013
0.0087
0.0102
0.0264
0.2653
0.7719
0.5157
0.0093
0.0
0.0
0.0
0.0
50.0000
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0000
0.0050
0.0087
0.0027
0.1032
0.5226
0.4951
0.0931
0.0
0.0
0.0
0.0
Table B.4. Predicted numbers of larvae per 1000 cubic meters in the first
eight 5.0 millimeter length classes for twenty time periods.
-------
TIMEs
TIME=
TIME=
TIflE=
TIME=
TIME*
TIME=
TIMBs
TIME=
TIME=
TIME=
TIME=
TIME*
TIME*
TIMEs
TISE*
THE*
TIME=
TIME=
TIME=
0.0
0.2500
0.5000
0.7500
1.0000
1.2500
1.5000
1.7500
2.0000
2.2500
2.5000
2.7500
3.0000
3.2500
3.5000
3.7500
4.0000
<*.2500
4.5000
U.7500
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
TOTAL
POPULATION
POPULATIONS
POPULATIONS
POPULATION^
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
POPULATIONS
0.0
0.7855
1.9765
1.6756
16.7127
95.1125
148.6581
113.6211
67.9819
40.650-f
24.3228
14.5276
8.6641
5.1604
3.0725
1.7931
0.8475
0.2036
0.0012
0.0
Table B.5. Predicted total population numbers per 1000 cubic meters for
twenty time periods.
87
-------
APPENDIX C: THE COMPUTER PROGRAM
88
-------
MAIN
DEC.
oo
1
2
3
H
5
6
7
3
9
10
11
12
13
11
15
16
17
18
19
20
21
22
23
2H
25
26
27
23
29
30
31
32
33
3*
35
36
37
38
39
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
IMPLICIT REAL*
-------
IAIN
DEC,
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
63
69
70
71
72
73
74
75
76
77
78
79
30
C ~
C —
C —
C —
C —
C
C
C
C
C
C
C
C
C
C
C
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
C -
NSIZEC
NCHNGE
DELCLA
DELRUN
DELSIZ
DELSZA
TO
TMAX
TBRT (I) ~
TIHEA(I) -
TIHEPL(T)
SIZECL(I)
SIZE (1) —
SIZE (I) —
ZNA(I)
ZIA(I)
ZA(I)
GA(I)
BA(I)
Nil
NUT
NUI3T
NINIT
S HUH IS (I)
DENDPA
NGROW
FRACA(I) -
FEACB(I) -
NUMBERS SPAWNED (OR RECRUITED INTO A SPECIFIED SIZE
CLASS) ARE GIVEN
NUMBER 0? INPUT SIZE CLASS DATA FOR ZA (S,T) AND GA(S,T)
NUMBER OF TIMES M WHICH ZA(S,T) AND GA(S,T) CHANGE OVER THE
PROJECTED SIMULATION RUN
LENGTH OF SIZE CLASSES FOR WHICH MORTALITY, ZA(S,T), AND
GROWTH RATE, GA(S,T), ARE GIVES AS INPUT DATA
LENGTH OF TIHE STEP
LENGTH OF SIZE CLASSES IN SIMULATION RUN
LENGTH OF DESIRED SIZE CLASS PRINTOUT
INITIAL TIME
END OF REPRODUCTIVE PERIOD
TIMES DURING SPAWNING PERIOD AT WHICH NUMBERS SPAWNED
(OR RECRUITED INTO A GIVEN SPECIDIED SIZE CLASS) PER
UNIT TIME (MONTH), BA (I), ABE GIVEN
TIMES AT WHICH 2A(SrT) AND G(S,T) VALUES ARE CHANGED (NORMALLY
EACH MO^TH
- TIMES FOR WHICH SIZE DISTRIBUTION IS CALCULATED
- SIZE CLASSES FOR FRICH INPUT DATA ZA(S,T) AND GA(S,T)
ARE GIVEN
SIZE AT BIRTH
AVERAGE LENGTH OF AN INDIVIDUAL IN SIZE CLASS I
SIZE-SPECIFIC NATURAL MORTALITY RATE IN SIZE CLASS I
SIZE-SPECIFIC IMPINGEMENT MORTALITY RATE IN SIZE CLASS I
SIZE-SPECIFIC MORTALITY RATE FOR SIZE CLASS I
SIZE-SPECIFIC GROWTH RATES
REPRODUCTION RATE DURING SPAWNING PERIOD (NUMBERS/MONTH)
THE INTEGER CONTROLS THE PRINTING OUT OF DETAILED COMPUTATIONS
AT GIVES TIME STEPS. FOR EXAMPLE, I? WE SET Nil =10, MASY
COMPUTATIONAL DETAILS WILL BE PRINTED OUT EVERY TIME THE
INDEX I IN DO-LOOP 1QQ IS A MULTIPLE OF SIX.
THE INTEGER CONTROLS THE PRINTING OUT OF DETEILED COMPUTATIONS
NUMBER OF INITIAL SIZE DISTRIBUTION POINTS
IF KINIT = 1, THERE IS AN INITIAL SIZE DISTRIBUTION
- INITIAL NUMBER OF FISH IN SIZE CLASS I
COEFFICIENT OF EFFECTS OF NUMBER DENSITY ON GROWTH RATE
NUMBER OF SUBCOHORTS HAVING DIFFERENT GROWTH RATES
FPACA(LI)*GAA(J) IS THE GROWTH RATE OF SUBCOHORT LI AT SIZE J
FRACB (LI)*SNUMIN (J) IS THE INITIAL POPULATION OF SUBCOHORT LI
IN SIZE CLASS J
-------
SAIN
DEC.
91
82
83
84
85
86
87
88
39
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
C
C
READ (5, 1000) NSIZES, NSIZEC, NBIRTH, NCHNGE, NRUNTM, Nil, NIII
1,»GBOW
1000 FORHAT(14I5)
READ (5, 1001) TO
1001 FORHAT(7E10.0)
READ (5, 1001) DELSIZ,
READ (5, 1001)
READ (5, 1001)
READ(5,1001)
BEAD (5, 1001)
IF(NBIRTH .EQ
READ (5, 1001)
READ(5,1001)
2 CONTINUE
TIMEA(1) = 0.
READ(5,1001)
READ (5, 1001)
READ (5, 1001)
READ (5, 1001)
READ (5, 1001)
DO 3 1=1,500
SNUHIN(I) = 0.0
SZNUHG(I) = 0.0
DO 3 J=1,40
SDISTQ(J,I) =0.0
3 CONTINUE
READ (5, 1000)
IF (NINIT .EQ
READ(5,1001)
7 CONTINUE
DO 10 1=1, NSIZEC
ZA(I) = ZNA (I) + ZIA(I)
GAA(I) = GA(I)
10 CONTINUE
TBRT{1) = TO
SIZECL(1) = SIZE(1)
DO 18 1=1, NSIZEC
SIZECL(H-I) = SI2ECL(I)
DELCLA, DELSZA
SIZE (1)
(ZNA (I) ,1=1, NSIZEC)
(ZI A (I) ,1=1, NSIZEC)
(GA (I) ,1=1 , NSIZEC)
. 0) GO TO 2
(BA (I) ,1=1 , NBIRTH)
(TBRT (I) ,1=1 .NBIRTH)
(TIBEA (I) ,I=2,N3HNGE)
(TIHEPL (I) ,1=1 ,NP.UNTH)
DENDPA
(FRACA(I) ,I=1,NGROH)
(FRAC3 (I) ,I=1,NGROW)
NINIT,
0) GO TO
(SHUHIN(I)
NQINT
7
I=1,NOINT)
DELCLA
-------
SAIN
121 18 CONTINUE
122 WRITE(6,2000)
123 2000 FORMAT(1H1,20X,'PARTIAL DIFFERENTIAL EQUATION MODEL OF A FISH COHO
124 1RT SIZE DISTRIBUTION',////)
125 WRITE(6,2006) NSIZES, NSIZEC, NBIRTH, NCHNGE
126 2006 FORMAT(1H ,5X,»NSIZES=',14, 4X,»NSIZEC=»,I4,4X,'NBIBTH=«,I4,4X,
127 1'NCHNGE=»,I4,//)
128 WRITE{6,2016) NRUNTM,NII, NIII
129 2016 FORMAT (1H ,5X,«NRUNTM=«,14,4X,•NII=«,I4,4X,'NIII=«,I4,//)
130 WRITE(6,2017) NINIT, NUINT
131 2017 FORMAT(1H ,5X, »NINIT=',I4,4X, • NUINT=» ,14
132 HRITE(6,2007) DELSIZ, DELCLA
133 2007 FORMAT(1H ,5X,'DELSIZ= •,E15. 8,2X,'DELCLA=
134 WRITE (6,2009) DELSZA
135 2009 FORMAT(1H ,5X,'DELSZA= »,F10.4,//)
136 WRITE(6,2014)
13"» 2014 FORMAT (1H ,///, 10X,'TIMES FOR WHICH SIZE DISTRIBUTION IS CALCULATE
138 1D',//)
139 DO 16 I=1,NRUNTM
140 WRITE(6,2099) I,TIMEPL(I)
jo 141 16 CONTINUE
142 WRITE(6,2004)
143 BRITE(6,2001)
144 2001 FORMAT(////,5X,«SIZE CLASS',10X,'MORTALITY RATE',6X,'GROWTH RATE',
145 1//)
146 DO 20 1=1,NSIZES
147 SIZE (1+1) = SIZE (I) + DELSIZ
148 20 CONTINUE
149 DO 30 1=1,NSIZEC
150 WRITE(6,2002) SIZECL(I), ZA (I) , GA (I)
151 2002 FORMAT(1H , 5X,6 (E15.8,3X) )
152 30 CONTINUE
153 WRITE(6,2011)
154 2011 FORMAT (1H ,////,20X,'PARAMETER CHANGES OCCUR AT',//)
155 DO 35 1=1,NCHNGE
156 HRITE{6,2012) TIMEA (I)
157 2012 FORMAT(1H ,5X,'TIMSA= ',E15.8)
158 35 CONTINUE
159 WRITE{6,2003)
160 2003 FORMAT(////,5X,'TIMS',16X,»BIRTH SATE',//)
DEC.
-------
MAIN
DEC.
vo
CO
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
19°
200
40
2008
2013
45
2022
2023
2024
C
C
C
C
DO 40 I=1,NBIRTH
WRITE(6,2002) TBET(I), BA (I)
CONTINUE
WRITF.(6,2008) DENDPA
FORSAT(1H ,/r5X,« DENDPA= »,E15.8r//)
SSIZE = DELSZA/DELSIZ
LSIZES = NSIZES/HSIZE
WRITE(6,2013) "SIZE, LSIZES
FORMAT(1H ,///,5X,«SSIZE= ' ,15,5X,•LSIZES= «,I5,//)
DO 45 1=1,40
DO 45 L=1,20
SDISTP(I,L) = 0.0
SDISTH(I,L) =0.0
CONTINUE
WRITE (6, 2022) NGROW
FORMAT (1H , 5X,« NUH3ER OF SDB30HOSTS, NGROW= «,I5,//)
WRITE (6,2023) (FRACA(I) rI=1,NGROW)
WRITE (6,2024) (FRACB(I) ,I=1,NSROW)
FORMAT{1H ,5X,«FRACA= • ,1 4 (F7. 3, 1X) )
FORHAT(1H r5X,»FRACB= f , 1 4 (F7. 3r 1X))
CALCULATE THE ANALYTICAL SOLUTION
ITERATE MODEL OVER TIME
DO 900 LI=1,flGPOW
DO 43 I=1,NSIZES
SDISTG(1,I) = 0.0
43 CONTINUE
JTIME = 0
T = TO
TLAST = TO
JLAST = 1
NRUNTP = NRUNTM * 1
TMAX = TBRT(NBIRTH)
DO 700 1=1,NRUNTM
T = TIMEPL(I)
STOREI =0.0
TEGI (1) = STOREI
STOREJ =0.0
2045 FORMAT(1H ,5X,«I= «
-------
MAIN DEC.
201 DO 46 J=1,NCHNGE
202 TF(T .LE. TIMEA(J+1)) GO TO 47
203 46 CONTINUE
204 47 CONTINUE
205 JTIME = J
206 TI = TIMEA(JTIME)
207 DO 50 J=1,NSTZES
208 IF (I .EQ. 1) GO TO 49
209 SZNUMG(J) =0.0
210 49 CONTINUE
211 SZNNEW(J) =0.0
212 50 CONTINUE
213 DO 51 J=1,NSIZEC
214 GA(J) = FRACA(LI) *GAA{J)
215 51 CONTINUE
216 IF(T .SE. TMAX) GO TO 60
217 DO 55 J=1,NBIRTH
218 IF(TBRT(J) .LE. T) GO TO 55
21Q GO TO 56
220 55 CONTINUE
221 GO TO 60
222 56 CONTINUE
223 BIRTH = BA (J-1) + ( (T-TBET(J-1))/DELBIH) * (BA (J) - BA(J-1»
224 GO TO 61
225 60 CONTINUE
226 BOTH = 0.0
227 61 CONTINUE
228 GHIN = GA(1)
229 SZNNEW(1) = BIRTH/GSIN
230 IF (I .NE. (I/NII) *NII) GO TO 63
231 WRITE(6,2004)
232 2004 FORMAT (1H1)
233 63 CONTINUE
234 C
225 C
236 C CALCULATION OF PORTION OF N(S.T) RESULTING FROM REPRODUCTION IN
237 C THE TIME INTERVAL
238 C
239 DO 100 J=2,NSIZES
240 CALL T3A? (SIZE (J-1) ,DELSIZ, 8, 1 . ,FACTRI,GNV)
-------
IAIN
DEC,
2*1
2*2
2*3
V£5
Ln
2*5
2*6
2*7
2*8
2*9
250
251
252
253
25*
255
256
257
258
259
260
261
262
263
26*
265
266
267
268
269
270
271
272
273
27*
275
276
277
278
279
280
= • ,E15.8, 3X,» J (S) = »,E15.8)
FACTRI = FACTPI * STOSEI
STOREI = FACTRI
TEST (J) = STOSEI
TB = T • - FACTRI
CALL TSAP(SIZE(J-1) ,DELSIZ, 8, 0. , FACTRJ, GNV)
IF(T .HE. (I/NIII)*NIII) GO TO 69
WRITE(6,2099) J,SIZE(J-1) ,DELSTZ, FACTRI, GNV, STOREJ,TB
2099 FORMAT(1H , 5X,I5, 2X,6 (E15.8, 1X) )
69 CONTINUE
FACTRJ = FACTRJ + STOSEJ
STOREJ = FACTRJ
IF(I .WE. (I/NII) *NII) GO TO 75
WRITE(6,2005) I, J, FACTRI, FACTRJ
2005 FORHAT(1H ,5Xrl3, 2X,I3,2X,« I
75 CONTINUE
TLOWER = TLAST
IF(T .GE. TMAX) GO TO 90
IF(TB .GT. TMAX -OR. TB .LT. TLOWER) GO TO 90
DO 80 K=1,NBIRTH
IF(TBRT(K) .LE. TB) GO TO 80
GO TO 81
80 CONTINUE
81 CONTINUE
IF(K .EQ. 1) GO TO 90
BIRTH = BA (K-1) + ( (TB-TBRT (K-1) ) / (TBRT (K)-TBRT (K-1) ) ) *
1 (BA (K) - BA(K-1))
GO TO 91
90 CONTINUE
BIRTH = 0.0
91 CONTINUE
IF (I .NE. (I/NII) *NII) GO TO 95
WRITE(6,2010) TB, BIRTH, GNV,TBPT(K) ,BA(K)
2010 FORHAT(1H ,5X,*TB= • ,E15. 8, 2X, «BIRTH= « ,E15.8,2X, • GNV= »,E15.8,
12X,fTBP.T= ',E15.8,2X,'BA= »,E15.8)
95 CONTINUE
IF (FACTRJ .GT. 30.) FACTRJ = 30.
SZNNEW (J) = BIRTH*GNV*EXP (-FACTRJ)
100 CONTINUE
105 CONTINUE
IF (I .NE. (I/NIII)*NIII) GO TO 111
-------
TAIN
DEC.
ON
281
282
283
28U
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
30*
305
306
307
303
309
310
311
312
313
31U
315
316
317
318
319
320
DO 110 J=1,NSIZES
WRITE(6,2015) SIZE(J), TEGI (J)
2015 FORMAT(1H ,5X,'SIZE= • ,E1 5. 8, 2X, • TEGI= «,E15.8)
110 CONTINUE
111 CONTINUE
C
C
C
150 CONTINUE
DO 155 J=1,NSIZES
IF(I ,EQ.
SZNUMG(J)
GO TO 15U
153 CONTINUE
SZNUMG(J)
SZNUHG(J)
15U CONTINUE
IF(I .HE.
WRITE(6,203H)
203U FORMAT (1H ,5?
155 CONTINUE
156 CONTINUE
WRITE(6,2018)
2018 FORMAT(///)
IF(JTIME . LE.
READ (5,1001)
READ (5,1001)
READ (5,1001)
»BITE(6,2019)
1) GO TO 153
= SDISTG(I-1, J)
= SNUMIN(J)
= SNU!«IN(J)*FRAC« (LI)
(I/NIII)*NIII) GO TO 155
J, SZNUMG(J)
•J= »,I5,5X,«SZNUMG= »,E15.8)
JLAST) GO TO 175
(ZNA (II) ,II=1,NSIZEC)
(ZIA(II) rII=1,NSIZEC)
(GA (II) rII=1,NS!ZEC)
T
2019 FORMAT(A5Xr«CHANGE IN GROWTH AND MORTALITY RATES AT TIME T =
iFio.a,/)
HRITE(6,2001)
DO 180 II=1rNSIZEC
ZA(II) = Z?»A(II) + ZIA(II)
WRITE(6,2002) SIZECL (II) , ZA(II), GA (II)
180 CONTINUE
HRTTE{6,2066) I,JTIM3,JLAST,T
2066 FORMAT(1H ,//,5X,3l5, U (E1 5. 8, 2X))
175 CONTINUE
-------
MAIN DEC.
321 C CALCULATION OF THAT PORTION OF N(S,T) RESULTING FROM INDIVIDUALS
322 C ALREADY PRESENT AT THE START OF THE INTERVAL
323 C
324 IF (I .NE. (I/NIII)*NIII) GO TO 339
325 WRTTE(6,2097)
326 2097 FORMAT (1H ,5X,«IN T .GT. TIMEA (I) • ,/)
327 339 CONTINUE
323 DO 500 J=1,NSIZES
329 ARG = TEGI(J) * TLAST - T
330 IF (ARG .LE. TEGI(1)) GO TO 499
331 DO 350 K=1,500
332 IF (ARG .GT. TEGI (K) ) GO TO 350
333 GO TO 351
334 350 CONTINUE
335 351 CONTINUE
336 SUBS = SIZE (K-1) + ( (ARG-TEGI (K-1) )/(TEGI (K)-TEGI (K-1) ))* (SIZE (K)
337 1 - SIZE (K-1))
338 IF (I .NE. (I/NII) *NII) GO TO 395
339 WRITE(6,2032) J,K,T
310 2032 FORMAT(1H ,2X,2I5,F10. 5)
341 395 CONTINUE
342 DO 400 K=1,500
343 IF (SUBS .GT. SIZE (K)) GO TO 400
344 GO TO 401
345 400 CONTINUE
346 401 CONTINUE
347 FACTRA = SZNUMG (K-1) + ( (SUBS-SIZE (K-1) )/DELSIZ) *
348 1(SZNUMG(K) - SZNUMG(K-I))
349 IF (FACTRA .LE. 0.0) FACTRA = 0.0
350 DO 407 K=1,50
351 IF (SUBS .GT. SIZECL(K)) GO TO 407
352 GO TO 406
353 407 CONTINUE
354 406 CONTINUE
355 G= GA(K-1) + ( (SUBS-SIZECL (K-1) )/DELCLA) *(GA (K)-GA (K-1) )
356 DO 412 K=1,50
357 IF(SIZE(J) .GT. SIZECL(K)) GO TO 412
358 GO TO 413
359 412 CONTINUE
360 413 CONTINUE
-------
SAIN DEC,
361 GD = GA(K-1) + ((SIZE(J) - SIZECL (K-1) ) /DELCLA) * (GA (K) -GA (K- 1) )
362 IF (I .HE. (I/NII) *NII) GO TO 415
363 WRITE(6,2033) ARG,TEGI (J) ,TE3I (K) ,SUBS,FACTRA,G
361 2033 FORHAT(1H ,4X,«ARG= • ,E15.8,1 X, »TEGJ= • ,E15. 8, 1 X, • TEGK= «,E15.8,
365 11X,«SA= »,E15.8r«FACTRA= • , E1 5.8, 1X, »G= »,E15.8)
366 415 CONTINUE
367 GNV = 1./G
368 C
369 C CALCULATION OF THE INTEGRAL R(TI,T)
370 C SEE EQUATION (44) IN TEXT
371 C
372 TOTT = TIHEPL(I) - TIKEPL (1-1)
373 JTEND = 50
374 DJTEND = JTEND
375 DELTAU = TOTT/DJTEND
376 TAU = TLAST + TEGI(J)
377 TAUTEG =0.0
378 DO 450 JT=1,JTEND
379 ARC = TAU - T
^ 380 IF(ARG -LT. TEGI(1)) GO TO 449
00 381 DO 430 K=1,100
382 IF(ARG ,GT. TEGI (K) ) GO TO 430
383 GO TO 431
384 430 CONTINUE
385 431 CONTINUE
386 SUBS = SIZE (K-1) * ( (ARG-TEGI (K-1) ) /(TEGI (K) -TEGI (K-1) ))* (SIZE (K)
387 1 - SIZE (K-1))
388 DO 440 K=1,100
389 IF (SUBS .GT. SIZECL (K)) GO T3 440
390 GO TO 441
391 440 CONTINUE
392 441 CONTINUE
393 Z = ZA(K-1) * ((SUBS-SIZECL(K-1U /DELCLA) *(ZA(K) - ZA(K-1))
394 DG = (GA(K) - GA (K-1))/DELCLA
395 DELL =1.0
396 IF(JT .EQ. 1 .OR. JT .EQ. NSIZES) DELL = 0.5
397 TAUTEG = TAUTEG + DELL*DELTAU*Z
398 449 CONTINUE
399 TAU = TAU + DELTAU
400 450 CONTINUE
-------
CONTINUE
CONTINUE
CONTINUE
TLAST = T
JLA-ST = JTIME
CONTINUE
OUTPUT RESULTS
IF(LI .NE. 4* (LI/4)) 50 TO 851
NSIZMX = SIZS(NSIZSS)
-------
•UIN DEC.
481 T = TISEPL(K)
482 WRITE{6,4102) T, (SDISTH (K,I) ,I=NM,NV)
483 87Q CONTINUE
484 IF(NV .LT. LSIZSS) GO TO 860
485 DO 895 I=1,NRUNTK
486 DO 880 J=1,NSIZES
487 SDISTQ(I,J) = SDISTQ(I,J) + SDISTG(I,J)
488 SIZE{J+1) = SIZE(J) * DELSIZ
489 880 CONTINUE
490 DO 890 J=1,LSIZES
491 SDISTP(I,J) = SDISTP(I,J) + SDISTH (I,J)
492 890 CONTINUE
493 895 CONTINUE
494 900 CONTINUE
495 NV = 0
496 950 NH = NV * 1
497 NV = NM * 9
498 T = 0.0
499 IF(NV .GT. NSIZES) NV = NSIZES
i- 500 WRITE(6,4104)
o 501 WBIT?(6,4103) (SIZE (I) ,I=NM,NV)
502 DO 960 K=1,NHDNTM
503 T = TIMEPL(K)
504 WRITE(6,4102) T, (SDISTQ (K,I) ,I=NM, NV)
505 960 CONTINUE
506 IF(NV .LT. NSIZES) GO TO 950
507 DO 975 L=1fLSIZES
508 SIZE(L-H) = SIZE(L) * DELSZA
509 975 CONTINUE
510 NV = 0
511 970 NM = NV + 1
512 NV = NM * 9
513 T = 0.0
514 IF(NV .GT, LSIZES) NV = LSIZES
515 WRITE(6,4104)
516 BRITE{6,4103) (SIZE (I) ,I=NM,NV)
517 DO 980 K=1,NRUNTM
518 T = TIMEPL(K)
519 BRITE(6,4102) T, (SDISTP (K, 1} , I=NM,NV)
520 980 CONTINUE
-------
MAIN
DEC.
441
442
800
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
4104
4103
4102
810
NRUNT = NRUNTM - 1
NV = 0
NM = NV + 1
NV = NM +9
T = 0.0000
IF ( NV.GT.NSIZES)
WRITE(6,4104)
NV = NSIZES
840
2021
850
851
855
860
E(6,4104)
FOR3AT(lHl,43Xr« ESTIMATED SIZE DISTRIBUTIONS1 ,//,57X,f SIZES1 ,//)
WRITE(6,4103) (SIZE (I) ,I=NS,NV)
FORHAT(1H ,7x,«TTME«,10(Fl1.4n
DO 810 K=1,NRONTH
T = TIMEPL(K)
WRITE (6,4102) T, (SDISTG (K,I) ,I=NM,HV)
FORMAT(1H ,11 (F11.4))
CONTINUE
IF(NV .LI. NSIZES) GO TO 800
T = 0.0
WRITE{6,2004)
DO 850 K=1, NRUNTM
TOTAL =0.0
DO 840 L=1, NSIZES
TOTAL = TOTAL + SDISIG (K,L)
CONTINUE
T = TIMEPL(K)
WRITE{6,2021) Tr TOTAL
FORMAT(1H ,5X,«TIHE= • , F8.4,4X, • TOTAL POPULATION= *,F11.4)
CONTINUE
CONTINUE
DO 855 L=1rLSIZES
SI-ZE{L*1) = SIZEC
CONTINUE
NV = 0
NH = NV * 1
NV = NM + 9
T = 0.0
IF(NV .GT. LSIZES) NV = LSIZES
WRITE(6,4104)
HRITE(6,4103) (SIZE CD ,I=NM,NV)
DO 87Q K=1,NRUNTM
DELSZA
-------
NMN
DEC,
521 IF(NV .LT. LSIZES) GO TO 970
522 STOP
523 END
-------
PLOTT DEC.
1 SUBROUTINE TRAP (SI, DS , M,X,F, 3NV)
2 C
3 C TRAPEZOID METHOD IS USED TO EVALUATE THE
4 C INTEGRALS I(S) (X=1) AND J(S) (X=0) .
5 C
6 C • SUBROUTINE TRAP SAMPLE CALL
7 C CALL TRAP (S,D2LS, HOET, STEPS,X,F,GROWTH)
8 C WHERE S = SIZE CLASS
9 C DELS = LENGTH OF SIZE CLASS
10 C MOP.T = SPECIFIC MORTALITY RATE
11 C STEPS = NO. OF SUBDIVISIONS OVER WHICH THE INTEGRAL IS
12 C APPROXIMATED
13 C X = 0. TO INTEGRATE Z/G
14 C 1. TO INTEGRATE 1/G
15 C F = VALUE OF THE INTEGRAL OVER SIZE CLASS I TO 1*1
16 C GROWTH = 1/GROWTH RATE
T* C
18 COHMON/RATEBK/GA(100) , ZA (103) , STZECL{100), DELCLA
19 DIMENSION 5N(20) ,S(21)
20 COMHON SMAX,V
21 S (1) = SI
22 H = DS / M
23 H2 = M + 1
24 DO 200 1=1,K2
25 DO 50 J=1,100
26 IF(SIZECL(J) .LE. SI) GO TO 50
27 GO TO 51
28 50 CONTINUE
29 51 CONTINUE
30 Z = ZA(J-1) * ((SI - SIZECL(J-I)) /DELCLA) * (ZA (J) - ZA(J-1))
31 G = GA(J-1) + ((SI - SIZECL(J-D)/DELCLA) *(GA(J) - GA(J-1))
32 GN(I) = 1./G
33 S (I+-1) = S (I) + H
34 200 CONTINUE
35 SUH = (Gfl(1) * GN(M2))/2.
36 DO 250 1=2, .1
37 SOf! = SUM * GN(I)
38 250 CONTINUE
3« Y = 1.0
40 IF(X.EQ.O.) Y = Z
-------
DEC.
PLOTT
(... F = Y * H * SOB
tt2 GHV = GN(«2)
43 RETOHN
EUD
-------
TECHNICAL REPORT DATA
(Please read Instructions on the reverse before completing)
. REPORT NO.
EPA-600/7-80-068
2.
3. RECIPIENT'S ACCESSION NO.
. TITLE AND SUBTITLE
A Partial Differential Equation Model of Fish Popula-
tion Dynamics and Its Application in Impingement
Impact Analysis
5. REPORT DATE
March 1980
6. PERFORMING ORGANIZATION CODE
. AUTHOR(S)
P.A.Hackney and T. A.McDonough (TVA,Norris, TN);
D. L. DeAngelis and M. E. Cochran (ORNL)
8. PERFORMING ORGANIZATION REPORT NO.
TVA EDT-101
. PERFORMING ORGANIZATION NAME AND ADDRESS
Tennessee Valley Authority
Division of Energy Demonstrations and Technology
hattanooga, Tennessee 37401
10. PROGRAM ELEMENT NO.
INE624A
11. CONTRACT/GRANT NO.
EPA Interagency Agreement
D8-E721-BE
2. SPONSORING AGENCY NAME AND ADDRESS
EPA, Office of Research and Development
Industrial Environmental Research Laboratory
Research Triangle Park, NC 27711
13. TYPE OF REPORT AND I
Final; 10A8-2/80
PERIOD COVERED
14. SPONSORING AGENCY CODE
EPA/600/13
15. SUPPLEMENTARY NOTES IERL-RTP project officer is Theodore G. Brna, Mail Drop 61,
919/541-2683. TVA project director is Hollis B. Flora TJ.
16. ABSTRACT The report gives results of a study to: (1) develop a mathematical model
describing fish populations as a function of life process dynamics and facilities that
impose additional mortality on fish populations; and (2) improve objective impinge-
ment impact prediction. The model accounts for hatching, growing, and mortality
as functions of time and permits computer simulation of impingement impact. It also
accounts for the genetic and environmental heterogeneity effects on the growth of a
cohort of fish. Gizzard shad data collected by TVA were used to corroborate the
model. Simulated impingement impacts for the steam-electric generating plant reser
voir studied were much less than could be measured in field studies. For a tenfold
increase over observed impingement losses , the model predicted that gizzard shad
stock levels would fall by < 10% for any age group. Similarly, the model with a
hundredfold increase over the observed losses predicted that age IV gizzard shad
stock levels were reduced about 65% from baseline values, with younger groups
showing less response. Model simulations revealed that intake-induced mortality
reduced the total numbers of gizzard shad in each age class by < 1%. These findings
show little effect for a species having high natural mortality, but they cannot be
generalized to other species with significantly different natural mortality patterns.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS
c. COSATI Field/Group
Pollution Impingement
Fishes Reservoirs
Population Electric Power Plants
Mathematical Models
Life Cycles
Mortality
Pollution Control
Stationary Sources
Life Process Dynamics
13B 14B
06C,08A 14G
12A 10B
06F
05K
13. DISTRIBUTION STATEMENT
Release to Public
19. SECURITY CLASS (This Report)
Unclassified
21. NO. OF PAGES
106
20. SECURITY CLASS (Thispage)
Unclassified
22. PRICE
EPA Form 2220-1 (9-73)
105
------- |