MODEL/REX WORKSHOP TO DEVELOP A MARINE ECOSYSTEMMODEL OF THE NORTH PACIFIC OCEAN INCLUDING PELAGICFISHES(Co-conveners: Bernard A. Megrey and Michio J. Kishi)

Summary
A 4-day MODEL/REX workshop made several
significant achievements:
1. Assembled an international team of marine
biologists, fisheries biologists, and physicaloceanographers who collectively achieved aconsensus on the structure and function of a
PICES Climate Change and Carrying Capacity(CCCC) prototype lower trophic level (LTL)
ecosystem model for the North Pacific Ocean
that included pelagic fishes, and named it“NEMURO.FISH”;
2. Developed a computer simulation model offish bioenergetics and growth;

3. Coupled the fish model to the NEMURO
lower trophic level model;
4. Adapted the fish bioenergetics model toPacific herring (Clupea harengus pallasi) in
the eastern North Pacific, and Pacific saury(Cololabis saira) in the western North Pacific;
5. Made recommendations for future modelingactivities.
The significance of these achievements will
ultimately be evaluated by how well the CCCC
Program effectively utilizes and embraces thesemodels as a basis of future modeling activity.

1.0 Workshop overview
Introduction
The North Pacific Marine Science Organization(PICES) organizes and promotes an international
science program, CCCC, in the temperate andsubarctic regions of the North Pacific Ocean.
Ecosystem modeling is one of five key research
activities defined by the CCCC Implementation
Panel. The PICES CCCC MODEL Task Team isgiven the role to encourage, facilitate andcoordinate modeling activities within the member
nations with respect to the goals and objectives of
the PICES-CCCC Program. At the 2000 Nemuro
Workshop, the MODEL Task Team developed
NEMURO, a lower trophic level marine
ecosystem model.
NEMURO has beeninternationally recognized, and recently, the Max
Planck Institute has adopted the use of NEMURO.
At PICES IX in Hokkodate, the REX and MODEL
Task Teams met and agreed it would be useful toextend NEMURO to include higher trophic levelcomponents. Based on some presentations there,
we agreed to try Pacific herring as a candidate
higher trophic level species and plans began for a
joint workshop. Dr. Michio Kishi prepared a
proposal to the Heiwa-Nakajima foundation ofJapan to help fund attendance to the workshop.
The proposal was successful and planning began
to hold the next workshop in Nemuro, Japan in
2002.
Goals and objectives of the Workshop
The goals of the 2002 Nemuro PICES workshop
were to (1) develop a bioenergetics-based fishmodel for Pacific herring (Clupea harengus
pallasi) and Pacific saury (Cololabis saira) and (2)to couple this model with output from the
NEMURO lower trophic level model developed atthe 2000 Nemuro PICES workshop.

77Organizing Committee, participants, sponsors
and venue
Drs. Michio J. Kishi, Bernard A. Megrey andFrancisco E. Werner organized the meeting. Drs.
Megrey and. Kishi served as workshop cochairmen. The Heiwa-Nakajima foundation ofJapan, PICES, and the city of Nemuro provided
financial support and access to excellent meetingrooms in the City Hall. The Nemuro Support
Committee supplied local logistical support.
The venue was set at the Multi Purpose Hall, alarge octagon-shaped room, in the Nemuro City
Cultural Center. The hall had a local area network
which included a server workstation, laser andcolor printers, and another personal computerconnected to the Internet. A classroom style table
was arranged in the center of the room for theplenary session. A set of LCD projectors andscreens and AC power outlets for participants’

laptop computers were available and werearranged in each work area to make group work
more effective.
Twenty six scientists from China, Korea, Russia,
Japan, Canada and the United States (Fig. 1)
convened in Nemuro, Japan, between January 25
and January 27, 2002, to participate in a modelingworkshop focused on developing a coupled lower
trophic level-higher trophic level model of themarine ecosystem. Most scientists arrived with
their own laptop computers.
Participants(Appendix 1) consisted of plankton scientists,
modelers, and individuals with biologicalknowledge of herring and saury. Key regional
data sets were also provided by many workshop
participants. The workshop was continued at the
Frontier Research System for Global Change(FRSGC) facilities in Yokohama on January 29,
2002.

Fig. 1 Nemuro Workshop participants. Left to right –Top Row: Douglas Hay, Tomokazu Aiki,
Masakatsu Inada, Daiki Mukai, Lan S. Smith, Vadim V. Navrotsky, Alexander V. Leonov, Francisco E.
Werner, Robert A. Klumb, Bernard A. Megrey, Toshio Katsukawa, Takeshi Okunishi, Yasuhiro
Yamanaka, Tomonori Azumaya. Bottom Row: Chul-hoon Hong, Sanae Chiba, Yuri I. Zuenko, Daji
Huang, Masahiko Fujii, Kazuaki Tadokoro, Shin-ichi Ito, Shoichi Hamaya (Nemuro City Supporter),Michio J. Kishi, Makoto B. Kashiwai.
78

Workshop schedule

13:00-15:30

Date: January 25th-29th, 2002
Venue: Nemuro City Culture Center* (25-27 Jan.
2002), FRSGC (28,29 Jan. 2002)
Conveners:
Michio J. Kishi (Hokkaido
University), Bernard A. Megrey (NOAA),Francisco E. Werner (University of North
Carolina)
Workshop Co-Chairmen:
M. Kishi and B.
Megrey15:40-16:00
16:30-

AgendaJanuary 25th, Friday
18:00
Opening ceremony
19:00
Welcome reception
January 26th, Saturday
09:00-09:10
Remarks by M. J. Kishi
09:10-09:30
Review of NEMURO (North
Pacific Ecosystem Model for
Understanding
Regional
Oceanography) developed by
PICES MODEL Task Team in2000 (Michio Kishi)
09:30-10:30
Review of NEMURO FORTRAN
code (Yasuhiro Yamanaka)
10:30-11:00
Fishbioenergetics/biomass
modeling:
an application to
Pacific herring (Bernard Megrey)
10:30-11:00
Review of NEMURO FORTRAN
code (Yasuhiro Yamanaka)
11:00-11:30
Review of Clupeid biology withemphasis on energetics (Robert
Klumb)
11:30-12:00
Analysis of change in Pacificherring distributions (Douglas
Hay)
12:00-13:00
Lunch
13:00-13:30
ReviewofPacific
saury(Cololabis saira) study under
VENFISH (Shin-ichi Ito)
13:30-17:00
Grouping of scientists (“teamherring” and “team saury”)
January 27th, Sunday
09:00-12:00
Continue working in teams
12:00-13:00
Lunch

18:30-20:30

Discussion on the results andmodification of model
Closing ceremony
Press conference (Megrey, Kishi,
Werner)Farewell party by Nemuro city (athotel)

January 28th, Monday
Move to Frontier Research System for Global
Climate Change
January 29th, Tuesday
09:00-12:00
Discussion on the results of newmodel and future strategy
12:00-13:00
Lunch
13:00-17:00
Seminar at FRSGC
13:00-13:30
Zuenko13:30-14:00
Navrotsky
14:00-14:30
Huang
14:30-15:00
Klumb
15:00-15:30
Hong
15:30-16:00
Tea break
16:00-16:30
5-minute speech of Japaneseparticipants
16:30-17:00
Discussion of future work
Workshop activity
After an opening ceremony with the people of
Nemuro and a welcome party held the day before,
the participants convened at the venue to start theworkshop.
On the first day, the workshop officially opened
with a welcome to all who had endured a long
journey to come back to Nemuro. In the morning
session, individual presentations were made on the
NEMURO LTL model, a review of the
FORTRAN program to execute NEMURO, theproposed fish bioenergetics model, andpresentations on herring and saury biology as
outlined in the agenda. During the afternoon
session, the workshop participants split into twogroups, to adapt the generalized fish bioenergetics
model for Pacific herring (“team herring”) and
Pacific saury (“team saury”).
The second day was taken up primarily with thetwo working groups dealing with their specifictasks. Results of the Pacific herring and saury
79

applications were presented for discussion in theafternoon. Also on the second day the coupledlower trophic level-higher trophic level model wasnamed NEMURO.FISH (North Pacific Ecosystem
Model
for
Understanding
Regional
Oceanography.For Including Saury and Herring).
Robert Klumb suggested the name.
The participants received closing remarks from thevice-chairman of the Nemuro Supporting
Committee where appreciation was extended tohave brought into being such a productive
workshop. These feelings were amplified during a

Sayonara Party, which was full of warm
hospitality by the people of Nemuro city.
The third session was held at the Frontier Research
System for Global Change in Yokohama. The
group discussed the structure and organization of
the final report, made writing assignments,
generated a list of workshop recommendations,discussed where the MODEL Task Team should
be going next, and the possibility of holding future
workshops. Several individual seminars were
presented by workshop participants dealing with
their personal research topics.

2.0 Workshop presentations
This section contains abstracts, extended abstracts,
or fully prepared reports and workshop summaries
given at the workshop. The reports that follow areorganized by authors, according to the schedule

provided in the agenda. The authors whose last
name is in underline and bold font made thepresentation. Model versions referenced in thesereports are described in Megrey et al. (2000).

2.1 A generalized fish bioenergetics/biomass model with an application to
Pacific herring
Bernard A. Megrey1, Kenny Rose2, Francisco E. Werner3, Robert A. Klumb4 and Douglas Hay5
1
2

3
4
5

National Marine Fisheries Service, Alaska Fisheries Science Center, 7600 Sand Point Way NE,
Seattle, WA 98115, U.S.A. E-mail: bern.megrey@noaa.gov
Coastal Fisheries Institute and Departm ent of Oceanography and Coastal Sciences, Wetlands
Resources Building, Louisiana State University, Baton Rouge, LA 70803, U.S.A. E-mail:karose@lsu.edu
Marine Sciences Department, CB # 3300, University of North Carolina, Chapel Hill, NC 27599-3300,
U.S.A. E-mail: cisco@unc.edu
Department of Natural Resource, Cornell Biological Field Station, Cornell University, 900
Shackelton Point Road, Bridgeport, NY 13030, U.S.A. E-mail: rak11@cornell.edu
Pacific Biological Station, Fisheries and Oceans Canada, 3190 Hammond Bay Rd, Pacific Biological
Station, Nanaimo, British Columbia, Canada V9R 5K6. E-mail: hayd@pac.dfo-mpo.gc.ca

We chose to use bioenergetics/biomass modelingto represent fish growth because (1) the theory isbased on the Law of Thermodynamics, (2) outputs
must equal inputs, ie., the energetic budget must
balance (Law of Conservation of Mass), (3) terms
in the equations are simple to biologicallyinterpret, (4) fish physiological terms are wellknown and in general can be directly measured,and (5) this modeling approach allows users tofocus on important external regulators such astemperature and diet composition.
Modelformulation and parameters for Pacific herringfollowed the approach used by Rudstam (1988) for
Atlantic herring (Clupea harengus).
The growth rate of an individual Pacific herring(non reproductive) is calculated as weightincrement per unit of weight per time and isdefined by:80

(2.1.1)

CAL z
dW= [C − ( R + S + F + E )] ⋅
⋅W
dt
CAL f

where C is consumption (g prey·g fish-1·d-1), E isexcretion or losses of nitrogenous excretory wastes(g prey·g fish-1·d-1), F is egestion or losses due tofeces (g prey·g fish-1·d-1), R is respiration or losses
through metabolism (g prey·g fish-1·d-1), S is
specific dynamic action or losses due to energy
costs of digesting food (g prey·g fish-1·d-1), W is
the weight of the fish (g wet weight), t is time(days) CALz is the caloric equivalent of
zooplankton (cal·g zooplankton-1), and CALf is thecaloric equivalent of fish (cal·g fish-1). Note that
(2.1.1) does not include energetic costs ofreproduction (spawning).
If we define CALz as calories·g zooplankton-1

CAL z =

2580 joules
1 cal⋅= 617.22
gram zoop 4.18 joules

and CALf as calories·g fish-1

CAL f =

5533 joules
1 cal⋅= 1323.68
gram fish 4.18 joules

then once the change in weight from 2.1.1 iscomputed in terms of g zooplankton·g fish-1·d-1, we
can multiply it by the weight of the fish (W, g) to
get g zooplankton·d-1, and finally convert gzooplankton·d-1 to g fish·d-1 by multiplying thechange in weight (dW/dt) by the ratio CALz/CALf.
In the simulations described in this report,
equation 2.1.1 was solved using an Euler
numerical integration routine using a dt=0.01.
The formulation of the individual processes
represented by the terms in equation 2.1.1 isdescribed individually below. Consumption andrespiration are nonlinear functions of fish weight
and water temperature.
In addition to the physiological parameters, themodel requires information about caloric content
of herring (which can change seasonally), caloric

content of the prey, diet composition, preydensities, and water temperatures.
Consumption
Consumption is estimated as the proportion of
maximum daily ration for herring at a particularmass and temperature.
Maximum daily
consumption rate (g of prey per g body mass of
herring per day) is estimated using an allometricfunction of mass from ad libitum feeding
experiments conducted at the optimum
temperature.
The basic form of the consumption function is(2.1.2) C = C MAX ⋅ p ⋅ f C (T )b

(2.1.3) C MAX = aC ⋅ W C
where C is the specific consumption rate (g prey·g
fish-1·d-1), CMAX is the maximum specific feeding
rate (g prey·g fish-1·d-1), p is the proportion of
maximum consumption, fC(T ) is a temperaturedependence function for consumption, T is watertemperature (ºC), W is herring mass (g wetweight), ac is the intercept of the allometric massfunction (for a 1 g fish at 0°C), and bc is the slopeof the allometric mass function. The subscript C
on the parameters refers to the consumption
process.
In equation (2.1.2), the maximum specific feeding
rate is modified by a water temperaturedependence function described below and an
additional proportionality constant that accounts
for ecological constraints on the maximum feedingrate. The p can range from 0 to 1, with 0
representing no feeding and 1 indicating the fish is
feeding at its maximum rate, based on its bodymass and water temperature. The lower panel of
figure 2.1.1 shows the relationship between fishweight and consumption from equation 2.1.3.
Temperature dependence for cool and cold
water species (Thornton and Lessem 1978)
The Thornton and Lessem description oftemperature dependence is essentially the product81

of two sigmoid curves: one curve is fit to theincreasing portion of the temperature dependence
function (gcta), and the other to the decreasing
portion (gctb). Four temperatures and percentages
are needed. We used two sets of parameters, one
for herring ages ≤1 year old, and one set for
herring > 1 years old.
As an example, parameters for the second set are
xk1=0.1, xk2=0.98, xk3=0.98, xk4=0.01, te1=1.0,
te2=13.0, te3=15.0, and te4=23.0.
For theincreasing part of the curve, te1 is the lower
temperature at which the temperature dependence
is a small fraction (xk1) of the maximum rate, andte2 is the water temperature corresponding to alarge fraction (xk2) of the maximum consumption
rate. For the decreasing portion of the curve, te3 is
the water temperature (≥ te2) at which dependence
is a fraction (xk3) of the maximum, and te4 is thetemperature at which dependence is some reduced
fraction (xk4) of the maximum rate.

function for a theoretical set of parameters. The
upper panel of Figure 2.1.1 shows the Thornton
and Lessem temperature adjustment function overa typical temperature range, and Figure 2.1.3
shows the flexibility of this curve by adjusting te2
for a range of temperatures. Finally, Figure 2.1.4
shows the multi-dimensional relationship between
consumption, body mass and water temperaturefrom equation 2.1.2.

The temperature dependence model is given by(2.1.4) f C (T ) = gcta ⋅ gctb
where T is water temperature (ºC)

1(te2 − te1)⎡
(1.0 − xk1) ⎤t 5 = tt 5 ⋅ ln ⎢ xk 2 ⋅(0.02 ⋅ xk1)⎥⎦⎣
t 4 = e [t 5⋅(T −te1)]
1
tt 7 =(te4 − te3)⎡
(1.0 − xk 4) ⎤t 7 = tt 7 ⋅ ln ⎢ xk 3 ⋅(0.02 ⋅ xk 4)⎥⎦⎣

tt 5 =

Fig. 2.1.1
Relationship between consumption
and temperature from equation 2.1.4 (upper panel)and consumption and weight from equation 2.1.3(lower panel).

(13,0 .98)
te2, xk2(22, 0.98)
te3, xk3

4 temp – F of T
pairs define thefunctional
relationshipt 6 = e [t 7⋅(te 4−T )](xk1 ⋅ t 4)
(1.0 + xk1 ⋅ (t 4 − 1.0))(xk 4 ⋅ t 6)
gctb =(1.0 + xk 4 ⋅ (t 6 − 1.0))
gcta =

Figure 2.1.2 shows an example of the Thornton
and Lessem (1978) temperature adjustment

te1, xk1(1, 0.1)

te4, xk4(27, 0.01)

Fig. 2.1.2
Example of the Thornton andLessem (1978) temperature adjustment curve for atheoretical set of parameters.

82

9
18
12
15

Fig. 2.1.3
Example of the Thornton andLessem (1978) temperature adjustment curve from
Figure 2.1.2 as a result of changing te2 from 9, 12,
15, 18.

in g, fR(T) is the temperature dependence function
for respiration, T is temperature in ºC, aR is theintercept of the allometric mass function andrepresents
theweightspecific
oxygenconsumption rate of a 1 g fish (g O2·g fish-1·d-1) at
0ºC and no activity, bR is the slope of theallometric mass function for standard metabolism,
activity is the activity multiplier, S is the specificdynamic action, SDA is the proportion of
assimilated energy lost to specific dynamic action,
Cisthe
specific
consumption
rate(g prey·g fish-1·d-1) and F is the specific egestion
rate (g prey·g fish-1·d-1). The subscript R on theparameters refers to the respiration process. The
coefficient 5.258 converts g O2·g fish-1·d-1 into
g prey·g fish-1·d-1 using the conversion(2.1.7)⎛ 13560 joules
1 cal ⎞ ⎛ 2580 joules
1 cal ⎞
⎟÷⎜⎟ = 5.258 g zoop / g (O2 )
⋅
⋅
⎜⎜
gO2
4.18 joules ⎟⎠ ⎜⎝ g zoop
4.18 joules ⎟⎠⎝The temperature dependence function for
respiration is a simple exponential relationshipgiven by(c ⋅T )

(2.1.8) f R (T ) = e R

where cR approximates the Q10 (the rate at which
the function increases over relatively low watertemperatures).

Fig. 2.1.4
Plot of the consumption, temperature
and weight relationships from equation 2.1.2.

Activity is a power function of body weight
conditioned on water temperature and is given by( d ⋅U )(2.1.9) activity = e R

Respiration
The respiration or metabolic rate is dependent on
body weight, ambient temperature and activity(swimming speed).
Total metabolic rate isestimated by adding the costs of respiration to thecosts of digestion, specific dynamic action (SDA).
Metabolism is modeled asb(2.1.5) R = a R ⋅ W R ⋅ f R (T ) ⋅ activity ⋅ 5.258(2.1.6) S = SDA ⋅ (C − F )where R is resting respiration (i.e. standard
metabolism) in (g O2·g fish-1·d-1), W is wet weight

where U is swimming speed in cm·s-1 and dR is acoefficient relating swimming speed tometabolism. Swimming speed is calculated as afunction of body weight and temperature usingb(c ⋅T )

(2.1.10) U = a A ⋅ W A ⋅ e A

where aA= 3.9, bA = 0.13 and cA=0.149 if T <9.0 °C
and aA= 15.0, bA = 0.13 and cA=0.0.
if T ≥9.0 °C
Figure 2.1.5 shows the three dimensional
relationshipbetweenrespiration,water
temperature and fish weight.

83

(Cj, g prey·g fish-1·d-1) using(2.1.13) C r =

n

∑C

j

j =1

C MAX ⋅

(2.1.14) C j =Fig. 2.1.5
Relationshipbetweenstandard
respiration, weight and temperature from equation
2.1.5.
Egestion and excretion
Egestion (F, fecal waste) and excretion (E,nitrogenous waste) can be computed as a constant
proportion of consumption.
(2.1.11) F = a F ⋅ C(2.1.12) E = a E ⋅ (C − F )where aF and aE are constant proportions ofconsumption for egestion and excretion
respectively. The subscript F and E on theparameters refers to the egestion and excretion
process.
Multispecies feeding functional response

In most cases realized consumption is calculated
by adjusting CMAX from equation 2.1.2 by p, andthis would be sufficient if there were only one
prey type by using a Type II functional response
equation (Fig. 2.1.6). When there are multiple
prey types, realized consumption depends on preydensities, vulnerability of each prey item toherring (the predator), and half-saturation
constants governing the rate of herring saturation.
A Type II functional response equation for
multiple prey types (after Rose et al. 1999) is used
to compute realized daily consumption of eachherring i (Cr, g prey·g fish-1·d-1) and theconsumption
of
each
prey
type
j

PDij ⋅ vij

K ij
PDik ⋅ vik
1+ ∑
K ik
k =1
n

where CMAX , which is dependent on the weight ofan individual fish and water temperature, is theconsumption rate (g prey·g fish-1·d-1) of individual
herring i from equation 2.1.3, PDij is the density ofprey type j (g wet weight/m3), vij is thevulnerability of prey type j to herring i(dimensionless), and Kij is the half saturation
constant (g wet weight/m3) for individual herring i
feeding on prey type k (k=1, 2, … j, … n).
Because the herring model is tracking one fish,
there is only one predator.
A total of three prey types are represented in thecurrent fish model, microzooplankton, copepods
and euphausids. The prey densities are read infrom the NEMURO model (µmole N/liter) andconverted to g wet weight/m3 using the conversion
14µg N 1.0e −6 g
1g dry weight
1g wet weight 1.0e 3 liters
⋅
⋅
⋅
⋅
µg
µmole N
0.07 g N dry weight 0.2 g dry weightm3

In Figure 2.1.7, the time-dependent solution to the
NEMURO model for the three prey groups at theStation P location is shown. These data were used
to drive herring consumption using the multiple
species functional response model.
In the situation when there are multiple prey types,
Figure 2.1.6 becomes more difficult to graphically
represent.
Figures 2.1.8 to 2.1.11 represent
equations 2.1.13 and 2.1.14 for various parametervalues.
In Figure 2.1.8, we represent fish consumption of
three prey types from the NEMURO LTL model(small zooplankton, large zooplankton andpredatory zooplankton) as stacked bars, where theheight of the bar is cumulative consumption fromequation 2.1.13, and the colored segments within abar represent the consumption of each prey type.

84

FR =

a⋅ x
b+ x

a = 10b = 100

Fig. 2.1.6
TypeIIfunctional
response
describing the theoretical relationship between
available food density and feeding rate when there
is just one prey type.

Figure 2.1.8, the vulnerability parameter for small
zooplankton was varied from 0.0 to 1.0, while
keeping the vulnerability parameter for large
zooplankton and predatory zooplankton equal to
1.0 and using the parameters at the left of the
figure in equations 2.1.13 and 2.1.14. In the
middle panel, just the vulnerability for large
zooplankton was varied from 0.0 to 1.0, while
holding the vulnerabilities for small zooplankton
and predatory zooplankton at 1.0. In the bottom
panel, only the vulnerability for predatory
zooplankton was varied from 0.0 to 1.0.
These results show that, for the prey whose
vulnerability is changing (let us call it the target
prey type), the contribution of the target prey typeto total consumption ranges from 0.0 at 0.0
vulnerability, gradually increases as vulnerability
increases, until a vulnerability of 1.0, where its
contribution to total consumption is exactly one
third. Also total consumption gradually increases
as the proportion of the target prey type increases
with increasing vulnerability to the predator.
Also note that the right-most bar in each panel isthe same (height and contribution of each preytype) when vulnerability is 1.0 for all prey types.
Using Figure 2.1.8 as a base case, Figure 2.1.9
shows the change when the half saturation
constant for large zooplankton (K2) is changed
from 100.0 to 10.0. Now each panel in Figure
2.1.9 is similar to the corresponding panel in
Figure 2.1.8 (the base case), except that large
zooplankton make up the bulk of total
consumption regardless of which prey types
vulnerability is changed.

Fig. 2.1.7
NEMURO model output showing
time-dependent dynamics of small, large andpredatory zooplankton.
The parameters at the left of the figure were usedin equations 2.1.13 and 2.1.14. For each panel
within a figure, the vulnerability of one prey type
was changed from 0 to 1, while keeping all otherparameters the same and assigning thevulnerability parameter for the remaining two preytype to 1.0. For example, in the top panel of

Now using Figure 2.1.9 as a base case, Figure
2.1.10 shows the change when the density of
predatory zooplankton (PD3) is changed from 2.0
to 4.8. Now each panel in Figure 2.1.10 is similar
to the corresponding panel in Figure 2.1.9, exceptthat the contribution of predatory zooplankton tototal consumption is higher in each case. Also, theheight of each bar (total consumption) is higher inFigure 2.1.10 compared to Figure 2.1.19.
Figure 2.1.11 shows the results of the mutispecies
feeding functional response for the parametervalues used in the herring application.

85PD1=2.0
PD2=2.0
PD3=2.0
CMAX=2.5

PD1=2.0
PD2=2.0
PD3=4.8
CMAX=2.5

K1=100.0
K2=100.0
K3=100.0

K1=100.0
K2=10.0
K3=100.0

Fig. 2.1.8
An example of the mutispecies
functional response formulation (equations 2.1.13
and 2.1.14) for three prey groups, varying thevulnerability of the target prey group one at a time.

PD1=2.0
PD2=2.0
PD3=2.0
CMAX=2.5
K1=100.0
K2=10.0
K3=100.0

Fig. 2.1.9 An example of the multispecies
functional response formulation (equations 2.1.13
and 2.1.14) for three prey groups, varying thevulnerability of the target prey group one at a time,and changing the half saturation constant for preygroup 2 (K2) from 100.0 to 10.0.
Linking a fish bioenergetics model to the
NEMURO LTL model
The NEMURO LTL model and the fishbioenergetics
model
were
developed
independently. Linking the two models involves
paying close attention and reconciling two
important differences: 1) the way the two models

Fig.2.1.10 An example of the multispecies
functional response formulation (equations 2.1.13
and 2.1.14) for three prey groups, varying thevulnerability of the target prey group one at a time,changing the half saturation constant for preygroup 2 (K2) from 100.0 to 10.0, and changing thedensity of prey group 3 (PD3) from 2.0 to 4.8.

PD1=1.8
PD2=0.8
PD3=4.8
CMAX=2.5
K1=100.0
K2=10.0
K3=100.0

Fig. 2.1.11 An example of the multispecies
functional response formulation (equations 2.1.13
and 2.1.14) for three prey groups, varying thevulnerability of the target prey group one at a timeusing the parameters in the herring model.
account for time, and 2) the way NEMURO
generates phytoplankton and zooplankton densities(mole N/liter), and the way the fishbioenergeticsmodel expects phytoplankton andzooplankton densities (µmole N/liter). These
differences are presented in Table 2.1.1.
Reconciling these differences requires the use of
several conversion coefficients, which can be seenin the code presented in Appendices 4 and 5.

86

Table 2.1.1 Ways in which NEMURO and thefish bioenergetics model account for time and LTL
densities.
Model
NEMUROFish Bioenergetics

Timesecondsday

LTL Density
mole N/literµ mole N/liter

Linking the fish bioenergetics model to NEMURO
can be done in two ways. In a static linkage (Fig
2.1.12), the NEMURO model is run and a timeseries of small, large and predatory zooplankton
abundances are stored in an output file and used as
an input file for the fish bioenergetics model

where they influence the consumption term of thebioenergetics governing equation 2.1.1.
The
models are run sequentially and there is nofeedback between the two models.
In the dynamic linkage (Fig 2.1.13), the models
are run simultaneously, the zooplankton preygroups contribute to the consumption term of the
fish bioenergetics governing equation 2.1.1, theZOOS, ZOOL, and ZOOP state variables of
NEMURO are reduced by the amount eaten byherring, fish excretion waste is added to thenitrogen pool of NEMURO, and fish egestion
waste is added to the DOM pool of NEMURO.

Fish Model

Fish Model

1 dW= [C − ( R + SDA + F + E )]
W dt

1 dW
= [C − (R + SDA+ F + E)]
W dt

Run NEMURO,
save output file of
Small, Lg and Pred
Zoop and use todrive fish model

NEMURO Small
Zooplankton
NEMURO Large
Zooplankton
NEMURO Predatory
Zooplankton

Fig. 2.1.12 Example of a static linkage between
NEMURO and the bioenergetics fish model.

NEMURO Small
Zooplankton

NEMURO NH4
NEMURO PON

NEMURO Large
Zooplankton
NEMURO Predatory
Zooplankton

Fig. 2.1.13 Example of a dynamic linkage
between NEMURO and the bioenergetics fishmodel..

Table 2.1.2 Summary of parameter values used in the generalized herring bioenergetics model from
Rudstam (1988).
SymbolParameter description
Consumption, CMAX
aC
Intercept for CMAX at (te1+te3)/2bC
coefficient for CMAX versus weight
te1
Temperature for xk1 (in ºC)
te2
Temperature for xk2 (in ºC)
te3
Temperature for xk3 (in ºC)
te4
Temperature for xk4 (in ºC)
xk1
Proportion of CMAX at te1
xk2
Proportion of CMAX at te2
xk3
Proportion of CMAX at te3
xk4
Proportion of CMAX at te4
Metabolism, R
aR
bR
cR
dR
S

Intercept for R
Coefficient for R versus weight
Coefficient for R versus temperature
Coefficient for R versus swimming speed
Coefficient for Specific Dynamic Action

Value
0.642
-0.256

1a 1b
15a 13b17a 15a25a 23b

0.10
0.98
0.98
0.01
0.0033-0.2270.0548
0.03
0.175

87

Table 2.1.2 (cont.)
SymbolParameter description
Swimming Speed, U
aAIntercept U (< 9 ºC) (in cm/s)aA
Intercept U (≥ 9 ºC) (in cm/s)
Coefficient U versus weightbA
Coefficient U versus temperature (< 9 ºC)
cA
cA
Coefficient U versus temperature (≥ 9 ºC)

Value
3.9
15.0
0.13
0.149
0.0

Egestion and Excretion, F and E
aF
Proportion of consumed food egested
Proportion of consumed food excreted
aE

0.16
0.10

a - values for age 0 and 1herringb - values for age 2 and older herring

2.2 Review of Clupeid biology with emphasis on energetics
Robert A. Klumb
Department of Natural Resources, Cornell Biological Field Station, Cornell University, 900 Shackelton
Point Road, Bridgeport, NY 13030, U.S.A. E-mail: rak11@cornell.edu
The general bioenergetics model based on the Lawof Thermodynamics balances all consumed energy
as follows: G = C – R – F – U, where G=growth,
C=consumption, R=metabolism (respiration),
F=egestion, and U=excretion. Consumed energy
is first allocated to costs of metabolism and waste
losses with the remainder available for somatic
growth. Energy lost by the gametes released
during spawning can also be included. Formulasand parameters provided below for the individualcomponents in the bioenergetics model follow theterminology and symbols used in Hansen et al.
(1997). Energy equivalent conversion factors foroxygen consumption, carbohydrates, fats, andprotein can be found in Elliott and Davison(1975), with additional comments on theoxycalorific coefficient found in Brett (1985).
Consumption
Consumption (C) = Cmax*P-value*f(T) andCmax=CA*WCBConsumption (g prey·g-1·d-1), is generally modeled
as an allometric (power) function of weight.

Maximum consumption rates are determined inlaboratory experiments by feeding fish a known(by weight) ad libitum ration and then subtracting
uneaten food after a specified time interval. For
adult alewife Alosa pseudoharengus, the specificslope for weight dependence on maximum
consumption was -0.3 (Stewart and Binkowski1986), a value intermediate to that found in studiesof larval and juvenile clupeids (De Silva andBalbontin 1974; Theilacker 1987). The specificweight-dependent slope (CB) for maximum
consumption of northern anchovy Engraulismordax larvae (wet weight < 0.001 g) recalculated from data in Theilacker (1987) was-0.367, while the slope for Atlantic herring (wetweight 8 – 15 g) was -0.256 (De Silva andBalbontin 1974). Rudstam (1988) used the slope
and intercept derived by De Silva and Balbontin(1974) in the bioenergetics model for adult
Atlantic herring Clupea harengus consumption.
Due to a lack of data for larval and juvenile fishes,
the same relations for maximum consumption of
adult herring and alewives were applied to age-0
fish by Arrhenius (1998a) and Klumb et al. (inreview), respectively.

88

