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

-------