PRODUCTION SURPLUS REVISITED
Extensions to the Graham-Schaefer Production Model
The standard Graham-Schaefer surplus production model is based on various simplifying
assumptions, including the following:
Fishing mortality is proportional to fishing effort, F = qf ;
the response in biomass to any changes is immediate; and
the system is entirely deterministic with r and K constant.
Now we will relax each of these assumptions in turn to examine how they influence the behavior of a
fishery whose dynamics are governed by logistic growth.
The Assumption that F = qf
Instead of assuming that the rate of fishing mortality F is directly proportional to the amount of
fishing effort f, suppose there is a handling time for each unit of biomass caught and that
dY
q f
=
B
dt
1 + q h B
We saw this model in Catch Process 3.
Parameter (tau) is the fraction of a trip spent fishing and handling; h is the handling time per unit
catch. The presence of a nonzero handling time can fundamentally alter the equilibrium
characteristics of the exploited fish stock. In the Graham-Schaefer model there is never more than
one nonzero equilibrium value for B(t). With the handling time model, however, there can be two
nonzero equilibria provided the slope of dY/dt versus B at the origin is greater than r and provided
h is large enough so that the horizontal asymptote for dY/dt is sufficiently small. (The slope of
dY/dt versus B at the origin is qf; the asymptote is f/h.)
a
b
dB/dt
dB/dt
<= Stable
<= Unstable
B
Stable
/
B
In the Graham-Schaefer model (the graph on the left) the equilibria are always stable. After any
small perturbation to B, the system will return to the original equilibrium (provided the parameters
remain constant) because fishery removals exceed natural growth when biomass is greater than
the equilibrium level, which causes the biomass to decrease, and vice versa when biomass is less
than the equilibrium level. In the model with handling time (the graph on the right), when there are
FW431/531
Copyright 2008 by David B. Sampson
Surplus4 - Page 169
two nonzero equilibria (as in the higher yield curve), the one closer to the origin is unstable. Any
small perturbation to B will move the system either to the stable equilibrium at the origin, or to the
larger equilibrium biomass. Condrey (1984), on the Supplemental Reading list, discusses this
model and its properties.
Even with a yield model that is linear with B we can get a system that has two nonzero equilibria,
one of which is unstable. For example, suppose that per capita growth is relatively low when stock
size is small (as in the graph on the left below). This will produce an inflection point in the
ascending portion of the population growth function (the graph on the right below). Such a growth
function is sometimes described as exhibiting depensation (as opposed to compensation ).
Per Capita Production
Population Growth
a
b
dB/dt
(dB/dt)/B
<= Stable
<= Unstable
B
B
The equilibria nearer the origin are
unstable. Any small perturbations to
B will force the system away from the
equilibrium levels.
If there is only one point of
equilibrium, at (MSY, Bmsy ),
the system is unstable.
dB/dt
We can even get a system with two nonzero equilibria, one of which is unstable, with a logistic
growth model and a linear yield model, for example, if the harvest is always a fixed yield regardless
of the size of the stock. The yield function in this case is a horizontal line, which intersects the
growth model at either zero, one, or two points.
\
/
Stable
Unstable
B
The Assumption that Response Occurs Instantaneously
Suppose the population productivity does not immediately respond to changes in population size
but instead does so after a time lag.
dN
N(t )
= r N ( t ) 1
F N ( t)
dt
K
This type of equation is sometimes described as a delay-differential equation . Natural growth
FW431/531
Copyright 2008 by David B. Sampson
Surplus4 - Page 170
delay-differential equation
here depends on the population size time units in the past. We can think of parameter (tau) as
a constant that represents the average age at recruitment. (If we did a similar model for population
biomass rather than abundance, would not have such a clear interpretation.) It is also possible to
construct the model so that follows some distribution (e.g., 50% recruit after a one year lag, 30%
after two years, 15% after three years, and 5% after four years.) Models with distributed delays
require much more complicated mathematics, however.
One method for examining the behavior of a delay-differential equation is to conduct a local
stability analysis. A good reference for this type of analysis and for other advanced topics in
population modeling is the book by Nisbet and Gurney (1982) entitled "Modelling Fluctuating
Populations."
Local Analysis Stability for Logistic Growth with a Time Delay
Suppose we have some population whose dynamics are determined by the delay- differential
equation given above. If the population is at equilibrium, then
N ( t) = N ( t ) = Ne = K 1
F
r
This follows directly from setting dN/dt equal to zero, and N(t) equal to N(t-). We can linearize the
DE by taking the first order Taylor series expansion around the equilibrium point Ne. We can then
analyze the linearized form of the DE to determine the approximate behavior of N(t) when it is near
the equilibrium point. First, let us simplify the notation by writing the differential equation as
N = N ( t)
N
dN
= G N , N = r N 1
F N
dt
K
(
)
where
and
N = N ( t )
This differential equation is a function of two variables, N and N.
The linearized differential equation is
dN
dt
where
( N Ne)
dG ( Ne)
dG ( Ne)
+ ( N Ne)
dN
dN
dG ( Ne)
is the partial derivative of G with respect to N, evaluated at N = Ne.
dN
To simplify notation, let
Notice that
n = N Ne
dn
dN dNe
dN
=
=
dt
dt
dt
dt
and
n = N Ne.
because
We will work with dn/dt rather than dN/dt . The equation is
where
FW431/531
a=
dG ( Ne)
=F
dN
and
dNe
dt
is a constant.
dn
= a n b n
dt
b=
Copyright 2008 by David B. Sampson
dG ( Ne)
2
= r 1 Ne
dN
K
Surplus4 - Page 171
b = r 1
W e can simplify coefficient b,
2
F
K 1 = r 2 F .
K
r
Ne
Observe that
b>0
==>
r
= F msy .
2
F<
It can be shown that the solution for n(t) is of the form
n ( t) = C exp ( t) exp ( i t)
where C is an arbitrary constant that is determined by the initial conditions, i is the unit imaginary
number and is equal to the positive square root of minus one, and parameters (mu) and
(omega) satisfy the following equations.
= a + b exp ( ) cos ( )
= b exp ( ) sin ( )
and
Depending on the values of a, b, and , the solution n(t) will exhibit one of the following types of
behavior
II. Stable and Underdamped
n(t)
n(t)
I. Stable and Overdamped
III. Unstable and Overdamped
IV. Unstable and Underdamped
n(t)
time, t
n(t)
time, t
time, t
time, t
Exponential decay (stable)
or growth (unstable)
FW431/531
Regular oscillations with decreasing (stable)
or increasing (unstable) amplitude.
Copyright 2008 by David B. Sampson
Surplus4 - Page 172
For any given values of a, b, and we can determine what kind of behavior the solution will exhibit
by seeing where the values for a and b lie on the following diagram, adapted from Fig. 2.8 of
Nisbet and Gurney (1982).
(1)
The horizontal axis is a and the vertical axis
is b. In our system a = F and is always
positive or zero, and b = r-2F.
4
The 45 line corresponds to a = b , which is
equivalent to
F=r/3
IV
3
2
(2)
II
==>
1
F=r/3.
The horizontal axis corresponds to b = 0,
which is equivalent to
btau
F = r - 2F
atau
0
1
2
3
1
F = r / 2 = F max
The -45 line corresponds to -a = b, which is
equivalent to
4
I
(3)
2
3
III
F=r
-F =r - 2F
==>
F=r.
4
(4)
The 45 diagonals together with the a axis divide the half-plane into four sections representing
systems with the following characteristics:
(1)
F < r/3;
(3) r / 2 < F < r ;
and
(2) r / 3 < F < r / 2 ;
(4)
r < F.
The two curved lines in the upper quadrant together with the -45 diagonal also divide the
half-plane into four regions. Depending on which region the particular values of a and b lie in,
the solution n(t) will exhibit that corresponding type of behavior.
I. Stable / Overdamped
III. Unstable / Overdamped
II. Stable / Underdamped
IV. Unstable / Underdamped
When F is greater than r/2 (the Fmax level of fishing), the behavior of the system is independent of
the value of the time lag t. If r is less than F (section 4, region I), then there is unstable growth of n
away from the equilibrium (n=0). If F is between r and r/2 (section 3, the lower portion of region I),
then there is stable overdamped return back to the equilibrium.
W hen F is less than r/2, the behavior of n(t) depends of the value of the time lag. If F is between
r/2 and r/3 (section 2), then the system is stable and overdamped for near zero, but becomes
underdamped (but stable) as increases. If F is less than r/3 (section 1), then the system goes
through two transitions as increases from zero. Initially the system is stable and overdamped, but
becomes stable and underdamped as increases, and finally becomes unstable with oscillations as
increases further.
Explore a Schaefer model with a time lag using the Mathcad demonstration.
Surplus production models with time lags have not received much attention in the fisheries
literature. The Supplemental Reading list includes two references, Walter (1973) and Beddington
and May (1975).
FW431/531
Copyright 2008 by David B. Sampson
Surplus4 - Page 173