The “P-value” in the bioenergetics model refers to
the proportion of maximum consumption. Thisvalue is used to fit the bioenergetics model toobserved growth or can be set constant to check
resultant growth potential in varied environments.
Temperature dependence of consumption isusually modeled as simple or modified exponentialfunctions (Hansen et al. 1997). For cool- andcold-water species, the temperature dependence ofconsumption is generally modeled using a curve
proposed by Thorton and Lessem (1978), which
modified the logistic equation. This function is
the product of two intersecting sigmoid curves(one ascending and one descending) forming a“humped” curve across the entire temperature
range inhabited by a given species.
Required parameters include the approximate
temperatures for optimum consumption, and thehigh and low temperatures where consumption isdramatically reduced (~98%) compared tomaximum consumption.
Any temperaturesderived from laboratory or field data showing
maximum or reduced consumption levels can beused. If specific data relating consumption andtemperature are lacking, the optimum of
consumption is generally equated to the fish’s
thermal optimum for growth (Beitinger andMagnuson 1979), and the temperatures where
consumption is dramatically reduced are derived
from the thermal tolerances (survival limits) of a
species.
Metabolism/respiration
Total metabolism = Respiration + specificdynamic action (SDA)
where
Respiration (R) = RA*WRB * f(T)*Activity andf(T) = eRQ*T
Metabolism of fishes is determined by measuringoxygen consumption at various temperatures overa known time period, and generally modeled as anallometric function of weight and an exponentialfunction of temperature. Brett and Groves (1979)
distinguished three types of metabolism in fishes:standard, routine and active.
By definition,

standard metabolism is the minimum energy
requirements needed by a fish at rest (also known
as basal metabolism), and it is this metabolic statethat is used in bioenergetics models. Measuring
standard metabolism is difficult and requires useof anesthetized fish or fish with movements
confined by small respirometers.
Routinemetabolism includes normal spontaneous activity,
while active metabolism includes the cost for
activity above the spontaneous activity level.
Winberg (1956) stated that active metabolism wasapproximately twice standard metabolism (i.e. the“Winberg multiplier” of 2). However, Ware(1975) indicated active metabolism could range
from 2 to 3 times standard rates. Bioenergetics
models generally use allometric function
parameters derived for standard metabolism
multiplied by a temperature function and an
activity factor to estimate total respiration costs.
Respiration (g oxygen·g-1·d-1) of adult fishes
generally scales negatively with weight (i.e.
negative slope), and ranges from -0.25 to -0.15 on
a weight specific basis (Winberg 1956). For
clupeids, slopes of the metabolism-weight
relations ranged from -0.19 to -0.28 for Atlantic
menhaden Brevoortia tyrannus (Hettler 1976),-0.215 for alewife (Stewart and Binkowski 1986),
and -0.227 for Atlantic herring (De Silva andBalbontin 1974). Rudstam (1988) used -0.227 in
the adult Atlantic herring bioenergetics model, andthis value was also applied to age-0 herring (Kerrand Dickie 1985; Arrhenius 1998a). The slope
for the metabolism-weight relation of Maurolicus
muelleri, a mesopelagic planktivore, was -0.15(Ikeda 1996).
The relation of respiration to weight of fishes has
been found to change ontogentically, withisometric (mass independent) relations for larvae
switching to negative allometries in adults (Postand Lee 1996). However, the variability of slopes
found in the review of 31 species by Post and Lee(1996) highlighted the need to derive weightmetabolism relations for the larvae of individual
species. The final weight-metabolism relation
derived likely depends on the range of fish sizes
used. Studies of larval fishes encompassing
greater than three orders of magnitude in weight
documented
isometric
relations
betweenmetabolism and weight for Clupeidae (Klumb et89

al., in review), Cyprinidae (Kamler 1972), andScombridae (Giguère et al. 1988).
Specific dynamic action (SDA)
SDA = SDA*(C – F)
Specific dynamic action, or more appropriately
termed “apparent specific dynamic action” andalso known as the “heat increment”, is the energy
allocated to the digestive processes of food,principally deamination of proteins but also
includes energy costs of absorption, transportationand deposition of food (Beamish 1974). Oxygen
consumption by fasting and fed fish in flowthrough respirometers (where the fish is subjected
to a known level of activity, i.e., forced to swim
against a known current) is required to measure
SDA (Beamish and Trippel 1990). Beamish andTrippel (1990) found that SDA increased withmeal size and body weight but declined withweight at fixed rations.
However, in mostbioenergetic models, SDA is considered a constant
proportion of ingested energy with values for adultfish ranging from 10-29% (reviewed by Beamish
and Trippel 1990). The SDA parameter inbioenergetics models is generally borrowed from
non-related species because proper measurement
requires strict laboratory experiments using
specialized equipment. For adult alewife (Stewart
and Binkowski 1986) and adult Atlantic herring(Rudstam 1988), SDA was assumed to be 17.5%
based on data for aholehole Kuhlia sandvicensis(Muir and Niimi 1972).
Arrhenius (1998a)
lowered SDA to 15% for age-0 Atlantic herring.
Larval clupeids have been found to assimilate foodmore efficiently than adults (Kiørboe et al. 1987).
In energetic terms, Kiørboe et al. (1987) estimated
SDA for larval Atlantic herring to be 10% of
assimilated rations, and Limburg (1994) calculated
the mean SDA for American shad Alosasapidissima juveniles to be 13%.
Activity
Activity = eRTO*VEL,
where VEL = RK1*WRK4 for T ≥ RTL
or VEL = ACT*WRK4*eBACT*T when T < RTL
The energetic cost of activity is generally
considered a multiple of standard metabolism. A

simple constant, i.e. the “activity multiplier = 2” ofWinberg (1956), can be used to accord increased(aerobic) metabolic costs due to swimming.
Exponential functions have been used to model
activity costs of adult alewife (Stewart andBinkowski 1986) and Atlantic herring (Rudstam
1988). The exponential model is composed ofthree components: 1) VEL which is the weight
dependence of swimming speed (cm/s), 2) thetemperature (T) dependence of swimming speed(BACT), and 3) the relation of respiration toswimming speed (RTO). The parameter ACT is
the intercept (cm/s) for a 1-g fish at 0°C.
Swimming speed can change from temperaturedependence to independence (at T = RTL).
Swimming speeds of Atlantic herring were only
dependent on weight at temperatures > 9°C(Rudstam 1988), and alewife swimming speeds
were independent at > 15°C. (Stewart andBinkowski 1986)
The coefficient for swimming speed dependenceof metabolism (RTO) used in the adult alewife
model was assumed constant (RTO = 0.03) andbased on data in Muir and Niimi (1972). Data for
Cape anchovy E. capensis found that coefficients
before, during, and after feeding ranged from 0.01
to 0.04 (James and Probyn 1989). A coefficient
relating respiration and swimming speed of 0.03
has also been reported for adult menhaden (Durbin
et al. 1981), and the coefficient for adultcoregonids was 0.02 (Dabrowski 1985). However,the coefficient relating respiration rate toswimming speed increased substantially in larval
coregonids (Dabrowski 1986) and cyprinids(Kaufmann 1990); therefore, a constant relating
metabolic cost to swimming speed is inappropriate
for early life stages. Using an exponential activity
function and a constant relating swimming speed
to oxygen consumption resulted in essentially noenergetic costs for the activity of YOY alewife(Klumb et al., in review).
How to best model the activity costs of larval fish
is uncertain since existing data from the few
studies relating metabolism and swimming speeds
at early life stages are equivocal. Because theslope of metabolism versus swimming speed
varied with body size, Dabrowski et al. (1988)
found active metabolic rates of coregonids to be
5 - 50 times standard metabolism. A size-effect on90

the slope for the metabolism-swimming speed
relation also existed for larvae of two cyprinid
species (Kaufmann 1990); however, ratios ofroutine metabolism to standard metabolism werelow (< 1.5) and essentially flat from 0.005 - 0.300
g (wet weight). In Kaufmann’s (1990) study,
ratios of active to routine metabolism (i.e., thefactorial scope) ranged from 2 - 4. These
contrasting results may lie in the function chosen
to describe the metabolism-swimming speed
relation, i.e., exponential (Dabrowski 1986;
Dabrowski et al. 1988) or allometric (Kaufmann1990). However, using an exponential model,
Wieser and Forstner (1986) found the ratios ofactive to routine metabolism for larvae of threecyprinid species ranged from 1 - 4 and wereindependent of weight (0.01 - 0.3 g wet) andtemperature (12 - 24°C). Activity rates of fishes
can also vary widely with growth rate and fooddensity(Ware1975),
whilelaboratory
measurements of metabolism during activity may
be higher than actual costs in the wild, since larvae
are also passively moved by water currents.
Klumb et al. (in review) used routine metabolism
parameters without an activity multiplier in abioenergetics model for age-0 alewife.
Clupeids have pronounced changes in activity
patterns possibly due to circadian rhythms (Katz
1978; Batty 1987). Clupeids do not swim inschools during darkness (Limburg 1994).
Accuracy of bioenergetic estimates of herringconsumption were improved when including dielfeeding cycles (Arrhenius 1998a).
Egestion
Egestion (F) = FA*C
Egestion is modeled as a constant proportion of
consumption. Assimilation efficiency (in terms ofenergy) of adult menhaden ranged from 86 to 92%(Durbin and Durbin 1981). In the adult alewife(Stewart and Binkowski 1986) and Atlantic
herring models (Rudstam 1988), egestion was
assumed to be 16% of consumption.
Data for egestion processes and models are notcommon; most extensive studies have been done
for brown trout Salmo trutta (Elliot 1976a, 1976b).
Egestion has been found to be a function oftemperature and ration (Elliott 1976a). However,
Stewart and Binkowski (1986) found smallchanges in estimated consumption when making
the simplified assumption of egestion being a
constant proportion of consumption in the
bioenergetics model for alewife.
The proportion of consumption egested has beenfound to be low in larval and juvenile clupeids(Kiørboe et al. 1987; Limburg 1994). Arrhenius(1998) used 16%, the value from the adult Atlantic
herring (Rudstam 1988) and alewife (Stewart andBinkowski 1986) models for the proportion of
assimilated ration egested by larval Atlantic
herring. Both Kiørboe et al. (1987) and Limburg(1994) found the percentage of food egested was
10% (by mass). However, Klumpp and vonWesternhagen (1996) found egestion for Atlantic
herring larvae age 8 - 33 days averaged 17.6%(range 13.4 - 25.6%) of ration (Artemia sp.
nauplii) energy content.
Based on the above three studies on larval andjuvenile clupeids (Kiørboe et al. 1987; Limburg
1994; Klumpp and von Westernhagen 1996),
Klumb et al. (in review) chose 0.125 as a first
approximation for the proportion of consumption
egested by larval and juvenile alewife.
Excretion
Excretion (U) = UA*(C – F)
Excretion is modeled as a constant proportion of
assimilation (consumption minus egestion). In theadult alewife (Stewart and Binkowski 1986) andAtlantic herring models (Rudstam 1988), excretion
was assumed to be 10% of assimilation based on
rates measured for brown trout (Elliott 1976b).
Few studies on larval fish excretion have been
conducted. For three species, Blennius pavo,
plaice Pleuronectes platesssa, and Atlantic
herring, Klumpp and von Westernhagen (1996)
found the mean percent of the assimilated ration
excreted was 6.0, 6.6 and 10.7%, respectively.
Due to high mortality for Atlantic herring larvae in
Klumpp and von Westernhagen’s study, Klumb etal. (in review) used the average value of 7.8% for
all three species as a first approximation of the91

percent of assimilation excreted by larval andjuvenile alewife.
Data requirements
There are four data requirements for the
bioenergetics model: 1) diet (in proportions ofprey types), 2) energy density of the predator fish,3) energy density of the prey, and 4) watertemperatures. The bioenergetics model is anindividual based model but can incorporate
populations by multiplying mean weight by
population number.
Diet
Diet information was summarized by Douglas E.
Hay (Pacific Biological Station, Fisheries andOceans Canada) based on recent observations(Hay and McCarter 2001, and older literature suchas Wailes 1936). Depending on the population,
herring diets can be simple or complicated. The
simple story is that herring eat mainly copepod
eggs and nauplii as larvae, copepod adults andnauplii as juveniles and euphausiids as adults.
This over-simplified` story gets messy when thesmaller, non-migratory marginal populations are
examined because they appear to eat a wider
variety of taxa. Herring feed intensely in thesummer months but they also eat during winter. In
all areas, winter diets, although small in relation tototal annual consumption, may be more variable
than summer diets. Perhaps the main point toemphasize, however, is that in southern British
Columbia, most herring feed in shelf waters where
the main food items are euphausiids. Atlantic
herring and sprat Sprattus sprattus diets consistedof 70 - 73% copepods, 12 - 14% Oikopleura, and 9- 12% cladocerans (De Silva 1973).
Adult clupeids feed by filtering or particulate
feeding (Blaxter and Hunter 1982; Janssen 1976).
Janssen (1976) found for alewives that the filter
feeding mode displayed by large alewives (total
lengths > 170 mm) was not size-selective, while
particulate feeders (total lengths 50 - 115 mm)selected zooplankton > 1.0 mm. Transition oflarvae to adult body morphology and feeding
modes occurs at metamorphosis (~35 mm) aftergill rakers and the upper and lower jaws become
developed (Blaxter and Hunter 1982). Although

activity of clupeids may be lower at night (Katz
1978), filter feeders can still feed in darkness(Hettler 1976; Janssen and Brandt 1980; Grabe1996).
Feeding activity of larval herring has been foundto be dependent on densities of copepod nauplii(Munk and Kiorboe 1985) with success a functionof prey size (Hunter and Blaxter 1982). Atlantic
herring (length 25 mm) larvae were able toconsume prey sizes ≥ 1.0 mm (Sherman andHoney 1971, cited in Hunter and Blaxter 1982).
Foraging behavior of Altantic herring larvaechanged with prey size and was related to larval
length by the equation: prey length = 0.027*larval
length (Munk 1992), and attack success was
directly related to relative prey size. Fiksen and
Folkford (1999) included the mouth size of herring
larvae, perception (visual field and reaction
distance), light intensity, and the length, width anddensity of plankton prey when modeling encounterrates and probabilities of successful strikes.
Energy density of predator and prey
Energy density, also called caloric content andenergy content, in bioenergetics models is used interms of wet weight.
Dry-weight data arecustomarily converted (approximated) to wetweight assuming dry weight is 10 - 20% of total
weight. Hartman and Brandt (1995) provided
many equations for estimating energy density fromthe percent dry weight of various marine andfreshwater fish species.
Assuming constantenergy densities or using values that are too highor low can greatly affect bioenergetics model
consumption estimates (Stewart and Binkowski1986).
Energy density (ED) of clupeids has been found tovary seasonally, peaking in fall and declining
through winter (Arrhenius and Hansson 1996;
Flath and Diana 1985; Paul et al. 1998). Age-0
EDs are lower than older fish (Arrhenius andHansson 1996; Flath and Diana 1985; Paul et al.
1998). For age-1 alewife in Lake Michigan, ED in
June and July was 4520 J·g-1, increased in Augustand September to 4729 and 5440, respectively,
then declined to 4729 in April, and 4436 by May(Flath and Diana 1985). For age-0 Atlantic
herring, Arrhenius and Hansson (1996) found ED92

increased from 2600 J·g-1 in mid-July to 4500 J·g-1
in October. The ED of age-0 Baltic Sea sprat
increased from 4000 J·g-1 in August to
approximately 5250 J·g-1 by December (Arrhenius
1998b). October ED of alewife was 5020 J·g-1(Flath and Diana 1985).
Higher energy densities were found for Pacificherring off Alaska (Paul et al. 1998; Foy and Paul
1999) compared to Great Lakes alewives andBaltic Sea clupeids. Age 2+ Pacific herring had
EDs in fall that ranged from 9400 to 10200 J·g-1
and declined over winter to 5200 to 6300 J·g-1 byspring. Females had higher energy densities inboth seasons than males by 200 - 400 J·g-1. Age-0
herring had EDs of 5700 J·g-1 in fall which
declined to 4400 J·g-1 by the following spring(Paul et al. 1998). Equations to predict ED fromstandard length of juveniles by month are provided
in Paul and Paul (1998a). The ED for age-0
captive fasting herring declined 23 J·g-1·d-1 from
Table 2.2.2

December to the end of January (Paul and Paul
1998b).
Energy densities of freshwater and marine
invertebrates can be found in Cummins andWuycheck (1971), while good tables of the caloric
content of marine invertebrates (with references)are presented in Foy and Norcross (1999) and Foy
and Paul (1999). Laurence (1976) provides energy
densities for marine calanoid copepods in the
Atlantic.
Like fish, the energy density forinvertebrates has been found to vary seasonally.
Table 2.2.1 Energy densities (jouls/gram) for
main food items of Pacific herring.
Food Item
Copepoda
Euphausiids (per gram wet weight)
Fish eggs (per gram wet weight)

J·g-1
2580
5020
4520

Existing bioenergetic models.

Reference
Winberg 1956
Kitchell et al. 1974
Elliott 1976b; 1979
Stewart et al. 1983
Rudstam 1988
Kerr and Dickie 1985
Arrhenius 1998
Fiksen and Folkford 1999
Stewart and Binkowski
1986
Hewett and Stewart 1989
Klumb et al. In review
Durbin and Durbin 1983

Comments
General models
extensive early work but reference not that accessible
results of International Biological Program (IBP) workshops and first
paper of the “Wisconsin” bioenergetics model – applied to bluegill sunfish(Lepomis macrochirus) in terms of mass balance
general review of energetics resulting from his extensive work with brown
trout (Salmo trutta)changed the Kitchell et al. 1974 model from mass to energy balance, forlake trout (Salvelinus namaycush)
Clupeid bioenergetics models
Adult Atlantic herring (Clupea harengus)
Age-0 Atlantic herring
Age-0 Atlantic herring
Larval Atlantic herring– Individual based model, which includes
metabolism, ingestion, prey encounter success, and multiple preyfunctional response
Adult alewife (Alosa pseudoharengus)
Age-0 alewife: (only temperatures for the consumption component
differed from the adult model)
Age-0 alewife
Adult menhaden (Breoortia tyrannus):– in terms of energy and Nitrogen932.3 Reflections of factors affecting size-at-age and strong year classes ofherring in the North PacificDouglas E. Hay
Pacific Biological Station, Fisheries and Oceans Canada, 3190 Hammond Bay Road, Nanaimo, British
Columbia, Canada V9R 5K6. E-mail: hayd@pac.dfo-mpo.gc.ca
One approach to the investigation of linkages
between oceanographic process and subsequent
impacts on marine fish populations, isretrospective analyses of age-specific growth rates(size-at-age) from archive collections of scales orotoliths.
This approach can be linked toindependent observations on (1) temporalvariation in abundance, (2) synchrony or
asynchrony of year-class strength among different
populations, or different species, and (3) habitat
requirements of life history stages (eggs, larvae,juveniles, adults) that have different spatial andtrophic characteristics.
Widespread geographic synchrony sometimes
occurs in Pacific herring (Hollowed and Wooster
1995, Hay et al. 2001). An exceptionally strong
year-class occurred in 1977 over a broad andgeographic range (Fig. 2.3.1). It was strong in
northern BC, parts of south-eastern and central
Alaska and the Bering Sea (Hollowed andWooster 1995). This 1977 year-class developed indifferent populations with different spawning
times, with a range of about 3 months from theearliest to the latest mean spawning time. Pacificherring spawn in shallow, inshore inter-and
subtidal waters. In many areas of the Pacific coastof North America, spawn deposition is monitored
and quantified annually. Spawn deposition was
not exceptional in 1977. Therefore, it follows thatin 1977, survival from eggs to the juvenile andrecruit stage, between 1977 and 1980, wasrelatively higher (or mortality was lower) thanmost other years. It also follows that the factor(s)that promoted the strong year-class were widelydistributed in space and time.
Retrospective analysis of archived herring scales(Fig. 2.3.2) from northern BC populations,indicates that individuals of the 1977 year-classwere of normal size, or slightly larger than normal,
relative to samples from other years (Fig. 2.3.3).

After age 4, the relative size-at-age of individuals
in the 1977 year-class declined, and was smaller
than normal, which indicates that growth rate
declined in older individuals. This retrospective
analysis of growth from scale analysis wascorroborated by analyses of catch-sampling data,collected routinely for the last 70 years. The sizeat-age of 3-year-old members of the 1977 yearclass was normal in most areas in 1980, but sizeat-age of older individuals (e.g. age 6 fish
collected in 1983) was smaller than normal (Fig.
2.3.4).
A strong 1977 year-class also occurred in several
other species, including blackcod and lingcod(Hollowed and Wooster 1995). Climate-relatedchanges, but not necessarily increases inabundance, also occurred in other marine species
including salmonids (Beamish et al. 1999) andpollock Theragra chalcogrammus (Ohtani and
Azumaya 1995). Further, there are periods when
there has been synchrony of strong year-classes
among different species in the North Pacific(Hollowed et al. 1987), which is evidence ofenvironmental influence on the production of yearclass strength.
The habitats occupied by age 1 and 2 herring aremainly inshore (Haegele 1997), whereas most of
the older age groups (age 3 and older) tend tooccupy shelf waters. During intensive summer
feeding periods, juvenile herring are found mainly
in shallow, nearshore waters of less that 50 m. In
general, age 1 juveniles occur in shallower waters,
closer to shore, than age 2 herring. In general,
herring form shoals of similar-sized individuals so
the two larger age groups do not mix, althoughboth age groups of juveniles occur in the same
vicinities, herring juveniles are widely dispersed
through all BC coastal waters.
Over the last 70 years in British Columbia (BC),herring stomachs have been examined by different

94people, in different years, at different places and atdifferent herring life history stages. Wailes (1936)
summarized the food of young herring mainly inthe first summer of life. At very young stages,
eggs (ova) and nauplii from various invertebrates
are most important. Copepod nauplii seem todominate the food but food composition variedwith location. The youngest juveniles (age 1) fed
mainly on copepods. Older, larger juveniles took
various zooplankton, with euphausiids being
common. More recent work examined gut datafrom herring juveniles in Georgia Strait, BC,
Hecate Strait and Prince William Sound Alaska(Haegele 1997, Foy and Norcross 1999, Hay and
McCarter 2001). In general the main food for
herring at ages 1 and 2 is copepods. Therefore if
the abundant 1977 year-class ate mainly copepods
at ages 1 and 2, then copepods must have beenabundant in nearshore northern waters, both in1977 and 1978. From our present understandingof herring life history, there is little opportunity fortrophic interaction (i.e. direct density-dependentcompetition for food) between age-classes: either
among juveniles (ages 1 versus age 2) or between
juveniles versus adults (age 3+ and older). In BC
waters, probably the first opportunity for directinteraction occurs during the third winter of life, atage 2+, when (BC) herring start to mature sexuallyand join the adult stock. At this time, however,
winter feeding is minimal and growth is slight.
The observations above can be summarized asfollows. In 1977, and some other years, we seethat strong year-classes can develop over broad
areas of time and space. They develop in years
when spawn deposition is normal, and sometimes
even lower than normal. Further, sometimes they
can be synchronous over broad areas of time andspace. Synchrony may develop in other species.
Retrospective analysis of herring scales indicatessuperior juvenile growth among strong cohorts,
but decreased growth during older adult stages (in1977). Strong year-classes can arise in years ofnormal or modest spawn deposition.
These
observations indicate that survival, between theegg and recruit stages, is enhanced.
Suchenhanced survival must occur during the juvenilestages that consume mainly copepods in nearshore
habitats.
Therefore strong year-classes may
develop as a consequence of changes in thesehabitats.

If lower mortality of early life history stages ispart of the explanation for the formation of the
1977 year-class - or other year classes, why did
this happen? Presumably it must reflect decreases
in mortality by starvation, disease or predation? In
1977, starvation seems unlikely, because juvenile
growth was enhanced compared to other years.
We have no evidence to suggest that disease
routinely limits survival. Rather, outbreaks seem
episodic, and this could explain years withexceptionally bad year-classes, but not the reverse.
A decrease in predation, between the egg/larvalstages and pre-recruit stage could occur if there
were (i) fewer predators, or (ii) if the predators‘switched’ or decreased predation on herring for adifferent prey species. Were predation rates onjuvenile herring lower in 1977 and 1978? We
have no data on this, but we observe that somecommon herring predators (lingcod and blackcod
and some piscivorous salmon) also had strong
1977 year-classes. Therefore it seems improbablethat there was a decrease in the potential
community of herring predators between 1977 and1980.
From the observations and reasoning above, we
conclude that the most parsimonious explanation
for the development of the strong 1977 year-class
was a general decrease in predation of juveniles
because the main herring predators had alternate
prey. Such a reduction in predation could occur
through predator switching during early life
history stages - specifically, predators of herringchose to feed on an alternate food source. If thisalternate food source was an unusually abundant
supply of copepods, available both to the juvenilesof herring and their predators, this could explain
our observations.
Specifically if predators
preferentially switched to copepods, instead ofherring juveniles, the consequence of a substantial
increase in copepod availability would be bothenhance survival and growth of juvenile herring.
If the cause(s) of the strong 1977 year-class was
similar in all geographic areas where it occurred,from northern BC to the Bering Sea, and if thecause was from decreased predation associatedwith availability of an alternate food source, thenclearly the factors which promoted this alternate
food source were widespread. There have been95

AlaskaBC

some suggestions (Hollowed and Wooster 1992;
Polovina et al. 1995) that there can be suchlinkages between offshore oceanographic changes
and changes in productivity or food abundance on
shelf and inshore waters, resulting from mid-gyrechanges, but these are not well understood. If
there were such a relationship, the impact of an

E dge

A nnulus 2

A nnulus 1
Focus

b

a

3

annulus 22annulus 1
1Year-class

Fig. 2.3.3
Retrospective analyses of scale
growth of 5-year-old-herring from archived
collections of scales from northern BC. Scale
growth, corresponding to juveniles at age 1 (bluerectangles) and age 2 (red rectangles), as estimatedfrom comparison of focus: annuli distances, wasnormal in the 1977 year-class. The 1977 yearclass is shown in yellow. The boxes and verticallines represent the range and 95% confidence
limits about the mean.
1977 year- class

250

Age 8
Age 7
Age 6
Age 5

200

Age 4
Age 3
Age 2

150

51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

c

Age 5 only

1977 year class

67 70 71 74 75 77 79 80 82 85 88 90 92

Median (mm) length at age

Fig. 2.3.1
Approximate locations of the strong
1977 year-class, indicated by arrows. Red arrows
with dark outlines show locations where the 1977
year-class made up 70% or more of the spawning
population as age 3 in 1980, or age 4 in 1981.
Plain red and yellow arrows show populations
where the 1977 year-class represented over 50%
and 40% of the populations, respectively.

Focus to annulus1 (mm)

4

Year

Fig. 2.3.2
A herring scale, showing the focus(start of growth) and the first and second annuli.
Retrospective indices of age-specific growth rates
during the first year (green bar a) and second year
(yellow bar b) was determined by direct
measurement of scales.

Fig. 2.3.4
Comparison of the size-at-age of the
1977 year-class with those of other years fromcatch sampling data collected in northern BC. The
1977 year-class (large dark circles) was normal (orslightly larger than normal at age 2). Thereafter,
the relative size-at-age, relative to previous yearclasses, decreased until age 8.

96

abundant production of zooplantion, specificallycopepods, could explain both enhanced growthand year-class survival in herring and otherspecies. There is a precedent for assuming that anabundant source of an alternate zooplankton preyspecies can reduce predation on herring. Wareand McFarlane (1995) showed that increased
euphausiid production resulted in a decreased hakepredation on adult herring off the west coast of

Vancouver Island.
Similar mechanisms
mightoperate at the juvenile stages, so factors
promoting a strong year-class of herring mightalso support strong year-classes of other species,
leading to synchrony between unrelated species
such as such as blackcod and lingcod. Again, theanswer is a tentative yes. Both of these species
have early life stages (first several years of life) innearshore waters.

2.4 Review for Pacific saury (Cololabis saira) study under the VENFISH projectShin-ichi Ito1, Yutaka Kurita1, Yoshioki Oozeki2, Satoshi Suyama3, Hiroya Sugisaki1, Yongjun
Tian2
1
2
3

Tohoku National Fisheries Research Institute, 3-27-5 Shinhamacho, Shiogama, Miyagi 985-0001,
Japan. E-mail: goito@affrc.go.jp, sugisaki@mgy.affrc.go.jp
National Research Institute of Fisheries Science, 2-12-4 Fukuura, Kanazawa-ku, Yokohama,
Kanagawa 236-8648, Japan. E-mail: oozeki@affrc.go.jp
Hachinohe Branch, Tohoku National Fisheries Research Institute, Same, Hachinohe, Aomori 0310841, Japan. E-mail: suyama@myg.affrc.go.jp

VENFISH (Comprehensive study of the Variationof the oceanic ENvironment and FISH populationsin the northwestern Pacific) project was started in
April 1997 and will end in March 2002. This
project has been supported by Japan Agriculture
Forest Fisheries Agency. The aim of this projectis clarification of bottom-up control process for
Pacific saury and walleye pollock in the
Northwestern Pacific. More than 20 scientists
from National Fisheries Research Centers at
Hokkaido, Tohoku, Yokohama and Shimizu, andHokkaido University and Tohoku University
joined this project.
The VENFISH team is composed of 5 teams andthere are primary production, zooplankton and fishteams. The fish team is composed of Pacific saury
and walleye pollock groups. Between these threeteams there is a plankton ecosystem model teamand a fish population model team. In this report
we will note our studies of saury, which is only
one portion of this project.
The main target area of the VENFISH project iseast of 160ºE in the northwestern Pacific, and inthat region there is a warm Kuroshio current and acold Oyashio current. Between these two westernboundary currents, there is a mixed water region,and in that area many eddies are detached from theKuroshio and Oyashio and make very complicatedenvironments. The saury spawning starts in themixed water region in autumn, moves to the
Kuroshio area in winter, and moves back to themixed water region in spring (Fig. 2.4.1) (Odate
1977; Watanabe and Lo 1989; Watanabe et al.
1997). Juveniles are advected to the Kuroshio
extension region, then grow and migrate to the
Oyashio region through the mixed water region for
feeding. After sufficient feeding they migrateback to the Kuroshio region for spawning. On thesouthward migration, they are fished in theJapanese coastal zone. We will briefly report the
new findings for Pacific saury in the later sections.
Feeding habitat
The feeding habitat of Pacific saury (Cololabis
saira) changes according to the life stage and thelocation. Larvae smaller than 15 mm mainly feedson Oncea and Oitona sp. (Nakata and Koyama
2002), whereas larvae and juvenile larger than 15mm prefer Calanus sp. Young saury which
migrate to the mixed water region mainly feed on
Euphausia pacifica. In the Oyashio region they
feed mainly on Euphausia pacifica and
Neocalanus cristatus and the ration becomes the97

maximum in this season. On the way of theirbackward migration, they feed Euphausia pacifica
and Sagitta elegans, but the ration decreases to theminimum. In the spawning area they feed on
calanoid copepods and the ration is higher than inautumn (Sugisaki and Kurita, in preparation;
Kurita and Sugisaki, in preparation).

45

Feeding groundin summer
Fishing groundin autumn

40

Oyashio Front

Mixed water region
Migration spawning groundin spring & autumn
35

30
130

Kuroshio Extension
spawning groundin winter

Migration

135

145

140

150

155

160

165

Fig. 2.4.1
Schematic picture of Pacific saury(Cololabis saira) life history. Spawning starts in
September and continues until June, shifting
location from the mixed water region andKuroshio region. The main spawning season is
winter. Juveniles are advected to the Kuroshio
extension region and migrate to the Oyashio
region through the mixed region for feeding.
After sufficient feeding they migrate back to the
Kuroshio region for spawning. On the southwardmigration, they are fished in the Japanese coastal
zone.
Spawning density
Kurita and Sugisaki (in preparation) surveyed theseasonal change of the saury distribution and theratio of mature stage in the three regions. In early
autumn, half of the saury occur in the Oyashio
region and they are immature. In winter almost allof the saury are in the Kuroshio and they are
mature. In spring, half of the saury exist in the
Kuroshio and most of them are mature. But theother half occurs in the mixed water region andonly about 70% of them are mature. These results
show that the most important area is the Kuroshio
region and the most important season is winter forthe saury spawning.
Kurita and Sugisaki (in preparation) estimated thespawning interval and batch fecundity. Using

these values and ratio of mature saury to the total,they estimated the spawning density for eachseason.
Their result showed that the mostimportant season for spawning is winter.
Larvae and juvenile
Many studies have been done about larval andjuvenile saury (Watanabe et al. 1997; Oozeki and
Watanabe 2000; Oozeki and Watanabe, in
preparation). Using widely sampled field data,
Watanabe et al. (1997) and Oozeki and Watanabe(2000) estimated the production of hatched larvaein each season since 1990 to 1997. The averagevalue for 8 years showed the highest value inautumn and the lowest value in spring. They also
estimated the growth rate and mortality of larvaeand production of juveniles. Growth rate showeda maximum in autumn and a minimum in spring.
Mortality was highest in autumn and lowest in
spring. As a result, the production of surviving
juveniles showed a maximum in spring and aminimum in autumn. But the fluctuation ofjuvenile production in spring is very high andstable in winter. Watanabe et al. (1997) suggested
that the stable winter juvenile might contribute tostable recruitment and middle size saury landings
in autumn. Also Watanabe and Lo (1989) pointed
out that winter was the most active spawningseason using larval catch data during 1973-1986.
Oozeki and Watanabe (2000) conductedlaboratory incubation experiments on saury eggs.
They reared same age larvae at three different
temperatures and observed growth rate. This wasdone for three different age larvae (9, 20, 30 days)
and the dependency of growth rate on age was also
tested. The result showed that the growth rateincreased linearly with temperature and also
increased with age. Analysis of the otolithincrement and the knob length of the larvae
showed the possibility of the estimation of growthrate of saury juveniles from the otolith field data.
Then, they estimated the instantaneous growth ratefrom otolith field data and analyzed therelationship between the recent growth rate andoceanic environments (Oozeki and Watanabe: inpreparation). Their result showed that the SST
and food density affected larval growth during theearly stages, and SST and chlorophyll become
more important in the later stage.

98

Growth rate of adults
Suyama et al. (in preparation) analyzed thepresence of a hyaline zone in the otoliths ofPacific saury. Usually the size decomposition isdone by knob length, but sometimes it is difficult
to divide them only from body length information.
On the other hand, the otoliths of large size saury
have the hyaline zone whereas the small andmiddle sizes do not. They analyzed the existence
the hyaline zone and found out that the large andmiddle size cohort can be decomposed by theboundary of 50% existence ration of hyaline zone.
Using this definition they decomposed the large
and middle cohort and analyzed the inter-annual
variability in the growth of each cohort. The
middle size fluctuated between 264 and 286 mm,
and the large size fluctuated between 303 and 314
mm, and the fluctuation was larger in the middle
size.
For example, the growth increment of the large
and middle cohort from July to November 1999,
was 11.3 and 19.3 mm respectively. On the other
hand, they increased to 12.5 and 31.3 mm
respectively in 2000. This result suggests that thegrowth rate of the large size cohort is more stably
estimated compared to the middle size cohort.

hyaline zone. According to his scenario, saury
which are born in the earlier season spawn in thefirst winter and also in the second winter. But thelater spawned saury do not spawn in the first year
and spawn in the second year.
Energy for migration and spawning
Kurita (personal communication) analyzed
seasonal variation of lipid and protein content in
30 cm knob length saury. The protein content did
not vary much but lipid variation showed verylarge variability. The average lipid content isabout 40 g in summer. In winter, which is theactive spawning season, mature saury contained
little neutral lipid. Moreover, protein seemed to be
utilized as energy sources because the sum ofprotein and water content was constant. From thisresult he concluded that saury need to feed inorder to spawn eggs in the Kuroshio region.
Thus, the environment may be very important forthe saury reproduction in the Kuroshio region.
From the energy balance between the food nutrientand egg production, he estimated that about 35.6%
of total assimilated energy was used for winter egg
production in the Kuroshio.
Population dynamics model for Pacific saury

Growth rate between juvenile and large size
Using the hyaline zone information from theotolith it is possible to estimate the growth rate ofyoung and adult saury, but the growth rate
between juvenile and young saury is very difficult
to estimate because of the existence of the hyalinezone. We cannot count the increment of theotolith because the increment is unclear in thehyaline zone. So, we cannot determine the age ofadult saury.
For
thisproblem,
Kurita(personal
communication) developed a new method toestimate the hatch date from the age at which theotolith increment width reached a second
maximum. It became possible to estimate the ageof saury using this method even if there is ahyaline zone. He estimated the hatch date of thelarge size saury and developed a new scenario of
the life history of Pacific saury combined with theinformation of the growth of the saury with noTian et al. (2002b) analyzed the interannual
variability of the saury stock using a population
dynamics model. In his model there are twocohorts. One is a cohort spawned during autumn winter and the other is spawned during winter spring. The life span of the saury was assumed to
be two years, and as a result the large size saury
included both cohorts. The governing equations
were growth rate, population, fishing effort andreproduction equations.
In the population
dynamics the mortality included environmental
effects. As environment factors they adapted SST
in the Kuroshio Extension zone (KE SST) and SOI(Southern Oscillation Index) according to theresult of Tian et al. (2000a).
The results showed that the effect of KE SST was
important to the longer-period variability, and the
SOI effect was important to both the longer-period
and inter-annual variability.

99Conclusion
Under the VENFISH project, much has beenlearned about Pacific saury and a new life historyof the saury was proposed. But information about
the time between the juvenile and small saury
stages are still limited. In the future more study isneeded on these stages.

A population dynamics model was constructedunder VENFISH and the effect of KE SST and
SOI was tested.
But in that model theenvironment influenced only mortality. In the
future we should include the environmental
influence on production and clarify the bottom-up
control mechanism of Pacific saury.

2.5 Formalization of interactions between chemical and biologicalcompartments in the mathematical model describing the transformation of
nitrogen, phosphorus, silicon and carbon compounds
Alexander V. Leonov1 and Gennady A. Kantakov2
1
2

Institute of Oceanology, Russian Academy of Sciences, 36 Nakhimovsky Ave., Moscow, 117851,
Russia. E-mail: leonov@sio.rssi.ru
Sakhalin Research Institute of Fisheries and Oceanography, 196 Komsomolskaya St., YuzhnoSakhalinsk, 693023, Russia. E-mail: okhotsk@sakhniro.ru

In the significant part of the ecological models
used for studying the joint dynamics of the
microorganism biomasses and biogenic substanceconcentrations in the natural waters, several mostimportant biological functions are formalized.
They are connected with the consumption ofbiogenic substances (UP) by microorganisms,
excretion of the metabolic products (L) by them,the microorganism mortality (S) and grazing (G)by microorganisms of higher trophic levels. Thechange of the microorganism biomass in thecourse of time (dB/dt) in the ecological models, asa rule, is represented by the following structural
equation:(2.5.1)dB / dt = (UP - L - S) * B - G * В*here B* is the biomass of microorganisms fromthe higher trophic level, and due to grazing they
have an influence on the development and activityof the considered microorganism group B; UP, L,
S, and G are specific rates of the biogenic
substance consumption, the metabolic productexcretion, the mortality of microorganisms B andtheir grazing by B*, respectively (day-1).
Biomasses B and B* are calculated in the units ofbiogenic elements (N, P, C or Si).

The simulation of processes of the substrateconsumption by microorganisms
For the simulation of processes of the substrateconsumption by microorganisms (bacterio-, phytoand zooplankton), the equation of MichaelisMenten-Monod is traditionally used:(2.5.2)UP = K (T, L) * Ci / (Km + Ci)

where UP is the growth rate of the microorganism
biomass (or the substrate uptake), day-1; Ci is theconcentration of concrete substratum, mg/l; Km is
the Michaelis constant, mg/l; K (T, L) is the
maximum growth rate of the microorganism
biomass (or the substrate uptake) corrected to thetemperature (T) and radiation (L) conditions in thewater environment, mg/(l day).
Thus, fordescription of the process of the substrate uptake
by one group of microorganism (by bacterio-,
phyto- or zooplankton) it is necessary to estimate
the values of two coefficients - K (T, L) and Km.
Using this equation form for the description of theconsumptionofseveralsubstratabymicroorganisms, means that the process of thesubstrate uptake is described independently ofeach other for any substrate, and in this case, thevalues of the rate constants for the consumption ofeach substrata should be evaluated. If the numberof such substrata reaches five (ammonia, nitrites,100nitrates, phosphates, silicates) then the number of
evaluated coefficients should be equal to ten.
The form of equation (2.5.2) with somemodifications is used for describing the processesof the substrate consumption by microorganisms
in the models developed by PICES MODEL Task
Team for the studying of chemical and biologicalcompartment dynamics in the marine environment.
However, in the marine environment, the substrateconcentrations are always little and therefore it isvery difficult to describe the dynamics of thebiomasses and substrate concentrations even for
one season. Frequently the very task of thesimulation of the chemical and biologicalcompartment dynamics in the marine systems is arather difficult labor-consuming or eveninsurmountable problem.

(2.5.4)

where BN is the biomass of the studied group of
microorganisms, in units of N, mg N/l; CN is theconcentration of N fractions consumed by thesemicroorganisms, mg N/l.
If there are several N-containing substrata in thewater environment (for example, ammonia NH4,
nitrites NO2, and nitrates NO3) and thesesubstrates are interchangeable and may beconsumed by the microorganism (let us mark it as
F and taking into account that the biomass isexpressed in units of N, it may be written as FN),
then the expression for CN can be represented in
the form of the pool on N (PoolFN) for the studiedgroup of microorganism:(2.5.5)Here we present the logic of the simulation of thesubstrate
consumption
processbythemicroorganisms that is used for development of
the model describing the transformation of N, P, C
and Si compounds in the polysubstrateenvironment (Leonov and Saposhnikov 1997).
First we shall transform the equation (2.5.2)
subdividing the terms of equation in the numerator
and the denominator on Ci. As a result, we obtainthe following equation:(2.5.3)UP = K (T, L) /(1 + Km / Ci)

The analysis of literature shows that the value of
Km in different examples of using equation (2.5.2)
for describing of processes in the natural waterschanges by 2-3 orders.
Consequently, theconvincing arguments of the application of theequation (2.5.2) for describing the substrateconsumption processes in the marine ecosystems
is clearly insufficient (large number of coefficientsfor the polysubstratal environment and the large
variability of the coefficient Km). The value of thecoefficient Km for the marine ecosystems may becompared with the values of the microorganism
biomasses. Therefore we have all reasons to use,
instead of the coefficient Km, the value of thebiomass in the units of biogenic element (N, P, C
or Si) from which biomass can be evaluated. If
the biomass is considered in N, then the equation(2.5.3) can be written as:

UP = K (T, L) /(1 + ВN / CN)

PoolFN = d(1) * NH4 + d(2) * NO2 +d(3) * NO3

Here the coefficients d(i) show preferences in the
consumption of each substrate by the
microorganism for this N-substrates (NH4, NO2,
and NO3). Value of the coefficients d(i) for each
substrate can change from 0 to 1, and their sum forthe selected set of substrata is 1.
How are the values of coefficients d(i) evaluated?
It is known from literature that the phytoplankton
consumes more preferably ammonium N thanother mineral forms. The nitrate N is in second
place. So, in the first approximation, we canassign the values of the preference coefficients inthe uptake of indicated substrates by the studiedgroup of phytoplankton: d(1) = 0.5; d(2) = 0.2;d(3) = 0.3. Inserting the equation (2.5.5) into theequation (2.5.4), we obtain:(2.5.6)UPFN = K (T, L) /(1 + ВN / PoolFN)

or(2.5.6a) UPFN = K (T, L) /(1 + ВN /(d(1)*NH4 +
d(2)*NO2 + d(3)*NO3))
The general rate of the N-containing substrata
consumption, UPFN, is composed of the rates of
the consumption of the individual substrates:(2.5.7)UPFN = UPNH4 + UPNO2 + UPNO3

101

Let us write down the equations, which describe
the consumptions of individual substrates (NH4,
NO2 and NO3) by phytoplankton taking into
account that in the water environment several
substrates are interchangeable on N, as indicatedby the equation (2.5.6a). Making elementary
algebraic conversions, we shall obtain theequations, which describe the consumption of eachstudied substrates by the phytoplankton:(2.5.8)UPNH4 = K (T, L) * d(1) * NH4 / (PoolFN+ ВN)(2.5.9)UPNO2 = K (T, L) * d(2) * NO2 / (PoolFN+ ВN)(2.5.10) UPNO3 = K (T, L) * d(3) * NO3 / (PoolFN+ ВN)
The suggested form of the description of theinterchangeable substrates by the microorganism
assumes that the rates of the consumption of each
substrate will be compared only in such a case,when the product of substrate concentrations totheir preference coefficient will be close. With the
maximum rate will be consumed that substratum,
for which the product of its preference coefficientto the concentration will be the greatest (in thiscase, from of three given substrates). This form of
the equations for the consumption of theinterchangeable substrates by microogranism (inparticular, by phytoplankton) gives the possibilityof switching for the intensive consumption by thehydrobionts only of those substrata whose
concentrations to these are greatest in thecomparison with other substrata. It gives apossibility for the water environment to restore thepool of those substrata, which in the process of thebiomass growth descend to the smallest(sometimes critically small) values.
Thisphenomenon in the description of the processes ofincreasing of the biomass and substrateconsumption is impossible by equations the
traditionally used for the simulation of marine
ecosystems, in which the substrate consumption
by different groups of microorganism is assigned
independently of each other.
Let us consider the case, when there are several
substrates as the interchangeable (as it was
examined above), so also not interchangeable, inthe water environment for the phytoplankton. If
we want correctly describe in the model thesubstrate uptake processes then we should
remember the basic Odum’s postulate thateverything is interrelated in the natural water
environment. The requirements of phytoplankton
in P cannot be compensated by N or Si
compounds, and vice versa.
Therefore thecompounds of different biogenic elements cannot
be considered as interchangeable for the formationof the microorganism biomass and the kinetics of
the uptake substrate processes should beformulated with point of view their mutual
influence on each other and not theirinterchangeability.
Taking into account these reasons, let us write
down the equation (2.5.6) for the rate of biomass
growth (or the different substrate uptake) for theconditions of the combined influence of N and P
compounds on the biomass of the considered
microorganism group (for example, thephytoplankton) keeping the logic of all foregoingreasons. Then we obtain, that
(2.5.11) UPF = K (T, L) /(1 + ВN / PoolFN + ВР /PoolFР)Here BP is the biomass in units of P, mg P/l;
PoolFp - the pool of the P substrates, mg P/l, that
may be consumed by the phytoplankton, and thesesubstrates are the dissolved mineral (DIP) andorganic (DOP) forms of P:(2.5.12) PoolFp = d(4)*DIP + d(5)*DOP
In this case the total rates of the uptake of N and P
compounds by the given group of microorganisms
are represented as:(2.5.13) UPFN = UPNH4 + UPNO2 + UPNO3(2.5.14) UPFP = UPDIP + UPDOP
Accordingly to the same logic, let us formulate
equations for describing the individual substrates
taking into account the influence of each of themon the kinetics of the formation of biomass and thesubstrate uptake being oriented toward generalequation (2.5.11):

102

(2.5.15) UPNH4 = K (T, L) * d(1) * NH4 / MF
(2.5.16) UPNO2 = K (T, L) * d(2) * NO2 / MF
(2.5.17) UPNO3 = K (T, L) * d(3) * NO3 / MF
(2.5.18) UPDIP = K (T, L) * d(4) * DIP / MF
(2.5.19) UPDOP = K (T, L) * d(5) * DOP / MF
where
(2.5.20) MF=PoolFN * PoolFP + BN * PoolFP + BP
* PoolFN.
When the joint consumption of N, P, and SI
compounds is considered for the same group of
microorganism, the equation (2.5.11) takes theform:(2.5.21) UPF = K (T, L) /(1 + ВN / PoolFN + ВР /PoolFР + Вsi / PoolFSi)
where(2.5.22) PoolFSi = d(6) * DISi
and DISI is the content of dissolved inorganic
silicon, mg Si/l.
In accordance to the accepted logic for theformulations of kinetic dependences, the equations
describing the individual substrate uptake and theirmutual influence on each other are written in thefollowing form:(2.5.23) UPNH4 = K (T, L) * d(1) * NH4 / MF1(2.5.24) UPNO2 = K (T, L) * d(2) * NO2 / MF1(2.5.25) UPNO3 = K (T, L) * d(3) * NO3 / MF1(2.5.26) UPDIP = K (T, L) * d(4) * DIP / MF1(2.5.27) UPDOP = K (T, L) * d(5) * DOP / MF1(2.5.28) UPDISi = K (T, L) * d(6) * DISi / MF1
where(2.5.29) MF1 = PoolFN * PoolFP * PoolFSi + BN *PoolFP * PoolFSi + BP * PoolFN * PoolFSi+ BSi * PoolFN * PoolFP

A similar form of equations may be used for any
functional group of microorganism taking into
account any the variety of the substrate assortment
including the components of the water
environment pollution (for example, oil products).
The substrate assortment for the organisms ofhigher trophic levels is higher than for theorganisms of lowest trophic status.
In thisassortment fall the dissolved and particulatedorganic compounds of biogenic elements,
including biomasses of certain microorganisms
and detritus.
The equation for the term G (the specific grazing
rate of the microorganism from the lowest trophic
levels by the organisms of higher levels) isconstructed, on the basis of the presented above
principles considering the high-constituent natureof water environment and the mutual influence of
the uptake of individual substrates on each other inthe process of the microorganism biomass growth.
The value of the maximum growth rate of the
microorganism biomass (or the substrateconsumption), K (T, L) should be corrected to theconditions on the temperature and for light (for theplanktonic organisms) in the water environment.
The analysis of ecological models existing atpresent shows that there are many methods ofcarrying out a similar correction.
In the model of the transformation of nitrogen,phosphorus, silicon and carbon compounds thetemperature dependence is considered by the
exponential function, which differs for thedifferent groups of microorganisms in the slope
and the optimum values of temperature. The
dependence of the plankton biomass growth as afunction of light conditions in water environment
is considered by the traditional functions that areused at the simulation of the processes ofphytoplankton photosynthesis and daily verticalmigration of zooplankton.
Formalization of the excretion processes ofmetabolic products by microorganisms
At first stages of mathematical simulation model
development as the independent scientificdirection in the studies of the natural aquasystems
state, this important biological function of103

microorganisms was not considered at all. At
present time, in the majority of the cases the
specific rate of the metabolic excretion by
microorganisms (L) is formulated in the ecologicalmodels by the simplest method and, as a rule, it isrepresented in the form of a constant quotas (α)
from the UP function:(2.5.30) L = α * UP
The experience of the experimental research of the
microorganism population dynamics and thesimulation of the conditions for the biomass
growth shows that the excretion fraction of the
products of metabolic exchange in different
microorganisms differs very substantially, and it
can change considerably in the process of thebiomass growth in each group of microorganisms.
This fact was taken into account, and during thedevelopment of the mathematical model of thetransformation of nitrogen, phosphorus, silicon,
and carbon compounds the different forms of theexpression of the excretion fraction of metabolic
products changing in the time were checked. It
was found the form of equation for α that mostcompletely consider the special features of the
microorganism biomass growth, and it isformulated as the dependence from the specificrate of the substrate uptake, UP:(2.5.31) α = a * UP / (1 + b * UP) + (1 - a / b)where a and b are constants (moreover a < b)whose values determine the nature of change inthe excretion fraction in the dependence on thevalues of the total substrate uptake by consideredgroup of microorganisms.
The first term of the equation (2.5.31) shows theforming quota of the metabolic excretion ofsubstance in the favorable on the nutrient
conditions of the water environment, when valuesof UP are significant.
The second term of the equation (2.5.31) shows
the quota of the metabolic excretion at thesubstrate deficite when values of UP become
minimum.

With the values of coefficients a and b can be
reproduced the significant spectrum of theconditions for the microorganism biomass growthwhich can be evaluated in the units of different
biogenic elements (N, P, C or Si).
Formalization of the
microorganism mortality

processesoftheThe processes of development and growth of the
microorganism biomass are continuous with theprocesses of the internal losses of biomasses (S).
It is possible to assume that the natural
physiological losses of the biomasses of any groupof microorganisms compose 5-10% of the total
biomass although this problem remains
insufficiently studied experimentally for all
microorganism groups. In the process of the
microorganism mortality, the detritus (or the deadsuspended matter) is formed in the water
environment. The biogenic substances containing
in it are actively included in turnover by bacteria
and zooplankton which transform detritus into thelabile nutrients well assimilated by othermicroorganisms. Under the conditions of reduced
temperatures, the detrital links become the mostimportant in the nutrition and growth of the
populations of fishes.
At the first ecological models, the microorganism
mortality S is not taken into account at all or only
natural physiological biomass losses are
considered.
The modern ecological models
include the natural internal biomass losses andtake into account losses inevitable at thestimulation of the biomass growth processes. It
may be differently formulated.
In themathematical model of the transformation of N, P,
Si, and C compounds this important biological
function is described by the equation:(2.5.32) SN = g (1) + g (2) * ВN / UPFN
where g (1) and g (2) are constants describing theprocesses of the natural biomass losses andmortality
depending on the conditions ofactivating the growth, respectively. If the biomassof the microorganism group is evaluated in theunits of different biogenic elements (N, P, C or Si)
then respectively for each case their specific ratesof the internal losses of biomasses are evaluated104with the use of specific values of coefficients g (i),values of biomasses B and rates of the substrateconsumption UP.
The set of model coefficients applied in two casestudies (for the Okhotsk Sea (Leonov andSapozhnikov 1997) and Caspian Sea (Leonov andSrygar 1999) for the simulation of microorganism
dynamics is presented in Table 2.5.1.
Thus, in the mathematical model of thetransformation of N, P, Si, and C compounds, theinteractions between chemical (concentrations ofbiogenic substances) and biological (biomasses the
microorganisms
bacteria,
phytoandzooplankton) compartments are considered andreproduced the most important biologicalprocesses of the substrate uptake, excretion of themetabolic products and mortality of the
microorganisms. As a result of these processes,
the turnover of chemical substances (organic andmineral) are performed in natural marine
ecosystems. The special feature of this model isthe formalization of the important biological
functions (the excretion of the metabolic products
and mortality of microorganisms) in a dependence
on the consumption of different biogenic

substances by microorganism. These biogenic
substrates are subdivided on interchangeable (onone biogenic element) and not interchangeable (onthe different biogenic elements). The used form ofequations for the description in this model of themost important biological functions serves as theexample for the formalization of the processes of
the internal regulation (self regulation) of the
microorganism activity within the ecosystems.
The account of a similar internal regulation
mechanism of the microorganism activity makes
this model sufficiently resistant and allow us toapply it without the significant correction of theparameters in the study the aqueous ecosystems
which essentially differ in the environmental
conditions (temperature, radiation, water regime,
transparency).
There are several positive
experiences in the application of this model tostudy the special features of the ecosystem
functioning of the Sea of Okhotsk (Leonov andSapozhnikov 1997) and Caspian Sea (Leonov andStygar 1999). The first results are also obtainedon the simulation of the intraannual dynamics ofbiogenic substances in the ecosystem in La
Perouse Strait and Aniva Bay (Sea of Okhotsk)(Pischalnik and Leonov 2002).

Table 2.5.1 Values of model parameters used for description of biological compartment dynamics in
the Sea of Okhotsk and the Caspian Sea.
Case study 1 - The Sea of Okhotsk

Case study 2 - The Caspian Sea

Heterotrophic bacteria (B)Maximum growth rate: K=1.0
Preference coefficients for substrate uptake:for C-containing substrate: dDOC=1;
for Si-containing substrate: dDOSi=0.6; dDISi=0.01
dSID=0.39;
for N-containing substrate: dDON=0.6; dND=0.4
for P-containing substrate: dDOP=0.4; dPD=0.6
Excretion activity:
for C substrate: aC=0.05; bC=0.09
for Si substrate: aSi=0.05; bSi=0.088’
for N substrate: aN=0.05; bN=0.087
for P substrate: aP=0.05; bP=0.09
Mortality coefficients:
for C substrate: g(1)C=0.04; g(2)C=0.04
for Si substrate: g(1)Si=0.045; g(2)Si=0.05
for N substrate: g(1)N=0.035; g(2)N=0.035
for P substrate: g(1)P=0.055; g(2)P=0.055

Heterotrophic bacteria (B)Maximum growth rate: K=0.75
Preference coefficients for substrate uptake:for C-containing substrate: dDOC=1;
for Si-containing substrate: dDOSi=0.6; dDISi=0.01
dSID=0.39;
for N-containing substrate: dDON=0.6; dND=0.4
for P-containing substrate: dDOP=0.4; dPD=0.6
Excretion activity:
for C substrate: aC=0.05; bC=0.088
for Si substrate: aSi=0.05; bSi=0.088
for N substrate: aN=0.05; bN=0.1
for P substrate: aP=0.05; bP=0.086
Mortality coefficients:
for C substrate: g(1)C=0.03; g(2)C=0.025
for Si substrate: g(1)Si=0.045; g(2)Si=0.05
for N substrate: g(1)N=0.028; g(2)N=0.03
for P substrate: g(1)P=0.045; g(2)P=0.05

105

First phytoplankton group (F1-diatom algae)
Maximum growth rate: K=2.5
Preference coefficients for substrate uptake:for Si-containing substrate: dDOSi=0.3; dDISi=0.7
for N-containing substrate:dNH4=0.025; NO2=0.025
dNO3=0.9; dUR=0.05
for P-containing substrate: dDOP=0.3; dDIP=0.7
Excretion activity:for Si substrate: aSi=0.051; bSi=0.052
for N substrate: aN=0.05; bN=0.053
for P substrate: aP=0.05; bP=0.065
Mortality coefficients:
for Si substrate: g(1)Si=0.0; g(2)Si=0.08
for N substrate: g(1)N=0.0; g(2)N=0.02
for P substrate: g(1)P=0.0; g(2)P=0.02

First phytoplankton group (F1-diatom algae)
Maximum growth rate: K=2.5
Preference coefficients for substrate uptake:for Si-containing substrate: dDOSi=0.3; dDISi=0.7
for N-containing substrate: dNH4=0.2; dNO2=0.05
dNO3=0.7; dUR=0.05
for P-containing substrate: dDOP=0.05; dDIP=0.95
Excretion activity:for Si substrate: aSi=0.051; bSi=0.052
for N substrate: aN=0.05; bN=0.052
for P substrate: aP=0.05; bP=0.055
Mortality coefficients:
for Si substrate: g(1)Si=0.04; g(2)Si=0.03
for N substrate: g(1)N=0.05; g(2)N=0.049
for P substrate: g(1)P=0.05; g(2)P=0.07

Second phytoplankton group (F2-peridinium
algae)
Maximum growth rate: K=1.8
Preference coefficients for substrate uptake:for N-containing substrate: dNH4=0.15; dNO2=0.05
dNO3=0.2;
dUR=0.6
for P-containing substrate: dDOP=0.4;
dDIP=0.6
Excretion activity:
for N substrate: aN=0.049; bN=0.0495
for P substrate: aP=0.049; bP=0.053
Mortality coefficients:
for N substrate: g(1)N=0.0; g(2)N=0.05
for P substrate: g(1)P=0.0; g(2)P=0.1

Second phytoplankton group (F2-green algae)
Maximum growth rate: K=2.5
Preference coefficients for substrate uptake:for N-containing substrate: dNH4=0.2;
dNO2=0.05
dNO3=0.7;
dUR=0.05
for P-containing substrate: dDOP=0.05; dDIP=0.95
Excretion activity:
for N substrate: aN=0.049; bN=0.0495
for P substrate: aP=0.049; bP=0.052
Mortality coefficients:
for N substrate: g(1)N=0.04; g(2)N=0.03
for P substrate: g(1)P=0.04; g(2)P=0.06

Third phytoplankton group (F3-blue-green
algae)
Maximum growth rate: K=2.5
Preference coefficients for substrate uptake:for N-containing substrate: dNH4=0.2;
dNO2=0.05
dNO3=0.7;
dUR=0.05
for P-containing substrate: dDOP=0.05; dDIP=0.95
Excretion activity:
for N substrate: aN=0.049; bN=0.0495
for P substrate: aP=0.049; bP=0.052
Mortality coefficients:
for N substrate: g(1)N=0.04; g(2)N=0.03
for P substrate: g(1)P=0.04; g(2)P=0.06
First zooplankton group (Z1-herbivorous)
First zooplankton group (Z1-herbivorous)
Maximum growth rate: K=1.5
Maximum growth rate: K=0.75
Preference coefficients for substrate uptake:
Preference coefficients for substrate uptake:for Si-containing substrate: dDOSi=0.15; dDISi=0.02 for Si-containing substrate: dDOSi=0.15; dDISi=0.02
dSiD=0.77; dBSi=0.01
dSiD=0.77; dBSi=0.01
dF1Si=0.05
dF1Si=0.05
for N-containing substrate: dND=0.48; dF1N=0.34
for N-containing substrate: dND=0.5;
dF1N=0.05
dF2N=0.05; dF3N=0.02
dF2N=0.25; dF3N=0.1

Third phytoplankton group (F3-green algae)
Maximum growth rate: K=1.8
Preference coefficients for substrate uptake:for N-containing substrate: dNH4=0.15; dNO2=0.05
dUR=0.6
dNO3=0.2;
for P-containing substrate: dDOP=0.4;
dDIP=0.6
Excretion activity:
for N substrate: aN=0.049; bN=0.0495
for P substrate: aP=0.049; bP=0.0523
Mortality coefficients:
for N substrate: g(1)N=0.0; g(2)N=0.05
for P substrate: g(1)P=0.0; g(2)P=0.1

106

dBN=0.11;
dBN=0.1;
for P-containing substrate: dPD=0.78; dF1P=0.15
for P-containing substrate: dPD=0.73; dF1P=0.1
dF2P=0.025; dF3P=0.025
dF2P=0.025; dF3P=0.025
dBP=0.02;
dBP=0.02; dDOP=0.1
Excretion activity:
Excretion activity:for Si substrate: aSi=0.048; bSi=0.052
for Si substrate: aSi=0.035; bSi=0.052
for N substrate: aN=0.041; bN=0.05
for N substrate: aN=0.041; bN=0.05
for P substrate: aP=0.035; bP=0.05
for P substrate: aP=0.035; bP=0.052
Mortality coefficients:
Mortality coefficients:
for Si substrate: g(1)Si=0.0; g(2)Si=0.2
for Si substrate: g(1)Si=0.05; g(2)Si=0.2
for N substrate: g(1)N=0.0; g(2)N=0.4
for N substrate: g(1)N=0.05; g(2)N=0.4
for P substrate: g(1)P=0.0; g(2)P=0.8
for P substrate: g(1)P=0.035; g(2)P=0.5
Second zooplankton group (Z2-predatory)
Second zooplankton group (Z2-predatory)
Maximum growth rate: K=0.75
Maximum growth rate: K=0.5
Preference coefficients for substrate uptake:
Preference coefficients for substrate uptake:for N-containing substrate: dND=0.55; dF1N=0.31 for N-containing substrate: dND=0.55; dF1N=0.2
dZ1N=0.1;
dF2N=0.02; dF3N=0.02; dZ1N=0.15; dBN=0.06
dBN=0.04
for P-containing substrate: dPD=0.8;
for P-containing substrate: dPD=0.75; dF1P=0.05
dF1P=0.1
dBP=0.05; dZ1P=0.05
dBP=0.05; dZ1P=0.05; dDOP=0.1
Excretion activity:
Excretion activity:
for N substrate: aN=0.0276; bN=0.0287
for N substrate: aN=0.0276; bN=0.03
for P substrate: aP=0.0276; bP=0.0287
for P substrate: aP=0.0276; bP=0.032
Mortality coefficients:
Mortality coefficients:
for N substrate: g(1)N=0.0; g(2)N=0.5
for N substrate: g(1)N=0.05; g(2)N=0.4
for P substrate: g(1)P=0.0; g(2)P=1.0
for P substrate: g(1)P=0.035; g(2)P=0.6-1
Note: the dimension of parameters: K - day , di, ai, bi - (undimension), g(1) - day-1, g(2)i [(mg Element/l)-1 (day-2)].

3.0 Herring group report and model results
Douglas E. Hay1, Robert A. Klumb2, Bernard A. Megrey3, S. Lan Smith4 and Francisco E.
Werner5 (authors listed alphabetically)
1
2
3
4
5

Pacific Biological Station, Fisheries and Oceans Canada, 3190 Hammond Bay Road, Nanaimo, B. C.,
Canada, V9R 5K6. E-mail: hayd@pac.dfo-mpo.gc.ca
Department of Natural Resources, Cornell Biological Field Station, Cornell University, 900
Shackelton Point Road, Bridgeport, NY 13030, U.S.A. E-mail: rak11@cornell.edu
National Marine Fisheries Service, Alaska Fisheries Science Center, 7600 Sand Point Way NE,
Seattle, WA 98115, U.S.A. E-mail: bern.megrey@noaa.gov
Frontier Research System for Global Change, Showa-machi 3173-25, Kanazawa-ku, Yokohama,
Kanagawa, 236-001, Japan. E-mail: lanimal@jamstec.go.jp
Marine Sciences Department, CB# 3300, University of North Carolina, Chapel Hill, NC 27599-3300,
U.S.A. E-mail: cisco@unc.edu

Summary report from the herring group
Specific data for most physiological parameters ofPacific herring are lacking. The first task of“Team Herring” towards linking the LTL

NEMURO model to pelagic fish required
modifications of the existing Atlantic herringbioenergetics model of Rudstam (1988). Threemain areas focused on at the workshop included:1) modifying the temperature dependence function107for consumption and cutoff temperature values,
where swimming speed changes in the Rudstam
model to describe the actual temperatures
inhabited by Pacific herring, 2) accounting for
known differences in larval and juvenile fish
physiology (age-0) compared to adults, and
3) incorporating known seasonal changes inenergy density of adult Pacific herring. Trends insize-at-age were discussed and potential
hypotheses to be tested after completion of themodel were proposed. In the application to Pacificherring our objectives were to model one fish,generate data to compare to observed size-at-age,
follow one cohort through time, and provide ameans to perform regional comparisons.

growth versus abundance and temperature in Haist
and Stocker (1985)).
As a preliminary approximation, Bernard Megrey
normalized the Baltic Sea temperatures to a zeroone scale, based on the maximum and minimum
temperatures observed for a location off the west
coast of Vancouver Island, Amphitrite lighthouse(Fig. 3.1).

Temperature-dependence of consumption andswimming speed
Douglas Hay provided diet data for Pacific herringfrom near Vancouver, British Columbia, fromwhich he and Robert Klumb tried to extract thetemperature-dependence terms for the herringconsumption equation.
The original herringbioenergetics model was formulated for the Baltic
Sea, but the Vancouver site has lower
temperatures and less seasonal variation oftemperature. Because temperature is one of themain process-mediating functions in the
bioenergetics model, we had to modify theparameters for temperature dependence on
consumption function to agree with thetemperature ranges inhabited by Pacific herring off
the coast of Vancouver.
Vadim Navrotsky
suggested that we formulate this temperaturedependence in terms of DT = T - Topt , where Toptis the optimal temperature for consumption(depending upon location). One could also use thetemperature of the waters in which the growth ofherring is maximized as a proxy for thetemperature at which their consumption rate is
maximum (e.g., 12°C, based on the data for peakAge 0Amphitrite
Baltic8.0
110.89715
11.31
1712.55223

Age 1
Amphitrite
8.0
10.897
11.31
12.552

Fig. 3.1 Location of B.C. Lighthouse stations
including Amphitrite Point Lighthouse, the sourceof the temperature data used in the model.
Values for the two temperature series were
rescaled using the formula:
(3.1)

TA =

(TB − TBmin ) ⋅ (TAmax − TAmin )+ TAmin(TBmax − TBmin )where TBmax=30.0, TBmin=1.0, TAmax=14.0,
TAmin=8.0, TA refers to temperatures fromAmphitrite lighthouse, and TB refers totemperatures from the Baltic Sea.
The re-scaled temperatures used for the Thornton
and Lessem (1978) temperature dependence
function for consumption were as follows:
Age > 1
Baltic1
15
17
23

Amphitrite
8.0
10.48310.897
12.553

Baltic
1
13
15
23

108

Temperatures in the respiration model, where
activity changed needed to be re-computed forage-0 and age-1 herring:
Old
15°C
9°C

►

New
10.897°C
9.655°C

Finally the equation describing the annual
temperature signal needed to be re-computed
based on observed mean monthly sea surface
temperature (SST) data from Amphitrite
lighthouse. The following equation
(3.2)

⎛
⎛
⎛ 2 ⋅ π ⋅ ( JDAY − 30 ) ⎞ ⎞ ⎞⎟T = 7.717 + ⎜⎜ 5.6796 ⋅ 0.5 ⋅ ⎜⎜1 − cos⎜⎟ ⎟⎟ ⎟
365
⎝
⎠⎠⎠
⎝
⎝

Fig. 3.3 Results of the Baltic Sea herring modelfrom Rudstam (1988). The solid line represents
model output and the open circles are weight-atage values from field observations.

was fit to the observed data (Fig 3.2) where JDAY
is Julian day and T is water temperature.

Julian Day

Fig. 3.2 Observed and predicted mean SST at
Amphitrite Lighthouse.
Validation to Pacific herring
To validate the bioenergetics model to herring, weused model structure and parameters after
Rudstam (1988) for Baltic Sea herring, but
included no young-of-the-year (YOY) dynamics,
no multispecies functional response, and no
spawning (Rudstam model has spawning).
Results of the model (Fig 3.4) can be compared to
the Rudstam results (Fig. 3.3) and good agreement
in dynamical behavior can be noted.

Fig. 3.4 Simulated growth from the herringbioenergetics model.
Separate age 0 and adult formulations
Describing the growth of a YOY fish involves
more than just rescaling process equations derived
for adult fishes. Often the process rates differ
substantially between different life stages (Postand Lee 1996). Cisco Werner and Rob Klumb
modified the Atlantic herring bioenergetics model
for age-0 herring. Rob’s parameters for respiration
were based on his laboratory measurements fromage-0 alewife, another clupeid, which used routine
metabolism without an activity multiplier.
Literature values for larval and juvenile clupeids

109

were also used that lowered SDA, egestion, andexcretion parameters compared to the adult
Atlantic herring parameters.

if age≥=1 and T ≤9.655 °C thenaR= 0.0033, bR = -0.227, cR= 0.0548, aA= 3.9,bA = 0.13, cA=0.149 and dR=0.03

The YOY formulation for herring respiration
proposed by Arrhenius (1998) along with ourconversion factor from wet weight (g) to energy(J) was:if age≥=1 and T >9.655 °C thenaR= 0.0033, bR = -0.227 cR= 0.0548, aA= 15.0,
bA = 0.13, cA=0.0 and dR=0.03(3.5) R = a R ⋅ W

bRIn the final set of simulations, the Arrhenius(1998) equations 3.5 and 3.6 were modified after
Klumb et al. (in press) to use the parameters.

⋅ f R (T ) ⋅ activity ⋅ 5.258

where the units are the same as in equation 2.1.5
and aR= 0.0033, bR = -0.227.
The temperature
respiration(3.6) f R (T ) = e

dependence

functionfor

if age=0 thenaR= 0.00528, bR = -0.007, cR= 0.0548, aA= 1.0,
bA = 0.0, cA=0.0 and dR=0.0.
In all simulations, equations for age 1 and older
Pacific herring were the same as described in
Arrhenius (1998).

(c R ⋅T )for a age-0 herring was similar to equation 2.1.8.
Activity is a power function of body weight
conditioned on water temperature and is given by( d ⋅U )The coefficients of SDA, egestion, and excretionin equations 2.1.6, 2.1.11, and 2.1.12 were madeage dependent with the parameters given in Table
3.1.

(3.7) activity = e R

Formulation for energy density

where U is swimming speed in cm·s-1 and dR is acoefficient relating swimming speed tometabolism. Swimming speed is calculated as afunction of body weight and temperature using

The energy density of clupeids varies seasonally.
Instead of using constant conversion factors, as inequation 2.1.1, we incorporated a simple energy
cycle based on data in Paul et al. (1998) for age-2
and greater herring. Paul et al. (1998) foundenergy density peaked at 9800 J/g wet wt. (range
9400 - 10200) in fall (October 1), and in spring(March 1) dropped to 5750 J/g wet wt. (range
5200 - 6300). For age-0 and age-1 herring weassumed a constant energy density of 4460 J/g wet
wt. (Foy and Paul 1999). Age-0 herring do exhibita seasonal energy cycle from 5000 J/g wet wt. in
November to 3900 J/g wet wt. in March, and could
be included in future modifications of the model.

b(c ⋅T )

(3.8) U = a A ⋅ W A ⋅ e A

Swimming speeds have been observed to switchfromtemperaturedependence(atlowtemperatures) to temperature independence (athigh temperatures). Formulations by life stage forchanges in swimming speeds versus the adjusted
temperatures from temperature-dependence ofconsumption and swimming speed section (given
earlier) were:if age=0 and T ≤10.897 °C thenaR= 0.0033, bR = -0.227, cR= 0.0548,
aA= 5.76, bA = 0.386, cA=0.238 and dR=0.03
if age=0 and T >10.897 °C thenaR= 0.0033, bR = -0.227, cR= 0.0548, aA= 8.6,
bA = 0.386, cA=0.0 and dR=0.03

The following code was used to implement a
straight-line approximation to a sinusoid that
described seasonal changes in energy density. The
period between March 1 and October 1 consistedof 214 days. The period prior to March 1 and theperiod after October 1, together summed to 151
days.

110

Table 3.1 Summary of final parameter values used in the herring bioenergetics model.
SymbolParameter description

Consumption, CMAX
aC
Intercept for CMAX atbC
coefficient for CMAX versus weight
te1
Temperature for xk1 (in ºC)
te2
Temperature for xk2 (in ºC)
te3
Temperature for xk3 (in ºC)
te4
Temperature for xk4 (in ºC)
xk1
Proportion of CMAX at te1
xk2
Proportion of CMAX at te2
xk3
Proportion of CMAX at te3
xk4
Proportion of CMAX at te4
Metabolism, R
aR
Intercept for RbR
Coefficient for R versus weightcR
Coefficient for R versus temperaturedR
Coefficient for R versus swimming speed
S
Coefficient for Specific Dynamic Action
Swimming Speed, U
aAIntercept U (< 9.655 ºC) (in cm/s)aA
Intercept U (≥9.655 ºC) (in cm/s)bA
Coefficient U versus weightcA
Coefficient U versus temperature (<9.655 ºC)
cA
Coefficient U versus temperature (≥9.655 ºC)

Value
0.642
-0.256

0.10
0.98
0.98
0.01

0.0a
0.125a
3.9bc
15.0bc
0.0bc

Egestion and Excretion, F and E
aF
Proportion of consumed food egested
aE
Proportion of consumed food excreted
Multispecies Functional Response
V11
Vulnerability of prey group 1 to predator 1
V12
Vulnerability of prey group 2 to predator 1
V13
Vulnerability of prey group 3 to predator 1
Half saturation constant for prey group 1 to predator
K11
1 (g wet weight/m3)
K12
Half saturation constant for prey group 2 to predator
1 (g wet weight/m3)
Half saturation constant for prey group 3 to predator
K13
1 (g wet weight/m3)

8.0a
8.0b
8.0c
10.897a 10.897b 10.483c
11.310a 11.310b 10.897c
12.552a 12.966b 12.552c

0.00528a 0.0033bc-0.007a -0.227bc
0.083a 0.0548bc
0.03bc0.175b 0.175c

0.13bc0.149bc

0.125a 0.16bc
0.078a 0.10bc
1.0
1.0
1.0
750.0
75.0
750.0

a - values for age-0 herring, b - values for age-1 herring, c - values for age-2 and older herring

111

-age predicted by the herring bioenergetics modelby adjusting the “p” parameter of equation 2.1.2.

if(iage.ge.2)then
enMar1=5750.
jdMar1=60
enOct1=9800.
jdOct1=274
if(jjday.lt.60)then
delen=(enMar1-enOct1)/151
en=enOct1+(90+jjday)*delenend ifif(jjday.ge.60.and.jjday.lt.274)then
delen=(enOct1-enMar1)/(jdOct1-jdMar1)
en=enMar1+(jjday-jdMar1)*delenend ifif(jjday.ge.274)then
delen=(enMar1-enOct1)/151
en=enOct1+(jjday-jdOct1)*delenend ifelse
en=4460.
end ifp=0.425

Fig. 3.10 Observed size-at-age of the 1973
herring year-class and size-at-age predicted fromthe herring bioenergetics model using p=0.425.

Figure 3.9 shows the straight line approximation toseasonal energy density. Forcing prey fields are
given in Figure 2.1.7.

Fig. 3.11 Observed size-at-age of the 1973
herring year-class and size-at-age predicted fromthe herring bioenergetics model using p=0.40.

Fig. 3.9 Straight line approximationseasonal energy density curve.

toaSimulation results and final modifications to
the herring model
The base model included YOY processes,
included no multi-species functional response andprovided a comparison of observed and predictedsize-at-age, and included no spawning (observed
data were taken after feeding but before
spawning). Figures 3.10 to 3.12 show the fit ofobserved size (weight)-at-age compared to size-at

p = 0 .4 7 5

Fig. 3.12 Observed size-at-age of the 1973
herring year-class and size-at-age predicted fromthe herring bioenergetics model using p=0.475.
112

Observed herring size-at-age data were taken fromthe Straight of Georgia herring data using the 1973
age class, seen as age-1 in 1973 and present in thefishery until age 12 in 1984. As can be seen fromthese figures the model predictions of size-at-age
were extremely sensitive to changing thisparameter, the best fit being when p=0.425. A
long-term simulation with these parameters isshown in Figure 3.13.
The base case was modified to include YOY
improvements, age specific rates, multispecies
functional response, location specific temperaturedescription and change of temperature curveparameters, re-adjustment of p and k’s totemperature change, and seasonal and age
dependent energy density for fish.

Figure 3.14 shows the base case. Figure 3.15
shows the base curve plotted against a run where
the SDA, E and F equations were made age
dependent. Age dependent parameters are given
in Table 3.1. Also plotted in Figure 3.15 aremodel predictions when modifying the respiration
equation to more accurately reflect the metabolic
requirements of an age-0 herring (R) (Klumb et al.
In review). The final curve in Figure 3.15 (All)demonstrates model output when all of thesefeatures were activated.

Fig. 3.15 YOY sensitivity. Observed andpredicted size-at-age due to implementing age
specific formulation for Specific Dynamic Action,
egestion, excretion, respiration one at a time. The“all” line represents the run where all processes
are age dependent and compared to the base run.

Fig. 3.13 Example of a long-term simulation ofherring growth using tuned model parameters.

Fig. 3.14 Comparison of observed and predictedsize-at-age, base case.

Fig. 3.16 Simulation run incorporating new
temperature dependent values (seasonal range 814ºC) and the seasonal herring energy density
algorithm. Comparisons are made of observed
size-at-age between adjusting little “p” in equation
2.1.2 and using the multispecies functional
dependence function (equations 2.1.13 and 2.1.14)
with the k values as described above.

113

Finally, Figure 3.16 demonstrates results of thecustomized herring model for a p value of 0.6375,
and for a run where the multispecies functional
response to three prey types was activated using
the parameters (k1, k2 and k3) shown in Figure
3.16, as well as the seasonal energy density
algoritm. Note that to implement the multispecies
functional
responsefeaturetheline“con=0.75*gcmax” in the FORTRAN code needs
to be commented out.
Trends in size-at-age:
hypothesis testing

some

ideas

for

Douglas Hay has data for size-at-age of Pacificherring over several decades. Over the last 20
years, the mean size-at-age has decreased atseveral locations for fish aged greater than 3 years.
However, the mean size-at-age for ages 1–3 years
did not show a significant decrease which may
result from difficulties in sampling small fish (i.e.
gear selectivity). In agreement to the observed

size-at-age data, measurements from scale annuli
collected over the same period from larger herring
also showed no consistent decrease in growth forfish during the first 3 years of life. This decreasein size-at-age first appears when herring can begin
to eat euphausiids in addition to copepods (age3+). Generally, when euphausiids are abundant,
the predation on herring by other piscivores that
also eat euphausiids is reduced. Given this double
benefit of more available food and less predation,
the growth of herring should be highly sensitive toeuphausiid
production.
The
predatory
zooplankton (ZP) compartment in the NEMURO
model was designed to represent euphausiids.
Thus the coupled NEMURO-herring bioenergetics
model could be used to examine the effects oftemperature and other physical forcings (e.g.,
Pacific Decadal Oscillation) on the production of
euphausiids and thereby on the size-at-age ofherring.

4.0 Saury group report and model results
Shin-ichi Ito1, Michio J. Kishi2, Yasuhiro Yamanaka3 and Masahiko Fujii4
1
2
3
4

Tohoku National Fisheries Research Institute, 3-27-5 Shinhamacho, Shiogama, Miyagi 985-0001,
Japan. E-mail: goito@affrc.go.jp
Hokkaido University, Minato-cho 3-1-1, Hakodate, Hokkaido 041-8611, Japan.
E-mail:kishi@salmon.fish.hokudai.ac.jp
Graduate School of Environmental Earth Science, Hokkaido University, North 10, West 5 Kita-ku,
Sapporo, Hokkaido 060-0810, Japan. E-mail: galapen.ees.hokudai.ac.jp
National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba, 305-8506, Japan. E-mail:fujii.masahiko@nies.go.jp

The members of “team saury” were D. Huang, C.
Hong, Y. I. Zuenko, T. Katukawa, T. Azumaya, S.
Chiba, M. Fujii, M. J. Kishi, K. Tadokoro, M. B.
Kashiwai, Y. Yamanaka, T. Okunishi, A. Tsuda,
D. Mukai, M. Inada, T. Aiki and S. Ito.
According to the life history of Pacific saury, S.
Ito proposed to have saury bioenergetics model
coupled with the ecosystem model composed of athree- ocean-box model which corresponds toKuroshio, Oyashio, and the mixed water region.
But the three-box model is a little complicated tostart with. As a first step we started from acoupled saury bioenergetics-ecosystem model with

one box, and adapted the same type of governing
equations for bioenergetics model as the ones for
Pacific herring.
Model parameters are discussed for applying themodel to Pacific saury. Here we report thediscussion summary and model results.
Life history stages
Pacific saury are spawned in the Kuroshio and themixed water region from autumn to spring. The
larvae are advected to the Kuroshio extension
region and juveniles migrate to the Oyashio region

114

through the mixed water region. After sufficientfeeding in the Oyashio region, they migrate backto the spawning region. The swimming activity,
feeding habitat and metabolism are different
according to the life history stages. Odate (1977)
and Kosaka (2000) divided the Pacific saury life
history stages according to knob length (KL)
(Table 4.1).
Table 4.1 Life stages of Pacific saury after
Odate (1977).

average ration of the Pacific saury are 5.0
gww/day/individual
for
20
cm,
7.2
gww/day/individual for 26 cm, and 10.2
gww/day/individual for 30 cm saury (Kurita andSugisaki; in preparation).
These data were
estimated in the Oyashio region. Comparing thiswith observational data, we adapted 0.6 for ac and-0.256 for bc parameters. Figure 4.1 shows the
CMAX curve and observational value of ration perunit wet weight of the Pacific saury.
2.5

Stage

Knob length

larvaejuvenile
earlier young
later young
small
adult

< 2.5 cm
2.5 - 5.9 cm
6.0 - 9.9 cm
10.0 - 14.9 cm
15.0 - 19.9 cm> 20.0 cm

Cmax

2
1.5

Cmax
feeding/W

1
0.5
0
0

About the earlier stage growth, Watanabe and Kuji
(1991) reared the saury larvae from hatching andthey showed that it takes 60 days to grow to 79
mm KL. Watanabe et al. (1988) analyzed thegrowth rate of Pacific saury and they showed thatit takes about 100 days to grow to 100 mm KL.
According to their result, it takes about 180 days
to become adult saury. Suyama et al. (1996)
showed lower growth rate and it takes about 200
days to become an adult. For simplicity, only
three life stages are assumed in the saurybioenergetics model (Table 4.2).
Table 4.2 Life stages of Pacific saury in thesaruy bioenergetics model.
Stagelarvae and juvenile
young and premature
adult

Age
< 60 days
60-180 days
> 180 days

Maximum consumption rate CMAX
Because adult Pacific saury are too difficult to rearin laboratories, there is no experimental estimationof consumption rate. Field data showed the

50

100

150

wet weightFig. 4.1 CMAX curve and observational value ofration per unit wet weight of the Pacific saury.
Temperature dependency for CMAX
Oozeki (in preparation) analyzed the relationshipbetween saury growth rate and environmental
factors using the same field data reported in
Watanabe et al. (1997). His result showed
positive contributions from surface temperatureand food density to growth rate. The SST range
was between 16-22°C. Oozeki and Watanabe(2000) reared Pacific saury in the laboratory withdifferent water temperatures and found a strong
dependence of growth rate on temperature. The
temperature range was between 12-24°C.
For adult saury we have no measures of growthrate at different temperatures. But the habitat
temperature is between 16 and 20°C. We adapted
the following values for the temperaturedependency parameters for CMAX of Pacific saury(Table 4.3). Figure 4.2 shows the temperaturedependence function for each stage.

115

On the other hand, the wet weight (g)-knob length(cm) relation was proposed by Kosaka (2000) as:

1.2
1stage1stage2

0.6

stage3

0.4
0.2
0
0

10

20

30

Temperature (degC)

Fig. 4.2 Temperature dependence function of
consumption rate of Pacific saury for each life
stage.
Table 4.3 Temperature dependency parameters
for CMAX.

te1
te2
te3
te4
xk1
xk2
xk3
xk4

Stage Stage
1
2
Temperature for 5
5
xk1 (in ºC)
Temperature for 20
16
xk2 (in ºC)
Temperature for 26
20
xk3 (in ºC)
Temperature for 30
30
xk4 (in ºC)
Proportion
of
0.10
CMAX at te1
Proportion
of
0.98
CMAX at te2
Proportion
of
0.98
CMAX at te3
Proportionof
0.5
CMAX at te4

Stage
3
5

larvae and juvenile
logW = -2.069+2.42439logL
earlier young
logW = -2.483+3.06174logL
later young
logW = -2.335+2.93760logL
small
logW = -2.688+3.22526logL
adult
logW = -2.685+3.21229logL
Figure 4.4 shows the wet weight-knob length
relation curves of Kosaka (2000). If we adapt thesimple one curve for all stages, it becomes
W = ( KL / 6.13 ) 3
and the curve will look like Figure 4.4. The
broken line in Figure 4.3 shows the same curve.

16
20
30

weight (g)

fc(T)

0.8

160

70

140

60

120

50

100

40

80

weightswimming speed

30

60
40

20

20

10

0

0
0

10

20

30

knob length (cm)Fig. 4.3 Swimming speed (cm/s) and wet weight(g) as a function of body length (cm).
180
160

Swimming speed
Although we do not have actual data on swimming
speed of Pacific saury, other small pelagic fish
swim at speeds of several times their body length
per second. We assumed the normal swimming
speed is two times of the knob length (nearly sameas body length) per weight (Fig. 4.3).
U = 2.0 KL

weight (g)

140
120
100

weight
obs

80
60
40
20
0
0

10

20

30

body length (cm)Fig. 4.4 Wet weight (g)-knob length (cm)relation curves of Kosaka (2000) (red) and fitting
curve (black).

116

The last equation could be rewritten as

70

0.33

60

KL = 6.13 W

50

U = 12.3 W0.33
and we adapted 12.3 as aA parameter when thetemperature is higher than 12°C and 0.33 for bA
value. For temperatures less than 12°C we
adapted 2.0 as aA. The weight - swimming speed
relation looks like Figure 4.5 when thetemperature is higher than 12°C.
The parameters we adapted for Pacific saurybioenergetics model are summarized on Table 4.4.

U (cm/s)and the swimming speed becomes

40
U

30
20
10
0
0

50

100

150

wet weight (g)Fig. 4.5 Wet weight (g)-swimming speed(U-cm/s) relation curve for higher temperature.

Table 4.4 Summary of parameter values used in the saury bioenergetics model.
SymbolParameter description

Value

Consumption, CMAX
aC
Intercept for CMAX at (te1+te3)/2bC
coefficient for CMAX versus weight
te1
Temperature for xk1 (in ºC)
te2
Temperature for xk2 (in ºC)
te3
Temperature for xk3 (in ºC)
te4
Temperature for xk4 (in ºC)
xk1
Proportion of CMAX at te1
xk2
Proportion of CMAX at te2
xk3
Proportion of CMAX at te3
xk4
Proportion of CMAX at te4

0.6-0.2565, 5, 520, 16, 16
26, 20, 20
30, 30, 30
0.10
0.98
0.98
0.5

Metabolism, R
aR
Intercept for RbR
Coefficient for R versus weightcR
Coefficient for R versus temperaturedR
Coefficient for R versus swimming speed
S
Coefficient for Specific Dynamic Action

0.0033-0.2270.05480.03
0.175

Swimming speed, U
aAIntercept U (< 12 ºC) (in cm/s)aA
Intercept U (≥ 12 ºC) (in cm/s)bA
Coefficient U versus weight
Coefficient U versus temperature (< 12 ºC)
cA
cA
Coefficient U versus temperature (≥ 12 ºC)

2.0
12.3
0.33
0.149
0.0

Egestion and excretion, F and E
Proportion of consumed food egested
aF
aE
Proportion of consumed food excreted

0.16
0.10

Multispecies functional response (by saury size groups)
117

V11
V12
V13
K11
K12
K13
V21
V22
V23
K21
K22
K23
V31
V32
V33
K31
K32
K33

Vulnerability of prey group 1 to predator 1
Vulnerability of prey group 2 to predator 1
Vulnerability of prey group 3 to predator 1
Half saturation constant for prey group 1 topredator 1 (g wet weight/m3)
Half saturation constant for prey group 2 to
predator 1 (g wet weight/m3)
Half saturation constant for prey group 3 topredator 1 (g wet weight/m3)
Vulnerability of prey group 1 to predator 2
Vulnerability of prey group 2 to predator 2
Vulnerability of prey group 3 to predator 2
Half saturation constant for prey group 1 topredator 2 (g wet weight/m3)
Half saturation constant for prey group 2 to
predator 2 (g wet weight/m3)
Half saturation constant for prey group 3 topredator 2 (g wet weight/m3)
Vulnerability of prey group 1 to predator 3
Vulnerability of prey group 2 to predator 3
Vulnerability of prey group 3 to predator 3
Half saturation constant for prey group 1 topredator 3 (g wet weight/m3)
Half saturation constant for prey group 2 to
predator 3 (g wet weight/m3)
Half saturation constant for prey group 3 topredator 3 (g wet weight/m3)

start day is February 1ststage 1 0-50mm
stage 2 50-200mm
stage 3 >200mm

1.0
0.0
0.0
100.0
100.0
100.0
1.0
1.0
0.0
100.0
100.0
100.0
0.0
1.0
1.0
100.0
100.0
100.0

0-30days
30-150days
150day-720days

Pacific saury

Model result

Figure 4.7 shows the interannual experiment of
ecosystem-saury coupled model with realistic

200
weight (g)The parameters which are revised for the Pacificsaury were used to integrate the bioenergetics
model coupled with the ecosystem model. Figure
4.6 shows the result of the integration, and showsthat the weight of saury reached 120 g after one
year. This seems reasonable for Pacific saury.
The model shows a high growth rate around 13ºC
water temperature. This corresponds to the habitat
temperature in the Oyashio region during thefeeding season.

25
20

150

15

100

10

50

weight
temp
light

5

0
0

200

400

600

0
800

age (day)Fig. 4.6 Result of Pacific saury bioenergetics
model. Light is in units of ly/min.

118

forcing of A7 (Akkeshi line St 7 off Hokkaido,
Japan). The model results show low growth ratein the third and fourth year cohort. The result
strongly depends on water temperature.

This model is not perfect and needs improvements
in several respects.

200
1
150
WeightFuture work

Pacific saury

2

100

3

50

5

4

0
DM J S DM J S DM J S DM J S DM J S DM J S
Month

•
•

•

The weight of the earliest stage is not
reproduced well. We should re-parameterizevalues for this stage.
More than half of the Pacific saury spawn in
the first year and all of them spawn in thesecond year (Kurita and Sugisaki; in
preparation). We should include the effect ofspawning in this model.
In this model only one ocean region is
included. But the saury migrate from thesubtropical to the subarctic region through themixed water region, each with its own
seasonal cycle of temperature and prey. We
should include at least three ocean regions in
the ecosystem-saury coupled model. We
suggest Figure 4.8 as a prototype three ocean
region model. This kind of model is very
useful for the analysis of interannual
variability of saury growth.

Fig. 4.7 Result of Pacific saury bioenergetics
model with realistic forcing.

Kuroshio

Mixed water
regionOyashioFig. 4.8 Schematic picture of a three-ocean-box
model. This model includes three ocean regions
but only one saury bioenergetic model.

5.0 Model experiments and hypotheses
Several model experiments were discussed to test
hypotheses regarding the effects of climatechange. The details of the experiments andhypotheses are described below.

Understanding the extent to which environmental
conditions account for these temperature
differences in herring growth is important forpredicting climate change effects and for effective
management of these fisheries in the future.

Space hypothesis
Geographic variation in fish growth: Differencesin environmental conditions, and resulting
differences in lower trophic conditions, canaccount for the differences in herring growth rates
among selected sites in the North Pacificecosystem. There exist long-term data sets on
size-at-age of herring from many locations in the
North Pacific. These data sets show that herring
growth rates over the past decades have varied
consistently among the different locations.

Key regions contributing to fish growth andbiomass variations: Pacific saury are spawned in
the subtropical and transition zone from autumn tospring, and migrate from the subtropical to thesubarctic ocean through a transition zone. The
environments of these three regions show different
interannual variability, and it is very difficult todistinguish which location (or season) is mostimportant to the interannual variability of fish
growth and biomass. We will tune up the NPZ
model coupled to the fish growth model with long-119

term climate and weather records by comparing
the model results of fish growth with interannual
variability in observed size composition. To
understand which region contributes most to theinterannual variability in fish growth and biomass,
sensitivity tests of this model will be very useful.
Understanding key regions in fish growth andbiomass variation is also important for predicting
climate change effects and for effective
management of these fisheries in the future.

environmental conditions.
For example, thedominant species changed from sardine to saury
across the regime shift in 1987.

Time hypothesis

Global climate change effects on energy
pathways and fish production: Climate change
may result in energy being diverted from thepelagic pathway and shunted through themicrobial pathway, resulting in less food for
pelagic fish and consequently slower fish growth
rates. We will use the coupled NPZ and fishmodels, the long-term datasets, and defined
climate change scenarios to predict how climatechange might affect energy cycling, shift thedominance among different phytoplankton andzooplankton groups, and affect fish growth andproduction in the North Pacific ecosystem. Model
simulations will be performed under present-day(baseline) environmental conditions, and for asuite of realistic climate change scenarios.
Comparing these linkages and pathways underbaseline and climate change scenarios for a varietyof locations that have different environmental
conditions (e.g., shallow coastal versus deep blue
water) will aid in the interpretation andgeneralization of our results.

Understanding regime shifts:
Synchronouschanges in herring growth rates across locations
may be accounted for by basin-wide decadal-scalechanges in environmental conditions. Preliminary
examination of herring growth rates at several
locations showed sudden shifts in growth rates
occurring in the same years across all locations.
We will combine the long-term datasets on herring
growth, where possible, with regional and local
long-term climate and weather records, and use the
NPZ model coupled to the fish growth model toexamine possible environmental regime changes.
Understanding how regime shifts cascade up thefood web may be our best chance for using pastconditions to infer future effects of climatechange.
Change of dominant species: Changes in thedominant small pelagic fish species seems tocoincide with basin-wide decadal-scale changes in

Comparing different fish bioenergetics models
with the same NPZ model is very useful tounderstand the climate change effects on ocean
ecosystems through bottom-up processes.
Climate change hypothesis

6.0 Recommendations
Results of the model work accomplished at theworkshop resulted in several recommendations.
They are listed here in priority. A description of
workplan scheduling is indicated after each item inparentheses and in italics:
1. Develop site-specific applications (to bescheduled);
2. Perform herring comparison between the Seaof Okhotsk and Vancouver Island (to bescheduled);

3. Incorporate data observations into NEMURO(to be programmed);a. Obtain physical parameters (radiation,
cloud cover, wind stress);b. Obtain realistic time series of SST andphotosynthetically active light;
c. Obtain physical observed time series;
d. Obtain observed zooplankton time series;
4. Executeadynamics
linkage
NEMURO.FISH (to be scheduled);

in

120

5. Revise physiological parameters (fish andLTL) (to be scheduled);
6. Public web distribution (PICES to support) (tobe scheduled);
7. Meet in Qingdao with (1) and (2) finished (tobe scheduled);
8. Consider explicit spatial (x, y, z) and temporal
structure (to be programmed);
9. Dissemination of NEMURO.FISH results in
GLOBEC
newsletterandscientificpublications (to be scheduled);
10. Develop ascheduled);

project

home

page

(to

be

11. Incorporate age structure, reproduction andearly life history into NEMURO.FISH (to be
programmed).
The last recommendation was discussed at length
because of a perceived need to provide a tool forthe management of fisheries.
It was noted that fish biomass at any time can be
represented as the product of fish weight (W) andfish numbers (N).
(6.1) Bt = N t ⋅ Wt
The rate of change of fish biomass can be written
as
(6.2)

dB
dW
dN
=N
+W
dt
dt
dt

where we know from the fish bioenergetics model
that
(6.3)

dW
= C − ( R + SDA + E + F ) ⋅ W
dt

and from fish population dynamics we know that
(6.4)

dN
= e −( F + M ) ⋅ N
dt

where F is the instantaneous fishing mortality rateand M is the instantaneous natural mortality rate.

Thus all of the components of equation 6.2 are
known or can be estimated.
Also note that with this approach we can compare
observed growth to von Bertalanffy growth, agrowth model very commonly used in fisheriesscience.
The von Bertalanffy empirical growth model (vonBertalanffy 1938) is written(6.5) W (t ) = W∞ ⋅ (1 − e− k ⋅( t −t 0 ) 3)

where W∞ is the asymptotic weight, W is theweight at time t, k is the growth parameter withunits t-1, and t0 is the theoretical age the fish would
have zero weight had they always grown as
described by (6.5).
The differential form of (6.5) is
(6.6)

2
1
dW
= 3k ⋅ (W 3 ⋅ W∞ 3 − W )dt

We know that the change in weight formulationfrom the fish bioenergetics model (6.3) can becollapsed into a simpler form because most terms
are dependent on consumption thus they areproportional to C and the C and R terms are weight
dependent. So (6.3) can be simplified to
(6.7)

dW
= p3 ⋅ C − p1 ⋅ W p 2
dt

Equating (6.7) and (6.6) we can calculate the rateof consumption required for von Bertalanffy
growth, C*, as
(6.8)

C* =

2
1
1 ⎡
⋅ 3k ⋅ W 3 ⋅ W∞ 3 − 3k ⋅ W + p1 ⋅ W p 2 ⎤
⎥⎦
p3 ⎢⎣

So it seems the theoretical foundation of extending
NEMURO.FISH to a population level model
useful for fisheries management is possible with afew minor modifications.

121

7.0 Achievements and future steps
The achievements of the Workshop can be listed
as follows:
1. Developed
theprototypemodel,NEMURO.FISH.
This was an extremely
important step because it translates our science
into something tangible and economically
relevant to human populations that rely on
fishes for food - obtaining food from the seas
on a sustainable basis;
2. Assembled an international team of marine
biologists, fisheries biologists, and physicaloceanographers who collectively achieved aconsensus on the structure and function of a
PICES Climate Change and Carrying Capacity(CCCC) prototype lower trophic level (LTL)
ecosystem model for the North Pacific Ocean
that included pelagic fishes, and named it“NEMURO.FISH”;

3. Developed a computer simulation model offish bioenergetics and growth;
4. Coupled the fish model to the NEMURO
lower trophic level model;
5. Adapted the fish bioenergetics model toPacific herring (Clupea harengus pallasi) in
the eastern North Pacific and Pacific saury(Cololabis saira) in the western North Pacific;
6. Made recommendations for future modelingactivities.
The significance of these achievements will
ultimately be evaluated by how well the CCCC
Program effectively utilizes and embraces thesemodels as a basis of future modeling activity.

8.0 Acknowledgements
This workshop was proposed and convened by
PICES, more precisely the PICES/CCCCIP/MODEL and REX Task Teams. On behalf of
the workshop participants, the co-conveners would
like to express sincere thanks for giving us atimely and valuable opportunity to participate in
the development of a linked higher trophic level lower trophic level marine ecosystem model
common among component programs of the
PICES GLOBEC Program. The Heiwa-Nakajima
Foundation of Japan contributed the major part offinancial support for this workshop. We note thisvaluable support, for without it, the workshop
would not have been possible. Nemuro wasselected as the venue due in large part to theinvitation from the Nemuro Supporting
Committee. The conveners express their very
deep appreciation for the warm welcome andperfect support arranged and given by the Nemuro
Supporting Committee, their staff, and the peopleof Nemuro.

9.0 References
Arrhenius, F. 1998a. Variable length of daily
feeding period in bioenergetics modelling: atest with 0-group Baltic herring. Journal of
Fish Biology 52:855-860.
Arrhenius, F. 1998b. Food intake and seasonalchanges in energy content of young Baltic Sea
sprat (Sprattus sprattus L.). ICES Journal of
Marine Science 55:319-324.

Arrhenius, F. and Hansson, S. 1996. Growth andseasonal changes in energy content of young
Baltic Sea herring (Clupea harengus L.).
ICES Journal of Marine Science 53:792-801.
Batty, R. S. 1987. Effect of light intensity on
activity and food-searching of larval herring,
Clupea harengus: a laboratory study. Marine
Biology 94:323-327.

122

Beamish, R.J., Noakes, D.J., Mcfarlane, G.A.,
Klyashtorin, L., Ivanov, V.V., and Kurashov,R. 1999. The regime concept and natural
trends in the production of Pacific salmon.
Can. J. Fish. Aquat. Sci. 56: 516-526.
Beamish, F. W. H. 1974. Apparent specificdynamic action of largemouth bass,
Micropterus salmoides. Journal of theFisheries Research Board of Canada 31:17631769.
Beamish, F. W. H. and Trippel, E. A. 1990. Heat
increment: A static or dynamic dimension inbioenergetic models? Transactions of theAmerican Fisheries Society 119:649-661.
Beitinger, T. L. and Magnuson, J. J. 1979.
Growth rates and temperature selection of
bluegill, Lepomis macrochirus. Transactionsof the American Fisheries Society 108:378382.
Bertalanffy, L. von. 1938. A quantitative theoryof organic growth (Inquiries on growth laws
II). Human Biology, 10(2): 181-213.
Blaxter, J. H. S. and Hunter, J. R. 1982. The
biology of clupeoid fishes. In J. H. S. Blaxter,
F. S. Russell, and M. Yonge (eds.) Advancesin Marine Biology Volume 20. Academic
Press, London, UK.
Brett, J. R.
1985.
Correction in use ofoxycalorific equivalent. Canadian Journal of
Fisheries and Aquatic Sciences 42:1326-1327.
Brett, J. R. and Groves, T. D. D.
1979.
Physiological energetics. In W. S. Hoar, D. J.
Randall, and J.R. Brett (eds.) Fish physiology
Vol. VIII bioenergetics and growth.
Academic Press, Inc., New York, pp. 279 –352.
Cummins, K. W. and Wuycheck, J. C. 1971.
Caloric equivalents for investigations inecologicalenergetics.
Mitteilungen
Internationale Vereinigung fur Theoretische
und Angewandte Limnologie 18.
Dabrowski, K. 1985. Energy budget of coregonid(Coregonus spp.) fish growth, metabolism andreproduction. Oikos 45:358-364.
Dabrowski, K. 1986. Energy utilization duringswimming and cost of locomotion in larval
and juvenile fish. Sonderdruck aus Journal ofApplied Ichthyology 2:110-117.
Dabrowski, K., Takashima, F. and Law, Y. K.
1988. Bioenergetic model of planktivorous
fish feeding, growth and metabolism:

theoretical optimum swimming speed of fishlarvae. Journal of Fish Biology 32:443-458.
De Silva, S. S. 1973. Food and feeding habits of
the herring Clupea harengus and the sprat C.
sprattus in inshore waters of the west coast of
Scotland. Marine Biology 20:282-290.
De Silva, C. D. and Tytler, P. 1973. The
influence of reduced environmental oxygen on
the metabolism and survival of herring andplaice larvae. Netherlands Journal of Sea
Research 7:345-362.
De Silva, S. S. and Balbontin, F.
1974.
Laboratory studies on food intake, growth andfood conversion of young herring, Clupea
harengus (L.). Journal of Fish Biology 6:645658.
Durbin, E. G. and Durbin, A. G.
1981.
Assimilation efficiency and nitrogen excretionof a filter-feeding planktivore, the Atlantic
menhaden, Brevoortia tyrannus (Pisces:
Clupeidae). Fishery Bulletin U.S. 79:601-616.
Durbin, E. G. and Durbin, A. G. 1983. Energy
and nitrogen budgets for the Atlantic
menhaden, Brevoortia tyrannus (Pisces:
Clupeidae), a filter-feeding planktivore.
Fishery Bulletin U.S. 81:177-199.
Durbin, A. G., Durbin, E. G., Verity, P. G., andSmayda, T. J. 1981. Voluntary swimming
speeds and respiration rates of a filter-feedingplanktivore,theAtlanticmenhaden,Brevoortia tyrannus (Pisces: Clupeidae).
Fishery Bulletin U.S. 78:877-886.
Elliott, J. M. 1976a. Energy losses in the waste
products of brown trout (Salmo trutta L.).
Journal of Animal Ecology 45:561-580.
Elliott, J. M. 1976b. The energetics of feeding,metabolism, and growth of brown trout (Salmo
trutta L.) in relation to body weight, watertemperature, and ration size. Journal ofAnimal Ecology 45:923-948.
Elliott, J. M. 1979. Energetics of freshwaterteleosts. Symposium Zoological Society of
London 44:29-61.
Elliott, J. M. and Davison, W. 1975. Energy
equivalents of oxygen consumption in animalenergetics. Oecologia 19:195-201.
Flath, L. E. and Diana, J. S. 1985. Seasonalenergy dynamics of the alewife insoutheastern Lake Michigan. Transactions of
the American Fisheries Society 114:328-337.

123

Fiksen, O. and Folkvord, A. 1999. Modelling
growth and ingestion processes in herring
Clupea harengus larvae. Marine Ecology
Progress Series 184:273-289.
Foy, R. J. and Norcross, B. L. 1999. Spatial andtemporal variability in the diet of juvenile
Pacific herring (Clupea pallasi) in Prince
William Sound, Alaska. Canadian Journal of
Zoology 77:697-706.
Foy, R. J. and Paul, A. J. 1999. Winter feeding
and changes in somatic energy content of age0 Pacific herring in Prince William Sound,
Alaska.
Transactions of the American
Fisheries Society 128:1193-1200.
Grabe, S. A. 1996. Feeding chronology andhabits of Alosa spp. (Clupeidae) juveniles
from the lower Hudson River estuary, NewYork.
Environmental Biology of Fishes
47:321-326.
Giguère, L. A., Côté, B., and St-Pierre, J.-F. 1988.
Metabolic rates scale isometrically in larval
fishes. Marine Ecology Progress Series 50:1319.
Haegele, C.W. 1997. The occurrence, abundance
and food of juvenile herring and salmon in theStrait of Georgia, British Columbia in 1990 to
1994. Can. Man. Rep. Fish. Aquat. Sci. 2390:
124 pp.
Haist, V. and Stocker, M. 1985. Growth andmaturation of Pacific herring (Clupea
harengus pallasi) in the Strait of Georgia.
Proceedings of The Symposium on the
Biological Characteristics of Herring and
Their Implication for Management. Can. J.
Aquat. Sci., vol. 42(l), pp. 138-146.
Hanson, P. C., Johnson, T. B., Schindler, D. E.,
Kitchell, J. F. 1997. Fish bioenergetics 3.0
for Windows. University of Wisconsin Sea
Grant Institute. Technical Report WISCU-T97-001, Madison, Wisconsin, USA.
Hartman, K. J. and Brandt, S. B.
1995.
Estimatingenergydensityoffish.
Transactions of the American Fisheries
Society 124:347-355.
Hay, D.E. and McCarter, P.B. 2001. Spatial,
temporal and life-stage variation in herring
diets in British Columbia. PICES Workshop
Reports. PICES Sci. Report 15: 95-98.
Hay, D.E., Thompson, M.J., and McCarter, P.B.
2001. Anatomy of a strong year class:Analysis of the 1977 year class of Pacificherring in British Columbia and Alaska.
Herring: Expectations for a new millennium.
University of Alaska Sea Grant, AK-SG-0104, Fairbanks, pp. 171-198.
Hettler, W. F. 1976. Influence of temperature andsalinity on routine metabolic rate and growthof young Atlantic menhaden. Journal of Fish
Biology 8:55-65.
Hollowed, A.B. and Wooster, W.S.
1992.
Variability of winter ocean conditions andstrong year-classes of Northeast Pacificgroundfish. ICES Mar. Sci. Symp. 195: 433444.
Hollowed, A.B. and Wooster, W.S.
1995.
Decadal scale variations in the eastern
subarctic Pacific: II. Response of Northeast
Pacific fish stocks. In R.J. Beamish. (ed.)
Climate change and northern fish populations.
Canadian Special Publication of Fisheries andAquatic Sciences 121: 375-386.
Ikeda, T. 1996. Metabolism, body composition,and energy budget of the mesopelagic fish
Maurolicus muelleri in the Sea of Japan.
Fishery Bulletin U.S. 94:49-58.
James, A. G. and Probyn, T. 1989. The
relationshipbetweenrespiration
rate,
swimming speed and feeding behaviour in theCape anchovy Engraulis capensis Gilchrist.
Journal of Experimental Marine Biology andEcology 131:81-100.
Janssen, J. 1976. Feeding modes and prey size
selectioninthealewife(Alosapseudoharengus). Journal of the Fisheries
Research Board of Canada 33:1972-1975.
Janssen, J. and Brandt, S. B. 1980. Feeding
ecology and vertical migration of adultalewives (Alosa pseudoharengus) in Lake
Michigan. Canadian Journal of Fisheries andAquatic Sciences 37:177-184.
Kamler, E. 1972. Respiration of carp in relation
to body size and temperature.
Polskie
Archiwum Hydrobiologii 19: 325-331.
Katz, H. M. 1978. Circadian rhythms in juvenile
American shad, Alosa sapidissima. Journal of
Fish Biology 12:609-614.
Kaufmann, R.
1990.
Respiratory cost ofswimming in larval and juvenile cyprinids.
Journal of Experimental Biology 150:343-366.
Kerr, S. R. and Dickie, L. M. 1985. Bioenergetics
of 0+ Atlantic herring (Clupea harengus

124harengus). Canadian Journal of Fisheries andAquatic Sciences 42 (1):105-110.
Kiørboe, T., Munk, P. and Støttrup, J. G. 1985.
First feeding by larval herring Clupea
harengus L. Dana 5:95-107.
Kiørboe, T., Munk, P. and Richardson, K. 1987.
Respiration and growth of larval herring
Clupea harengus: relation between specificdynamic action and growth efficiency. Marine
Ecological Progress Series 40:1-10.
Kitchell, J. F., Koonce, J. F., O'Neill, R. V.,
Shugart, H. H., Jr., Magnuson, J. J. and Booth,
R. S. 1974. Model of fish biomass dynamics.
Transactions of the American Fisheries
Society 103:786-798.
Klumb, R. A., Rudstam, L. G., and Mills, E. L.
Respiration rates and swimming speeds oflarvalandjuvenile
alewives
Alosapseudoharengus:
implicationsforbioenergetics models. Transactions of theAmerican Fisheries Society (in review).
Klumpp, D. W. and Westernhagen, H. von. 1986.
Nitrogen balance in marine fish larvae:
influence of developmental stage and preydensity. Marine Biology 93:189-199.
Kosaka S. 2000. Life history of the Pacific saury
Cololabis saira in the northwest Pacific andconsiderations on resources fluctuations basedon it. Bull. Tohoku Natl. Fish. Res. Inst. 63:
1-96.
Kurita Y. and Sugisaki, H. 2002. Seasonalchanges in daily ration of Pacific saury,
Cololabis saria. to be submitted to Fish.
Oceanogr.
Laurence, G. C. 1976. Caloric content of somenorth Atlantic calanoid copepods. Fishery
Bulletin U.S. 78: 218-220.
Leonov A.V and Sapozhnikov V.V.
1997.
Biohydrochemical model of organic substancetransformations and its applying for the
calculation of primary production in the
Okhots Sea ecosystem. In Complex studiesof the Okhotsk Sea ecosystem. Publ. House
VNIRO, Moscow, pp. 143-166. (In Russian).
Leonov A.V. and Stygar O.V. 1999. Seasonalvariationsofbiogenic
substance
concentrations and bioproductivity of waters
in northern part of the Caspian Sea. Vodnye
Resoursy V. 26 (6), pp. 743-756. (In Russian).
Limburg, K. E. 1994. Ecological constraints on
growth and migration of juvenile American

shad (Alosa sapidissima Wilson) in the
Hudson River Estuary, New York. Doctoral
dissertation. Cornell University, Ithaca, NewYork.
Muir, B. S. and Niimi, A. J. 1972. Oxygen
consumption of the euryhaline fish aholehole(Kuhlia sandvicensis) with reference tosalinity, swimming, and food consumption.
Journal of the Fisheries Research Board of
Canada 29:67-77.
Munk, P. 1992. Foraging behaviour and prey size
spectra of larval herring Clupea harengus.
Marine Ecology Progress Series 80:149-158.
Munk, P. and Kiorboe, T. 1985. Feeding
behaviour and swimming activity of larval
herring (Clupea harengus) in relation todensity of copepod nauplii. Marine Ecology Progress Series 24:15-21.
Nakata, K., Ito, H., Ichikawa, T., and Sasaki, K.
2002. Reproduction of three Oncaea speices
in the Kuroshio Extention in spring. Fish.
Oceanogr. (submitted ).
Odate, S. 1977. On distribution of Pacific sauryin the North Pacific Ocean. Res. Inst. North
Pac. Fish. Sp. Vol. 10 Faculty of Fisheries,
Hokkaido University, Hakodate, Japan, pp.
353-382.
Ohtani, K. and Azumaya, T. 1995. Influence ofinterannual changes in ocean conditions on theabundance of walleye pollock (Theragra
chalcogramma) in the eastern Bering Sea. In
R.J. Beamish. (ed.) Climate change andnorthern fish populations. Canadian Special
Publication of Fisheries and Aquat. Sci. 121:
87-95.
Oozeki, Y. and Watanabe, Y. 2000. Comparisonof somatic growth and otolith increment
growth in laboratory-reared larvae of Pacificsaury, Cololabis saira, under different
temperature conditions. Marine Biology 136:
349-359.
Oozeki, Y. and Watanabe, Y. 2002. Environment
factor affecting the larval growth of Pacificsaury, Cololabis saira, in the northwestern
Pacific Ocean. Fish. Oceanogr. (submitted ).
Paul, A. J. and Paul, J. M. 1998a. Spring andsummer whole-body energy content of
Alaskan juvenile Pacific herring. Alaska
Fishery Research Bulletin 5:131-136.
Paul, A. J. and Paul, J. M. 1998b. Comparisonsof whole body energy content of captive

125

fasting age zero Alaskan Pacific herring(Clupea pallasi Valenciennes) and cohorts
over-wintering in nature.
Journal of
Experimental Marine Biology and Ecology
226:75-86.
Paul, A. J., Paul, J. M., and Brown, E. D. 1998.
Fall and spring somatic energy content for
Alaskan Pacific herring (Clupea pallasi
Valenciennes 1847) relative to age, size andsex. Journal of Experimental Marine Biology
and Ecology 223:133-142.
Pischalnik V.M. and Leonov A.V.
2001.
Experience of combined application ofelectronic atlas of oceanographical data forshelf zone of Sakhalin Island and simulation
mathematical model for the study ofbiotransformation processes in marine
environment (examples for La Perouse Strait
and Aniva Bay). (in press, in Russian).
Polovina, J.J., Mitchum, G.T. and Evans, G.T.
1995. Decadal and basin-scale variation in
mixed layer depth and the impact on
biological production in the Central and North
Pacific 1960-1988.
Deep Sea Research
42:1701-1716.
Post, J. R. and Lee, J. A. 1996. Metabolic
ontogeny of teleost fishes. Canadian Journalof Fisheries and Aquatic Sciences 53:910-923.
Rose, K.A., Rutherford, E.S., McDermot, D.S.,
Forney, J.L., and Miles, E.L.
1999.
Individual-based model of yellow perch andwalleye populations in Oneida Lake.
Ecological Monographs 69(2): 127-154.
Rudstam, L. G. 1988. Exploring the dynamics ofherring consumption in the Baltic: applicationsof an energetic model of fish growth. Kieler
Meeresforschung Sonderheft 6:312-322.
Sherman, K. and Honey, K. A. 1971. Seasonalvariations in the food of larval herring incoastal waters of central Maine. Rapports et
Proces-Verbaux des Reunions.
ConseilInternational pour l’Exploration de la Mer
160:121-124.
Stewart, D. J. and Binkowski, F. P. 1986.
Dynamics of consumption and foodconversion by Lake Michigan alewives: anenergetics-modeling synthesis. Transactionsof the American Fisheries Society 115:643661.
Stewart, D. J., Weininger, D., Rottiers, D. V., and
Edsall, T. A. 1983. An energetics model forlake trout, Salvelinus namaycush: application
to the Lake Michigan population. Canadian
Journal of Fisheries and Aquatic Sciences 40:
681-698.
Sugisaki H. and Kurita, Y. 2002. Daily rhythm
and seasonal variation of feeding habit ofsaury. Fish. Oceanogr. (submitted ).
Suyama S., Sakurai, Y., and Shimazaki, K. 1996
Age and growth of pacific saury Cololabis
saira (Brevoort) in the Western North Pacific
Ocean estimated from daily otolith growthincrements. Fish. Science 62: 1-7.
Suyama S., Kurita, Y., Kamei, Y., Kajiwara, Y.,
and Ueno, Y. 2002. Annual change of thesize in each year class of Pacific saury(Cololabis saira) estimated based on thehyaline zone in the otolith. Fish. Oceanogr.
(submitted ).
Theilacker, G. H. 1987. Feeding ecology andgrowth energetics of larval northern anchovy,
Engraulis mordax. Fishery Bulletin U.S.
85:213-228.
Thornton, K. W. and Lessem, A. S. 1978. A
temperaturealgorithm
formodifyingbiological rates. Transactions of the American
Fisheries Society 107:284-287.
Tian Y., Akamine, T., and Suda, M. 2002a.
Long-term variability in the abundance ofPacific saury in the northwestern Pacific oceanand climate changes during the last century.
Bull. Jpn. Soc. Fish. Oceanogr. 66; 16-25.
Tian Y., Ueno, Y., Akamine, T., and Suda, M.
2002b. Climate-ocean variability and the
response of Pacific saury (Cololabis saira) in
the northwestern Pacific during the lastcentury. Fish. Sci. (submitted).
Wailes, G. H. 1936. Food of Clupea pallasii insouthern British Columbia waters. Journal of
the Biology Board of Canada 1:477-?
Watanabe. Y., and Butler, J. L., and Mori, T.
1988. Growth of Pacific saury, Cololabis
saira, in the northeastern and northwestern
Pacific Ocean. Fish. Bull. U.S. 86: 489-498.
Watanabe. Y., and Kuji, Y. 1991. Verification ofdaily growth increment formation in saury
otolith by rearing larvae from hatching. J.
Ichthyol. 38: 11-15.
Watanabe, Y., Oozeki, Y., and Kitagawa, D.
1997.
Larval parameters determiningpreschooling juvenile production of Pacificsaury (Cololabis Saira) in the northwestern

126

Pacific. Can. J. Fish. Aquat. Sci., 54: 10671076.
Watanabe. Y., and Lo, N. C. H. 1989. Larval
production and mortality of Pacific saury,
Cololabis saira, in the northwestern Pacific
Ocean. Fish. Bull. U.S. 86: 601-613.
Ware, D.M.and McFarlane, G.A. 1995. Climateinduced changes in Pacific hake (Merluccius
productus) abundance in the Vancouver Island
upwelling system. In R.J. Beamish. (ed.)
Climate change and northern fish populations.
Canadian Special Publication of Fisheries andAquat. Sci. 121: 509-521.

Ware, D. M. 1975. Growth, metabolism, andoptimal swimming speed of a pelagic fish.
Journal of the Fisheries Research Board of
Canada 32:33-41.
Wieser, W. and Forstner, H. 1986. Effects oftemperature and size on the routine rate ofoxygen consumption and on the relative scopefor activity in larval cyprinids. Journal ofComparative Physiology B 156:791-796.
Winberg, G.G. 1956. Rate of metabolism andfood requirements of fishes. Belorussian
University, Minsk. Translated from Russian:Fisheries Research Board of Canada
Translation Service 194, 1960, Ottawa.

127

128

Appendix 1

List of Nemuro 2002 workshop participants.

Canada
HAY, Douglas E.
Fisheries and Oceans Canada
Pacific Biological Station
3190 Hammond Bay Road
Nanaimo, BC.
Canada. V9R 5K6
hayd@pac.dfo-mpo.gc.ca

Japan
AZUMAYA, Tomomori
Hokkaido National Fisheries Research InstituteKatsurakoi-116,
Kushiro, 085-0802
Japanazumaya@fra.affrc.go.jp
CHIBA, Sanae
Frontier Research System for Global Change
Showa-machi 3173-25, Kanazawaku, Yokohama
Kanagawa, 236-001
Japanchibas@jamstec.go.jp
FUJII, Masahiko
National Institute for Environmental Studies
16-2 Onogawa, TsukubaJapanfujii.masahiko@nies.go.jp.ac.jp
INADA, Masakatsu
Hokkaido University
Minato-cho 3-1-1
Hakodate, 041-8611
Japan
f990061@ec.hokudai.ac.jp
ITO, Shin-ichi
Tohoku National Fisheries Research InstituteShinhama-cho 3-27-5, Shiogama
Miyagi, 985-0001
Japan
goito@affrc.go.jp
KATUKAWA, Toshio
Ocean Research Institute, University of Tokyo
Minamidai-1-15-1, Shinjuku-ku
Tokyo, 164-8639
Japankatukawa@ori.u-tokyo.ac.jp

KASHIWAI, Makoto B.
Hokkaido National Fisheries Research InstituteKatsurakoi-116,
Kushiro, 085-0802
Japankashiwai@fra.affrc.go.jp
KISHI, Michio J.
Hokkaido University
Minato-cho 3-1-1
Hakodate, 041-8611
Japankishi@salmon.fish.hokudai.ac.jp
KOBAYASHI, Tokimasa
Tohoku National Fisheries Research Institute
Samemachi Shimomekurakubo
Hachinohe, Aomori, 031-0841
Japan
tokikoba@fra.affrc.go.jp
MUKAI, Daiki
Hokkaido University
Minato-cho 3-1-1
Hakodate, 041-8611
Japan
f980102@ec.hokidai.ac.jp
OKUNISHI, Takeshi
Hokkaido University
Techo-park 1-2-14 Shinonoppro Atshetsu-ku
Sapporo, Hokkaido, 005-8601
Japan
t-okunishi@econixe.co.jp
SMITH, Lan S.
Frontier Research System for Global Change
Showa-machi 3173-25, Kanazawaku, Yokohama
Kanagawa, 236-001
Japan
lanimal@jamstec.go.jp
TADOKORO, Kazuaki
Frontier Research System for Global Change
Showa-machi 3173-25, Kanazawaku, Yokohama
Kanagawa, 236-001
Japandenden@jamstec.go.jp

129

TOMOKAZU, Aiki
Hokkaido University
Minato-cho 3-1-1
Hakodate, 041-8611
Japan
f990001@ec.hokudai.ac.jp
TSUDA, Atsushi
Hokkaido National Fisheries Research InstituteKatsurakoi-116,
Kushiro, 085-0802
Japan
tsuda@fra.affrc.go.jp
YAMANAKA, Yasuhiro
Graduate School of Environmental Earth Science
Hokkaido University
North 10, West 5 Kita-ku
Sapporo, 060-0810
Japan
galapen@ees.hokudai.ac.jp
YOSHIE, Naoki
Graduate School of Environmental Earth Science
Hokkaido University
North 10, West 5 Kita-ku
Sapporo, 060-0810
Japannaoki@ees.hokudai.ac.jp

People`s Republic of China
HUANG, Daji
2nd Institute of Oceanography
P. O. Box 1207, Hangzhou
Zhejiang, 310012
People`s Republic of China
dajih2001@yahoo.com

Russian FederationZUENKO, Yury I.
Pacific Fisheries Research Center
4 Shevchenko Alley
Vladivostok, 690950
Russiakheng@tinro.ru
LEONOV, Alexander
Shirshov Institute of Oceanology
36 Nakhimovsky Prospekt
Moscow, 117997
Russialeonov@sio.rssi.ru
NAVROTSKY, Vadim V.
Pacific Oceanological Institute
43 Baltiyskaya Street
Vladivostok, 690041
Russianavr@online.vladivostok.ru

U.S.A.
KLUMB, Robert A.
Cornell University
Cornell Biological Field Station
900 Shackelton Point Road
Bridgeport, NY 13030
U.S.A.
rak11@cornell.edu
MEGREY, Bernard A.
National Marine Fisheries Service
Alaska Fisheries Science Center
7600 Sand Point Way N.E.
Seattle, WA 98115-0070
U.S.A.
bern.megrey@noaa.gov

Republic of Korea
HONG, Chul-Hoon
Pukyong National University599-1 Daeyeon3-Dong Nam-Gu
Pusan, 608-737
Republic of Korea
hongch@pknu.ac.kr

WERNER, Francisco E.
Department of Marine Science
University of North Carolina
Chapel Hill
NC 27599-3300
U.S.A.
cisco@marine.unc.edu

130Appendix 2

Herring bioenergetic model FORTRAN code for the base case.

C ------------------------------------------------------------------------C Bioenergetic herring model based on the paper of Rudstam (1988).
C Exploring the dynamics of herring consumption in the Baltic:
C Aplications of an energetic model of fish growth.
C Kieler Meeresforsch Sonderth 6:312-322.
C
C Originally coded as difference equation in FORTRAN by Kenny Rose 26 Dec 01
C
C Corrected and changed to differential equation with Euler andC Runge Kutta numerical integration schemeCC
01/03/02 Bernard A. MegreyC
C Added observed and predicted size at age dataC
01/12/02 Bernard A. MegreyC
C Added YOY formulations per Arrhenius (1998)
C
01/22/02 Bernard A. Megrey
C
C All relic code removed for general distribution at
C Nemuro 2002 workshop
C
01/25/02 Bernard A. MegreyCC
C ------------------------------------------------------------------------program NemuroHerring
include 'stuff.cmn'
include 'state.cmn'
include 'sizeaa.cmn'
REAL NYEARS, NSTEPS, NSTEPOPEN(UNIT=11,FILE='nemuro.txt',STATUS='unknown')
OPEN(UNIT=8,FILE='compareEuler.out',STATUS='unknown')
OPEN(UNIT=9,FILE='sizeatage.out',STATUS='unknown')C
C-----read in the 3 zoop groups from Nemuro output (7th, 9th and 11th columns)C
do 45 ii=1,731
READ(11,999)id(ii),zop1(ii),zop2(ii),zop3(ii)
999 FORMAT(1x,i3,1x,5(13x,1x),2(e13.6,1x,13x,1x),e13.6)
CC----- take the first year for now
C
IF(ii.le.365)then
zoop1(ii)=zop1(ii)zoop2(ii)=zop2(ii)zoop3(ii)=zop3(ii)
endif

131

45

continueC------convert Nemuro zoop in uM N/L to g ww/m3
C------ tt1 is conversion from uM N/liter to g ww/m3
C------ 14 ug N/uM * 1.0e-6 g/ug * 1 g dw/0.07 g N dw
C------ * 1 g ww/0.2 g dw *C------ 1.e3 liters/m3
do 55 i=1,365
tt1=14.0*1.0e-6*(1.0/0.07)*(1.0/0.2)*1.0e3
zoop1(i)=zoop1(i)*tt1
zoop2(i)=zoop2(i)*tt1
zoop3(i)=zoop3(i)*tt1
C----initial weight and age of newly metamphosed herringx(1)=0.2
iage=0maxage=055
continueC
C --- number of state variablesC
nstate=1
C
C - time initializationC
TZERO = 0.0
NYEARS = 11.0
NSTEPS = 100.0
TEND = 365.0*NYEARS
dt= 1/NSTEPSTPRINT = 1.0
TEPS = 1.0E-05
C
C --- MAIN TIME LOOPCNSTEP=0.0
DO WHILE (time .LT. TEND)NSTEP=NSTEP+1.0
time= TZERO + NSTEP*dt
CALL EULER(x,xdot,nstate,dt,time)C
CALL KUTTA(x,xdot,nstate,dt,time)
iday=int(amod(time,365.0))+1
iyr=int(time/365.)+1C
C----- update age every time iday 365 goes by andC----- after NSTEP have gone by
C----- the NSTEPS test is to avoid incrementing iage
C----- every dt
C
IF(iday.eq.365 .and. amod(NSTEP,NSTEPS).eq.0) then132

iage=iage+1
jpage(iage)=iage
psizeaa(iage)=x(1)if(iage .gt. maxpage) maxpage=iage
endifCC --- check for time to printC
IF(ABS(AMOD(time,TPRINT)) .LT. TEPS) THEN
write(8,1001) time, x(1), wtemp, gcmax
1001 format(1x,4(f9.3,1x))
C
ENDIF
ENDIF
END DO
kage=min(maxoage,maxpage)C
C --- write predicted and observed size-at-age output
C
do 73 i=1,kage
write(9,1002) i, psizeaa(i), osizeaa(i)73 continue
1002 format(1x,i4,2(1x,f9.3))
close(8)close(9)STOP
END
SUBROUTINE DER(x,xdot,time)C-----------------------------------------------------C
C Herring bioenergetics differential equation process.
C Prey base are in units of micromoles N /m^3 and are
C converted
C via converson factors
C
C programmed by BAM 01/05/02
C
C Added YOY formulations per Arrhenius (1998)
C
01/22/02 Bernard A. MegreyC
C-----------------------------------------------------include 'stuff.cmn'
include 'state.cmn'
C

WRITE(*,*) 'IN DER', time
iday=int(amod(time,365.0))+1
iyr=int(time/365.)+1

133

C
C zero out xdot
C
DO 15 i=1,nstatexdot(i)=0.0
15 CONTINUE
C---- start age-0 on day 200
C---- jday is julian day (1,..., 400) but goes past 365
C---- iday is counter for day in model simulation
C----- jjday is julian day (i.e., jday reset for >365)
C---- I start on day 200; i you want different then change 200
C--- below and 165, which 365 minus the start day.
jday=iday+200
IF(jday.le.365)then
jjday=jdayelse
jjday=iday-165
endifCC----- generate daily temperatures for a year -- made up
C
t1=float(jjday)
t2=12.75-10.99*cos(0.0172*t1)-6.63*sin(0.0172*t1)
wtemp=t2-5.0
IF(wtemp.le.1.0)wtemp=1.0
C
write(*,*) "wtemp",wtemp
C50 continue
C
C Herring = x(1)C
C
C----- set vulnerabilities and k values for 3 zoop groupsC
vul(1)=1.0
vul(2)=1.0
vul(3)=1.0
k(1)=100.0k(2)=10.0k(3)=100.0
CC-----if using constant p rather than functional response, set pC-----hereC
p=0.6
C----- loop over yearsC
do 100 iyr=1,9134

C---- loop over days for each yearC
do 200 iday=1,365
Cwrite(*,*) 'iday jday jjday', iday,jday,jjday
C-----iiday is running value of days in simulation (1,...., 2000)
C
iiday=(iyr-1)*365+idayC

tt1=1.0/x(1)t1=0.0033*tt1**0.227
write(*,*) 't1 tt1', x(1), t1, tt1c----***this is the new stuff from Arrhenius (1998) for YOY only***c----- The 5.258 puts resp is in units of g zoop/g fish/dayc----- [13560 joules/gram oxygen]/4.18 joules/cal = 3244 cal/gO2c------[2580 joules/gram zoop]/4.18 joules/cal = 617 cal/g zoopc------ 3244/617 = 5.258
IF(iage.eq.0)thenIF(wtemp.le.15.0)then
v=5.76*EXP(0.0238*wtemp)*x(1)**0.386
endif
IF(wtemp.gt.15.0)then
v=8.6*x(1)**0.386
endif
a=EXP((0.03-0.0*wtemp)*v)
resp=t1*EXP(0.0548*wtemp)*a*5.258
endifc------***back to the old equations for respiration for age-1 andC------older ***IF(iage.ge.1)thenIF(wtemp.le.9.0)then
u=3.9*x(1)**0.13*EXP(0.149*wtemp)
else
u=15.0*x(1)**0.13
endif
resp=t1*EXP(0.0548*wtemp)*EXP(0.03*u)*5.258
endifC
write(*,*) 'after new stuff'
C
IF(wtemp.lt.9.0)thenC
u=3.9*x(1)**0.13* exp(0.149*wtemp)
CelseC
u=15.0*x(1)**0.13
Cendifc-- --- 13,560 joules/g O2 1 cal/4.18 joules 1 g ww/5533 calC
resp=t1*EXP(0.0548*wtemp)*EXP(0.03*u)*0.59
C
C----- Thornton and Lessem (1978)temperature effectc----**Arrhenius (1998) for age-0 changed te4 from 25 to 23 degrees***IF(iage.eq.0)then
xk1=0.1135xk2=0.98
xk3=0.98
xk4=0.01
te1=1.0
te2=15.0
te3=17.0
te4=23.0
endifCIF(iage.eq.1)then
xk1=0.1
xk2=0.98
xk3=0.98
xk4=0.01
te1=1.0
te2=15.0
te3=17.0
te4=25.0
endif
IF(iage.gt.1)then
xk1=0.1
xk2=0.98
xk3=0.98
xk4=0.01
te1=1.0
te2=13.0
te3=15.0
te4=23.0
endifCC- non age dependent temperature effect on consumptionC
xk1=0.1C
xk2=0.98
C
xk3=0.98
C
xk4=0.01C
te1=1.0C
te2=13.0C
te3=15.0C
te4=23.0
tt5=(1.0/(te2-te1))
t5=tt5 * alog(0.98*(1.0-xk1)/(0.02*xk1))
t4=exp(t5*(wtemp-te1))
tt7 = 1.0/(te4-te3)t7=tt7*alog(0.98*(1.0-xk4)/(0.02*xk4))
t6=exp(t7*(te4-wtemp))

136

gcta=(xk1*t4)/(1.0+xk1*(t4-1.0))
gctb=xk4*t6/(1.0+xk4*(t6-1.0))
gctemp=gcta * gctb
gcmax=0.642*tt1**0.256*gctemp
C
C----- no tempeature effectC
gcmax=0.642*tt1**0.256*1.0
C----- either use fixed p or call functional responseC
con=p*gcmax
Cwrite(*,*) 'der time iday iyr jjday zoop123',time,iday,
Ciyr,C
jjday, zoop1(jjday),zoop2(jjday), zoop3(jjday)
cnum=zoop1(jjday)*vul(1)/k(1)+zoop2(jjday)*vul(2)/k(2)$+zoop3(jjday)*vul(3)/k(3)c1=gcmax*zoop1(jjday)*vul(1)/k(1)c2=gcmax*zoop2(jjday)*vul(2)/k(2)c3=gcmax*zoop3(jjday)*vul(3)/k(3)con1=c1/(1.0+cnum)
con2=c2/(1.0+cnum)
con3=c3/(1.0+cnum)
con=con1+con2+con3
C
C --- for comparison to Kenny's version
CC
con=0.75*gcmax
CC --- to tune to observed size at age data
con=0.48*gcmax
CC --- egestionC
f=0.16*con
C
C --- excretion
e=0.1*(con-f)C
C --- Specific Dynamic ActionCc------ *******Arrhenius (1998) changed SDA from 17.5% to 15% ****IF(iage.eq.0)sda=0.15*(con-f)
IF(iage.ge.1)sda=0.175*(con-f)c----- J/g ww 1 cal=4.18 JC
con1=con*2580.0/5533.0
C
write(*,*) 'der',con1,resp,xdot(1)C
C --- bioenergetics differential equationC
xdot(1)=(con-resp-f-e-sda)*x(1)*2580./5533.

137

IF(wtemp.le.1.0)xdot(1)=0.0
t1=float(jjday)
if(amod(t1,365.0).ge.152.0.and.amod(t1,365.0).le.156.0) thenwrite(*,*) 'in spawn'
xdot(1)=(con-resp-f-e-sda-0.20)*x(1)*2580./5533.
endifc----- update age every time day 365 goes bycIF(iday.eq.365) thenciage=iage+1cjage(iage)=iagec
psizeaa(iage)=x(1)cwrite(*,*)'iday, iyr, jjday, iage, maxage',iday, iyr,c $jjday, iage, maxagecif(iage .gt. maxage) maxage=iagec
c
endif
C

WRITE(*,*) 'OUT OF DER'
RETURNEND
SUBROUTINE EULER(x,xdot,nstate,dt,time)C-----------------------------------------------------C
C USE THE EULER METHOD TO SOLVE A SYSTEM OF NONLINEARC DIFFERENTIAL EQUATIONS. A SUBROUTINE DER IS NEEDEDC TO COMPUTE THE DERIVATIVS OF THE STATE VARIBALESCC X - STATE VARIABLE ARRAYC XDOT - ARRAY OF DERIVATIVES OF STATE ARIABLESC NSTATE - NUMBER OF STATE VARIABLESC DT - TIME STEPC TIME - CURRENT TIMEC
C-----------------------------------------------------include 'stuff.cmn'
include 'state.cmn'
INTEGER ICALL DER(x,xdot,time)
DO 10 i=1, nstatex(i) =x(i) + dt * xdot(i)10 CONTINUE
RETURNEND
SUBROUTINE KUTTA(x,xdot,nstate,dt,time)138

C-----------------------------------------------------C
C USE THE 4TH ORDER RUNGE KUTTA METHOD TO SOLVE A SYSTEM OF NONLINEARC DIFFERENTIAL EQUATIONS. A SUBROUTINE DER IS NEEDEDC TO COMPUTE THE DERIVATIVS OF THE STATE VARIBALESCC X - STATE VARIABLE ARRAYC XDOT - ARRAY OF DERVATIVES OF STATE VARIABLESC NSTATE - NUMBER OF STATE VARIABLESC DT - TIME STEPC TIME - CURRENT TIMEC
C programmed by Bernard A. Megrey 01/06/02
C
C
C-----------------------------------------------------include 'stuff.cmn'
include 'state.cmn'
INTEGER IREAL SUMDX(16), DTO2, XPLUS(16)C
WRITE(*,*) 'IN KUTTA'
DTO2 = dt/2.0
CALL DER(x,xdot,time)
DO 10 I=1, nstate
XPLUS(I) = x(I) + DT02 * xdot(I)
SUMDX(I)= xdot(I)10 CONTINUECALL DER(XPLUS,xdot,time)
DO 20 I=1, nstate
XPLUS(I) = x(I) + DTO2 * xdot(I)
SUMDX(I) = SUMDX(I) + 2.0 * xdot(I)20 CONTINUECALL DER(XPLUS,xdot,time)
DO 30 I=1, nstate
XPLUS(I) = x(I) + dt * xdot(I)
SUMDX(I) = SUMDX(I) + 2.0 * xdot(I)30 CONTINUECALL DER(XPLUS,xdot,time)
DO 40 I=1, nstate
SUMDX(I) = SUMDX(I) + xdot(I)
x(I) = x(I) + dt * SUMDX(I) / 6.0

139

C

40 CONTINUEWRITE(*,*) 'OUT OF KUTTA'
RETURN
END

C==========================================================CC Include file: state.cmn
C state variable declarationsC
C BAM 1/02/02
C---------------------------------------------------------REAL*4 dt, time, x(1), xdot(1)INTEGER*4 nstate
C==========================================================CC Include file: stuff.cmn
C miscellaneous common block variablesC
C BAM 1/2/02
C---------------------------------------------------------COMMON /ZOO/ zoop1(365), zoop2(365), zoop3(365)COMMON /ISV/ wtemp, con, gcmax,f, e ,sda, con1, resp
COMMON /TIMER/ iday, jday, iyr, iage
REAL*4 zoop1,zoop2,zoop3,k(3),vul(3)REAL*4 zop1(731),zop2(731),zop3(731)
INTEGER*4 id(731)
C==========================================================CC Include file: sizeaa.cmn
C Size at age common blockC
C BAM 1/06/02
C
C---------------------------------------------------------COMMON /SIZEAA/joage(25),jpage(25),psizeaa(25),osizeaa(25),$maxoage,maxpage
C --- joage Observed age
C --- jpage Predicted age
C --- psizeaa Predicted size at age
C --- osizeaa Observed size at age
INTEGER*4 joage, jpage, maxoage, maxpage
REAL*4 psizeaa, osizeaa
DATA maxoage /11/
DATA joage /1,2,3,4,5,6,7,8,9,10,11,14*0/

140

CC --- observed size at age Pacific herring data (wt) from the C
C --- 1973 year class
C --- as supplied by Doug Hay
C
DATA osizeaa /8.75,63.0,87.66,124.13,139.26,152.98,177.43,$
185.78,187.64,195.0,208.73,14*0.0/

141

142

Appendix 3 Fully customized herring subroutine. See Appendix 2 for main program andcommon block including files.
SUBROUTINE DER(x,xdot,time)C-----------------------------------------------------C
C Herring bioenergetics differential equation process.
C Prey base are in units of micromoles N /m^3 and are converted
C via converson factors
C
C programmed by B.A. Megrey 01/05/02
C
C Added YOY formulations per Arrhenius 1998
C
01/22/02 B.A. Megrey
C modifications and customizations by F.E.Werner and R.A.Klumb 1/26/02
C while at the Nemuro workshop
C
C-----------------------------------------------------include 'stuff.cmn'
include 'state.cmn'
CC- calculate date and yearC
iday=int(amod(time,365.0))+1
iyr=int(time/365.)+1C
C zeroout xdot
C
DO 15 i=1,nstatexdot(i)=0.0
15 CONTINUE
cc---- start age-0 on day 200c---jday is julian day (1,..., 400) but goes past 365
c---iday is counter for day in model simulationc---jjday is julian day (i.e., jday reset for >365)c---I start on day 200; i you want different then change 200
c---below and 165, which 365 minus the start day.
jday=iday+200
IF(jday.le.365)then
jjday=jdayelse
jjday=iday-165
endifCC----- generate daily temperatures for a year -- made up
C143

t1=float(jjday)! BaseCase T
t2=12.75-10.99*cos(0.0172*t1)-6.63*sin(0.0172*t1)! BaseCase T
wtemp=t2-5.0! BaseCase T
IF(wtemp.le.1.0)wtemp=1.0! BaseCase T
pi=acos(-1.)
wtemp=7.717+(5.6796*0.5*(1.-cos(2.*pi*(t1-30.)/365.)))C
C--- Herring weight state variable = x(1)C--- weight affect on respirationC
tt1=1.0/x(1)t1=0.0033*tt1**0.227
CC --- *********this is the new stuff from Arrhenius (1998) for YOY only*********
C --- The 5.258 puts resp (g oxygen/fish) into units of g zoop/g fish/day
C --- [13560 joules/gram oxygen]/4.18 joules/cal = 3244 cal/gO2
C --- [2580 joules/gram zoop]/4.18 joules/cal = 617 cal/g zoop
C --- so respiration in grams/oxygen/g fish/day is multiplied
C --- by 3244/617 = 5.258
C --- to get food energy equivalents of a gram of oxygen respiredCc
IF(iage.eq.0)then! BASEc
IF(wtemp.le.15.0)then! BASEcv=5.76*EXP(0.0238*wtemp)*x(1)**0.386! BASEcendif! BASEc
IF(wtemp.gt.15.0)then! BASEcv=8.6*x(1)**0.386! BASEcendif! BASEc
a=EXP((0.03-0.0*wtemp)*v)! BASEc
resp=t1*EXP(0.0548*wtemp)*a*5.258! BASEcendif! BASEC
IF(iage.eq.0)then! R.A. Klumb (26 Jan 2002)
t1=0.00528*tt1**0.007! R.A. Klumb (26 Jan 2002)
resp=t1*EXP(0.083*wtemp)*5.258! R.A. Klumb (26 Jan 2002)
end if! R.A. Klumb (26 Jan 2002)
C
C --- *********back to the old equations for respiration for age-1 andC --- older********
C
IF(iage.ge.1)then
C BaseIF(wtemp.le.9.0)thenIF(wtemp.le.9.655)then
u=3.9*x(1)**0.13*EXP(0.149*wtemp)
else
u=15.0*x(1)**0.13
endif
resp=t1*EXP(0.0548*wtemp)*EXP(0.03*u)*5.258
endifCC --- Thornton and Lessem temperature effect144

C --- age dependent values
C --- *******Arrhenius (1998)for age-0 changed te4 from 25 to 23 degrees******
C
IF(iage.eq.0)then
xk1=0.1
xk2=0.98
xk3=0.98
xk4=0.01
C BaseC BaseC BaseC Basete1=1.0
te2=15.0
te3=17.0
te4=23.0

te1=8.0
te2=10.897
te3=11.310
te4=12.552
endifCIF(iage.eq.1)then
xk1=0.1
xk2=0.98
xk3=0.98
xk4=0.01
C BaseC BaseC BaseC Basete1=1.0
te2=15.0
te3=17.0
te4=25.0

te1=8.0
te2=10.897
te3=11.310
te4=12.966
endif
IF(iage.gt.1)then
xk1=0.1
xk2=0.98
xk3=0.98
xk4=0.01
C BaseC BaseC BaseC Basete1=1.0
te2=13.0
te3=15.0
te4=23.0
te1=8.0
te2=10.483
te3=10.897

145te4=12.552
endif
tt5=(1.0/(te2-te1))
t5=tt5 * alog(0.98*(1.0-xk1)/(0.02*xk1))
t4=exp(t5*(wtemp-te1))
tt7 = 1.0/(te4-te3)t7=tt7*alog(0.98*(1.0-xk4)/(0.02*xk4))
t6=exp(t7*(te4-wtemp))
gcta=(xk1*t4)/(1.0+xk1*(t4-1.0))
gctb=xk4*t6/(1.0+xk4*(t6-1.0))
gctemp=gcta * gctb
gcmax=0.642*tt1**0.256*gctemp
C
C --- multispecies functional response
C --- usse either this or adjust little pC
C----- set vulnerabilities and k values for 3 zoop groupsC
vul(1)=1.0
vul(2)=1.0
vul(3)=1.0
k(1)=0.3638
k(2)=0.0364
k(3)=0.3638ccck(1)=2000.0k(2)=200.0k(3)=2000.0
cnum=zoop1(jjday)*vul(1)/k(1)+zoop2(jjday)*vul(2)/k(2)$+zoop3(jjday)*vul(3)/k(3)c1=gcmax*zoop1(jjday)*vul(1)/k(1)c2=gcmax*zoop2(jjday)*vul(2)/k(2)c3=gcmax*zoop3(jjday)*vul(3)/k(3)con1=c1/(1.0+cnum)
con2=c2/(1.0+cnum)
con3=c3/(1.0+cnum)
con=con1+con2+con3

C
C-----if using constant p rather than functional response, set p here
C --- to tune to observed size at age data. If using functional response
C --- comment the next line outC
con=0.6375*gcmax
CC --- egestionC146

f=0.16*con! Base CaseIF(iage.eq.0)f=0.125*con! Age-dependent – R.A. Klumb (26 Jan 2002)
C
C --- excretion
e=0.1*(con-f)! Base CaseIF(iage.eq.0)e=0.078*con! Age-dependent – R.A. Klumb (26 Jan 2002)
C
C --- Specific Dynamic ActionCc------ *******Arrhenius (1998) age dependent SDA from 17.5% to 15% ****IF(iage.eq.0)sda=0.15*(con-f)! Base CaseIF(iage.eq.0)sda=0.125*(con-f) ! Age-dependent – R.A. Klumb (26 Jan 2002)
IF(iage.ge.1)sda=0.175*(con-f)C
C --- use the ratio of calories/g of zoop (2580) to calories/g of fish (5533)
C
C --- bioenergetics differential equation - constant energy density for herringCC
xdot(1)=(con-resp-f-e-sda)*x(1)*2580./5533.
! Base CaseC
C include seasonal variation of energy density for .ge. 2 yr oldsC
if(iage.ge.2)then
enMar1=5750.
jdMar1=60
enOct1=9800.
jdOct1=274
if(jjday.lt.60)then
delen=(enMar1-enOct1)/151
en=enOct1+(90+jjday)*delenend ifif(jjday.ge.60.and.jjday.lt.274)then
delen=(enOct1-enMar1)/(jdOct1-jdMar1)
en=enMar1+(jjday-jdMar1)*delenend ifif(jjday.ge.274)then
delen=(enMar1-enOct1)/151
en=enOct1+(jjday-jdOct1)*delenend ifelse
en=4460.
end ifxdot(1)=(con-resp-f-e-sda)*x(1)*2580./en
IF(wtemp.le.1.0)xdot(1)=0.0
CC --- Spawning section. Assume loose 20% of body weight/day
C
t1=float(jjday)cif(amod(t1,365.0) .ge. 152.0 .and.
c &
amod(t1,365.0) .le. 156.0) thencxdot(1)=(con-resp-f-e-sda-0.20)*x(1)*2580./5533.

147

c

endif
RETURN
END

148

Appendix 4 NEMURO.FISH (NEMURO FORTRAN code supplied by Yasuhiro Yamanaka)with the herring bioenergetic model (base case) (supplied by Bernard Megrey and Ken Rose).
The herring model is linked to NEMURO in a one-way static link.
!***********************************************************************! NEMURO model Jun 13, 2002 written by Yasuhiro Yamanaka!***********************************************************************
program NEMURO.FISH
implicit none
! ..... Control for Time Ingration .....
character(19) :: Cstart = '0001/07/20 00:00:00' ! Starting datecharacter(19) :: Cend = '0011/07/21 00:00:00' ! Ending datecharacter(19) :: Cstep = '0000/00/00 01:00:00' ! Time step
character(19) :: Cmon = '0000/00/01 00:00:00' ! Monitor Interval
character(19) :: CTime
real(8):: dt, TTime, Tbefore, Season, Tmon
integer:: Iyr, Imon, Iday, Ihour, Imin, Isec
! ..... scale conversion .....
real(8),parameter :: d2s= 86400.0d0 ! day ---> secreal(8),parameter :: mcr= 1.0d-6! micro! ..... Prognostic Variables (with initial conditions) and Thier Source Term .....
real(8):: TPS = 0.1D-6, QPS ! Small Phytoplankton [molN/l]real(8):: TPL = 0.1D-6, QPL ! Large Phytoplanktonreal(8):: TZS = 0.1D-6, QZS ! Small Zooplanktonreal(8):: TZL = 0.1D-6, QZL ! Large Zooplanktonreal(8):: TZP = 0.1D-6, QZP ! Pradatory Zooplanktonreal(8):: TNO3 = 5.0D-6, QNO3 ! Nitratereal(8):: TNH4 = 0.1D-6, QNH4 ! Ammmonium
real(8):: TPON = 0.1D-6, QPON ! Particulate Organic Nitrogenreal(8):: TDON = 0.1D-6, QDON ! dissolved Organic Nitrogenreal(8):: TSiOH4 = 10.0D-6, QSiOH4 ! Silicate
real(8):: TOpal = 0.1D-7, QOpal ! Particulate Opal
! ..... Prognostic Variables (with initial conditions) and Thier Source Term .....
!
real(8):: THrr = 0.2D0, QHrrl ! Particulate Opal
! ...... Light Condition Parameters .......
real(8),parameter :: alpha1 = 4.0D-2 ! Light Dissipation coefficient of sea water[/m]real(8),parameter :: alpha2 = 4.0D4! PS+PL Selfshading coefficientS+PL[l/molN/m]real(8),parameter :: IoptS = 0.15D0 ! PS Optimum Light Intensity S[ly/min]real(8),parameter :: IoptL = 0.15D0 ! PL Optimum Light Intensity[ly/min]integer,parameter :: LLN= 10! Number of sublayer for calculating of Lfc
real(8):: LfcS! Light factor for PS
real(8):: LfcL! Light factor for PLreal(8):: kappa, Lint, dLint, LfcUS, LfcUL, LfcDS, LfcDL
integer:: L! ...... biological Parameters ......
real(8),parameter :: VmaxS = 0.4D0/d2s ! PS Maximum Photosynthetic rate @0degC [/s]real(8),parameter :: KNO3S = 1.0D-6 ! PS Half satuation constant for Nitrate [molN/l]real(8),parameter :: KNH4S = 0.1D-6 ! PS Half satuation constant for Ammonium [molN/l]real(8),parameter :: PusaiS = 1.5D6! PS Ammonium Inhibition Coefficient[l/molN]real(8),parameter :: KGppS = 6.93D-2 ! PS Temp. Coeff. for Photosynthetic Rate [/degC]real(8),parameter :: MorPS0 = 5.85D4/d2s ! PS Mortality Rate @0degC
[/s]

149

real(8),parameter :: KMorPS = 6.93D-2 ! PS Temp. Coeff. for Mortality[/degC]real(8),parameter :: ResPS0 = 0.03D0/d2s ! PS Respiration Rate at @0degC[/s]real(8),parameter :: KResPS = 0.0519D0 ! PS Temp. Coeff. for Respiration[/degC]real(8),parameter :: GammaS = 0.135D0 ! PS Ratio of Extracell. Excret. to Photo.[(nodim)]real(8),parameter :: VmaxL = 0.8D0/d2s ! PL Maximum Photosynthetic rate @0degC [/s]real(8),parameter :: KNO3L = 3.00D-6 ! PL Half satuation constant for Nitrate [molN/l]real(8),parameter :: KNH4L = 0.30D-6 ! PL Half satuation constant for Ammonium [molN/l]real(8),parameter :: KSiL = 6.00D-6 ! PL Half satuation constant for Silicate [molSi/l]real(8),parameter :: PusaiL = 1.50D6 ! PL Ammonium Inhibition Coefficient[l/molN]real(8),parameter :: KGppL = 6.93D-2 ! PL Temp. Coeff. for Photosynthetic Rate [/degC]real(8),parameter :: MorPL0 = 2.90D4/d2s ! PL Mortality Rate @0degC[/s]real(8),parameter :: KMorPL = 6.93D-2 ! PL Temp. Coeff. for Mortality[/degC]real(8),parameter :: ResPL0 = 0.03D0/d2s ! PL Respiration Rate at @0degC[/s]real(8),parameter :: KResPL = 0.0519D0 ! PL Temp. Coeff. for Respiration[/degC]real(8),parameter :: GammaL = 0.135D0 ! PL Ratio of Extracell. Excret. to Photo.[(nodim)]real(8),parameter :: GRmaxS = 0.40D0/d2s ! ZS Maximum Rate of Grazing PS @0degC [/s]real(8),parameter :: KGraS = 6.93D-2 ! ZS Temp. Coeff. for Grazing[/degC]real(8),parameter :: LamS = 1.40D6 ! ZS Ivlev constant[l/molN]real(8),parameter :: PS2ZSstar= 0.043D-6 ! ZS Threshold Value for Grazing PS[molN/l]real(8),parameter :: AlphaZS = 0.70D0 ! ZS Assimilation Efficiency[(nodim)]real(8),parameter :: BetaZS = 0.30D0 ! ZS Growth Efficiency[(nodim)]real(8),parameter :: MorZS0 = 5.85D4/d2s ! ZS Mortality Rate @0degC[/s]real(8),parameter :: KMorZS = 0.0693D0 ! ZS Temp. Coeff. for Mortality[/degC]real(8),parameter :: GRmaxLps = 0.10D0/d2s ! ZL Maximum Rate of Grazing PS @0degC [/s]real(8),parameter :: GRmaxLpl = 0.40D0/d2s ! ZL Maximum Rate of Grazing PL @0degC [/s]real(8),parameter :: GRmaxLzs = 0.40D0/d2s ! ZL Maximum Rate of Grazing ZS @0degC [/s]real(8),parameter :: KGraL = 6.93D-2 ! ZL Temp. Coeff. for Grazing[/degC]real(8),parameter :: LamL = 1.4000D6 ! ZL Ivlev constant[l/molN]real(8),parameter :: PS2ZLstar= 4.00D-8 ! ZL Threshold Value for Grazing PS[molN/l]real(8),parameter :: PL2ZLstar= 4.00D-8 ! ZL Threshold Value for Grazing PL[molN/l]real(8),parameter :: ZS2ZLstar= 4.00D-8 ! ZL Threshold Value for Grazing ZS[molN/l]real(8),parameter :: AlphaZL = 0.70D0 ! ZL Assimilation Efficiency[(nodim)]real(8),parameter :: BetaZL = 0.30D0 ! ZL Growth Efficiency[(nodim)]real(8),parameter :: MorZL0 = 5.85D4/d2s ! ZL Mortality Rate @0degC[/s]real(8),parameter :: KMorZL = 0.0693D0 ! ZL Temp. Coeff. for Mortality[/degC]real(8),parameter :: GRmaxPpl = 0.20D0/d2s ! ZP Maximum rate of grazing PL @0degC [/s]real(8),parameter :: GRmaxPzs = 0.20D0/d2s ! ZP Maximum rate of grazing ZS @0degC [/s]real(8),parameter :: GRmaxPzl = 0.20D0/d2s ! ZP Maximum rate of grazing ZL @0degC [/s]real(8),parameter :: KGraP = 6.93D-2 ! ZP Temp. Coeff. for grazing[/degC]real(8),parameter :: LamP = 1.4000D6 ! ZP Ivlev constant[l/molN]real(8),parameter :: PL2ZPstar= 4.00D-8 ! ZP Threshold Value for Grazing PL[molN/l]real(8),parameter :: ZS2ZPstar= 4.00D-8 ! ZP Threshold Value for Grazing ZS[molN/l]real(8),parameter :: ZL2ZPstar= 4.00D-8 ! ZP Threshold Value for Grazing ZL[molN/l]real(8),parameter :: PusaiPL = 4.605D6 ! ZP Preference Coeff. for PL[l/molN]real(8),parameter :: PusaiZS = 3.010D6 ! ZP Preference Coeff. for ZS[l/molN]real(8),parameter :: AlphaZP = 0.70D0 ! ZP Assimilation Efficiency[(nodim)]real(8),parameter :: BetaZP = 0.30D0 ! ZP Growth Efficiency[(nodim)]real(8),parameter :: MorZP0 = 5.85D4/d2s ! ZP Mortality Rate @0degC[/s]real(8),parameter :: KMorZP = 0.0693D0 ! ZP Temp. Coeff. for Mortality[/degC]real(8),parameter :: Nit0 = 0.03D0/d2s ! NH4 Nitrification Rate @0degC[/s]real(8),parameter :: KNit = 0.0693D0 ! NH4 Temp. coefficient for Nitrification [/degC]

150

!

!

real(8),parameter :: VP2N0 = 0.10D0/d2s ! PON Decomp. Rate to Ammonium @0degC [/s]real(8),parameter :: KP2N = 6.93D-2 ! PON Temp. Coeff. for Decomp. to Ammon. [/degC]real(8),parameter :: VP2D0 = 0.10D0/d2s ! PON Decomp. Rate to DON @0degC[/s]real(8),parameter :: KP2D = 6.93D-2 ! PON Temp. Coeff. for Decomp. to DON [/degC]real(8),parameter :: VD2N0 = 0.20D0/d2s ! DON Decomp. Rate to Ammonium @0degC [/s]real(8),parameter :: KD2N = 6.93D-2 ! DON Temp. Coeff. for Decomp. to Ammon. [/degC]real(8),parameter :: VO2S0 = 0.10D0/d2s ! Opal Decomp. Rate to Silicate @0degC [/s]real(8),parameter :: KO2S = 6.93D-2 ! Opal Temp. Coeff. for Decomp.to Silicate[/degC]real(8),parameter :: RSiN = 2.0D0!Si/N ratio[molSi/molN]
real(8),parameter :: RCN= 106.0D0/16.0D0 !C/N ratio[molC/molN]..... bottom boundary Condition .....
real(8),parameter :: setVPON = 40.0D0/d2s ! Settling velocity of PON[m/s]real(8),parameter :: setVOpal = 40.0D0/d2s ! Settling velocity of Opal[m/s]real(8),parameter :: TNO3d = 25.0d-6 ! Nitrate Concentraion in the Deep Layer [molN/l]real(8),parameter :: TSiOH4d = 35.0d-6 ! Silicate Concentraion in the Deep Layer [molSi/l]...... Paramters of ZL Vertical Migration .....
character(19) :: CZup ='0000/04/01 00:00:00' ! Date Coming up to the Euphotic Layer
character(19) :: CZdwn='0000/09/01 00:00:00' ! Date Returning to the Deep Layer
real(8):: TZup, TZdwn, TZLd, SVRate=0.2D0
integer:: IyrU, ImonU, IdayU ,IhourU, IminU, IsecU
integer:: IyrD, ImonD, IdayD ,IhourD, IminD, IsecD!
real(8)real(8)real(8)real(8)real(8)real(8)real(8)real(8)real(8)integer
!
!

!

!
!
!

:: GppPSn, GppNPSn, GppAPSn, RnewS, ResPSn, MorPSn, ExcPSn:: GppPLn, GppNPLn, GppAPLn, RnewL, ResPLn, MorPLn, ExcPLn:: GppPLsi, GppSiPLsi,
ResPLsi, MorPLsi, ExcPLsi:: GraPS2ZSn, GraPS2ZLn, GraPL2ZLn, GraPL2ZLsi, GraZS2ZLn:: GraPL2ZPn, GraPL2ZPsi, GraZS2ZPn, GraZL2ZPn:: EgeZSn, MorZSn, ExcZSn, EgeZLn, EgeZLsi, MorZLn, ExcZLn::
EgeZPn, EgeZPsi, MorZPn, ExcZPn:: DecP2N, DecP2D, DecD2N, DecO2S, Nit:: ExpPON, ExpOpal,ExcNO3, ExcSiOH4:: lt=0 , nt..... Environmental Condition .....
real(8):: Temp! Temperature [degC]real(8):: Lint0! Light Intencity at sea surface [ly/min]real(8):: MLD = 30.0d0! Mixed Layer Depth [m]real(8):: ExcTime = 1.0d0 / (100.0d0*d2s) ! Exch. Coeff. between Sur-Deep [/s]..... statement function & def. type of functions .....
real(8):: cd2tt, nd2tt
character(19) :: tt2cd
real(8):: Td, GraF, Mich, a, b, c
Td (a,b) = a * exp(b*Temp)
GraF(a,b,c) = MAX( 0.0D0, 1.0 - exp(a * (b - c)))Mich(a,b) = b / ( a + b )***** Initial Setting *****..... for time control .....
TTime = cd2tt(Cstart)! Starting Date
CTime = TT2CD(TTime)! present time (charactor form)
dt = cd2tt(Cstep) - cd2tt('0000/00/00 00:00:00') ! Time Step (real8 form)

151

!

!

!
!
!

!
!

!
!

Tmon = cd2tt(Cmon) - cd2tt('0000/00/00 00:00:00') ! Monitor Interval (real8 form)
nt = NINT( ( cd2tt(Cend) - cd2tt(Cstart) ) / dt ) ! Total Time Steps..... for Vertical Migration ......
TZup = CD2TT( CZup )TZdwn = CD2TT( CZdwn )call TT2ND(IyrU, ImonU, IdayU ,IhourU, IminU, IsecU ,TZup )call TT2ND(IyrD, ImonD, IdayD ,IhourD, IminD, IsecD ,TZdwn)
TZLd = TZL ! ZL living in the deep layer at the initial condition
TZL = 0.0..... File Open for monitoring output .....
open( 10, file='Results.csv', form='FORMATTED' )write(10,'(A,13(",", A))') 'Time(day)', &'NO3' , 'NH4' , 'PS' , 'PL',&
'ZS' , 'ZL'
, 'ZP' , 'PON' , &'DON' , 'SiOH4' , 'Opal' , 'TotalN' , 'TotalSi'
write(10,'(A,11(",", F8.4))') CTime, &TNO3/mcr, TNH4 /mcr, TPS /mcr, TPL /mcr, &TZS /mcr, TZL /mcr, TZP /mcr, TPON/mcr, &TDON/mcr, TSiOH4/mcr, TOpal/mcr
open( 11, file='Forcing.csv', form='FORMATTED' )write(11,'(A,13(",", A))') 'Time(day)','Lint0', 'TMP', 'MLD', 'ExcTime'write(11,'(A,13(",", 1PE10.4))') CTime, Lint0, Temp, MLD, ExcTime*d2s****** Main Loop *****
do lt = 1, nt..... time control (Season : 0 to 1, parcentage in a year).....
Tbefore = TTime! one step before present time
TTime = TTime + dt
! present time (real8 form)
CTime = TT2CD(TTime)! present time (charactor form)
CALL TT2ND(Iyr, Imon, Iday ,Ihour, Imin, Isec ,TTime)
Season = ( TTime- ND2TT(Iyr ,1,1,0,0,0) )/ &( ND2TT(Iyr+1,1,1,0,0,0) - ND2TT(Iyr ,1,1,0,0,0) )..... Example of Boundray condition .....
Lint0 = 0.1d0 * ( 1.0D0 + 0.3d0 * cos( 2.0d0*3.1415926536d0*(Season - 0.50D0) ) )Temp =
6.0D0 + 4.0d0 * cos( 2.0d0*3.1415926536d0*(Season - 0.65D0) )if (Temp .lt. 4.0 ) then
MLD = MLD + dt * (150.0d0 - MLD ) / ( 100.0d0 * d2s )ExcTime = ExcTime + dt * ( 1.0d0/( 40.0d0*d2s) - ExcTime ) / (100.0d0*d2s)else
MLD = MLD + dt * ( 30.0d0 - MLD ) / ( 5.0d0 * d2s )ExcTime = ExcTime + dt * ( 1.0d0/(100.0d0*d2s) - ExcTime ) / ( 5.0d0*d2s)
end if..... Light Factors (LfcS, LfcL).....
Lint = Lint0
LfcDS = Lint/IoptS * exp(1.0D0 - Lint/IoptS)
LfcDL = Lint/IoptL * exp(1.0D0 - Lint/IoptL)
LfcS = 0.0D0
LfcL = 0.0D0
Kappa = alpha1 + alpha2 * ( TPS + TPL )

152

!

!

!

!
!
!

!

dLint = exp( -Kappa * (MLD/LLN) )do L = 1, LLN
LfcUS = LfcDS
LfcUL = LfcDL
Lint = Lint
* dLint
LfcDS = Lint/IoptS * exp( 1.0D0 - Lint/IoptS )LfcDL = Lint/IoptL * exp( 1.0D0 - Lint/IoptL )LfcS = LfcS + ( LfcUS + LfcDS ) * 0.5D0 / LLN
LfcL = LfcL + ( LfcUL + LfcDL ) * 0.5D0 / LLN
end do..... Photosynthesis of PS .....
GppNPSn = Mich( KNO3S, TNO3 ) * exp( - PusaiS * TNH4 )GppAPSn = Mich( KNH4S, TNH4 )GppPSn = Td(VmaxS, KGppS) * LfcS * TPS * ( GppNPSn + GppAPSn )RnewS = GppNPSn / ( GppNPSn + GppAPSn )ResPSn = Td( ResPS0, KResPS ) * TPSMorPSn = Td( MorPS0, KMorPS ) * TPS * TPS
ExcPSn = GammaS* GppPSn..... Photosynthesis of PL .....
GppNPLn = Mich( KNO3L, TNO3 ) * exp( - PusaiL * TNH4 )GppAPLn = Mich( KNH4L, TNH4 )GppSiPLsi = Mich( KSiL , TSiOH4 )GppPLn = Td(VmaxL, KGppL) * LfcL * TPL * min( ( GppNPLn + GppAPLn ), GppSiPLsi )RnewL = GppNPLn / ( GppNPLn + GppAPLn )ResPLn = Td( ResPL0, KResPL ) * TPLMorPLn = Td( MorPL0, KMorPL ) * TPL * TPL
ExcPLn = GammaL* GppPLn..... Grazing PS, PL, ZS, ZL --> ZS, ZL, ZP .....
GraPS2ZSn = Td(GRmaxS, KGraS) *GraF(LamS,PS2ZSstar,TPS) *TZSGraPS2ZLn = Td(GRmaxLps,KGraL) *GraF(LamL,PS2ZLstar,TPS) *TZLGraPL2ZLn = Td(GRmaxLpl,KGraL) *GraF(LamL,PL2ZLstar,TPL) *TZLGraZS2ZLn = Td(GRmaxLzs,KGraL) *GraF(LamL,ZS2ZLstar,TZS) *TZLGraPL2ZPn = Td(GRmaxPpl,KGraP) *GraF(LamP,PL2ZPstar,TPL) *TZP * exp( -PusaiPL *(TZL+ TZS))GraZS2ZPn = Td(GRmaxPzs,KGraP) *GraF(LamP,ZS2ZPstar,TZS) *TZP * exp( -PusaiZS * TZL )GraZL2ZPn = Td(GRmaxPzl,KGraP) *GraF(LamP,ZL2ZPstar,TZL) *TZP..... Mortality, Excration, Egestion for Zooplanktons..... Commented out after Saito-san Meeting at 19 Jun, 2000 .....
BetaZS = 0.3 ** ( 1.0 + Mich( TPL, TPS ) )ExcZSn = (AlphaZS- BetaZS) * GraPS2ZSn
EgeZSn = (1.0 - AlphaZS) * GraPS2ZSn
MorZSn = Td( MorZS0, KMorZS ) * TZS * TZSExcZLn = (AlphaZL- BetaZL) * (GraPS2ZLn+GraPL2ZLn+GraZS2ZLn)
EgeZLn = (1.0 - AlphaZL) * (GraPS2ZLn+GraPL2ZLn+GraZS2ZLn)
MorZLn = Td( MorZL0, KMorZL ) * TZL * TZLExcZPn = (AlphaZP- BetaZP) * (GraPL2ZPn+GraZS2ZPn+GraZL2ZPn)
EgeZPn = (1.0 - AlphaZP) * (GraPL2ZPn+GraZS2ZPn+GraZL2ZPn)
MorZPn = Td( MorZP0, KMorZP ) * TZP * TZP..... Decomposition PON, DON, Opal ---> NH4, DON, SiOH4 .....
DecP2N = Td( VP2N0 , KP2N ) * TPON
DecP2D = Td( VP2D0 , KP2D ) * TPON

153

!

!
!

!
!

!
!

DecD2N = Td( VD2N0 , KD2N ) * TDON
DecO2S = Td( VO2S0 , KO2S ) * TOpal
Nit= Td( Nit0 , KNit ) * TNH4..... silica fluxes .....
GppPLsi = GppPLn * RSiNResPLsi = ResPLn * RSiNMorPLsi = MorPLn * RSiNExcPLsi = ExcPLn * RSiNGraPL2ZLsi = GraPL2ZLn * RSiNGraPL2ZPsi = GraPL2ZPn * RSiNEgeZLsi = GraPL2ZLsi
EgeZPsi = GraPL2ZPsi..... Tendency Terms for biological processes .....
QNO3 = -( GppPSn - ResPSn ) * RnewS &-( GppPLn - ResPLn ) * RnewL + Nit
QNH4 = -( GppPSn - ResPSn ) * (1.0 - RnewS) &-( GppPLn - ResPLn ) * (1.0 - RnewL) &- Nit + DecP2N + DecD2N + ExcZSn + ExcZLn + ExcZPn
QPS = GppPSn - ResPSn - MorPSn - ExcPSn - GraPS2ZSn - GraPS2ZLn
QPL = GppPLn - ResPLn - MorPLn - ExcPLn - GraPL2ZLn - GraPL2ZPn
QZS = GraPS2ZSn - GraZS2ZLn - MorZSn - ExcZSn - EgeZSn - GraZS2ZPn
QZL = GraPL2ZLn + GraZS2ZLn - MorZLn - ExcZLn - EgeZLn + GraPS2ZLn - GraZL2ZPn
QZP = GraPL2ZPn + GraZS2ZPn - MorZPn - ExcZPn - EgeZPn + GraZL2ZPn
QPON = MorPSn + MorPLn + MorZSn + MorZLn + MorZPn &+ EgeZPn + EgeZSn + EgeZLn - DecP2N - DecP2D
QDON = ExcPSn + ExcPLn + DecP2D - DecD2N
QSiOH4 =-GppPLsi + ResPLsi + ExcPLsi + DecO2S
QOpal = MorPLsi + EgeZLsi + EgeZPsi - DecO2S..... Exchange Fluxes between the Surface and Deep Layers .....
ExpPON = setVPON / MLD * TPON
ExpOpal = setVOpal / MLD * TOpal
ExcNO3 = ExcTime * ( TNO3d - TNO3 )ExcSiOH4 = ExcTime * ( TSiOH4d - TSiOH4 )QNO3 = QNO3 + ExcNO3
QSiOH4 = QSiOH4 + ExcSiOH4QPON = QPON - ExpPON
QOpal = QOpal - ExpOpal...... Time Integration with Forward Scheme .....
TNO3 = TNO3 + dt * QNO3
TNH4 = TNH4 + dt * QNH4TPS = TPS + dt * QPS
TPL = TPL + dt * QPLTZS = TZS + dt * QZSTZL = TZL + dt * QZLTZP = TZP + dt * QZPTPON = TPON + dt * QPONTDON = TDON + dt * QDONTSiOH4 = TSiOH4 + dt * QSiOH4

154

TOpal = TOpal + dt * QOpal!
!

..... Vertical Migration of ZL .....
TZdwn = ND2TT(Iyr, ImonD, IdayD ,IhourD, IminD, IsecD )TZup = ND2TT(Iyr, ImonU, IdayU ,IhourU, IminU, IsecU )if ( (Tbefore .lt. TZdwn).and.(TTime .ge. TZdwn) ) then
TZLd = TZL
TZL = 0.0
write(*,*) '*** Down ***', CTimeend ifif ( (Tbefore .lt. TZup).and.(TTime .ge. TZup) ) then
TZL = SVRate * TZLd
write(*,*) '*** UP ***', CTimeend if!
call Herring(TTime, Tbefore, TZS, TZL, TZP, Temp)!
!
!

..... Monitor .....
if ( int(TTime/Tmon).ne. int(Tbefore/Tmon) ) thenwrite(*,'(A,13(",", F8.4))') CTime, Season
write(10,'(A,11(",", F8.4))') CTime, &TNO3/mcr, TNH4 /mcr, TPS /mcr, TPL /mcr, &TZS /mcr, TZL /mcr, TZP /mcr, TPON/mcr, &TDON/mcr, TSiOH4/mcr, TOpal/mcr
write(11,'(A,13(",", 1PE10.4))') CTime, Lint0, Temp, MLD, ExcTime*d2s
end ifend do!
close(10); close(11)
!
stop
end
!************************************************************************************
Subroutine Herring(TTime, Tbefore, TZS, TZL, TZP, Temp)!
implicit nonereal(8):: TTime, Tbefore, TZS, TZL, TZP, Temp
real(8),parameter :: d2s= 86400.0d0 ! day ---> sec
integer:: Iyr, Imon, Iday, Ihour, Imin, Iseccharacter(19) :: CAge ='0000/07/19 00:00:00' ! Date of Aging ( JJday = 200 )character(19) :: CTime
real(8):: TAge
integer:: iage = 0
integer, save :: IyrA, ImonA, IdayA ,IhourA, IminA, IsecA
integer:: JJdayreal(8):: ZooP1, ZooP2, ZooP3, tt1
real(8):: t1,t2,wtemp
real(8):: x(1) =0.2d0, xdot(1)real(8):: cd2tt, nd2tt
character(19) :: tt2cd
real(8):: vul(3), k(3)155integer(4):: id(365)real(8):: zop1(365), zop2(365), zop3(365)real(8):: v, a, u, resp
real(8) :: xk1,xk2,xk3,xk4,te1,te2,te3,te4,tt5,t5,t4,tt7,t7,t6, gcta,gctb,gctemp,gcmax
real(8) :: cnum,c1,c2,c3,con1,con2,con3,conreal(8) :: f,e, sda!
integer, save!
!

:: First = 1==================
if ( First .eq. 1 ) then; First = 0
TAge = CD2TT( CAge )call TT2ND(IyrA, ImonA, IdayA ,IhourA, IminA, IsecA ,TAge )open( 20, file='Herring.csv', form='FORMATTED' )!
!!!!
OPEN(UNIT=111,FILE='nemuro.txt',STATUS='unknown')!!!!
-----read in the 3 zoop groups from Nemuro output last 3 columns!!!!
do JJday=1,365!!!!
READ(111,999)id(JJday),zop1(JJday),zop2(JJday),zop3(JJday)!!!! 999
FORMAT(1x,i3,1x,3(e13.6,1x))!!!!
end doend if! ==================!
CTime = TT2CD(TTime)! present time (charactor form)
CALL TT2ND(Iyr, Imon, Iday ,Ihour, Imin, Isec ,TTime)
JJday = 1 + ( TTime - ND2TT(Iyr ,1,1,0,0,0) ) / d2s!
!------convert Nemuro zoop in uM N/L to g ww/m3!------ tt1 is conversion from uM N/liter to g ww/m3!------ 14 ug N/uM * 1.0e-6 g/ug * 1 g dw/0.07 g N dw * 1 g ww/0.2 g dw *!-----1.0e3 liters/m3!
tt1=14.0*1.0e-6*(1.0/0.07)*(1.0/0.2)*1.0e3
zoop1 = TZS*tt1 *1.0d6
zoop2 = TZL*tt1 *1.0d6
zoop3 = TZP*tt1 *1.0d6!!!! zoop1 = zop1(JJday) * tt1!!!! zoop2 = zop2(JJday) * tt1!!!! zoop3 = zop3(JJday) * tt1!
! ..... Temperature Seting .....
!
t1=float(jjday)
t2=12.75-10.99*cos(0.0172*t1)-6.63*sin(0.0172*t1)
wtemp=t2-5.0
IF(wtemp.le.1.0)wtemp=1.0!
!
!

write(*,*) TT2CD(cd2tt('0002/01/01 00:00:00')+200.0*86400.0)
stop..... Aging of Herring ......

156

TAge = ND2TT(Iyr, ImonA, IdayA ,IhourA, IminA, IsecA )if ( (Tbefore .lt. TAge).and.(TTime .ge. TAge) ) thenwrite(*,*) '*** Aging +1 of Herring ***', CTime
iage = iage + 1
end if!
!--- Herring weight state variable = x(1)!
!----- set vulnerabilities and k values for 3 zoop groups!
vul(1) = 1.0; vul(2) = 1.0; vul(3) = 1.0
k (1) = 0.3638; k (2) = 0.0364; k (3) = 0.3638!
! --- weight affect on respiration!
tt1 = 1.0 / x(1)t1 = 0.0033 * tt1**0.227! --- *********this is the new stuff from Ahhrenius for YOY only*********! --- The 5.258 puts resp (g oxygen/fish) into units of g zoop/g fish/day! --- [13560 joules/gram oxygen]/4.18 joules/cal = 3244 cal/gO2! --- [2580 joules/gram zoop]/4.18 joules/cal = 617 cal/g zoop
! --- so respiration in grams/oxygen/g fish/day is multiplied by 3244/617 = 5.258! --- to get food energy equivalents of a gram of oxygen respired!
IF (iage .eq. 0 )thenIF(wtemp.le.15.0)thenv = 5.76 * exp( 0.0238 * wtemp ) * x(1)**0.386
else
v = 8.6 * x(1)**0.386
endif
a=EXP((0.03-0.0*wtemp)*v)
resp=t1*EXP(0.0548*wtemp)*a*5.258! --- *********back to the old equations for respiration for age-1 and older********
else ! (iage .gt. 0)
IF (wtemp.le.9.0)then
u=3.9*x(1)**0.13*EXP(0.149*wtemp)
else
u=15.0*x(1)**0.13
endif
resp=t1*EXP(0.0548*wtemp)*EXP(0.03*u)*5.258
endif!C!C --- Thornton and Lessem temperature effect!C --- age dependent values!C --- *******Arrhenius for age-0 he changed te4 from 25 to 23 degrees******!C
if ( iage .eq. 0 ) thenxk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.01
te1 = 1.0; te2 = 15.0; te3 = 17.0; te4 = 23.0
else if ( iage .eq. 1 ) thenxk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.01

157

te1 = 1.0; te2 = 15.0; te3 = 17.0; te4 = 25.0
else if( iage .gt. 1 ) thenxk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.01
te1 = 1.0; te2 = 13.0; te3 = 15.0; te4 = 23.0
endif!
tt5 = ( 1.0 / ( te2 - te1 ) )t5 = tt5 * log( 0.98 * ( 1.0 - xk1 ) / ( 0.02 * xk1 ) )t4 = exp( t5 * ( wtemp - te1 ) )!
tt7 = 1.0 /( te4 - te3 )t7 = tt7 * log( 0.98 * ( 1.0 - xk4 ) / ( 0.02 * xk4 ) )t6 = exp( t7 * ( te4 - wtemp ) )!
gcta = ( xk1 * t4 ) / ( 1.0 + xk1 * ( t4 - 1.0 ) )gctb = xk4 * t6 / ( 1.0 + xk4 * ( t6 - 1.0 ) )gctemp= gcta * gctb
gcmax = 0.642 * tt1**0.256 * gctemp!
! --- multispecies functional response! --- usse either this or adjust little p!
cnum=zoop1 * vul(1)/k(1) + zoop2*vul(2)/k(2) +zoop3 * vul(3)/k(3)c1=gcmax*zoop1*vul(1)/k(1)c2=gcmax*zoop2*vul(2)/k(2)c3=gcmax*zoop3*vul(3)/k(3)con1=c1/(1.0+cnum)
con2=c2/(1.0+cnum)
con3=c3/(1.0+cnum)
con= con1+con2+con3!
!-----if using constant p rather than functional response, set p here
! --- to tune to observed size at age data
! con=0.425*gcmax!
! --- egestion!
f=0.16*con!
! --- excretion
e=0.1*(con-f)!
!
! --- Specific Dynamic Action!
!c------ *******Arrhenius age dependent SDA from 17.5% to 15% ****IF ( iage .eq. 0 ) thensda=0.15*(con-f)else
sda=0.175*(con-f)
end if158

!C!C --- use the ratio of calories/g of zoop (2580) to calories/g of fish (5533)!C!C --- bioenergetics differential equation!C
xdot(1)=(con-resp-f-e-sda)*x(1)*2580./5533.
!
IF(wtemp.le.1.0)xdot(1)=0.0!C!C --- Spawning section. Assume loose 20% of bosy weight/day!C
t1=float(jjday)!
if( mod(JJday,365) .ge. 152.0 .and. mod(JJday,365) .le. 156.0) then!
xdot(1)=(con-resp-f-e-sda-0.20)*x(1)*2580./5533.
!
write(*,*) '### Spawning ###'
!
endif
!
! if (iage .eq. 1 ) then!
write(*,*) JJday, wtemp, x(1), xdot(1)!
stop!
end if!
write(*,'(A,I4,3(1PE14.5))') Ctime, JJday, wtemp, x(1), xdot(1)!
! Time Integration!
x(1) = x(1) + 3600.0d0 /d2s * xdot(1)!
! ..... for Check .....
if ( int(TTime/d2s) .ne. int(Tbefore/d2s) ) then!!
write(*,'(A,I4,3(1PE14.5))') Ctime, JJday, wtemp, x(1), xdot(1)
!!
stop
!!
write(*,*) TZS, zop1(JJday), TZL,zop2(JJday), TZP,zop3(JJday)!!
write(*,*) TZP*1.0d6, zop3(JJday)
write(20,'(A,11(",", F12.4))') CTime, x(1), wtemp, gcmax
end if!
return!
stop
end!************************************************************************************!* Utilities for Date Control Writtien by Yasuhiro Yamanaka (galapen@ees.hokudai.ac.jp) *!************************************************************************************! exp. 1997/12/31 23:59:59 --> 6.223158719900000E+10! exp. 0000/01/01 00:00:00 --> 0.000000000000000E+00!************************************************************************************
real(8) function CD2TT( Cdate )!
integer:: Iyr, Imon, Iday , Ihour, Imin, Isecreal(8):: ND2TT
character(19) :: Cdate!

159

if ( len( Cdate ) .ne. 19 ) thenwrite(*,*) '### Length of date is no good ###'
stop
end ifread (Cdate( 1: 4),*) Iyr
read (Cdate( 6: 7),*) Imon
read (Cdate( 9:10),*) Iday
read (Cdate(12:13),*) Ihour
read (Cdate(15:16),*) Iminread (Cdate(18:19),*) Isec!
CD2TT = ND2TT(Iyr, Imon, Iday , Ihour, Imin, Isec)!
returnend function!***********************************************************************! exp. 6.223158719900000E+10 --> 1997/12/31 23:59:59!***********************************************************************
character(19) function TT2CD(tt)!
integer :: Iyr, Imon, Iday , Ihour, Imin, Isecreal(8) :: tt!
call TT2ND( Iyr, Imon, Iday, Ihour, Imin, Isec , tt )!
write(TT2CD,'(I4.4,5(A,I2.2))') Iyr, '/', Imon, '/', Iday, &' ', Ihour, ':', Imin, ':', Isec!
returnend function!***********************************************************************! exp. 1997,12,31,23,59,59 --> 6.223158719900000E+10!***********************************************************************
real(8) function ND2TT(Iyr, Imon, Iday, Ihour, Imin, Isec)!
integer :: IM2D(12,0:1) = &reshape( (/ 0,31,59,90,120,151,181,212,243,273,304,334, &0,31,60,91,121,152,182,213,244,274,305,335 /), (/12,2/) )integer :: Iyr, Imon, Iday, Ihour, Imin, Isec
integer :: Iy4, Iy1, Ileap, Im, Itt!
!
Iy4 = 1461 * ( Iyr / 4 )Iy1 = 365 * mod( Iyr, 4 )!
if ( mod( Iyr, 4 ) .ne. 0 ) thenIleap = 0else
Ileap = 1
end if
Im = IM2D( Imon, Ileap)160

!
Itt = Iy4 + Iy1 + Im + Iday - Ileap!
ND2TT = Ihour * 3600 + Imin * 60 + IsecND2TT = ND2TT + Itt * 86400.0D0!
returnend function!***********************************************************************! exp. 6.223158719900000E+10 --> 1997,12,31,23,59,59!***********************************************************************
subroutine TT2ND(&
Iyr , Imon , Iday , Ihour, Imin, Isec, & !O & I
tt )!
integer :: Iyr, Imon, Iday , Ihour, Imin, Isec
integer :: Itt, Iy, Iy4, Iyd, Iy1, Ileap, Imd, Im, Its
integer :: IM2D(12,0:1) = &reshape( (/ 0,31,59,90,120,151,181,212,243,273,304,334, &0,31,60,91,121,152,182,213,244,274,305,335 /), (/12,2/) )integer :: IY2D(4) = (/0,366,731,1096/)
real(8) :: tt, tt0, ND2TT!
!
! ..... ITT [day] .....
Itt = 1 + tt / 86400.0D0!
Iy4 = (Itt-1) / 1461
Iyd = Itt - Iy4 * 1461
do IY = 1, 4
if ( IY2D(Iy) + 1 .le. Iyd ) thenIy1 = Iy
end ifend do!
Iyr = Iy4 * 4 + Iy1 - 1
if ( mod(Iyr,4) .ne. 0 ) thenIleap = 0else
Ileap = 1
end if
IMD = IYD - IY2D(IY1)!
do IM = 1, 12
if ( IM2D(IM,ILEAP)+1 .le. IMD ) then
IMON = IM
end ifend doIDAY = IMD - IM2D(IMON,ILEAP)!
TT0 = ND2TT(IYR, IMON, IDAY ,0,0,0)161

ITS = nint( TT - TT0 )Ihour = ITS / 3600
Imin = ( ITS - Ihour * 3600 ) / 60
Isec = ITS - Ihour * 3600 - Imin * 60!
returnend subroutine!

162

Appendix 5 NEMURO.FISH (NEMURO FORTRAN code supplied by Yasuhiro Yamanaka)with the saury bioenergetic model (base case) (supplied by Bernard Megrey and Ken Rose andmodified by “team saury”). The saury model is linked to NEMURO in a one-way static link.
!************************************************************************************! NEMURO model Jun 13, 2002 written by Yasuhiro Yamanaka!
modified by Masahiko Fujii!
Shin-ichi Ito!************************************************************************************
program NEMURO
implicit none
! ..... Control for Time Ingration .....
character(19) :: Cstart = '0001/02/01 00:00:00' ! Starting datecharacter(19) :: Cend = '0003/02/01 00:00:00' ! Ending datecharacter(19) :: Cstep = '0000/00/00 01:00:00' ! Time step
character(19) :: Cmon = '0000/00/01 00:00:00' ! Monitor Interval
character(19) :: CTime
real(8):: dt, TTime, Tbefore, Season, Tmon
integer:: Iyr, Imon, Iday, Ihour, Imin, Isec
! ..... scale conversion .....
real(8),parameter :: d2s= 86400.0d0 ! day ---> secreal(8),parameter :: mcr= 1.0d-6! micro! ..... Prognostic Variables (with initial conditions) and Thier Source Term .....
real(8):: TPS = 0.1D-6, QPS ! Small Phytoplankton [molN/l]real(8):: TPL = 0.1D-6, QPL ! Large Phytoplanktonreal(8):: TZS = 0.1D-6, QZS ! Small Zooplanktonreal(8):: TZL = 0.1D-6, QZL ! Large Zooplanktonreal(8):: TZP = 0.1D-6, QZP ! Pradatory Zooplanktonreal(8):: TNO3 = 5.0D-6, QNO3 ! Nitratereal(8):: TNH4 = 0.1D-6, QNH4 ! Ammmonium
real(8):: TPON = 0.1D-6, QPON ! Particulate Organic Nitrogenreal(8):: TDON = 0.1D-6, QDON ! dissolved Organic Nitrogenreal(8):: TSiOH4 = 10.0D-6, QSiOH4 ! Silicate
real(8):: TOpal = 0.1D-7, QOpal ! Particulate Opal
! ..... Prognostic Variables (with initial conditions) and Thier Source Term .....
!
real(8):: THrr = 0.2D0, QHrrl ! Particulate Opal
! ...... Light Condition Parameters .......
real(8),parameter :: alpha1 = 4.0D-2 ! Light Dissipation coefficient of sea water[/m]real(8),parameter :: alpha2 = 4.0D4! PS+PL Selfshading coefficientS+PL[l/molN/m]real(8),parameter :: IoptS = 0.15D0 ! PS Optimum Light Intensity S[ly/min]real(8),parameter :: IoptL = 0.15D0 ! PL Optimum Light Intensity[ly/min]integer,parameter :: LLN= 10! Number of sublayer for calculating of Lfc
real(8):: LfcS! Light factor for PS
real(8):: LfcL! Light factor for PLreal(8):: kappa, Lint, dLint, LfcUS, LfcUL, LfcDS, LfcDL
integer:: L! ...... biological Parameters ......
real(8),parameter :: VmaxS = 0.4D0/d2s ! PS Maximum Photosynthetic rate @0degC [/s]real(8),parameter :: KNO3S = 1.0D-6 ! PS Half satuation constant for Nitrate [molN/l]real(8),parameter :: KNH4S = 0.1D-6 ! PS Half satuation constant for Ammonium [molN/l]163

real(8),parameter :: PusaiS = 1.5D6! PS Ammonium Inhibition Coefficient[l/molN]real(8),parameter :: KGppS = 6.93D-2 ! PS Temp. Coeff. for Photosynthetic Rate [/degC]real(8),parameter :: MorPS0 = 5.85D4/d2s ! PS Mortality Rate @0degC[/s]real(8),parameter :: KMorPS = 6.93D-2 ! PS Temp. Coeff. for Mortality[/degC]real(8),parameter :: ResPS0 = 0.03D0/d2s ! PS Respiration Rate at @0degC[/s]real(8),parameter :: KResPS = 0.0519D0 ! PS Temp. Coeff. for Respiration[/degC]real(8),parameter :: GammaS = 0.135D0 ! PS Ratio of Extracell. Excret. to Photo.[(nodim)]real(8),parameter :: VmaxL = 0.8D0/d2s ! PL Maximum Photosynthetic rate @0degC [/s]real(8),parameter :: KNO3L = 3.00D-6 ! PL Half satuation constant for Nitrate [molN/l]real(8),parameter :: KNH4L = 0.30D-6 ! PL Half satuation constant for Ammonium [molN/l]real(8),parameter :: KSiL = 6.00D-6 ! PL Half satuation constant for Silicate [molSi/l]real(8),parameter :: PusaiL = 1.50D6 ! PL Ammonium Inhibition Coefficient[l/molN]real(8),parameter :: KGppL = 6.93D-2 ! PL Temp. Coeff. for Photosynthetic Rate [/degC]real(8),parameter :: MorPL0 = 2.90D4/d2s ! PL Mortality Rate @0degC[/s]real(8),parameter :: KMorPL = 6.93D-2 ! PL Temp. Coeff. for Mortality[/degC]real(8),parameter :: ResPL0 = 0.03D0/d2s ! PL Respiration Rate at @0degC[/s]real(8),parameter :: KResPL = 0.0519D0 ! PL Temp. Coeff. for Respiration[/degC]real(8),parameter :: GammaL = 0.135D0 ! PL Ratio of Extracell. Excret. to Photo.[(nodim)]real(8),parameter :: GRmaxS = 0.40D0/d2s ! ZS Maximum Rate of Grazing PS @0degC [/s]real(8),parameter :: KGraS = 6.93D-2 ! ZS Temp. Coeff. for Grazing[/degC]real(8),parameter :: LamS = 1.40D6 ! ZS Ivlev constant[l/molN]real(8),parameter :: PS2ZSstar= 0.043D-6 ! ZS Threshold Value for Grazing PS[molN/l]real(8),parameter :: AlphaZS = 0.70D0 ! ZS Assimilation Efficiency[(nodim)]real(8),parameter :: BetaZS = 0.30D0 ! ZS Growth Efficiency[(nodim)]real(8),parameter :: MorZS0 = 5.85D4/d2s ! ZS Mortality Rate @0degC[/s]real(8),parameter :: KMorZS = 0.0693D0 ! ZS Temp. Coeff. for Mortality[/degC]real(8),parameter :: GRmaxLps = 0.10D0/d2s ! ZL Maximum Rate of Grazing PS @0degC [/s]real(8),parameter :: GRmaxLpl = 0.40D0/d2s ! ZL Maximum Rate of Grazing PL @0degC [/s]real(8),parameter :: GRmaxLzs = 0.40D0/d2s ! ZL Maximum Rate of Grazing ZS @0degC [/s]real(8),parameter :: KGraL = 6.93D-2 ! ZL Temp. Coeff. for Grazing[/degC]real(8),parameter :: LamL = 1.4000D6 ! ZL Ivlev constant[l/molN]real(8),parameter :: PS2ZLstar= 4.00D-8 ! ZL Threshold Value for Grazing PS[molN/l]real(8),parameter :: PL2ZLstar= 4.00D-8 ! ZL Threshold Value for Grazing PL[molN/l]real(8),parameter :: ZS2ZLstar= 4.00D-8 ! ZL Threshold Value for Grazing ZS[molN/l]real(8),parameter :: AlphaZL = 0.70D0 ! ZL Assimilation Efficiency[(nodim)]real(8),parameter :: BetaZL = 0.30D0 ! ZL Growth Efficiency[(nodim)]real(8),parameter :: MorZL0 = 5.85D4/d2s ! ZL Mortality Rate @0degC[/s]real(8),parameter :: KMorZL = 0.0693D0 ! ZL Temp. Coeff. for Mortality[/degC]real(8),parameter :: GRmaxPpl = 0.20D0/d2s ! ZP Maximum rate of grazing PL @0degC [/s]real(8),parameter :: GRmaxPzs = 0.20D0/d2s ! ZP Maximum rate of grazing ZS @0degC [/s]real(8),parameter :: GRmaxPzl = 0.20D0/d2s ! ZP Maximum rate of grazing ZL @0degC [/s]real(8),parameter :: KGraP = 6.93D-2 ! ZP Temp. Coeff. for grazing[/degC]real(8),parameter :: LamP = 1.4000D6 ! ZP Ivlev constant[l/molN]real(8),parameter :: PL2ZPstar= 4.00D-8 ! ZP Threshold Value for Grazing PL[molN/l]real(8),parameter :: ZS2ZPstar= 4.00D-8 ! ZP Threshold Value for Grazing ZS[molN/l]real(8),parameter :: ZL2ZPstar= 4.00D-8 ! ZP Threshold Value for Grazing ZL[molN/l]real(8),parameter :: PusaiPL = 4.605D6 ! ZP Preference Coeff. for PL[l/molN]real(8),parameter :: PusaiZS = 3.010D6 ! ZP Preference Coeff. for ZS[l/molN]real(8),parameter :: AlphaZP = 0.70D0 ! ZP Assimilation Efficiency[(nodim)]real(8),parameter :: BetaZP = 0.30D0 ! ZP Growth Efficiency[(nodim)]real(8),parameter :: MorZP0 = 5.85D4/d2s ! ZP Mortality Rate @0degC
[/s]

164

!

!

real(8),parameter :: KMorZP = 0.0693D0 ! ZP Temp. Coeff. for Mortality[/degC]real(8),parameter :: Nit0 = 0.03D0/d2s ! NH4 Nitrification Rate @0degC[/s]real(8),parameter :: KNit = 0.0693D0 ! NH4 Temp. coefficient for Nitrification [/degC]real(8),parameter :: VP2N0 = 0.10D0/d2s ! PON Decomp. Rate to Ammonium @0degC [/s]real(8),parameter :: KP2N = 6.93D-2 ! PON Temp. Coeff. for Decomp. to Ammon. [/degC]real(8),parameter :: VP2D0 = 0.10D0/d2s ! PON Decomp. Rate to DON @0degC[/s]real(8),parameter :: KP2D = 6.93D-2 ! PON Temp. Coeff. for Decomp. to DON [/degC]real(8),parameter :: VD2N0 = 0.20D0/d2s ! DON Decomp. Rate to Ammonium @0degC [/s]real(8),parameter :: KD2N = 6.93D-2 ! DON Temp. Coeff. for Decomp. to Ammon. [/degC]real(8),parameter :: VO2S0 = 0.10D0/d2s ! Opal Decomp. Rate to Silicate @0degC [/s]real(8),parameter :: KO2S = 6.93D-2 ! Opal Temp. Coeff. for Decomp.to Silicate[/degC]real(8),parameter :: RSiN = 2.0D0!Si/N ratio[molSi/molN]
real(8),parameter :: RCN= 106.0D0/16.0D0 !C/N ratio[molC/molN]..... bottom boundary Condition .....
real(8),parameter :: setVPON = 40.0D0/d2s ! Settling velocity of PON[m/s]real(8),parameter :: setVOpal = 40.0D0/d2s ! Settling velocity of Opal[m/s]real(8),parameter :: TNO3d = 25.0d-6 ! Nitrate Concentraion in the Deep Layer [molN/l]real(8),parameter :: TSiOH4d = 35.0d-6 ! Silicate Concentraion in the Deep Layer [molSi/l]...... Paramters of ZL Vertical Migration .....
character(19) :: CZup ='0000/04/01 00:00:00' ! Date Coming up to the Euphotic Layer
character(19) :: CZdwn='0000/09/01 00:00:00' ! Date Returning to the Deep Layer
real(8):: TZup, TZdwn, TZLd, SVRate=0.2D0
integer:: IyrU, ImonU, IdayU ,IhourU, IminU, IsecU
integer:: IyrD, ImonD, IdayD ,IhourD, IminD, IsecD!
real(8)real(8)real(8)real(8)real(8)real(8)real(8)real(8)real(8)integer
!
!

!

!
!
!

:: GppPSn, GppNPSn, GppAPSn, RnewS, ResPSn, MorPSn, ExcPSn:: GppPLn, GppNPLn, GppAPLn, RnewL, ResPLn, MorPLn, ExcPLn:: GppPLsi, GppSiPLsi,
ResPLsi, MorPLsi, ExcPLsi:: GraPS2ZSn, GraPS2ZLn, GraPL2ZLn, GraPL2ZLsi, GraZS2ZLn:: GraPL2ZPn, GraPL2ZPsi, GraZS2ZPn, GraZL2ZPn:: EgeZSn, MorZSn, ExcZSn, EgeZLn, EgeZLsi, MorZLn, ExcZLn::
EgeZPn, EgeZPsi, MorZPn, ExcZPn:: DecP2N, DecP2D, DecD2N, DecO2S, Nit:: ExpPON, ExpOpal,ExcNO3, ExcSiOH4:: lt=0 , nt..... Environmental Condition .....
real(8):: Temp! Temperature [degC]real(8):: Lint0! Light Intencity at sea surface [ly/min]real(8):: MLD = 30.0d0! Mixed Layer Depth [m]real(8):: ExcTime = 1.0d0 / (100.0d0*d2s) ! Exch. Coeff. between Sur-Deep [/s]..... statement function & def. type of functions .....
real(8):: cd2tt, nd2tt
character(19) :: tt2cd
real(8):: Td, GraF, Mich, a, b, c
Td (a,b) = a * exp(b*Temp)
GraF(a,b,c) = MAX( 0.0D0, 1.0 - exp(a * (b - c)))Mich(a,b) = b / ( a + b )***** Initial Setting *****..... for time control .....

165

!

!

!
!
!

!
!

!
!

TTime = cd2tt(Cstart)! Starting Date
CTime = TT2CD(TTime)! present time (charactor form)
dt = cd2tt(Cstep) - cd2tt('0000/00/00 00:00:00') ! Time Step (real8 form)
Tmon = cd2tt(Cmon) - cd2tt('0000/00/00 00:00:00') ! Monitor Interval (real8 form)
nt = NINT( ( cd2tt(Cend) - cd2tt(Cstart) ) / dt ) ! Total Time Steps..... for Vertical Migration ......
TZup = CD2TT( CZup )TZdwn = CD2TT( CZdwn )call TT2ND(IyrU, ImonU, IdayU ,IhourU, IminU, IsecU ,TZup )call TT2ND(IyrD, ImonD, IdayD ,IhourD, IminD, IsecD ,TZdwn)
TZLd = TZL ! ZL living in the deep layer at the initial condition
TZL = 0.0..... File Open for monitoring output .....
open( 10, file='Results.csv', form='FORMATTED' )write(10,'(A,13(",", A))') 'Time(day)', &'NO3' , 'NH4' , 'PS' , 'PL',&
'ZS' , 'ZL'
, 'ZP' , 'PON' , &'DON' , 'SiOH4' , 'Opal' , 'TotalN' , 'TotalSi'
write(10,'(A,11(",", F8.4))') CTime, &TNO3/mcr, TNH4 /mcr, TPS /mcr, TPL /mcr, &TZS /mcr, TZL /mcr, TZP /mcr, TPON/mcr, &TDON/mcr, TSiOH4/mcr, TOpal/mcr
open( 11, file='Forcing.csv', form='FORMATTED' )write(11,'(A,13(",", A))') 'Time(day)','Lint0', 'TMP', 'MLD', 'ExcTime'write(11,'(A,13(",", 1PE10.4))') CTime, Lint0, Temp, MLD, ExcTime*d2s****** Main Loop *****
do lt = 1, nt..... time control (Season : 0 to 1, parcentage in a year).....
Tbefore = TTime! one step before present time
TTime = TTime + dt
! present time (real8 form)
CTime = TT2CD(TTime)! present time (charactor form)
CALL TT2ND(Iyr, Imon, Iday ,Ihour, Imin, Isec ,TTime)
Season = ( TTime- ND2TT(Iyr ,1,1,0,0,0) )/ &( ND2TT(Iyr+1,1,1,0,0,0) - ND2TT(Iyr ,1,1,0,0,0) )..... Example of Boundray condition .....
Lint0 = 0.1d0 * ( 1.0D0 + 0.3d0 * cos( 2.0d0*3.1415926536d0*(Season - 0.50D0) ) )Temp =
6.0D0 + 4.0d0 * cos( 2.0d0*3.1415926536d0*(Season - 0.65D0) )if (Temp .lt. 4.0 ) then
MLD = MLD + dt * (150.0d0 - MLD ) / ( 100.0d0 * d2s )ExcTime = ExcTime + dt * ( 1.0d0/( 40.0d0*d2s) - ExcTime ) / (100.0d0*d2s)else
MLD = MLD + dt * ( 30.0d0 - MLD ) / ( 5.0d0 * d2s )ExcTime = ExcTime + dt * ( 1.0d0/(100.0d0*d2s) - ExcTime ) / ( 5.0d0*d2s)
end if..... Light Factors (LfcS, LfcL).....
Lint = Lint0
LfcDS = Lint/IoptS * exp(1.0D0 - Lint/IoptS)
LfcDL = Lint/IoptL * exp(1.0D0 - Lint/IoptL)

166

!

!

!

!
!
!

LfcS = 0.0D0
LfcL = 0.0D0
Kappa = alpha1 + alpha2 * ( TPS + TPL )dLint = exp( -Kappa * (MLD/LLN) )do L = 1, LLN
LfcUS = LfcDS
LfcUL = LfcDL
Lint = Lint
* dLint
LfcDS = Lint/IoptS * exp( 1.0D0 - Lint/IoptS )LfcDL = Lint/IoptL * exp( 1.0D0 - Lint/IoptL )LfcS = LfcS + ( LfcUS + LfcDS ) * 0.5D0 / LLN
LfcL = LfcL + ( LfcUL + LfcDL ) * 0.5D0 / LLN
end do..... Photosynthesis of PS .....
GppNPSn = Mich( KNO3S, TNO3 ) * exp( - PusaiS * TNH4 )GppAPSn = Mich( KNH4S, TNH4 )GppPSn = Td(VmaxS, KGppS) * LfcS * TPS * ( GppNPSn + GppAPSn )RnewS = GppNPSn / ( GppNPSn + GppAPSn )ResPSn = Td( ResPS0, KResPS ) * TPSMorPSn = Td( MorPS0, KMorPS ) * TPS * TPS
ExcPSn = GammaS* GppPSn..... Photosynthesis of PL .....
GppNPLn = Mich( KNO3L, TNO3 ) * exp( - PusaiL * TNH4 )GppAPLn = Mich( KNH4L, TNH4 )GppSiPLsi = Mich( KSiL , TSiOH4 )GppPLn = Td(VmaxL, KGppL) * LfcL * TPL * min( ( GppNPLn + GppAPLn ), GppSiPLsi )RnewL = GppNPLn / ( GppNPLn + GppAPLn )ResPLn = Td( ResPL0, KResPL ) * TPLMorPLn = Td( MorPL0, KMorPL ) * TPL * TPL
ExcPLn = GammaL* GppPLn..... Grazing PS, PL, ZS, ZL --> ZS, ZL, ZP .....
GraPS2ZSn = Td(GRmaxS, KGraS) *GraF(LamS,PS2ZSstar,TPS) *TZSGraPS2ZLn = Td(GRmaxLps,KGraL) *GraF(LamL,PS2ZLstar,TPS) *TZLGraPL2ZLn = Td(GRmaxLpl,KGraL) *GraF(LamL,PL2ZLstar,TPL) *TZLGraZS2ZLn = Td(GRmaxLzs,KGraL) *GraF(LamL,ZS2ZLstar,TZS) *TZLGraPL2ZPn = Td(GRmaxPpl,KGraP) *GraF(LamP,PL2ZPstar,TPL) *TZP * exp( -PusaiPL *(TZL+ TZS))GraZS2ZPn = Td(GRmaxPzs,KGraP) *GraF(LamP,ZS2ZPstar,TZS) *TZP * exp( -PusaiZS * TZL )GraZL2ZPn = Td(GRmaxPzl,KGraP) *GraF(LamP,ZL2ZPstar,TZL) *TZP..... Mortality, Excration, Egestion for Zooplanktons..... Commented out after Saito-san Meeting at 19 Jun, 2000 .....
BetaZS = 0.3 ** ( 1.0 + Mich( TPL, TPS ) )ExcZSn = (AlphaZS- BetaZS) * GraPS2ZSn
EgeZSn = (1.0 - AlphaZS) * GraPS2ZSn
MorZSn = Td( MorZS0, KMorZS ) * TZS * TZSExcZLn = (AlphaZL- BetaZL) * (GraPS2ZLn+GraPL2ZLn+GraZS2ZLn)
EgeZLn = (1.0 - AlphaZL) * (GraPS2ZLn+GraPL2ZLn+GraZS2ZLn)
MorZLn = Td( MorZL0, KMorZL ) * TZL * TZLExcZPn = (AlphaZP- BetaZP) * (GraPL2ZPn+GraZS2ZPn+GraZL2ZPn)
EgeZPn = (1.0 - AlphaZP) * (GraPL2ZPn+GraZS2ZPn+GraZL2ZPn)

167

!

!

!
!

!
!

!
!

MorZPn = Td( MorZP0, KMorZP ) * TZP * TZP..... Decomposition PON, DON, Opal ---> NH4, DON, SiOH4 .....
DecP2N = Td( VP2N0 , KP2N ) * TPON
DecP2D = Td( VP2D0 , KP2D ) * TPON
DecD2N = Td( VD2N0 , KD2N ) * TDON
DecO2S = Td( VO2S0 , KO2S ) * TOpal
Nit= Td( Nit0 , KNit ) * TNH4..... silica fluxes .....
GppPLsi = GppPLn * RSiNResPLsi = ResPLn * RSiNMorPLsi = MorPLn * RSiNExcPLsi = ExcPLn * RSiNGraPL2ZLsi = GraPL2ZLn * RSiNGraPL2ZPsi = GraPL2ZPn * RSiNEgeZLsi = GraPL2ZLsi
EgeZPsi = GraPL2ZPsi..... Tendency Terms for biological processes .....
QNO3 = -( GppPSn - ResPSn ) * RnewS &-( GppPLn - ResPLn ) * RnewL + Nit
QNH4 = -( GppPSn - ResPSn ) * (1.0 - RnewS) &-( GppPLn - ResPLn ) * (1.0 - RnewL) &- Nit + DecP2N + DecD2N + ExcZSn + ExcZLn + ExcZPn
QPS = GppPSn - ResPSn - MorPSn - ExcPSn - GraPS2ZSn - GraPS2ZLn
QPL = GppPLn - ResPLn - MorPLn - ExcPLn - GraPL2ZLn - GraPL2ZPn
QZS = GraPS2ZSn - GraZS2ZLn - MorZSn - ExcZSn - EgeZSn - GraZS2ZPn
QZL = GraPL2ZLn + GraZS2ZLn - MorZLn - ExcZLn - EgeZLn + GraPS2ZLn - GraZL2ZPn
QZP = GraPL2ZPn + GraZS2ZPn - MorZPn - ExcZPn - EgeZPn + GraZL2ZPn
QPON = MorPSn + MorPLn + MorZSn + MorZLn + MorZPn &+ EgeZPn + EgeZSn + EgeZLn - DecP2N - DecP2D
QDON = ExcPSn + ExcPLn + DecP2D - DecD2N
QSiOH4 =-GppPLsi + ResPLsi + ExcPLsi + DecO2S
QOpal = MorPLsi + EgeZLsi + EgeZPsi - DecO2S..... Exchange Fluxes between the Surface and Deep Layers .....
ExpPON = setVPON / MLD * TPON
ExpOpal = setVOpal / MLD * TOpal
ExcNO3 = ExcTime * ( TNO3d - TNO3 )ExcSiOH4 = ExcTime * ( TSiOH4d - TSiOH4 )QNO3 = QNO3 + ExcNO3
QSiOH4 = QSiOH4 + ExcSiOH4QPON = QPON - ExpPON
QOpal = QOpal - ExpOpal...... Time Integration with Forward Scheme .....
TNO3 = TNO3 + dt * QNO3
TNH4 = TNH4 + dt * QNH4TPS = TPS + dt * QPS
TPL = TPL + dt * QPLTZS = TZS + dt * QZSTZL = TZL + dt * QZL168

TZP = TZP + dt * QZPTPON = TPON + dt * QPONTDON = TDON + dt * QDONTSiOH4 = TSiOH4 + dt * QSiOH4
TOpal = TOpal + dt * QOpal!
!

..... Vertical Migration of ZL .....
TZdwn = ND2TT(Iyr, ImonD, IdayD ,IhourD, IminD, IsecD )TZup = ND2TT(Iyr, ImonU, IdayU ,IhourU, IminU, IsecU )if ( (Tbefore .lt. TZdwn).and.(TTime .ge. TZdwn) ) then
TZLd = TZL
TZL = 0.0
write(*,*) '*** Down ***', CTimeend ifif ( (Tbefore .lt. TZup).and.(TTime .ge. TZup) ) then
TZL = SVRate * TZLd
write(*,*) '*** UP ***', CTimeend if!
call Saury(TTime, Tbefore, TZS, TZL, TZP, Temp)!
!
!

..... Monitor .....
if ( int(TTime/Tmon).ne. int(Tbefore/Tmon) ) thenwrite(*,'(A,13(",", F8.4))') CTime, Season
write(10,'(A,11(",", F8.4))') CTime, &TNO3/mcr, TNH4 /mcr, TPS /mcr, TPL /mcr, &TZS /mcr, TZL /mcr, TZP /mcr, TPON/mcr, &TDON/mcr, TSiOH4/mcr, TOpal/mcr
write(11,'(A,13(",", 1PE10.4))') CTime, Lint0, Temp, MLD, ExcTime*d2s
end ifend do!
close(10); close(11)
!
stop
end
!************************************************************************************
Subroutine Saury(TTime, Tbefore, TZS, TZL, TZP, Temp)!
implicit nonereal(8):: TTime, Tbefore, TZS, TZL, TZP, Temp
real(8),parameter :: d2s= 86400.0d0 ! day ---> sec
integer:: Iyr, Imon, Iday, Ihour, Imin, Iseccharacter(19) :: CAge ='0000/03/01 00:00:00' ! Date of Aging ( JJday = 200 )character(19) :: CTime
character(19) :: CAge2 ='0000/07/01 00:00:00' ! Date of Aging ( JJday = 200 )character(19) :: CTime2real(8):: TAgereal(8):: TAge2
integer:: iage = 0
integer, save :: IyrA, ImonA, IdayA ,IhourA, IminA, IsecA

169

integer, save :: IyrB, ImonB, IdayB ,IhourB, IminB, IsecB
integer:: JJdayreal(8):: ZooP1, ZooP2, ZooP3, tt1
real(8):: t1,t2,wtemp
real(8):: x(1) =0.2d0, xdot(1)real(8):: cd2tt, nd2tt
character(19) :: tt2cd
real(8):: vul(3), k(3)integer(4):: id(365)real(8):: zop1(365), zop2(365), zop3(365)real(8):: v, a, u, resp
real(8) :: xk1,xk2,xk3,xk4,te1,te2,te3,te4,tt5,t5,t4,tt7,t7,t6, gcta,gctb,gctemp,gcmax
real(8) :: cnum,c1,c2,c3,con1,con2,con3,conreal(8) :: f,e, sda
real(8) :: phalf!
integer, save:: First = 1!
PHalf=0.100!

==================
if ( First .eq. 1 ) then; First = 0
TAge = CD2TT( CAge )call TT2ND(IyrA, ImonA, IdayA ,IhourA, IminA, IsecA ,TAge )TAge2 = CD2TT( CAge2 )call TT2ND(IyrB, ImonB, IdayB ,IhourB, IminB, IsecB ,TAge2 )open( 20, file='saury.csv', form='FORMATTED' )!
!!!!
OPEN(UNIT=111,FILE='nemuro.txt',STATUS='unknown')!!!!
-----read in the 3 zoop groups from Nemuro output last 3 columns!!!!
do JJday=1,365!!!!
READ(111,999)id(JJday),zop1(JJday),zop2(JJday),zop3(JJday)!!!! 999
FORMAT(1x,i3,1x,3(e13.6,1x))!!!!
end doend if! ==================!
CTime = TT2CD(TTime)! present time (charactor form)
CALL TT2ND(Iyr, Imon, Iday ,Ihour, Imin, Isec ,TTime)
JJday = 1 + ( TTime - ND2TT(Iyr ,1,1,0,0,0) ) / d2s!
!------convert Nemuro zoop in uM N/L to g ww/m3!------ tt1 is conversion from uM N/liter to g ww/m3!------ 14 ug N/uM * 1.0e-6 g/ug * 1 g dw/0.07 g N dw * 1 g ww/0.2 g dw!-----1.0e3 liters/m3!
tt1=14.0*1.0e-6*(1.0/0.07)*(1.0/0.2)*1.0e3
zoop1 = TZS*tt1 *1.0d6
zoop2 = TZL*tt1 *1.0d6

170

zoop3 = TZP*tt1 *1.0d6!!!! zoop1 = zop1(JJday) * tt1!!!! zoop2 = zop2(JJday) * tt1!!!! zoop3 = zop3(JJday) * tt1!
! ..... Temperature Seting .....
!
t1=float(jjday)
t2=12.75-10.99*cos(0.0172*t1)-6.63*sin(0.0172*t1)
wtemp=t2-5.0
IF(wtemp.le.1.0)wtemp=1.0!
!
!

write(*,*) TT2CD(cd2tt('0002/01/01 00:00:00')+200.0*86400.0)
stop..... Aging of saury ......
TAge = ND2TT(Iyr, ImonA, IdayA ,IhourA, IminA, IsecA )if ( (Tbefore .lt. TAge).and.(TTime .ge. TAge) .and. (iage.eq.0)) thenwrite(*,*) '*** Aging +1 of saury ***', CTime
iage = iage + 1
end ifTAge2 = ND2TT(Iyr, ImonB, IdayB ,IhourB, IminB, IsecB )if ( (Tbefore .lt. TAge2).and.(TTime .ge. TAge2) .and. (iage.eq.1) ) thenwrite(*,*) '*** Aging +1 of saury ***', CTime
iage = iage + 1
end if!
!--- saury weight state variable = x(1)!
!----- set vulnerabilities and k values for 3 zoop groups!
if ( iage .eq. 0 ) thenvul(1) = 1.0; vul(2) = 0.0; vul(3) = 0.0
k (1) = phalf; k (2) = phalf; k (3) = phalf
else if ( iage .eq. 1 ) thenvul(1) = 1.0; vul(2) = 1.0; vul(3) = 0.0
k (1) = phalf; k (2) = phalf; k (3) = phalfelsevul(1) = 0.0; vul(2) = 1.0; vul(3) = 1.0
k (1) = phalf; k (2) = phalf; k (3) = phalfendif!
! --- weight affect on respiration!
tt1 = 1.0 / x(1)t1 = 0.0033 * tt1**0.227! --- *********this is the new stuff from Ahhrenius for YOY only*********! --- The 5.258 puts resp (g oxygen/fish) into units of g zoop/g fish/day! --- [13560 joules/gram oxygen]/4.18 joules/cal = 3244 cal/gO2! --- [2580 joules/gram zoop]/4.18 joules/cal = 617 cal/g zoop
! --- so respiration in grams/oxygen/g fish/day is multiplied by 3244/617 = 5.258! --- to get food energy equivalents of a gram of oxygen respired

171

!
!ccc
IF (iage .eq. 0 )then!ccc
IF(wtemp.le.15.0)then!ccc
v = 5.76 * exp( 0.0238 * wtemp ) * x(1)**0.386
!ccc
else
!ccc
v = 8.6 * x(1)**0.386
!ccc
endif
!ccc
a=EXP((0.03-0.0*wtemp)*v)!ccc
resp=t1*EXP(0.0548*wtemp)*a*5.258! --- *********back to the old equations for respiration for age-1 and older********!ccc
else ! (iage .gt. 0)!c
IF (wtemp.le.9.0)then
IF (wtemp.le.12.0)then!C
u=3.9*x(1)**0.13*EXP(0.149*wtemp)
u=2.0*x(1)**0.33*EXP(0.149*wtemp)else!c
u=15.0*x(1)**0.13
u=11.7*x(1)**0.33
endif
resp=t1*EXP(0.0548*wtemp)*EXP(0.03*u)*5.258
!ccc
endif
!C
!C --- Thornton and Lessem temperature effect!C --- age dependent values!C --- *******Arrhenius for age-0 he changed te4 from 25 to 23 degrees******!C
if ( iage .eq. 0 ) then!cxk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.01!c
te1 = 1.0; te2 = 15.0; te3 = 17.0; te4 = 23.0
xk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.5
te1 = 5.0; te2 = 20.0; te3 = 26.0; te4 = 30.0
else if ( iage .eq. 1 ) then!cxk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.01!c
te1 = 1.0; te2 = 15.0; te3 = 17.0; te4 = 25.0
xk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.5
te1 = 5.0; te2 = 16.0; te3 = 20.0; te4 = 30.0
else if( iage .gt. 1 ) then!cxk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.01!c
te1 = 1.0; te2 = 13.0; te3 = 15.0; te4 = 23.0
xk1 = 0.1; xk2 = 0.98; xk3 = 0.98; xk4 = 0.5
te1 = 5.0; te2 = 16.0; te3 = 20.0; te4 = 30.0
endif!
tt5 = ( 1.0 / ( te2 - te1 ) )t5 = tt5 * log( xk2 * ( 1.0 - xk1 ) / ( (1.0-xk2) * xk1 ) )t4 = exp( t5 * ( wtemp - te1 ) )!
tt7 = 1.0 /( te4 - te3 )t7 = tt7 * log( xk3 * ( 1.0 - xk4 ) / ( (1.0-xk3) * xk4 ) )t6 = exp( t7 * ( te4 - wtemp ) )!

172

gcta = ( xk1 * t4 ) / ( 1.0 + xk1 * ( t4 - 1.0 ) )gctb = xk4 * t6 / ( 1.0 + xk4 * ( t6 - 1.0 ) )gctemp= gcta * gctb!c
gcmax = 0.642 * tt1**0.256 * gctemp
gcmax = 0.6 * tt1**0.256 * gctemp!
! --- multispecies functional response! --- usse either this or adjust little p!
cnum=zoop1 * vul(1)/k(1) + zoop2*vul(2)/k(2) +zoop3 * vul(3)/k(3)c1=gcmax*zoop1*vul(1)/k(1)c2=gcmax*zoop2*vul(2)/k(2)c3=gcmax*zoop3*vul(3)/k(3)con1=c1/(1.0+cnum)
con2=c2/(1.0+cnum)
con3=c3/(1.0+cnum)
con= con1+con2+con3!
!-----if using constant p rather than functional response, set p here
! --- to tune to observed size at age data
! con=0.425*gcmax!
! --- egestion!
f=0.16*con!
! --- excretion
e=0.1*(con-f)!
!
! --- Specific Dynamic Action!
!c------ *******Arrhenius age dependent SDA from 17.5% to 15% ****IF ( iage .eq. 0 ) thensda=0.15*(con-f)else
sda=0.175*(con-f)
end if!C!C --- use the ratio of calories/g of zoop (2580) to calories/g of fish (5533)!C!C --- bioenergetics differential equation!C
xdot(1)=(con-resp-f-e-sda)*x(1)*2580./5533.
!
IF(wtemp.le.1.0)xdot(1)=0.0!C!C --- Spawning section. Assume loose 20% of bosy weight/day!C
t1=float(jjday)!
if( mod(JJday,365) .ge. 152.0 .and. mod(JJday,365) .le. 156.0) then!
xdot(1)=(con-resp-f-e-sda-0.20)*x(1)*2580./5533.

173

!
!
!
!
!
!
!
!
!
!
!

write(*,*) '### Spawning ###'endif
if (iage .eq. 1 ) thenwrite(*,*) JJday, wtemp, x(1), xdot(1)stop
end ifwrite(*,'(A,I4,3(1PE14.5))') Ctime, JJday, wtemp, x(1), xdot(1)Time Integrationx(1) = x(1) + 3600.0d0 /d2s * xdot(1)!
!

..... for Check .....
if ( int(TTime/d2s) .ne. int(Tbefore/d2s) ) then!!
write(*,'(A,I4,3(1PE14.5))') Ctime, JJday, wtemp, x(1), xdot(1)
!!
stop
!!
write(*,*) TZS, zop1(JJday), TZL,zop2(JJday), TZP,zop3(JJday)!!
write(*,*) TZP*1.0d6, zop3(JJday)
write(20,'(A,11(",", F12.4))') CTime, x(1), wtemp, gcmax
end if!
return!
stop
end!************************************************************************************!* Utilities for Date Control Writtien by Yasuhiro Yamanaka (galapen@ees.hokudai.ac.jp) *!************************************************************************************! exp. 1997/12/31 23:59:59 --> 6.223158719900000E+10! exp. 0000/01/01 00:00:00 --> 0.000000000000000E+00!************************************************************************************
real(8) function CD2TT( Cdate )!
integer:: Iyr, Imon, Iday , Ihour, Imin, Isecreal(8):: ND2TT
character(19) :: Cdate!
if ( len( Cdate ) .ne. 19 ) thenwrite(*,*) '### Length of date is no good ###'
stop
end ifread (Cdate( 1: 4),*) Iyr
read (Cdate( 6: 7),*) Imon
read (Cdate( 9:10),*) Iday
read (Cdate(12:13),*) Ihour
read (Cdate(15:16),*) Iminread (Cdate(18:19),*) Isec!
CD2TT = ND2TT(Iyr, Imon, Iday , Ihour, Imin, Isec)!

174

returnend function!***********************************************************************! exp. 6.223158719900000E+10 --> 1997/12/31 23:59:59!***********************************************************************
character(19) function TT2CD(tt)!
integer :: Iyr, Imon, Iday , Ihour, Imin, Isecreal(8) :: tt!
call TT2ND( Iyr, Imon, Iday, Ihour, Imin, Isec , tt )!
write(TT2CD,'(I4.4,5(A,I2.2))') Iyr, '/', Imon, '/', Iday, &' ', Ihour, ':', Imin, ':', Isec!
returnend function!***********************************************************************! exp. 1997,12,31,23,59,59 --> 6.223158719900000E+10!***********************************************************************
real(8) function ND2TT(Iyr, Imon, Iday, Ihour, Imin, Isec)!
integer :: IM2D(12,0:1) = &reshape( (/ 0,31,59,90,120,151,181,212,243,273,304,334, &0,31,60,91,121,152,182,213,244,274,305,335 /), (/12,2/) )integer :: Iyr, Imon, Iday, Ihour, Imin, Isec
integer :: Iy4, Iy1, Ileap, Im, Itt!
!
Iy4 = 1461 * ( Iyr / 4 )Iy1 = 365 * mod( Iyr, 4 )!
if ( mod( Iyr, 4 ) .ne. 0 ) thenIleap = 0else
Ileap = 1
end if
Im = IM2D( Imon, Ileap)!
Itt = Iy4 + Iy1 + Im + Iday - Ileap!
ND2TT = Ihour * 3600 + Imin * 60 + IsecND2TT = ND2TT + Itt * 86400.0D0!
returnend function!***********************************************************************! exp. 6.223158719900000E+10 --> 1997,12,31,23,59,59!***********************************************************************
subroutine TT2ND(&
Iyr , Imon , Iday , Ihour, Imin, Isec, & !O & I

175

tt

)

!
integer :: Iyr, Imon, Iday , Ihour, Imin, Isec
integer :: Itt, Iy, Iy4, Iyd, Iy1, Ileap, Imd, Im, Its
integer :: IM2D(12,0:1) = &reshape( (/ 0,31,59,90,120,151,181,212,243,273,304,334, &0,31,60,91,121,152,182,213,244,274,305,335 /), (/12,2/) )integer :: IY2D(4) = (/0,366,731,1096/)
real(8) :: tt, tt0, ND2TT!
!
!

..... ITT [day] .....
Itt = 1 + tt / 86400.0D0!
Iy4 = (Itt-1) / 1461
Iyd = Itt - Iy4 * 1461
do IY = 1, 4
if ( IY2D(Iy) + 1 .le. Iyd ) thenIy1 = Iy
end ifend do!
Iyr = Iy4 * 4 + Iy1 - 1
if ( mod(Iyr,4) .ne. 0 ) thenIleap = 0else
Ileap = 1
end if
IMD = IYD - IY2D(IY1)!
do IM = 1, 12
if ( IM2D(IM,ILEAP)+1 .le. IMD ) then
IMON = IM
end ifend doIDAY = IMD - IM2D(IMON,ILEAP)!
TT0 = ND2TT(IYR, IMON, IDAY ,0,0,0)
ITS = nint( TT - TT0 )Ihour = ITS / 3600
Imin = ( ITS - Ihour * 3600 ) / 60
Isec = ITS - Ihour * 3600 - Imin * 60!
returnend subroutine!

176