ON THE ASSESSMENT AND SUSTAINABILITY OF THUNDER BAY CISCO
By
Nicholas C. Fisch

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Fisheries and Wildlife – Master of Science
2018

ABSTRACT
ON THE ASSESSMENT AND SUSTAINABILITY OF THUNDER BAY CISCO
By
Nicholas C. Fisch
Cisco, Coregonus artedi, have been reduced to a small fraction of their historical
abundance throughout the Laurentian Great Lakes. Today, remnant spawning stocks
that continue to support commercial fisheries are confined to western Lake Superior.
Although these fish are of particular economic and ecologic importance in the region,
formal stock assessment models have yet to be developed for the species. In addition,
effects of current exploitation rates on these remnant stocks have yet to be evaluated
using quantitative methods. In this thesis, we first develop and compare multiple stateof-the-art stock assessment models for a spawning stock of cisco in Thunder Bay,
Ontario, in an effort to determine an appropriate assessment framework to model the
remaining cisco stocks in western Lake Superior. Results strongly suggest statistical
catch-at-age assessment (SCAA) methods are most appropriate for modeling cisco in
Thunder Bay, and should be applied to stocks in Minnesota and Wisconsin waters of
Lake Superior. We then perform a simplified management strategy evaluation of the
Thunder Bay cisco stock based on the SCAA in an effort to determine both the
sustainability of the current harvest control rule, and the performance of alternate
harvest control rules in managing cisco in Thunder Bay. Results suggest current
exploitation rates are sustainable in Thunder Bay; however, yield, long-term spawning
biomass, and risk of collapse can be improved by implementing control rules involving
biomass based thresholds that decrease exploitation rate at low stock sizes.

ACKNOWLEDGEMENTS

First and foremost, I’d like to thank my advisor Dr. James Bence for his patience,
guidance, and support throughout the development of this thesis. I would like to thank
the Great Lakes Fishery Commission for funding support on this project. I am also
indebted to the Ontario Ministry of Natural Resources and Forestry (OMNRF), who
were generous enough to lend us their outstanding dataset on cisco. Specifically within
OMNRF I’d like to thank Eric Berglund and Anthony Chiodo, who were available to
answer any question I had about the dataset. I would like to thank Ian Harding, who
lent us his dataset on otolith increments of individual cisco, which was crucial to our
project. I would like to thank Dan Yule, who helped me tremendously in understanding
and processing hydroacoustic data. I would also like to thank Jared Myers, Sam
Truesdell, and Rick Clark who aided in various aspects throughout the project and who,
along with Dr. Bence and Dan Yule, developed the proposal that led to the creation and
funding of this project. Lastly, I would like to thank all of my colleagues at the
Quantitative Fisheries Center at Michigan State University for creating such a great
environment to learn about quantitative methods in fisheries science.
The collective work of Jim Bence, Dan Yule, Jared Myers, Sam Truesdell, Rick
Clark, and myself on this project is reflected through the use of the plural “we”
throughout this thesis.

iii

TABLE OF CONTENTS

LIST OF TABLES........................................................................................................ vi
LIST OF FIGURES ..................................................................................................... vii
INTRODUCTION ....................................................................................................... 1
Current Management ........................................................................................ 2
Contents ............................................................................................................. 4
REFERENCES ................................................................................................... 6
CHAPTER 1. A Comparison of Age- and Size-Structured Assessment Models Applied to
a Stock of Cisco in Thunder Bay, Ontario .................................................................. 10
Abstract.............................................................................................................. 10
Introduction ...................................................................................................... 11
Methods ............................................................................................................. 14
Study Species ............................................................................................ 14
Stock Area ................................................................................................. 14
Data ........................................................................................................... 15
Collection & Processing ............................................................................ 16
Fishery Independent Data .......................................................... 16
Commercial Data ........................................................................ 20
Growth Increments ..................................................................... 21
Process Model ........................................................................................... 22
Observation Model ................................................................................... 28
Aging Error ............................................................................................... 33
Likelihood ................................................................................................. 34
Data Weighting ......................................................................................... 38
Model Fitting/Calibration/Troubleshooting ........................................... 40
SCAA ........................................................................................... 40
SCSA ........................................................................................... 41
Comparison .............................................................................................. 42
Results ............................................................................................................... 43
Biomass..................................................................................................... 43
Exploitation rate ....................................................................................... 44
Recruitment .............................................................................................. 46
Abundance ................................................................................................ 48
Aging Error ............................................................................................... 48
Growth ...................................................................................................... 48
Retrospective Analysis ............................................................................. 50
Model Fit to Data ...................................................................................... 52
Computational Intensity .......................................................................... 55
Discussion .......................................................................................................... 55
Conclusions and Recommendations ........................................................ 61
iv

APPENDIX ........................................................................................................ 64
REFERENCES ................................................................................................... 66
CHAPTER 2. Evaluating the Sustainability of a Cisco Fishery in Thunder Bay, Ontario
under Alternate Harvest Control Rules ..................................................................... 71
Abstract ............................................................................................................. 71
Introduction ...................................................................................................... 73
Methods ............................................................................................................. 76
Harvest Control Rules and Policy Parameters ......................................... 76
Performance Statistics .............................................................................. 79
Model ........................................................................................................ 80
Recruitment .............................................................................................. 81
Fishing Mortality ...................................................................................... 84
Assessment Error ..................................................................................... 86
Sensitivity Analyses .................................................................................. 87
Results ............................................................................................................... 88
Recruitment Scenario ............................................................................... 88
Average Yield ............................................................................................ 88
Risk (% of years SB < 20% unfished level) .............................................. 91
Absolute Annual Variation in Yield (AAV) .............................................. 94
Spawning Biomass at the end of the SPM (Final SB) .............................. 96
Performance metrics under similar levels of yield .................................. 99
Sensitivity ................................................................................................. 103
Discussion .......................................................................................................... 104
APPENDIX ........................................................................................................ 117
REFERENCES ................................................................................................... 120

v

LIST OF TABLES

Table 1.1. Data source years for each assessment. Composition refers to age and length
composition data ........................................................................................................ 16
Table 1.2. Description of survey sampling of cisco in Thunder Bay .......................... 19
Table 1.3. Description of how the OMNRF sampled commercial fishery data ......... 20
Table 1.4. Negative log-likelihood and SDNR values for common data sources ...... 52
Table 1.5. Final year (2015) relative differences [(SCAA-SCSA)/SCAA]................... 56
Table 2.1. Performance statistics are presented as results for the 4yr and 7yr recruitment
scenarios (4yr | 7yr). Values are presented as medians over simulations. Average yield
(kg) is mean yield over years. Risk is calculated as the percentage of years SB is below
20% of the unfished condition. AAV is defined in methods. Final spawning biomass is
the median SB of the last 5 years in a simulation. Catch levels for constant catch control
rules are presented in 100,000 kg (i.e. 100k=100,000 kg)....................................... 101

vi

LIST OF FIGURES

Figure 1.1. OMNRF quota management areas pre-2016. We are characterizing Thunder
Bay stock area as zones 1-4 ........................................................................................ 15
Figure 1.2. Sex-specific proportions of weight in samples of the fishery catch ........ 23
Figure 1.3. Age composition of female cisco caught in Thunder Bay, Ontario, from 1999
to 2015 ........................................................................................................................ 36
Figure 1.4. Length composition of female cisco caught in Thunder Bay, Ontario, from
1999 to 2015. Y axis labels denote length bin endpoints, i.e. “[150,160)” refers to a length
bin starting at 150 mm and ending at 160 mm .......................................................... 37
Figure 1.5. Biomass in millions of kg, spawning biomass in millions of kg of mature
females (>250 mm), exploitation rate (yield/biomass of fish > 250 mm), and abundance
in millions of fish for the SCAA and the SCSA. Shaded regions denote 95% HPD
intervals and dashed lines are medians of the posterior distribution. Squares denote
SCAA output while triangles denote SCSA output. Light shading denotes the HPD for
the SCAA, darker shading denotes the HPD for the SCSA, and the darkest shading is
where the two intervals overlap ................................................................................. 45
Figure 1.6. Recruitment for the SCAA and SCSA. Where both represent the number of
fish entering the model in a given year, they are not defined the same way in that for the
SCAA it is the number of age 2s entering the population and in the SCSA it is the
number of fish greater that 170 mm that are entering the population. Points denote
medians of the posterior distribution and error bars are the 95% HPD intervals.... 47
Figure 1.7. Upper left: Fit to growth increment data. Red line depicts the median of the
posterior distribution and shaded region represents 95% HPD intervals. Upper right:
Residuals, medians of the posterior distribution, from fit to growth increment data.
Lower panel: Growth transition matrix at the posterior medians for growth parameters.
Note that the area of the circles represent the probability of growing into a length bin
given a starting length bin. Length bins are represented on axes as midpoints. Plus
group is length bin 410-420 mm................................................................................ 49
Figure 1.8. Retrospective analyses for spawning biomass and exploitation rate. Shown
are times series estimates of spawning biomass and exploitation rate when five terminal
years of data are sequentially dropped from each assessment. Spawning biomass in
millions of kg of mature females and exploitation rate as yield/ biomass of fish > 250
mm .............................................................................................................................. 51
Figure 1.9. Fit to hydroacoustic estimates of spawning abundance. No data from 19992004 and in 2006. Spawning fish in millions of fish. Points denote medians of the
posterior distribution and error bars are 95% HPD intervals ................................... 53
vii

Figure 1.10. Fit to landings data. Harvest (first column) is in 10,000 kg. Note observed
data points in the first column cannot be seen as they are covered by model point
estimates ..................................................................................................................... 54
Figure 1.11. Comparison of models when natural mortality is estimated within SCAA and
when it is assumed known. First column are the original models, while the second
column natural mortality is assumed known for SCAA. Shaded regions denote 95% HPD
intervals and dashed lines are medians of the posterior distribution. Squares denote
SCAA output while triangles denote SCSA output. Light shading denotes the HPD for
the SCAA, darker shading denotes the HPD for the SCSA, and the darkest shading is
where the two intervals overlap ................................................................................. 56
Figure 1.12. Retrospective analyses where both SCAA and SCSA assumed the same
known natural mortality values. Shown are times series estimates of spawning biomass
and exploitation rate when five terminal years of data are sequentially dropped from
each assessment. Spawning biomass in millions of kg of mature females and
exploitation rate as yield/ biomass of fish > 250 mm ............................................... 58
Figure 2.1. Harvest control rules considered in this analysis and associated policy
parameters .................................................................................................................. 78
Figure 2.2. SR curves for each recruitment scenario. “Data”, medians of the posterior
distribution of the SCAA, are plotted as points. The 7yr scenario SR curve uses all “data”
points while the 4yr scenario solely uses the points highlighted in green ................ 84
Figure 2.3. Summary of the distributions of average harvest over the simulation period
for each respective control rule. Shown are medians (horizontal bar) and 25-75 quantiles
(box). Labels specify policy parameters that make up each control rule (CU = “U”; CUT1
= “U SBT”; CUT2 = “U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”).
Exploitation rates are presented as decimals and biomass thresholds as percentages. For
CUT2 control rules, a label of “0.10 50-20%” describes a control rule that has an
exploitation rate of 0.10 above 50% of the unfished spawning biomass, while that rate
linearly declines below that threshold to 0 at 20% of the unfished spawning biomass.
Catch levels are described in 100,000 kg (i.e. 100k = 100,000 kg) .......................... 90
Figure 2.4. Summary of the distributions of risk level for each respective control rule.
Risk is defined as the percentage of years in each simulation where SB is below 20% of
the unfished level. Shown are medians (horizontal bar) and 25-75 quantiles (box).
Labels specify policy parameters that make up each control rule (CU = “U”; CUT1 = “U
SBT”; CUT2 = “U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”). Exploitation
rates are presented as decimals and biomass thresholds as percentages. For CUT2
control rules, a label of “0.10 50-20%” describes a control rule that has an exploitation
rate of 0.10 above 50% of the unfished spawning biomass, while that rate linearly
declines below that threshold to 0 at 20% of the unfished spawning biomass. Catch
levels are described in 100,000 kg (i.e. 100k = 100,000 kg) .................................... 93

viii

Figure 2.5. Summary of the distributions of absolute annual variation in yield for each
respective harvest control rule Shown are medians (horizontal bar) and 25-75 quantiles
(box). Labels specify policy parameters that make up each control rule (CU = “U”; CUT1
= “U SBT”; CUT2 = “U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”).
Exploitation rates are presented as decimals and biomass thresholds as percentages. For
CUT2 control rules, a label of “0.10 50-20%” describes a control rule that has an
exploitation rate of 0.10 above 50% of the unfished spawning biomass, while that rate
linearly declines below that threshold to 0 at 20% of the unfished spawning biomass.
Catch levels are described in 100,000 kg (i.e. 100k = 100,000 kg) .......................... 95
Figure 2.6. Summary of the distributions of final spawning biomass for each respective
control rule, with final spawning biomass defined as the median of the last 5 years
spawning biomass in each simulation, characterized as a percentage of the unfished
level. Shown are medians (horizontal bar) and 25-75 quantiles (box). Labels specify
policy parameters that make up each control rule (CU = “U”; CUT1 = “U SBT”; CUT2 =
“U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”). Exploitation rates are
presented as decimals and biomass thresholds as percentages. For CUT2 control rules, a
label of “0.10 50-20%” describes a control rule that has an exploitation rate of 0.10
above 50% of the unfished spawning biomass, while that rate linearly declines below
that threshold to 0 at 20% of the unfished spawning biomass. Catch levels are described
in 100,000 kg (i.e. 100k = 100,000 kg) ..................................................................... 98
Figure 2.7. Comparison of performance metrics for CU rule with an exploitation rate of
0.05 (black points) and CC rule with a catch level of 250,000 kg (red points). Upper
panel plots risk (percentage of years SB < 20% unfished) versus average harvest (in
millions of kg) obtained in the same individual simulations. Middle panel compares
final spawning biomass (millions of kg) and average harvest. Lower panel compares
absolute annual variation in yield and average harvest. Only results from simulations
that produced an average harvest within the 25-75% quantiles (i.e., inside the
interquartile range) are plotted ................................................................................. 100
Figure 2.8. Spawning biomass for the projection of the current harvest control rule, a
10% exploitation rate. Shown are medians (horizontal bar) and 25-75 quantiles (box)
.................................................................................................................................... 105

ix

INTRODUCTION
Cisco, Coregonus artedi, also known as lake herring, are an ecologically and
economically important fish species of the Laurentian Great Lakes ecosystem. They
historically served as the primary prey for many dominant predators throughout the
region prior to the invasion of rainbow smelt, Osmerus mordax, and alewife, Alosa
pseudoharengus (Dryer et al. 1965; Smith 1968; Leach and Nepszy 1976; Christie et al.
1987; Jude et al. 1987; Diana 1990; Wolfert and Bur 1992; Conner et al. 1993; O’Gorman
and Stewart 1999; Ray et al. 2007). This has led some to postulate that restoration of a
key native predator, lake trout, Salvelinus namaycush, in the Great Lakes ecosystem is
dependent on the replacement of exotic alewives and rainbow smelt by the native cisco
(Bronte et al. 2008; Markham et al. 2008). In addition to their ecological significance, a
commercial fishery for cisco has operated in the Great Lakes for over a century primarily
targeting spawning females to supply foreign roe markets (Stockwell et al. 2009).
During the first half of the 20th century, this fishery accounted for a greater amount of
yield taken out of the Great Lakes than that of any other species (Baldwin et al. 2002).
However, since that time, cisco fishery yield has drastically decreased, due to the
sequential collapse of populations in each of the Great Lakes, largely attributed to
overexploitation; Lake Erie in the 1920s (Hartman 1973; Selgeby 1982), Lakes Ontario
and Huron in the mid-1950s (Berst and Spangler 1973; Christie 1973), Lake Michigan in
1960 (Wells and McLain 1973), and Lake Superior in the early 1960s (Selgeby 1982).
Today the fishery operates on a much smaller scale mostly focusing on stocks in western
Lake Superior that have rebounded during the post-collapse period. Due to their
importance in the region, rehabilitation and protection of self-sustaining stocks of cisco

1

to support a stable production of predators and a sustainable commercial fishery has
become a priority in Lake Superior (Schreiner et al. 2006).
Current Management
The cisco fishery in Ontario and Minnesota waters of Lake Superior is currently
managed based on a fixed exploitation rate control rule. A fixed exploitation rate control
rule aims to set catch quotas to some constant proportion of stock size (Walters and
Martell 2004). This builds in an inherent feedback system; as the stock declines, the
quotas do the same, and vice versa. For this type of control rule to be effective,
management agencies need estimates of population abundance each year, or each time
period they would like to enact catch quotas (Walters and Martell 2004). In Ontario and
Minnesota, this is facilitated by hydroacoustic surveys of spawning cisco abundance,
where total allowable catch (TAC) is set at 10% of cisco>250mm in Ontario, and 10% of
cisco>305mm in Minnesota, both calculated from the hydroacoustic surveys.
Management of cisco in the state of Wisconsin does not currently involve setting a TAC;
however numerous restrictions on gill net mesh sizes are employed in different areas, in
addition to the prohibition of fishing in a number of refuges and restricted-use areas
(Stockwell et al 2009). There is currently no targeted fishery for cisco in Michigan
waters of Lake Superior as harvest has been restricted to incidental catches in the chub
fishery.
Given these current management strategies in western Lake Superior, there are
several key concerns related to the long term sustainability of the fishery. Primarily,
there is large concern for overfishing in Wisconsin waters of Lake Superior, due to the
lack of quota allocations or TAC. This concern stems from the lucrative nature of the

2

fishery, where catch will largely be determined by market demands. It also stems from
the highly efficient nature of the fishery in western Lake Superior, which targets
aggregations of fish in November as they congregate to spawn. Under this scenario,
fishers can move from dense patch to dense patch of spawning fish maintaining high
catch rates, potentially depleting the resource as their catch-per-unit-effort does not
necessarily decrease (hyperstability, Hilborn and Walters 1992). Another concern
relates to the reliance on hydroacoustic surveys to set a TAC each year, as these surveys
are not necessarily done every year due to weather restrictions and budget constraints,
or in some years could produce anomalous results due to factors such as fluctuations in
spatial distributions. The integration of other sources of data such as fishery catch-atage into a formal stock assessment model can allow the calculation of abundance even in
years when no hydroacoustic survey is performed, in addition to improved estimates in
years when surveys are done. Lastly, while fixed exploitation rate control rules can
sometimes effectively achieve objectives (Walters and Martell 2004, Deroba and Bence
2008), the specific exploitation rate of 10% put into place in Ontario and Minnesota
waters has not been evaluated using quantitative methods. In a more formal harvest
policy analysis, a stock may be forecasted into the future under a variety of different
management strategies using a population model to determine which strategy is optimal
in promoting fishery objectives, such as high or consistent yield, and/or sustainability of
the stock. These 10% exploitation rates were actually put into place based on
exploitation rates seen as sustainable for other Lake Superior fish stocks such as lake
trout, lake whitefish, Coregonus clupeaformis, and lake sturgeon, Acipenser fulvescens
(Ebener et al. 2008, Stockwell et al. 2009). Where precautionary approaches to
management such as these are an important first step, a harvest policy tailored to what

3

is known about cisco stocks in western Lake Superior is needed to ensure sustainability
of the cisco fishery into the future.
Schreiner et al. (2006) described an objective to rehabilitate and protect
Coregonine stocks by exploring “the use of a population model to compliment the use of
an acoustics-based model” in determining TAC for cisco. In addition, the fisheries
management plan for Minnesota waters of Lake Superior from 2016-2025 (Goldsworthy
et al., 2015) states an objective to “Continue to support establishing a lake-wide stock
assessment model to complement existing acoustics-based quota calculations”. Ebener
et al. (2008) stated that “there is an overwhelming need to develop an overarching
research framework to better understand cisco population dynamics and ecology within
the context of managing them in a sustainable fashion”. These statements describe a
shared view that integrating additional sources of data in characterizing cisco
population dynamics is essential in promoting the sustainable management of cisco
stocks in Lake Superior.
Contents
The objectives of this thesis were to (1) develop an integrated assessment model
that is effective in estimating abundance and characterizing population dynamics of
cisco, and (2) to evaluate the sustainability of the current control rule, and the
performance of alternate harvest control rules in meeting cisco fishery objectives. Given
the superior quality and breadth of data available in Ontario waters of Lake Superior,
this thesis focused on the Thunder Bay cisco stock. The assessment and analysis of
Thunder Bay cisco can serve as a framework to later model other cisco stocks in Lake
Superior (Minnesota, Wisconsin). This thesis consists of two chapters, the first of which

4

is a comparison of two different stock assessment methods, catch-at-age and catch-atsize, applied to cisco in Thunder Bay in an effort to determine which assessment method
best characterizes population dynamics of cisco. The second chapter uses the best
assessment model developed in the first chapter to implement a simplified management
strategy evaluation of Thunder Bay cisco involving simulations of the stock under
multiple sources of uncertainty and harvest control rules to evaluate sustainability and
performance of specific harvest control rules in meeting fishery objectives.

5

REFERENCES

6

REFERENCES

Baldwin, N. A., R. W. Saalfeld, M. R. Dochoda, H. J. Buettner, and R. L. Eshenroder.
2002. Commercial fish production in the Great Lakes 1867–2000. Great Lakes
Fishery Commission, Ann Arbor, Michigan.
Berst, A. H., and G. R. Spangler. 1973. Lake Huron: the ecology of the fish community
and man’s effects on it. Great Lakes Fishery Commission Technical Report 21.
Bronte, C. R., C. C. Krueger, M. E. Holey, M. L. Toneys, R. L. Eshenroder, and J. L.
Jones. 2008. A guide for the rehabilitation of lake trout in Lake Michigan. Great
Lakes Fishery Comm ,;]ission, Miscellaneous Publication 2008-01, Ann Arbor,
Michigan.
Christie, W. J. 1973. A review of the changes in the fish species composition of Lake
Ontario. Great Lakes Fishery Commission Technical Report 23.
Christie, W. J., K. A. Scott, P. G. Sly, and R. H. Strus. 1987. Recent changes in the
aquatic food web of eastern Lake Ontario. Canadian Journal of Fisheries and
Aquatic Sciences 44(Supplement 2):37–52.
Conner, D. J., C. R. Bronte, J. H. Selgeby, and H. L. Collins. 1993. Food of salmonine
predators in Lake Superior, 1981–1987. Great Lakes Fishery Commission
Technical Report 59.
Deroba, J.J., and Bence, J.R. 2008. A review of harvest policies: understanding relative
performance of control rules. Fisheries Research 94: 210–223.
Diana, J. S. 1990. Food habits of angler-caught salmonines in western Lake Huron.
Journal of Great Lakes Research 16:271–278.
Dryer, W. R., L. F. Erkkila, and C. L. Tetzloff. 1965. Food of lake trout in Lake Superior.
Transactions of the American Fisheries Society 94:169–176.
Ebener, M.P., Stockwell, J.D., Yule, D.L., Gorman, O.T., Hrabik, T.R., Kinnunen, R.E.,
Mattes, W.P., Oyadomari, J.K., Schreiner, D.R., and Geving, S. 2008. Status of
cisco (Coregonus artedi) in Lake Superior during 1970-2006 and management
and research considerations. Ann Arbor, Michigan: Great Lakes Fishery
Commission, Lake Superior Technical Report 1.
Goldsworthy, C.A., K.A. Reeves, J.E. Blankenheim, and N.R. Peterson. 2015. Fisheries
Management Plan for the Minnesota Waters of Lake Superior. Minnesota
Department of Natural Resources, Division of Fish and Wildlife, Fisheries
Management Section.

7

Hartman, W. L. 1973. Effects of exploitation, environmental changes, and new species
on the fish habitats and resources of Lake Erie. Great Lakes Fishery Commission
Technical Report 22.
Hilborn, R. and C. J. Walters. 1992. Quantitative Fisheries Stock Assessment: Choice,
Dynamics and Uncertainty. London: Chapman and Hall. 175-191.
Jude, D. J., F. J. Tesar, S. F. Deboe, and T. J. Miller. 1987. Diet and selection of major
prey species by Lake Michigan salmonines, 1973–1982. Transactions of the
American Fisheries Society 116:677–691.
Leach, J. H., and S. J. Nepszy. 1976. The fish community in Lake Erie. Journal of the
Fisheries Research Board of Canada 33:622–638.
Markham, J., A. Cook, T. MacDougall, L. Witzel, K. Kayle, C. Murray, M. Fodale, E.
Trometer, F. Neave, J. Fitzsimons, J. Francis, and M. Stapanian. 2008. A
strategic plan for the rehabilitation of lake trout in Lake Erie, 2008–2020. Great
Lakes Fisheries Commission, Ann Arbor, Michigan.
O’Gorman, R., and T. J. Stewart. 1999. Ascent, dominance, and decline of the alewife in
the Great Lakes: food web interactions and management strategies. Pages 489–
513 in W. W. Taylor and C. P. Ferreri, editors. Great Lakes fisheries policy and
management. Michigan State University Press, East Lansing.
Ray, B. A., T. R. Hrabik, M. P. Ebener, O. T. Gorman, D. R. Schreiner, S. T. Schram, S. P.
Sitar, W. P. Mattes, and C. R. Bronte. 2007. Diet and prey selection by Lake
Superior lake trout during spring, 1986–1991. Journal of Great Lakes Research
33:104–113.
Schreiner, D.R., J.J. Ostazeski, T.N. Halpern, and S.A. Geving. 2006. Fisheries
management plan for the Minnesota waters of Lake Superior. Minnesota
Department of Natural Resources, Division of Fish and Wildlife, Fisheries
Management Section.
Selgeby, J. H. 1982. Decline of lake herring (Coregonus artedii) in Lake Superior: an
analysis of the Wisconsin herring fishery, 1936–1978. Canadian Journal of
Fisheries and Aquatic Sciences 39:554–563.
Smith, S. H. 1968. Species succession and fishery exploitation in the Great Lakes.
Journal of the Fisheries Research Board of Canada 25:667–693.
Stockwell, J.D., Ebener, M.P., Black, J.A., Gorman, O.T., Hrabik, T.R., Kinnunen, R.E.,
Mattes, W.P., Oyadomari, J.K., Schram, S.T., and Schreiner, D.R. 2009. A
synthesis of cisco recovery in Lake Superior: implications for native fish
rehabilitation in the Laurentian Great Lakes. North American Journal of
Fisheries Management 29(3): 626–652.

8

Walters, C.J. and S.J. Martell, 2004. Fisheries ecology and management. Princeton
University Press.
Wells, L., and A. L. McLain. 1973. Lake Michigan: man’s effects on native fish stocks and
other biota. Great Lakes Fishery Commission Technical Report 20.
Wolfert, D., and M. T. Bur. 1992. Selection of prey by walleyes in the Ohio waters of the
central basin of Lake Erie. U.S. Fish and Wildlife Service, Research Publication
182.

9

Chapter 1
A Comparison of Age- and Size-Structured Assessment Models Applied to a
Stock of Cisco in Thunder Bay, Ontario
Abstract
Stock assessment is a critical component of the fisheries management process,
involving the calculation of key variables used in making management decisions.
However; there is still considerable uncertainty in assessment science as to which class
of models is appropriate to use depending on circumstances. A common class of models
used when age data are available are statistical catch-at-age assessment (SCAA) models,
which track cohorts through time. When age data are unavailable, as is often the case in
invertebrate fisheries where the lack of a bony structure such as otoliths makes aging
difficult, statistical catch-at-size assessment (SCSA) models are more often employed,
tracking fish or invertebrates through time by size classes rather than ages. Do SCAA
models actually perform better than SCSA models when age data are available, or is this
just an assumption we make in fisheries research and management? We examined this
question as it relates to a specific case study by evaluating the effectiveness of both
SCAA and SCSA models in characterizing cisco, Coregonus artedi, population dynamics
in Thunder Bay, Ontario. Both models were fit using an integrated framework with
multiple sources of data including hydroacoustic estimates of spawning stock,
commercial and fishery independent age/length compositions, and landings data. Our
results suggest that for cisco in Thunder Bay, parameter confounding resulting in the
inability to estimate natural mortality hampered the utility of a SCSA model in
comparison with a SCAA model when age composition data were available.

10

Introduction
Stock assessment is a critical aspect of fisheries research and management,
supporting the calculation of key quantities such as spawning biomass, abundance,
exploitation rate, recruitment, and their associated uncertainties. Most assessments
done in the United States are based on age-structured assessment methods (Punt et al.,
2017), which, when statistically fit, can be referred to as statistical catch-at-age
assessment (SCAA) models. These models track cohorts of fish through time, using
observations of catch-at-age and auxiliary information to estimate population
parameters (Fournier and Archibald, 1982; Deriso et al., 1985). When catch-at-age data
are unavailable for a species of interest, as is the case in many invertebrate fisheries
where lack of a bony structure such as an otolith makes aging difficult, size-structured
assessment methods are often employed (Punt et al., 2013). Similarly, when statistically
fit these types of models can be referred to as statistical catch-at-size assessment (SCSA)
models. Sullivan et al. (1990) developed and applied a framework for SCSA, which differ
from SCAA in that they utilize observations of catch-at-size and track fish in size bins
rather than age classes through time, often making use of a growth model that
determines transition probabilities of size bins in subsequent time steps. While agestructured models can be fit using harvest size composition data (and a model to convert
predicted age compositions to size compositions; Fournier et al., 1990; Fournier et al.,
1998), often together with age composition data (Methot and Wetzel 2013, Punt et al.
2013), contemporarily the use of SCSA is preferred when the sole or primary harvest
composition data is for sizes rather than ages (Punt et al., 2013).

11

Each method offers distinct advantages and disadvantages. For size based
methods, the model can directly account for the size structure of removals from a
population (Punt et al., 2017), it can more appropriately model some fishery processes
such as selectivity as size-based (although in some cases both age and size may be
involved), and importantly, size composition data is almost always more abundant and
it is both easier and cheaper to collect. SCSA models can considerably decrease the
number of fish that need to be aged, as age compositions of the catch are not required.
SCSA is not without its challenges. Primary among them, as previously mentioned, an
SCSA needs a method, often a growth model, to determine transition probabilities of
fish or invertebrates through size bins for each time step, where additional aspects such
as time-varying or density dependent growth can add complexity. This is not so in SCAA
models which benefit from the fact that a fish must be a year older in the next (yearly)
time step; A caveat being that our ability to observe ages is not perfect, as there is
measurement error involved in aging organisms, and ignoring this error can result in
biased assessment output (Coggins and Quinn, 1998; Reeves 2003; Bertignac and
Pontual, 2007). Although aging error is generally overlooked in SCAA models, it can be
accounted for both within a model (Thompson et al., 2011; Methot and Wetzel, 2013)
and using quality control in aging techniques (Campana, 2001).
Most likely due to the ideal transition of fish through age bins (and the properties
that come along with it, i.e., recruitment into the model), very seldom are SCSA models
developed for species when age data are available. Additionally, few studies have
compared the two methods. One such study, Punt et al., (2017), used simulation analysis
to compare the performance of age-, size-, and age- and size-structured assessment

12

methods, concluding that based on an age- and size-structured operating model, sizestructured and age- and size-structured assessment methods performed best, while agestructured methods performed poorest. This study was done, as are most simulation
studies, based on known population dynamics pre-specified by researchers. The
advantage of simulation studies is the ability to compare assessment results to what is
pre-specified in the operating model as the true population dynamics of the stock. This
specification of the operating model can also be a disadvantage, if the researchers
conception on the dynamics of the stock and fishery (e.g. survey selectivity as age-based
process in Punt et al., 2017), do not actually reflect underlying processes. Fitting
alternative models to real data can be highly useful in helping to better define plausible
processes. Thus comparing alternative models fit to real data and simulation studies are
not alternative approaches, but rather synergistic ways to advance assessment methods.
In this paper we develop and fit both integrated SCAA and SCSA models for a
stock of cisco, Coregonus artedi, in Thunder Bay (Lake Superior), Ontario. Our
objective was to compare and contrast performance of the different assessment methods
when applied to a real stock to provide recommendations on which type of model may
be preferred under different scenarios. We were specifically interested in the overall
question: “when age data are present, are we, in fisheries research and management,
assuming age-structured models will outperform size-structured models because they
are more directly intuitive and easier to fit, or can size models be used instead to better
model fishery processes and perhaps decrease the burden on natural resource agencies
to age so many fish each year?” To our knowledge, only one study has performed a
comparison between age- and size-structured models on a real stock with true dynamics

13

unknown (Akselrud et al. 2017, concluding that age-structured fit data best). In a time of
shrinking natural resource agency budgets, it seems these comparisons are overdue, and
could provide important information for future simulation studies.
Methods
Study Species
Cisco are a planktivorous species native to the Laurentian Great Lakes. They are
largely pelagic, however form annual spawning aggregations during the month of
November in nearshore bays and areas of western Lake Superior (remaining stocks are
largely confined to western Lake Superior). These aggregations have supported a
lucrative commercial roe fishery for decades, as fisherman generally target spawning
fish during late November using suspended gillnets. Additionally, since 2005 these
aggregations have been surveyed annually using hydroacoustics in Thunder Bay.
Current management relies on a fixed exploitation rate control rule where 10% of the
spawning fish estimated from the hydroacoustic surveys are allocated as quota among a
limited number of fishers. No formal assessment models have been developed for this
stock or others in western Lake Superior.
Stock Area
We treated Ontario Ministry of Natural Resources and Forestry (OMNRF)
management zones 1-4 (Figure 1.1) as the stock area for Thunder Bay cisco. This stock
has been hypothesized to be a discrete spawning stock, primarily on the basis that an
adjacent depleted stock in Black Bay has not shown any sign of recovery since a collapse
in the 1980s, which would be expected if cisco from Thunder Bay were spilling over into

14

Black Bay under a non-discrete spawning stock scenario (Ebener et al. 2008).
Additionally this area was chosen based on coverage of the hydroacoustic surveys, which
generally sample over zones 1-4 in Thunder Bay.

Figure 1.1. OMNRF quota management areas pre-2016. We are characterizing Thunder
Bay stock area as zones 1-4.

Data
The SCAA and SCSA models made use of six main sources of observed data in the
fitting process (Table 1.1): (1) Number of cisco > 250 mm in Thunder Bay estimated
from hydroacoustic surveys (2005, 2007-2015), age or size composition of cisco caught

15

in (2) mid-water trawls (2005, 2007-2010, 2015) and (3) multi-mesh gillnets (2009,
2013-2015) that accompany the hydroacoustic surveys, (4) age or size composition of
the fishery catch subsamples (1999-2015), and (5) male and (6) female biomass in the
fishery yield each year (1999-2015). The SCSA made use of one additional source of
data; (7) individual growth increments of cisco.
Table 1.1. Data source years for each assessment. Composition refers to age and length
composition data.
Year
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015

Hydroacoustic
Survey

X
X
X
X
X
X
X
X
X
X

Fishery
Harvest
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

Fishery
Composition
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

MWT Survey
Composition

Gillnet Survey
Composition

X
X
X
X
X

X

X

X
X
X

Collection & Processing
Fishery Independent Data
The United States Geological Survey (USGS) and the OMNRF collectively
sampled spawner abundance in Thunder Bay, Ontario, using acoustics in 2005 and from
2007-2015. Acoustic data was collected and processed using techniques according to the
Great Lakes standard operating procedure for fishery acoustic surveys (Rudstam et. al
2009). The USGS performed the acoustic surveys in 2005 and from 2007-2010 while

16

the OMNRF performed the surveys from 2011-2015. Each year the acoustic gear
specifications differed slightly, utilizing different frequency split beam transducers with
a half-power beam width varying from 5.3-6.7 degrees. The transducer was generally
deployed from 1-5m below the surface depending on the year. The transducer emitted 35 pings/s, with pulse duration usually set at around 0.4 ms. Vessel position was
measured with a differentially corrected Global Positioning System unit (accurate to
1m), and survey path information was imbedded in the acoustic data files. Acoustic data
was processed with Echoview Software (Sonar Data, Tasmania, Australia). Thresholds
for the Sv (scattering volume) and single target echograms were set at -65 dB and -60
dB respectively. A line to exclude surface noise was set at 2-7m depth (depending on sea
state) and a line to exclude the bottom was set 0.5 m above the bottom signal. Softwaredefined bottom lines were examined carefully to ensure all bottom echoes were properly
excluded, and all segments of echograms containing electrical or other noise (i.e., echoes
obviously not from fish) were eliminated before estimating fish densities.
Total fish densities (number/ha) were calculated for 15-m vertical cells over 3km
horizontal intervals using echo integration methods described in Yule et al. (2006). To
obtain densities of spawning ciscoes (>250 mm TL), we multiplied total fish density in
each cell by the proportion of single targets larger than -35.6 dB, a threshold that is
consistent with the predicted target strength of a 250 mm cisco (Yule et al. 2006). We
summed all the vertical cells down to the lake bed for each 3 km interval. In order to
separate density of large cisco from the density of large fish belonging to other species
we took the density estimate of large fish (>250 mm) for each interval and matched it up
to the nearest mid-water trawl or multi-mesh gillnet sample by Euclidean distance. We

17

then multiplied the large fish acoustic density estimate for each cell by the proportion of
large (>250 mm) cisco caught out of all large fish caught in the nearest bio-data sample
(trawl or gillnet). For example, if there was a density estimate of 100 large fish/ha for an
interval and the nearest bio data sample caught 80% cisco out of fish caught greater
than 250 mm, then the large cisco density estimate would be 80 spawning cisco/ha for
that interval. To estimate abundance of spawning cisco in Thunder Bay each year, we
averaged the densities of large cisco for all 3 km intervals and multiplied this by the
surface area of the bay (77,749 ha, zones 1-4).
As mentioned above, fishery independent bio-data were collected during the
acoustic surveys to separate cisco densities from other species densities. Two different
types of gear were used in collecting biological samples, mid-water trawls and multimesh gillnets. The USGS used mid-water trawls to collect cisco biological samples in
2005 and from 2007-2010, while the OMNRF used multi-mesh gillnets to collect
biological samples from 2013-2015. No biological samples were collected in 2011 and
2012, for these years we decided to take the average over the time series of the
percentage of large cisco caught out of all large fish caught each year (≈95%) to
substitute for direct samples of the percentage of large targets attributed to cisco. In
2009 and 2015 each agency collected biological samples using both gears, trawls and
gillnets, in order to observe how sampling gear affected spawning cisco density
estimates (sampling gear had minimal impact on acoustic density estimates of large
cisco). In 2009 we decided to use the acoustic density estimates based on the mid-water
trawl apportionment due to the fact that only 4 gillnets were deployed (biological
samples from these gillnets were still used for composition data). In 2015 we averaged

18

the estimates from the mid-water trawl apportioned densities and gillnet apportioned
densities to calculate one spawning cisco density estimate for the year. In 2007 the
USGS performed two acoustic surveys, one in mid-November and one in late November,
to observe if spawner abundance changed as the season progressed (Acoustic densities
did not differ significantly between the two surveys). We computed an overall mean
acoustic density for 2007 based on a weighted average using the number of 3km
intervals that were sampled for each survey to obtain one spawning cisco estimate for
2007.
During the collection of biological samples, cisco were usually subsampled for age
although in some years all cisco caught were aged. A complete description of how cisco
caught in mid-water trawls and multi-mesh gillnets were subsampled for aging each
year can be found in Table 1.2. We applied sex-specific age-length keys (ALK) using
10mm bins to aging data each year to calculate the age composition of male and female
cisco caught using mid-water trawls and multi-mesh gillnets in Thunder Bay. For size
composition data, we allocated all cisco caught in survey gear each year into sex-specific
10mm length bins.
Table 1.2. Description of survey sampling of cisco in Thunder Bay.
Year
2005
2007
2008
2009
2010
2011
2012
2013
2014
2015

How ages were sampled (Fishery
Independent Data)
Targeted 20 fish per 50mm length bin
Targeted 40 otolith pairs for both males and females
per 50mm length bin
Targeted 40 otolith pairs for both males and females
per 50mm length bin
Uncertain for MWT, all gillnet fish were aged
All aged
No bio-data collected with AC survey
No bio-data collected with AC survey
All aged
All aged
All aged except 33 in gillnets

19

Number sampled
(MWT-Gillnet)
794 - 0
1845 - 0

Number Aged
(MWT-Gillnet)
196 - 0
441 - 0

559 - 0

278 - 0

994 - 302
520 - 0
0-0
0-0
0 - 678
0 - 135
478 - 824

297 - 302
520 - 0
0-0
0-0
0 - 677
0 - 134
473 - 795

Commercial Data
Licensed operators were required to report daily catch to the Ontario Commercial
Fisheries Association, this allowed us to calculate aggregate yield for management zones
1-4 each year. To obtain biological data from the fishery, the OMNRF also collects the
first 10 cisco caught in each gillnet on each day. These fish were measured to the nearest
millimeter (total length and fork length), sexed, and weighed to the nearest gram. Ages
were subsampled in many years based on ALK fixed allocation bin sampling where 10
fish belonging to a 10 mm length bin of each sex and zone were selected at random to be
aged. If less than 10 fish were in a 10 mm bin, then all fish in that bin were aged. Some
years all cisco sampled from commercial gillnets were aged. Annual subsampling of
cisco caught in commercial gillnets for aging is summarized in Table 1.3.
Table 1.3. Description of how the OMNRF sampled commercial fishery data.
Year
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015

Subsampled
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
No
No
Yes
No
No
No
Yes
Yes

Fish sampled – Fish aged
860 - 402
3241 - 594
1221 - 574
1147 - 676
1208 - 704
1091 - 527
661 - 280
644 - 378
839 - 330
654 - 654
638 - 637
500 - 219
563 - 562
478 - 477
429 - 427
733 - 517
705 - 457

We combined aging data for management zones 1-4 each year and applied sexspecific ALKs using 10mm bins to calculate the age composition of male and female
cisco caught in commercial gillnets in Thunder Bay. For size compositions, we allocated
20

sampled commercial cisco each year to sex-specific 10 mm length bins. If fish had a fork
length (FL) but no total length (TL) record, total length was estimated based on a TL-FL
linear regression derived from the rest of the commercial data that had both fork length
and total length. A sex-specific weight-length power relationship assuming a
multiplicative error structure

W  aLb
was fit by a linear regression of log weight on log length using the commercial data. This
was used to estimate missing weights of harvested fish (through a bias-corrected WL
formula)

W  aL e
b

(

2
2

)

where a is the exponential regression intercept estimate, b is the regression slope
estimate, and  is the regression estimate of the residual error variance. Aggregate
2

yield of management zones 1-4 in Thunder Bay each year was separated into male and
female yield by multiplying the total harvest by the proportion of male and female
weight in the commercial database each year, respectively (Figure 1.2).
Growth Increments
Sagittal otoliths were removed from cisco collected in the western arm of Lake
Superior in both Minnesota and Wisconsin waters during the period 1988-2015 by the
Wisconsin Department of Natural Resources. These otoliths were sectioned and data on
increments in the radii of these otoliths were collected as part of a Master’s thesis
(Harding 2017). Additional details on otolith preparation and radii measurement can be

21

found in appendix 1. We back-calculated growth of individual cisco using the scale
proportional hypotheses (Francis 1990), where the length of the cisco at each otolith
annuli Li , is calculated as

Li 

Si
SC

a a

 LC   
b b


where S i is the radius of the otolith at annuli i , S C is the radius of the full otolith at
capture, and LC is the length of the fish at capture. a and

b

are obtained by

regressing length of fish on the length of that fish’s otolith radius.

E ( SC | LC )  a  bLC
For the final data on individual cisco growth increments, we only used growth during
the last full annulus for each fish in an attempt to avoid a substantial Lee’s
phenomenon, giving us a total of 169 individual cisco growth increments.
Process Model
Predicted quantities needed to compare to the observed data listed above were
calculated using a variety of equations describing the stock and fishery. The models ran
from 1999 to 2015. The SCAA model ages began at age 2 and ended at a plus group age
of 15 (denoting all cisco older than 14) while the SCSA model size bins were divided in
10 mm increments beginning at 170 mm and ending at a plus group of 410 mm
(denoting all cisco ≥ 410 mm). Age or size bins are referenced throughout the
manuscript with subscript

j . Given the fishery operates primarily as a roe fishery, it

captures a disproportionate number of females in Thunder Bay each year (Figure 1.2),

22

we decided to make the assessment models sex specific, essentially tracking the number
of male and female fish separately through time. This is presented in subsequent
equations with the subscript s , denoting sex.

Figure 1.2. Sex-specific proportions of weight in samples of the fishery catch.

To begin, for the SCAA we freely estimated (as parameters) cisco numbers at age
for each sex in the first year of the assessment model. For the SCSA, cisco abundance in
the first year for each sex was estimated through a combination of size and abundance
components; one gamma distribution denoting the initial size composition multiplied
by two freely estimated values of abundance, one for each sex.

23

x


1
g( x | ,  ) 
x 1e 

( ) 

Ij  

j*5

j*5

g ( x |  ,  )dx

where  and  are the estimated shape and scale of the gamma distribution, I j is the
initial abundance proportion for size bin j , and

j*

is the size bin midpoint. We

attempted alternative parameterizations including assuming a separate initial size
composition for each sex and estimating one value of abundance, in addition to
estimating two size compositions and two abundance values, each of which did not
outperform in terms of model comparison criterion PSIS-LOO (Pareto smoothed
importance sampling leave-one-out; Vehtari et al., 2017). PSIS-LOO is an efficient
approximation of the exact cross-validation model comparison criterion and has been
shown to be asymptotically equal to WAIC however more robust in cases with weak
priors or influential observations (Vehtari et al., 2017), both of which occur in
assessment modeling.
Sex-specific numbers of cisco at each age or size in each year was obtained from
the exponential survival equation

SCAA

N j 1, y 1,s  N j , y ,s e
j'

SCSA

( Fj , y ,s  M s )

N j ', y 1,s   Pj , j ' N j , y ,s e
j 1

24

( F j , y , s  M s )

 R j ', y 1,s

where N j , y ,s denotes the number of cisco in age or size bin

j alive in year

y of sex s .

Pj , j ' is the size transition matrix, denoting the probability that a fish in length bin j will

grow into length bin
length bin

j' .

j'

in the next time step. R j ' is the number of fish that recruit into

F j , y ,s represents the sex specific instantaneous fishing mortality in year

y for a given age or size bin. M s is the instantaneous natural mortality for each sex,

assumed constant over time and ages or sizes. We decided to allow natural mortality to
vary by sex based on information from previous studies which indicated that male cisco
may experience higher natural mortality than their female counterparts (TeWinkel et al.
2002, Yule et al. 2008). We evaluated different scenarios where natural mortality was
assumed equal for both sexes. Ultimately, and primarily in the interest of numerical
stability, we chose to add a prior on each natural mortality variable based on the
updated Hoenig linear model surrogate equation from Then et al. (2014).
We calculated the size transition matrix, Pj , j ' , using growth parameters L and

K to model an average growth increment and parameters a and b to model the
variance in growth increment as a function of the expected growth increment
E (  j )  ( L  j  )(1  e K )

Var (  j )  a  b * E (  j )

where j  is the length bin midpoint. The probability of remaining in the same length bin
in the next time step,

Pj , j (the diagonal of the matrix), was calculated by integrating from

5 mm to 0 (assumes fish are at length bin midpoint). Assuming no negative growth, the

25

rest of the transition matrix,

Pj , j ' , was calculated by plugging in the expected growth

increment and variance in growth increment into the cumulative distribution function
(CDF) of the gamma distribution, integrating from length bin midpoint +5 mm ( j2' ) to
length bin midpoint -5 mm ( j1' ).
j2'

Pj , j '   ' g (  j |  j ,  j )d j
j1

where  j and  j can be calculated from E (  j ) and V (  j ) . A full description of the
gamma likelihood can be found in the likelihood section of this chapter.
For the SCAA, we modeled recruitment in each year as multiplicative deviations
about a mean recruitment level:

R y   y
The log of the deviations, log(  y ) , assumed a normal distribution with mean 0 and
variance  .
2

log(  y ) ~ N (0, 2 )
We assumed equal sex ratios at recruitment, apportioning 50% of the recruitment each
year to the model starting age of each sex.

N 2, y ,s  0.5R y
where N 2, y ,s denotes the number of cisco age 2 (model starting age) in year y of sex s .
Note that this model does not assume any prior relationship between the magnitude of

26

recruitment and stock size. We modeled recruitment in the SCSA model with an added
length based component p j , representing the proportion of recruits going into each size
class.
R j , y  Ry p j

where R y is modeled in the same fashion as in the SCAA. We calculated the proportion
of recruits going into each length class, p j , using a weighted average of the size
transition probabilities of size bins smaller than the model starting bin (i.e. what
percentage of fish from bins 15-165 mm would grow into bins > 170 mm in the next time
step). To calculate the weights, we used the mean and variance in length at age of cisco
ages 0-2 to generate a sample distribution of the size distribution of cisco less than 170
mm (model starting length). To account for depletion in abundance, we weighted
different ages using a natural mortality value of 0.3 (i.e. 1 age 0, 0.74 age 1, and 0.55 age
2). These composition weights were then multiplied by the size bin (15 mm-165 mm)
transition probabilities into bins greater than 170 mm. Once again we assumed equal
sex ratios at recruitment, apportioning 50% of the recruitment each year to the model
length bins of each sex, R j , y ,s  0.5 * R j , y .
We calculated instantaneous fishing mortality for both models using
Fj , y ,s  s j f y ,s

where s j is the fishery selectivity for cisco over ages or size bins and f y ,s is the fishing
intensity in a given y year for sex s . During preliminary analysis and in previous work
done in Lake Superior (Clark 2012), cisco gill net catch per unit effort (CPUE) displayed

27

a non-linear relationship with estimates of abundance from hydroacoustic surveys, and
over most of the observed range in hydroacoustic densities there was no clear
relationship between gill net catch and hydroacoustic density. Due to this, fishery effort
was not used when fitting the model. Instead we directly estimated fishing intensity in
each year.
We modeled fishery selectivity using a two parameter gamma function as in
Deriso et al. (1985)

j  e  j
sj 
sref
where  and  are gear selectivity coefficients and the denominator denotes the value
that would be obtained for the numerator for a reference category, made age 7 and size
bin 380-390 mm for each respective model. We initially estimated fishery selectivity
parameters independently for each sex however, given near identical estimates, in the
interest of parsimony decided to assume equal fishery selectivity for each sex.
Observation Model
Predictions of data source 1, the hydroacoustic estimates of the number of cisco
greater than 250 mm, Ĥ y , were modeled using

SCAA

Hˆ y  e  P( Fish j  250mm) N j , y ,s e
s

SCSA

j

Hˆ y  e 
s

L

N

j 250

28

j , y ,s

e

Z j , y ,s

5
6

Z j , y ,s

5
6

where P( Fish j  250mm) is the probability that a cisco in age bin

j is greater than 250

mm, L represents the terminal size bin, and  is the logarithm of the hydroacoustic


survey calibration coefficient (Hulson et al., 2008), which when presented as e , can be
referred to as hydroacoustic catchability. Given that selectivity of the hydroacoustic
survey is assumed to be knife edged at 250 mm, where all fish become fully selected to
the gear, hydroacoustic estimates of spawning stock are in theory absolute estimates of
spawning stock, so this coefficient was expected to be at or very near 0. We applied
mortality to numbers of fish at age or length for the first 10/12ths of the year (

N j , y ,s e

Z j , y ,s

5
6

) given the hydroacoustic surveys are performed during the spawning

season in November. The probability that a fish of a given age is greater than 250 mm,
P( Fish j  250mm) , was calculated outside of the model by characterizing the size

distribution of each age of cisco. This was done by first obtaining the mean and variance
in length of cisco at each age from 0-27 (or model ages, i.e. 2-15), then using the
cumulative distribution function of the normal distribution to calculate the probability
that a cisco age 𝑎 chosen at random would be greater than 250 mm. To reduce gear bias
we used a multiple gear approach similar to Wilson et al., (2015), combining data from
multi-mesh gillnets, mid-water trawls, and commercial gillnets from the assessment
model years 1999-2015. Data from both sexes was also combined based on the
observation that cisco length at age in Thunder Bay did not appear to differ by sex. To
try and reduce bias in mean length at age associated with fixed bin allocation ALK
sampling (i.e. 10 fish aged per 10 mm length bin, Quinn and Deriso 1999), we calculated
mean length at age by:

29

 (n * L )

n
i, j

Lj

i

i

i, j

i

Where L j is the mean length at age 𝑗, ni , j is the number of fish in length bin 𝑖 that are
age 𝑗, and Li is the midpoint length of length bin 𝑖. Variance in length at age was
calculated using

Var ( L j )  

ni , j * ( Li  L j )2

n

i

i, j

i

The bias-adjusted mean length at each age and variance in length at each age was then
plugged into the cumulative distribution function of the normal distribution to obtain
the probability that a cisco age 𝑗 selected at random would be greater than 250 mm. We
then fit a logistic function through these probabilities to obtain the final probability that
a cisco age 𝑗 would be greater than 250 mm, P( Fish j  250mm) .
Predictions of data sources 2 and 3, the age or size composition of Thunder Bay
cisco each year obtained from mid-water trawls and multi-mesh gillnets, was modeled
using

Pˆji, y ,s 

i
j

s N j , y ,s e

 s N
i
j

s

where

Z j , y ,s

j , y ,s

e

5
6

Z j , y ,s

5
6

j

s ij is the survey selectivity of gear i (midwater trawls or multi-mesh gillnets) for

age or size

j . Pˆji, y ,s is the sex specific predicted proportion at age or size caught from
30

each survey gear in a given year. Once again we applied mortality to numbers of fish at

age or length for the first 10/12ths of the year ( N j , y ,s e

Z j , y ,s

5
6

) given the survey bio-data

are collected during the spawning season in November. Survey selectivity was modeled
to adequately characterize the selective nature of mid-water trawl surveys and multimesh gillnets using a two parameter gamma function identical to the formula described
for the fishery, however without a denominator to standardize the selectivity, given we
were calculating relative values.
Predictions of data source 4, the age or size composition of the fishery catch, was
modeled using the Baranov catch equation

Cˆ j , y ,s 

F j , y ,s

N j , y ,s (1  e

F j , y ,s  M s

Pˆj f, y ,s 

Cˆ j , y ,s
 Cˆ
s

( F j , y , s  M s )

)

j , y ,s

j

where Cˆ j , y ,s is the predicted number of commercially caught cisco age or size
y of sex s , and

j in year

Pˆj f, y ,s is the predicted sex specific age or size composition of the fishery

catch each year.
Predictions of data sources 5 and 6, the yield of each sex in each year, was
modeled using

Yˆy ,s   Cˆ j , y ,s w j ,s
j

31

where Yˆy ,s is the predicted sex specific fishery yield each year and w j ,s is the mean
weight of a commercially caught cisco age or size

j . For the SCSA this term was

obtained from the weight-length regression previously described. For the SCAA, to
account for ALK fixed allocation bin sampling bias in mean weight-at-age (Quinn and
Deriso 1999), we calculated the adjusted average weight of a commercially caught cisco
similar to P( Fish j  250mm) , using

 (n *W )

n
i , j ,s

W j ,s

i

i

i , j ,s

i

Where W j ,s represents the mean weight of a commercially caught cisco age
ni , j ,s is the number of fish in bin i that are age

j of sex s ,

j of sex s , and Wi is the average of the

length bin endpoints converted to weight using the same weight-length relationship
described previously. It should be noted we used 10 mm bins, the same bin sizes used in
the ALK sampling procedure. Variance of the mean weight of a commercially caught
cisco age

j of sex s , Var (W j ,s ) , was then modeled using
Var (W j ,s )  
i

ni , j ,s * (Wi  W j ,s )2
(  ni , j ,s )2
i

A von Bertalanffy growth function (in weight) was then fit using the mean weights of
commercially caught cisco, W j ,s , weighted by the inverse of the variance,

1
, to
Var (W j ,s )

obtain w j ,s , the adjusted average weight of a commercially caught cisco age

32

j of sex s .

w j ,s  W,s (1  e

(  K s ( j t0 , s ) bs

)

The growth function used the beta parameter, bs , from the weight-length relationship
previously described.
Predictions of data source 7, the individual cisco growth increments, which were
only used in the SCSA, were calculated using
E ( l )  ( L  l )(1  e K )

Var (l )  a  b * E (l )
where l represents the starting length of a cisco (length at start of annulus) and E ( l )
denotes the expected growth increment of a cisco given starting length.
Aging Error
Aging error was included in the SCAA model by multiplying the model predicted
catch at age (true catch at age) and the predicted relative catch at age from survey gear
by an aging error matrix. The aging error matrix was estimated by characterizing the
expected coded age ( C j ) given true age ( T j ) and the coefficient of variation of coded age
given true age as linear functions.
(C j )  a  bT j
CV (C j )  c  dT j

Where a ,

b

,

c , and

d

are estimated parameters. For ease of computation, given

preliminary fits suggested a and d were ≈ 0, further runs of the model assumed both

33

parameters fixed at zero. The aging error matrix was then computed using the
cumulative distribution function of the lognormal distribution, integrating over coded
ages j-0.5 to j+0.5. The plus group was calculated by integrating from j-0.5 to infinity.
The true catch at age matrix, output by the model, was multiplied by the aging error
matrix to obtain the predicted catch at age matrix used in calculations of the predicted
age composition of the catch. The same was done for the relative catch at age from the
survey, which was multiplied by the aging error matrix prior to calculating predicted
survey compositions.
Likelihood
We calculated the log likelihood components, Li , for data sources 1, 5, and 6 in
each model through a lognormal likelihood

Li  

1
2 i

2

 log( x

)  log( xˆ i , y )   n log(  i )
2

i, y

y

where  i is the externally specified standard deviation for likelihood component i , xi , y
and xˆi , y are the observed and predicted values for year y , and n is the number of
observations.
The log likelihood components for data sources 2, 3, and 4 for the SCAA assumed
a robust multivariate normal likelihood equation as in Starr (1999)


Li   0.5 log (1  pi , j , y ) pi , j , y

y
j







2


ˆ

(
p

p
)
0.1 
i, j, y
i, j, y




log
exp

0
.
01





0
.
1
Nb 
 ~ 
 2 (1  p ) p
 Ni, y

i, j, y
i, j, y 

 

Nb 



34

~
where N
denotes the effective sample size from data source i in year y , pi , j , y and
i, y

pˆ i , j , y are the observed and predicted proportions of cisco in year y that are age

j from

data source i , and Nb represents the number of age bins. Robust likelihoods aid in
keeping a small number of outlier composition data points from unduly influencing
model fit (Fournier et al. 1990, Francis 2011). This is especially important given the
nature of cisco year classes in Thunder Bay, which exhibit a “boom-or-bust” pattern
where there may be a very large cohort moving through the time series followed by
many years of almost no recruitment (Figure 1.3). For the SCSA, given we expect less
outlier composition data points as disparity in year class strength is reduced due to
recruitment into size bins (Figure 1.4), we utilized a regular multinomial likelihood

~
Li   N i , y  pi , j , y log( pˆ i , j , y )
y

j

where pi , j , y and pˆ i , j , y are the observed and predicted proportions of cisco in year y that
are in size bin j from data source i . Data sources i for these likelihoods include each
composition dataset (age or size; fishery catch, MWT survey, and multi-mesh gillnet
survey). It should be noted that both sexes went into one likelihood for each i , meaning
~
only one value of N
was used for each data source. This results in double the number
i, y

of size bins ( Nb * 2 ) for each i , to account for both males and females.

35

Figure 1.3. Age composition of female cisco caught in Thunder Bay, Ontario, from 1999
to 2015.

36

Figure 1.4. Length composition of female cisco caught in Thunder Bay, Ontario, from
1999 to 2015. Y axis labels denote length bin endpoints, i.e. “[150,160)” refers to a length
bin starting at 150 mm and ending at 160 mm.

37

Data source 7, the individual growth increment data, which was only used in the
SCSA, assumed a gamma log likelihood

L7  l log( l )  log( (l ))  (l  1) log( l )  l l
l

where and  l and  l are the shape and rate parameters of the gamma distribution, and

 l denotes an observed growth increment of a cisco with starting length l . The expected
value is given by E (  l ) 

l

and similarly the variance is defined as Var (  l )  l2 .
l
l

Log prior components for natural mortality and recruitment deviations that were
not compared to data but rather to expectations specified as informative priors also
contributed to the objective function through a normal prior distribution

Li  

1
2 i

2

 (log( xˆ

i, j

))2  n log(  i )

j

where xˆ i , j represents the model predicted deviations.
The objective function was then the negative sum of the log likelihood and log
prior components

L   Li
i

Data Weighting
Standard deviations,  i , in likelihood equations for data sources 1, 5 and 6, and
for recruitment deviations, were modeled as one estimated parameter and two assumed
variance ratios, denoting what we might expect the multiplicative difference in standard
38

deviations to be. For example, we expected the yield (data sources 5 and 6, sharing a  )
to have the smallest standard deviation, since catch is usually assumed known, followed
by hydroacoustic abundance estimates, and finally recruitment deviations, which we
expected to vary considerably given the drastic difference in cohort strength between
years (Figure 1.3). We estimated the standard deviation for harvest, while assuming
variance ratios ( Vr ) of 0.04 and 0.0004 for hydroacoustic estimates of abundance and
recruitment deviations, respectively.

i 

1 2
f
Vr

where  f denotes the standard deviation for fishery harvest.
~
Effective sample sizes, N
, for the composition datasets were calculated using
i, y

iterative reweighting procedure T3.4 of Francis (2011)
~
N i , y  N i , y wi

where N i , y denotes the previous iterations effective sample size. wi was calculated using
TA1.8 of Francis (2011)
wi 

where Oi , y 

 j *O

i, j, y

1
 (Oi , y  Ei , y ) 
Var 
0.5 
 ( vi , y N i , y ) 

with similar formula for Ei , y , and vi , y 

j

( j E
2

i, j, y

)  Ei2,y . i in

j

these formulas denotes a composition data set (age or size; commercial fishery, MWT
survey, or multi-mesh gillnet survey) and

j denotes age or size bin. The initial effective

39

sample sizes were a maximum of 200 cisco for the commercial fishery, and a maximum
of 100 cisco for each of the survey gears. The iterative process was stopped when the
effective sample sizes appeared to converge. This process was robust to starting sample
~
sizes, with the N
converging on similar estimates regardless of starting values. We
i, y

utilized maximum likelihood estimation to obtain effective sample size convergence
prior to running the models using Bayesian methods in Automatic Differentiation Model
Builder (ADMB). Effective sample sizes converged on 62 for fishery compositions, 45 for
MWT compositions, and 50 for multi-mesh gillnet compositions for the SCAA and 58,
22, and 11 for the SCSA.
Model Fitting/Calibration/Troubleshooting
SCAA
The SCAA model was unable to converge on an estimate of  , which denotes the
logarithm of acoustic catchability. Essentially this parameter scales our entire
population by representing what proportion of spawning cisco the acoustic survey is
actually detecting. Due to this, we decided to assume a conservative scenario where

  0 , which assumes the acoustic survey is an absolute index of spawner abundance.
Here, by conservative, we mean that actual catchability is likely lower, abundance is
likely higher, and thus yield targets set based on the abundance estimates will likely be
lower than if the target exploitation rate were applied to an unbiased abundance
estimate. The hydroacoustic surveys are generally thought to be a conservative estimate
of abundance as all areas of the water column are not detectable with the gear, i.e.
surveys are likely missing some fish (Yule et al., 2012). The model was run for 10 million

40

iterations each saving every 500th, burning in the starting 2500 values from the final
chains. Chain burn in was assessed visually and convergence determined using
Geweke’s convergence diagnostic (Geweke, 1991).
SCSA
It became clear at the beginning of model calibration for the SCSA that the model
was going to be unable to output plausible estimates of natural mortality. The model
would confound estimates of recruitment, selectivity, and natural mortality. It was
unable to converge on plausible estimates of natural mortality even when given assumed
known growth parameters at levels previously estimated using fixed natural mortality at
prior point estimates. The model would increase natural mortality to an implausibly
high value, inflate recruitment, and make larger fish more selected. What the model was
doing was creating many fish through recruitment, killing them off at high rates through
natural mortality in order to have enough fish at spawning sizes to fit the hydroacoustic
data. Few large fish were predicted to survive, but fishery selectivity was highest for the
largest fish in order to fit the fishery composition data. Similar to the SCAA, the SCSA
was also unable to converge on an estimate of acoustic catchability. Given these issues,
we decided to fix natural mortality at its prior point estimates (0.283 for males, 0.256
for females), fix acoustic catchability at 1 (   0 ), and estimate growth. We ran the
model for double the number of samples (20 million) and also doubled chain burn in
(first 5000 iterations).

41

Comparison
Given different data used in each assessment, it was not possible to compare the
final models in terms of predictive accuracy/information theoretic measures such as
PSIS-LOO, WAIC, or DIC. Instead, final models were compared using a variety of
criteria. First, we considered what assumptions we had to make to fit each model. We
also looked at retrospective patterns, parameter/output uncertainty and model
fit/residuals. Retrospective analyses primarily focused on spawning biomass and
exploitation rate. Mohns rho (Mohn, 1999) was calculated for spawning biomass and
exploitation rate as the mean relative error for the last year of each peel compared to the
corresponding year in the full assessment.

 YF , p  YF ,ref 


Y
F ,ref



  

Where Y is the assessment output quantity, either spawning biomass or exploitation
rate, ref refers to the full assessment, and F refers to the final year of a given
assessment peel, p . Five years were removed from the assessment. We also calculated a
mean final year absolute bias for the retrospective analyses, as the mean absolute value
of the relative error for the last year of each peel compared to the corresponding year in
the full assessment. This statistic considers the bias in the final year of each peel as
opposed to whether or not there is a consistent pattern.
 YF , p  YF ,ref

YF ,ref


 

42






Both  and  were calculated using medians of the posterior distribution as point
estimates.
Residuals for common data sources were compared both visually in addition to
using the standard deviation of the normalized residuals (SDNR, Breen et al., 2003;
Francis, 2011; Carvalho et al., 2016). These were calculated as the standard deviation of
the normalized residual for each data point (formulas in Table B1 in Francis 2011). A
relatively good model fit is characterized by smaller residuals and a SDNR near 1
(Carvalho et al., 2016), although Francis (2011) notes that a value much less than 1 is not
a cause for concern but rather means that the data set is fitted better than was expected.
Due to their correlative nature, composition data points cannot be compared using this
metric (Francis, 2011), so these were compared visually.
Results
Point estimates of quantities output from the models are reported as medians of
the posterior distribution, with 95% highest posterior density (HPD) intervals reported
in parentheses. All parameters in each model indicated convergence based on Geweke’s
diagnostic at an alpha level of 0.01. The SCAA estimated a total of 89 parameters while
the SCSA estimated 67 parameters. Natural mortality estimates within the SCAA for
males and females were 0.284 (0.214-0.367) and 0.252 (0.182-0.336) respectively.
Biomass
SCAA biomass fluctuated throughout the time series beginning at 5.92 (2.5012.48) million kg in 1999, reaching a peak in 2006 of 8.41 (4.94-13.10) million kg, and
ending at 1.70 (1.07-2.47) million kg in 2015 (Figure 1.5). Spawning biomass, defined as

43

the mature female biomass (>250 mm), began at 4.96 (2.10-10.29) million kg, initially
declined then rose to an estimate of 4.88 (2.91-7.65) million kg in 2006 and ended the
time series at 0.96 (0.59-1.40) million kg.
SCSA biomass exhibited a similar trend to the SCAA, with an increase during the
beginning of the time series from 2.62 (1.58-4.54) million kg in 1999 to a peak in 2007
of 5.97 (4.35-7.84) million kg, followed by a decrease to a final year estimate of 2.85
(2.01-3.97) million kg in 2015. Spawning biomass began at 2.04 (1.30-3.48) million kg
and ended at 1.44 (0.98-2.05) million kg, with a peak in 2008 of 3.01 (2.17-4.01) million
kg.
Exploitation Rate
Exploitation in the SCAA was modest throughout the time series, hovering
around 3%, although in 2010 began to increase resulting in a final year estimate of 9%
(6%-14%, Figure 1.5). This resulted in fully selected fishing mortality rate estimates of
0.08 (0.04-0.12) and 0.20 (0.11-0.32) for males and females in 2015, respectively.
Exploitation rate was defined as yield divided by the biomass of fish larger than 250
mm. This is similar to the control rule allocation in Thunder Bay, where total allowable
catch is calculated as 10% of cisco biomass greater than 250 mm estimated from the
hydroacoustic survey. For the SCSA, exploitation rate decreased from 8% (4%-13%) at
the beginning of the time series to 3% (2%-4%) in 2007 and increased throughout the
rest of the time series to a final year estimate of 6% (4%-8%). Final year fully selected
fishing mortality rates were 0.05 (0.03-0.07) and 0.17 (0.10-0.25) for males and
females.

44

Figure 1.5. Biomass in millions of kg, spawning biomass in millions of kg of mature
females (>250 mm), exploitation rate (yield/biomass of fish > 250 mm), and abundance
in millions of fish for the SCAA and the SCSA. Shaded regions denote 95% HPD
intervals and dashed lines are medians of the posterior distribution. Squares denote
SCAA output while triangles denote SCSA output. Light shading denotes the HPD for
the SCAA, darker shading denotes the HPD for the SCSA, and the darkest shading is
where the two intervals overlap.

45

Recruitment
As expected, recruitment was highly variable throughout the time series, with
evidence of ~4 “boom” recruitment years in the SCAA, belonging to 1998, 2003, 2005,
and 2009 year classes (Figure 1.6). Estimates of recruitment (age-2 fish) for these years
(y+2; 2000, 2005, 2007, and 2011) were 19.19 (7.46-39.43), 36.69 (18.68-61.22), 1.86
(0.91-3.22), and 4.38 (2.30-7.11) million fish, respectively. Recruitment was estimated
to be around 15,000 fish in 10 years, and the 3 remaining years were estimated at
modestly low values, with estimates ranging from 0.85 (0 -1.78) million in 2004 to 0.21
(0-0.98) million fish in 1999.
Recruitment in the SCSA, where not defined in the same way as the SCAA, in that
it describes the number of fish greater than 170 mm that recruit into the model, showed
a similar trend to SCAA recruitment with 3-4 clear modes most likely attributed to the
1998, 2003, 2005, and 2009 “boom” year classes (Figure 1.6).

46

Figure 1.6. Recruitment for the SCAA and SCSA. Where both represent the number of
fish entering the model in a given year, they are not defined the same way in that for the
SCAA it is the number of age 2s entering the population and in the SCSA it is the
number of fish greater that 170 mm that are entering the population. Points denote
medians of the posterior distribution and error bars are the 95% HPD intervals.

47

Abundance
SCAA abundance echoed biomass results, with intermittent spikes due to “boom”
recruitment years and an overall declining trend over time. In 1999 the model predicted
there were around 12.62 (5.24-26.63) million cisco, a high of 44.07 (23.68-73.08)
million estimated in 2005, and in the final year 3.90 (2.46-5.69) million (Figure 1.5).
SCSA predicted abundance began the time series at 7.57 (3.71-13.89) million fish,
31.03 (21.35-41.31) million fish at its peak in 2005, followed by a decrease to around
5.61 (3.63-8.47) million fish in 2015.
Aging Error
Very little aging error was estimated within the SCAA model. Approximately no
bias was estimated in aging as true age increased ( b ≈ 1.00), in addition to a very low
estimated CV ( c = 0.02).
Growth

L and K were estimated at 428 mm (419-439) and 0.28 (0.25-0.30),
respectively (Figure 1.7).

48

Figure 1.7. Upper left: Fit to growth increment data. Red line depicts the median of the
posterior distribution and shaded region represents 95% HPD intervals.

49

Figure 1.7 (cont’d) Upper right: Residuals, medians of the posterior distribution, from fit
to growth increment data. Lower panel: Growth transition matrix at the posterior
medians for growth parameters. Note that the area of the circles represent the
probability of growing into a length bin given a starting length bin. Length bins are
represented on axes as midpoints. Plus group is length bin 410-420 mm.

Retrospective Analyses
Retrospective patterns for each model were very similar (Figure 1.8). Mohn’s rho
estimates for spawning biomass and exploitation rate were -0.01 and 0.01 for the SCAA
and 0.18 and -0.13 for the SCSA, respectively. All of these rho values are within a range
of values deemed “not a cause for concern” in retrospective analyses (Hurtado-Ferro et.
al, 2015). Absolute bias estimates for spawning biomass and exploitation rate were 0.14
and 0.16 for the SCAA and 0.21 and 0.15 for the SCSA, respectively.

50

Figure 1.8. Retrospective analyses for spawning biomass and exploitation rate. Shown
are times series estimates of spawning biomass and exploitation rate when five terminal
years of data are sequentially dropped from each assessment. Spawning biomass in
millions of kg of mature females and exploitation rate as yield/biomass of fish > 250
mm.

51

Model Fit to Data
Assessment model fits to the hydroacoustic data were very similar (Figure 1.9).
Both assessments treated the observed hydroacoustic spawning abundance estimate in
2011 as an outlier. Outside of that outlier data point, both models predicted a near linear
decline in spawning abundance since 2005, in accordance with the observed data
points. The median of the negative log-likelihood for the fit to acoustic data was lower
for the SCAA (Table 1.4). SDNR values for the acoustic data were also closer to 1 for the
SCAA, indicating better fit.
Fit to harvest data was nearly identical between the assessments (Figure 1.10). In
fact both models fit the data so closely that the observed data is no longer visible on
Figure 1.10. HPD credible intervals were slightly smaller for the SCSA model. The
medians of the negative log-likelihoods for male and female harvest were lower for the
SCSA (Table 1.4), and the SDNR values were smaller for the SCSA. Both model SDNR
yield values were well below 1, indicating better model fit than expected.
Both models fit the fishery composition data points well (see Supplemental
Figures 1.1-1.8). These two fits cannot be directly compared as they used different data.
Table 1.4. Negative log-likelihood and SDNR values for common data sources.

SCAA – NLL
SCSA – NLL
SCAA – SDNR
SCSA – SDNR

Male
Yield
-32.97
-35.30
0.09
0.06

Female
Yield
-32.94
-35.32
0.12
0.04

Hydroacoustic
Data
0.31
1.94
1.26
1.50

52

Figure 1.9. Fit to hydroacoustic estimates of spawning abundance. No data from 19992004 and in 2006. Spawning fish in millions of fish. Points denote medians of the
posterior distribution and error bars are 95% HPD intervals.

53

Figure 1.10. Fit to landings data. Harvest (first column) is in 10,000 kg. Note observed
data points in the first column cannot be seen as they are covered by model point
estimates.

54

Computational Intensity
The SCSA was considerably more computationally intensive than the SCAA,
necessitating about 3x the run time for the same number of iterations (where the SCAA
took ~5 hours for 10 million, SCSA took ~15 hours).
Discussion
Overall both models showed similar trends in outputs and modest differences in
final year estimates (Table 1.5, Figure 1.5). However, the SCAA model showed a larger
degree of uncertainty during the first half of the time series, which decreased
throughout the second half of the time series to actually end up being less than SCSA
uncertainty in the final year (Figure 1.5). This disparity in uncertainty at the beginning
of the time series is likely due to differences in the initial parameterization of each
model, where much more flexibility was afforded to the SCAA by estimating 26 initial
abundance parameters (1 for each age-sex combination above recruitment age).
Conversely the SCSA estimated only 4 initial abundance parameters; 2 for the initial size
composition, and 1 for each sex as abundance multipliers. Increased certainty in the
SCSA during most of the time series may also be driven by the assumption of known
natural mortality values. In fact, when we re-ran the SCAA with assumed known natural
mortalities at their prior point estimates, uncertainty in model output decreased
significantly (Figure 1.11), indicating that the certainty in output expressed by the SCSA
is likely in some part due to assuming known natural mortality values. In terms of
uncertainty in the final year of each assessment, one aspect that might be driving the
smaller credible interval for abundance outputs (biomass, spawning biomass,
abundance) in the final year for the SCAA is lack of a recent “boom” year class. The last

55

large recruitment year was in 2009 (Figure 1.6 left panel, and by boom standards it was
pretty weak), so effectively, without replenishment the model is consistently decreasing
the uncertainty in biomass and abundance from 2009-2015 through depletion of fish.
Table 1.5. Final year (2015) relative differences [(SCAA-SCSA)/SCAA].
Output
Abundance
Biomass
Spawning Biomass
Exploitation rate

Relative
Difference
-44%
-68%
-51%
40%

Figure 1.11. Comparison of models when natural mortality is estimated within SCAA and
when it is assumed known. First column are the original models, while the second
column natural mortality is assumed known for SCAA. Shaded regions denote 95% HPD
intervals and dashed lines are medians of the posterior distribution. Squares denote
SCAA output while triangles denote SCSA output. Light shading denotes the HPD for
the SCAA, darker shading denotes the HPD for the SCSA, and the darkest shading is
where the two intervals overlap.

56

In terms of model fit, the SCAA had better fit to the acoustic data while the SCSA
had better fit to the landings data [although the slightly better fit of the SCSA to
landings may not be all that significant, as all landings SNDRs were well below 1,
indicating better fit than expected by each model (Francis 2011)]. For the retrospective
analysis, where  estimates were larger for the SCSA, upon visual inspection the
patterns appear comparable between the two assessments, if not worse for the SCAA
(Figure 1.8). This result of smaller  estimates for the SCAA even though patterns may
appear more severe if not equal to those in the SCSA is driven by equally large bias in
the terminal years of peels for the SCAA in opposite directions (i.e. not a consistent
pattern but equal numbers of over and under estimates). Given no  estimates are at
values considered “cause for concern” (Hurtado-Ferro et al., 2015), it may be more
prudent to consider  in comparing retrospective analyses, which was larger for
exploitation rate and smaller for spawning biomass for the SCAA compared to the SCSA.
Once again it is likely assuming known natural mortality values at their prior point
estimates led to both a smaller  estimate for exploitation rate and the appearance of
less severe retrospective patterns within the SCSA. When we re-ran the SCAA
retrospective analysis with assumed known natural mortality values at their prior point
estimates, absolute bias estimates (  ) decreased. In fact, all retrospective statistics were
lower (closer to 0 for rho) for the SCAA with assumed known natural mortality than the
SCSA counterparts (Figure 1.12; SCAA fixed M  : Spawning Biomass – 12%,
Exploitation rate – 14%).

57

Figure 1.12. Retrospective analyses where both SCAA and SCSA assumed the same
known natural mortality values. Shown are times series estimates of spawning biomass
and exploitation rate when five terminal years of data are sequentially dropped from
each assessment. Spawning biomass in millions of kg of mature females, and
exploitation rate as yield/biomass > 250 mm.

58

Perhaps the most important result of our study is the inability to estimate natural
mortality within the SCSA. Given natural mortality is one of the most influential
quantities in stock assessment and its estimation within an assessment can be difficult
(Lee et al., 2011; Brodziak at al., 2011; Sippel et al., 2017), its estimability in the SCAA
certainly favors the SCAA as an assessment model choice. An interesting note is the
remarkable similarity of the prior natural mortality point estimates (the fixed, assumed
known SCSA M values; 0.283 and 0.256) to the estimated natural mortality point
estimates for the SCAA. In fact, when we ran the SCAA model without specifying
informative priors on natural mortality, a similar result occurred (Male M = 0.282,
Female M = 0.250), indicating this similarity is not due to the specification of an
informative prior in the SCAA but rather that the age composition data are providing
crucial information on natural mortality. The similarity between assumed known
natural mortality in the SCSA and estimated natural mortality in the SCAA in addition
to the utilization of the same hydroacoustic and landings data likely led to similar output
between the two assessments.
The inability to estimate natural mortality within the SCSA due to its
confounding with estimates of recruitment and selectivity is not a new finding, as
parameter confounding has been noted to be potentially more serious in size-structured
assessments (Punt et al., 2013). Where parameter confounding did not change growth
parameters much (mainly influenced selectivity, recruitment, and natural mortality), its
underlying cause may have been variation in growth, such that variation in size-at-age
makes it hard for size-structured models to discern cohorts from length composition
data (Punt et al., 2013). In fact, even if we assumed growth was known within the SCSA,

59

the model still confounded selectivity, recruitment, and natural mortality. Another
aspect that may have led to the inability to estimate natural mortality within the SCSA is
the range of vulnerability to the fishery for cisco in Thunder Bay, where by the time cisco
start to show up in the fishery length compositions they are at or very near asymptotic
size, by this time having substantially slowed their somatic growth (Figure 1.4). This
results in similarity in fishery size composition data between years making it difficult to
observe strong year classes pulse through the fishery composition data (Figure 1.4).
Where the fishery independent survey gear does select smaller fish and is, to some
extent, able to discern cohorts from its length composition data (likely why recruitment
in SCSA for 2003 and 2009 cohorts were approximated well), our survey composition
data was limited, only having started in 2005 and missing critical years in 2006 and
2011-2012. The missing survey data pre-2005 likely resulted in recruitment of the 1998
cohort being spread over ~5 years in the SCSA (Figure 1.6). Fishery independent survey
size composition data throughout the full time series would have likely resulted in a
better approximation of year-class strength and possibly allowed estimation of natural
mortality within the SCSA. Alternatively in the SCAA model, likely due to the boom-orbust recruitment pattern, the model was clearly able to distinguish 3-4 large year classes
moving through the fishery and estimate their associated depletion. Estimation of
natural mortality within the SCAA may also have been made possible by relatively light
exploitation, effectively making the major source of mortality and transition through the
population-age matrix one of natural depletion.
Estimation of natural mortality within size structured assessments is possible, as
within Punt et al. (2013), a review of integrated size-structured assessment methods,

60

three out of nine assessments that were reviewed in depth estimated natural mortality.
Two of these three assessments modeled selectivity as logistic, the other modeled it as a
double normal, while all three modeled recruitment as lognormal deviates entering the
population through a specified size distribution. These selectivity functions are less
flexible than a gamma function, which may indicate a reason they did not experience
parameter confounding to the extent we did with regard to natural mortality, selectivity
and recruitment. Although interesting to note that in our study even if we fixed
selectivity and growth at values estimated using assumed known natural mortalities,
and then tried to estimate natural mortality, the SCSA model would still inflate
recruitment and estimate implausibly high natural mortalities. Where some sizestructured models may indeed be able to estimate natural mortality, our study indicates
that this may be an even taller task than it is in SCAA models, and is dependent on a
multitude of factors from variability/patterns in recruitment, variability in growth and
size at age, and vulnerability range of organisms within size composition data.
Conclusions and Recommendations
Although more uncertain, primarily due to the estimability of natural mortality,
we conclude that the SCAA is more appropriate for modeling population dynamics of
cisco in the Laurentian Great Lakes. Where size based assessment models can
considerably decrease the amount of fish that need to be aged, as this study shows, age
composition data can be crucial to the ability to estimate natural mortality within a
model. We prefer to avoid reliance on assumed known scale parameters (acoustic
catchability) and natural mortalities, if this can be avoided. Where the assessments both
resulted in similar natural mortality estimates, in other case studies this may not occur,

61

and bias in using a surrogate equation for natural mortality may result in biased
assessment output. In addition, assuming known natural mortality may artificially
decrease model uncertainty. We do not necessarily expect this conclusion to apply for
all, or even most species. In fact, this result is likely largely driven both by the unique
life history of cisco in exhibiting boom-or-bust recruitment, and the fact that most
growth occurs before cisco are vulnerable to the fishery, allowing for the estimation of
natural mortality within the SCAA and not within the SCSA. For species with less
variable recruitment, less variable growth, and size composition data throughout the
growth period of their life span, size based assessment methods may perform equally
well, or better, than age-structured methods.
Our conclusion, that the SCAA was more appropriate than the SCSA when
applied to cisco, is largely driven by our desire to estimate natural mortality. While it is
tempting to contrast our conclusion with other comparisons of size- and age-based
assessment models (Akselrud et al., 2017; Punt et al., 2017), those studies did not
attempt to estimate natural mortality within the assessment models. Where Punt et al.,
(2017) concluded that age-structured methods performed poorest and Akselrud et al.,
(2017) concluded that age-structured methods fit the data best, we believe that the
conclusions of these studies might depend on their assumption that natural mortality
was known. Akselrud et al., (2017) and Punt et al., (2017) also considered assessment
models that take into account both age- and size-based processes in their analyses.
Where these models may improve assessment accuracy (Gilbert et al., 2006; McGarvey
et al., 2007; Punt et al., 2017), they are also very data- and computationally intensive.
We did not consider them in our analyses. It is possible that an age- and size-structured

62

model could outperform both SCAA and SCSA in application to cisco in Thunder Bay.
Additionally, while we believe the comparisons we made and conclusion we reached in
preferring the age-based model is valid, we cannot be sure that the estimated population
sizes and mortality rates are closer to true values than those generated by a size-based
model, given the truth is not known. Further, our analysis cannot define the conditions
under which the natural mortality is estimable and produces useable assessment results,
as we had only one data set resulting from one set of conditions. This is a potential
advantage of simulations like those of Punt et al. (2017) over empirical comparisons of
alternative models as shown here. Our empirical comparisons highlighted some aspects
of the performance of size- and age-based models contrast in real world applications
and thus can point the way for future simulations. More work is needed that directly
investigates the estimability of natural mortality and catchability within size-structured
assessment models both from a simulation perspective and in real world assessments.

63

APPENDIX

64

Excerpt from Harding (2017) on preparation, aging, and measurement of otolith
radii: “Transverse sections of the sagittal otoliths were taken as per Schreiner and
Schram (2001), and otolith sections were briefly etched in 1% acetic acid solution for 15
minutes to improve the visibility of annuli. Fish were aged using transmitted light at
10X magnification on a Nikon SMZ 1500 stereoscope. All fish were aged twice by the
same reader with greater than two weeks between readings. Fish with discrepancies
between the first and second readings were excluded from further analyses…. Digital
images of otoliths were captured at 4X magnification using a Nikon SMZ 1500
stereoscope and a Nikon DXM 1200 scope mounted camera. All images were calibrated
using a stage micrometer. The origin of the otolith was identified as the centroid of the
age-one annulus, and the measurement axis was defined as a 30-degree angle from the
longest axis of the age-one annulus through the origin to the otolith margin in the
ventral-medial direction (Figure 1-2). Measurements were made using ImageJ imaging
software (version 1.48, Research Services Branch, National Institute of Health).
Measurements were made along the measurement axis from the origin to the outer edge
of the discontinuous zone.”

65

REFERENCES

66

REFERENCES

Akselrud, C. I. A., Punt, A. E., & Cronin-Fine, L. 2017. Exploring model structure
uncertainty using a general stock assessment framework: The case of Pacific cod
in the Eastern Bering Sea. Fisheries Research, 193, 104-120.
Bertignac, M., and De Pontual, H. 2007. Consequences of bias in age estimation on
assessment of the northern stock of European hake (Merluccius merluccius) and
on management advice. ICES Journal of Marine Science, 64(5), 981-988.
Breen, P.A., Kim, S.W., and Andrew, N.L. 2003. A length-based Bayesian stock
assessment model for the New Zealand abalone Haliotis iris. Mar. Freshw. Res.
54(5): 619–634.
Brodziak, J., Ianelli, J., Lorenzen, K., & Methot Jr, R. D. 2011. Estimating natural
mortality in stock assessment applications. NOAA Technical Memorandum
NMFS-F/SPO, 119, 38p.
Campana, S. E. 2001. Accuracy, precision and quality control in age determination,
including a review of the use and abuse of age validation methods. Journal of fish
biology, 59(2), 197-242.
Carvalho, F., Punt, A. E., Chang, Y. J., Maunder, M. N., & Piner, K. R. 2016. Can
diagnostic tests help identify model misspecification in integrated stock
assessments?. Fisheries Research.
Clark, R. 2012. Direct estimates of the catchability of gill nets on spawning
congregations of cisco in Lake Superior. Study final report to Chippewa-Ottawa
Resource Authority.
Coggins, L. G., & Quinn, T. J. 1998. A simulation study of the effects of aging error and
sample size on sustained yield estimates. Fishery stock assessment models, 955975.
Deriso, R. B., Quinn Ii, T. J., & Neal, P. R. 1985. Catch-age analysis with auxiliary
information. Canadian Journal of Fisheries and Aquatic Sciences, 42(4), 815824.
Ebener, M. P., Stockwell, J. D., Yule, D. L., Gorman, O. T., Hrabik, T. R., Kinnunen, R.
E., Mattes, W. P., Oyadomari, J. K., Schreiner, D. R., Geving, S., Scribner, K.,
Schram, S. T., Seider, M. J., Sitar, S. P. 2008. Status of cisco (Coregonus artedi)
in Lake Superior during 1970-2006 and management and research
considerations. Gt. Lakes Fish. Comm., Lake Superior Tech. Rep, (1), 126.

67

Fournier, D., & Archibald, C. P. 1982. A general theory for analyzing catch at age data.
Canadian Journal of Fisheries and Aquatic Sciences, 39(8), 1195-1207.
Fournier, D.A., Sibert, J.R., Majkowski, J., and Hampton, J. 1990. MULTIFAN a
likelihood-based method for estimating growth parameters and age composition
from multiple length frequency data sets illustrated using data for southern
bluefin tuna (Thunnus maccoyii). Can. J. Fish. Aquat. Sci. 47(2): 301–317.
Fournier, D. A., Hampton, J., & Sibert, J. R. 1998. MULTIFAN-CL: a length-based, agestructured model for fisheries stock assessment, with application to South Pacific
albacore, Thunnus alalunga. Can. J. Fish. Aquat. Sci. 55(9), 2105-2116.
Francis, R. 1990. Back-calculation of fish length: a critical review. Journal of Fish
Biology 36:883–902.
Francis, R.I.C.C., 2011. Data weighting in statistical fisheries stock assessment models.
Can. J. Fish. Aquat. Sci. 68, 1124–1138.
Geweke, J. 1991. Evaluating the accuracy of sampling-based approaches to calculating
posterior moments. In Bayesian Statistics 4 (ed JM Bernado, JO Berger, AP
Dawid and AFM Smith). Clarendon Press, Oxford, UK
Gilbert, D. J., Davies, N. M., & McKenzie, J. R. 2006. Development of an age–length
structured model of the Hauraki Gulf–Bay of Plenty snapper (Pagrus auratus)
population. Marine and freshwater research, 57(5), 553-568.
Harding, I. C. 2017. Growth and Predatory Demand of Cisco (Coregonus artedi) in
Western Lake Superior (Masters Thesis, University of Minnesota).
Hulson, P. J. F., Miller, S. E., Quinn, T. J., Marty, G. D., Moffitt, S. D., & Funk, F. 2007.
Data conflicts in fishery models: incorporating hydroacoustic data into the Prince
William Sound Pacific herring assessment model. ICES Journal of Marine
Science, 65(1), 25-43.
Hurtado-Ferro, F., Szuwalski, C. S., Valero, J. L., Anderson, S. C., Cunningham, C. J.,
Johnson, K. F., Licandeo, R., McGilliard, C. R., Monnahan, C. C., Muradian, M.
L., Ono, K., Vert-Pre, K. A., Whitten, A. R., Punt, A. E. 2014. Looking in the rearview mirror: bias and retrospective patterns in integrated, age-structured stock
assessment models. ICES Journal of Marine Science, 72(1), 99-110.
Lee, H. H., Maunder, M. N., Piner, K. R., & Methot, R. D. 2011. Estimating natural
mortality within a fisheries stock assessment model: an evaluation using
simulation analysis based on twelve stock assessments. Fisheries Research,
109(1), 89-94.

68

McGarvey, R., Feenstra, J. E., & Ye, Q. 2007. Modeling fish numbers dynamically by age
and length: partitioning cohorts into" slices". Canadian Journal of Fisheries and
Aquatic Sciences, 64(9), 1157-1173.
Methot, R. D., & Wetzel, C. R. 2013. Stock synthesis: a biological and statistical
framework for fish stock assessment and fishery management. Fisheries
Research, 142, 86-99.
Mohn, R. 1999. The retrospective problem in sequential population analysis: An
investigation using cod fishery and simulated data. ICES Journal of Marine
Science, 56(4), 473-488.
Punt, A. E., Huang, T., & Maunder, M. N. 2013. Review of integrated size-structured
models for stock assessment of hard-to-age crustacean and mollusc species. ICES
Journal of Marine Science, 70(1), 16-33.
Punt, A.E., Allen Akselrud, C., Cronin-Fine, L., 2017. The effects of applying misspecified age- and size-structured models. Fish. Res. 188, 58–73.
Quinn, T. J., & Deriso, R. B. 1999. Quantitative fish dynamics. Oxford University Press.
Reeves, S. A. 2003. A simulation study of the implications of age-reading errors for
stock assessment and management advice. ICES Journal of Marine Science,
60(2), 314-328.
Rudstam, L. G., Parker-Stetter, S. L., Sullivan, P. J., & Warner, D. M. 2009. Towards a
standard operating procedure for fishery acoustic surveys in the Laurentian Great
Lakes, North America. ICES Journal of Marine Science, 66(6), 1391-1397.
Sippel, T., Lee, H. H., Piner, K., & Teo, S. L. 2016. Searching for M: Is there more
information about natural mortality in stock assessments than we realize?
Fisheries Research.
Starr, P.J., Bentley, N., and Maunder, M.N. 1999. Assessment of the NSN and NSS
stocks of red rock lobster (Jasus edwardsii) for 1998. New Zealand Fisheries
Assessment Research Document 99/34. Ministry of Fisheries, Wellington, New
Zealand.
Sullivan, P. J., Lai, H. L., & Gallucci, V. F. 1990. A catch-at-length analysis that
incorporates a stochastic model of growth. Canadian Journal of Fisheries and
Aquatic Sciences, 47(1), 184-198.
Then, A. Y., Hoenig, J. M., Hall, N. G., & Hewitt, D. A. 2014. Evaluating the predictive
performance of empirical estimators of natural mortality rate using information
on over 200 fish species. ICES Journal of Marine Science, 72(1), 82-92.

69

TeWinkel, L. M., T. Kroeff, G. W. Fleischer, and M. Toneys. 2002. Population dynamics
of bloaters (Coregonus hoyi) in Lake Michigan, 1973–1998. Archiv fuer
Hydrobiologie Special Issues 57:307–320.
Thompson, G.G., A’Mar, Z.T., Palsson, W.A., 2011. Assessment of the Pacific cod stock
in the Gulf of Alaska. In: Plan Team for Groundfish Fisheries of the Gulf of Alaska
(compiler), Stock Assessment and Fishery Evaluation Report for the Groundfish
Resources of the Gulf of Alaska. North Pacific Fishery Management Council, 605
W. 4th Avenue Suite 306, Anchorage, AK 99501, pp.161–306.
Vehtari, A., Gelman, A., & Gabry, J. 2017. Practical Bayesian model evaluation using
leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 14131432.
Wilson, K. L., Matthias, B. G., Barbour, A. B., Ahrens, R. N., Tuten, T., & Allen, M. S.
2015. Combining samples from multiple gears helps to avoid fishy growth curves.
North American Journal of Fisheries Management, 35(6), 1121-1131.
Yule, D. L., Stockwell, J. D., Cholwek, G. A., Evrard, L. M., Schram, S., Seider, M., &
Symbal, M. 2006. Evaluation of methods to estimate lake herring spawner
abundance in Lake Superior. Transactions of the American Fisheries Society,
135(3), 680-694.
Yule, D. L., J. D. Stockwell, J. A. Black, K. I. Cullis, G. A. Cholwek, and J. T. Myers.
2008. How systematic age underestimation can impede understanding of fish
population dynamics: lessons learned from a Lake Superior cisco stock.
Transactions of the American Fisheries Society 137:481–495.
Yule, D. L., Schreiner, D. R., Addison, P. A., Seider, M. J., Evrard, L. M., Geving, S. A., &
Quinlan, H. R. 2012. Repeat surveys of spawning cisco (Coregonus artedi) in
western Lake Superior: timing, distribution and composition of spawning stocks.
Advances in Limnology, 63, 65-87.

70

Chapter 2
Evaluating the Sustainability of a Cisco Fishery in Thunder Bay,
Ontario under Alternate Harvest Control Rules
Abstract
Sustainable management of fish stocks is promoted through the application of
Management Strategy Evaluations, providing advice to managers on the relative
performance of alternative management approaches (strategies) while accounting for
uncertainty. In this study, we developed a simplified management strategy evaluation of
a cisco, Coregonus artedi, stock in Thunder Bay, Ontario, in an effort to determine both
the sustainability of the current harvest control rule, a fixed exploitation rate of 10%,
and the performance of alternative harvest control rules in meeting fishery objectives.
Our simulations explicitly accounted for uncertainty in the frequency of “boom” year
classes being produced by cisco, the shape of the stock-recruit function, stock
abundance, and the sex-specific nature of the roe harvest. Assuming future productivity
is similar to that observed over the past 30 years, results suggest the current
exploitation rate of 10% is sustainable in terms of maintaining spawning biomass above
20% of the unfished level. Alternate control rules involving biomass thresholds defining
when exploitation rate is to decrease as a function of spawning stock size increased
yield, decreased risk, and increased the magnitude of spawning biomass at the end of
the simulation period, while resulting in more inter-annual variation in yield. Constant
catch control rules greatly underperformed constant exploitation rate control rules in
terms of magnitude in yield; however constant catch control rules did reduce interannual variation in yield compared to constant exploitation rate control rules, and use of

71

conditional versions of constant catch control rules (i.e., threshold stock sizes below
which catch was reduced) mitigated risks of staying at low stock size.

72

Introduction
Rational management of fish stocks to promote sustainable and economically
viable yields requires clearly defined objectives and quantitative analyses on the effect of
alternative harvest policies in achieving said objectives. This is often facilitated through
a process known as Management Strategy Evaluation (MSE), or the evaluation of
management strategies using simulation (Punt et al. 2008). A central tenet of these
simulations is attempting to account for uncertainty in key processes, such as the
assessment process, the stock-recruit relationship, or the implementation of a harvest
control rule, as accounting for these uncertainties has been shown to affect the outcome
of evaluations (Deroba and Bence, 2008). This can be done by including several possible
scenarios within an operating model that encompass the realistic range of key
uncertainties underlying the true dynamics of the fishery (Deroba and Bence, 2012).
MSEs can allow for tailoring specific harvest control rules to meet given fishery
objectives. Alternatively, due to limited information or analytical capacity, many
fisheries are managed through the calculation of biological reference points (Goodyear
1993) used in defining targets or limits (Quinn and Deriso 1999; Caddy and Mahon
1995) based on generalizable rules that have been proposed and applied across fisheries
with different life histories and harvest dynamics (i.e., fishing mortality should be lower
than F0.1, SPR40%). Time and data permitting, MSEs are preferred for fisheries
management.
Loosely defined, harvest policies are guidelines on how harvest levels should be
set in each season, whereas control rules refer to the specific formulae used for a given
harvest policy that specify the target amount of harvest given policy parameters and the

73

estimated state of the system (i.e. spawning stock biomass). Control rules generally fall
into three separate categories; constant exploitation rate, constant catch, and constant
escapement rules, in addition to derivatives of each aimed to correct perceived
weaknesses (Deroba and Bence 2008). Constant exploitation rate rules aim to set catch
quotas to a constant proportion of stock size (Walters and Martell 2004). This builds in
an inherent feedback system; as the stock declines, the quotas do the same, and vice
versa. Constant catch rules set catch quotas at some constant level regardless of stock
size, valuing the stability in catch. Constant escapement rules set catch at all biomass
over some predetermined level, that level generally being chosen to ensure sufficient
levels of spawning stock remain in the population to provide for adequate replacement.
Derivatives of these control rules can include the addition of thresholds, either biomassbased or exploitation rate-based, that aim to decrease exploitation rate or catch at low
stock sizes. Tuning or policy parameters refer to the specific exploitation rate, catch
level, or escapement level used to define a given harvest control rule and dictate the level
of harvest given the estimated state of the system. Policy parameters can also include
biomass or exploitation rate thresholds involved in derivatives of the three types of
harvest control rules. Previous work has not led to general conclusions regarding what
harvest control rule is best for given objectives and fishery dynamics (Deroba and
Bence, 2008), so it is important to consider a suite of different harvest control rules and
policy specific parameters of interest to stakeholders within the MSE.
Cisco, Coregonus artedi, currently support a roe fishery in Thunder Bay, Ontario,
and are managed via a constant exploitation rate control rule, where the total allowable
catch (TAC) is set to 10% of the estimated spawning stock biomass. The full harvest

74

policy includes estimation of the spawning biomass through hydroacoustic surveys, and
allocation of the TAC among a fixed set of license holders. While constant exploitation
rate control rules can sometimes effectively achieve objectives (Deroba and Bence 2008,
Walters and Martell 2004), the specific exploitation rate of 10% put into place in
Thunder Bay has not been evaluated using MSE, rather it was chosen based on
exploitation rates seen as sustainable for other Lake Superior fish stocks such as lake
trout, Salvelinus namaycush, lake whitefish, Coregonus clupeaformis, and lake
sturgeon, Acipenser fulvescens (Ebener et al. 2008, Stockwell et al. 2009). Whereas
precautionary approaches to management are an important first step, such as setting
conservative exploitation rates based on similar species, use of a harvest control rule
tailored to cisco, obtained through a MSE that explicitly accounts for uncertainties
related to cisco recruitment and assessment, could allow managers to better achieve
objectives.
We conducted a simplified MSE of the Thunder Bay cisco stock, projecting the
stock into the future under a variety of different harvest control rules using a stochastic
simulation model to determine which type of control rule and associated policy
parameters performs best in achieving fishery objectives. Our objectives for this analysis
were twofold: 1) determine whether the current exploitation rate of 10% promotes
sustainability of Thunder Bay cisco, and 2) evaluate the performance of alternative
harvest control rules at meeting cisco fishery objectives. In this paper we present a
stochastic simulation model that attempts to account for uncertainty in the recruitment
process, the assessment process, and the sex-specific nature of cisco harvest while
evaluating alternative harvest control rules and tuning parameters.

75

Methods
Harvest Control Rules and Policy Parameters
In preparation for this study, we presented our proposal and solicited input at the
Lake Superior Technical Committee (LSTC) meeting in Sault Ste. Marie, Ontario, in July
2016. The LSTC consists of fishery biologists from agencies around Lake Superior, their
purpose being to advise the Lake Superior Committee on technical information
regarding the status of stocks including management alternatives and guidelines in
making and evaluating fisheries management decisions. Specifically at this meeting we
inquired which type of harvest control rules the technical committee would like us to
consider and also which performance statistics were most important to them (i.e., what
are the objectives for the fishery). Based on input from the committee, we considered
two main types of harvest control rules and their derivatives; constant exploitation rate
and constant catch rules. We explicitly considered two derivatives within each control
rule in addition to their standard formulation (Figure 2.1). For constant exploitation
rate, we considered 1) Constant U (CU), a simple constant exploitation rate control rule
where catch will be proportional to stock size (Figure 2.1A). 2) Constant U threshold 1
(CUT1), defined as a constant exploitation rate until a threshold spawning stock biomass
(SBT) is reached, at which point the exploitation rate linearly declines as a function of
spawning stock size until both are zero (Figure 2.1B). 3) Constant U threshold 2 (CUT2),
defined as a constant exploitation rate until an upper threshold spawning stock biomass
(SBUT) is reached, at which point exploitation rate linearly declines as a function of
spawning stock size and becomes zero at some lower threshold of spawning stock size
(SBLT; Figure 2.1C). These derivatives of the CU rule aim to produce a compensatory

76

response by gradually decreasing fishing mortality below a threshold. For constant catch
control rules, we considered 1) regular constant catch (CC), where catch is constant
regardless of spawning stock size (Figure 2.1D). 2) Conditional Constant Catch 1 (CCC1),
defined as constant catch until some threshold exploitation rate (UT) is reached, a point
at which the control rule reverts to constant exploitation rate at the predetermined
threshold (Figure 2.1E; Clark and Hare 2004, Deroba and Bence 2008). 3) Conditional
Constant Catch 2 (CCC2), defined as constant catch until a threshold spawning stock
biomass (SBT) is reached at which point the catch is reduced to a new lower level of
constant catch (CL, Figure 2.1F). Each of these derivatives of the CC rule aim to keep
catch relatively stable while attempting to avoid high fishing mortality rates at low
spawning stock sizes.

77

U

SBT

SBLT

SBUT

UT
C
CL

SBT

Figure 2.1. Harvest control rules considered in this analysis and associated policy
parameters.

Spawning stock biomass thresholds (SBT, SBUT) were defined as 20, 30, 40, and
50% of unfished spawning stock biomass. Lower spawning stock biomass thresholds for
CUT2 (SBLT) were defined as 20 and 30% of unfished spawning stock biomass. We
decided to not go lower than 20% of unfished spawning stock biomass as a threshold for

78

CUT1 and CUT2 as this is a level at which it has been suggested that fishing should stop
(Thompson, 1993). In addition, numerous studies have determined spawning biomass
should be maintained between 20-50% of unfished spawning biomass (Clark, 1991;
Fujioka et al., 1997; Quinn et al., 1990). Exploitation rates for CU, CUT1, and CUT2 were
defined as 0.05, 0.10, 0.15, 0.20, and 0.25. Constant catch levels (C) were defined as
100,000 kg, 150,000kg and 200,000kg, 250,000kg, and 300,000kg. Exploitation rates
and catch levels were chosen based on their proximity to the current constant
exploitation rate (0.10) and to mean harvest levels over the past 17 years (~160,000kg,
sd ~ 25,000), respectively. Low catch levels may not be economically viable for
fisherman, and very high catch levels may exceed the current fishery capacity, as might
high exploitation rates. Threshold exploitation rates at which CCC1 would revert to CU
(UT) were defined as 0.15, 0.20, 0.25. CCC2 lower catch levels (CL), to be enacted when
spawning stock biomass is estimated below thresholds mentioned above, were defined
as half of the catch level (e.g. quotas of 100,000kg a year would have a CL of 50,000kg).
In total we compared 51 different harvest control rule combinations (Table 1).
Performance Statistics
Performance statistics the LSTC wanted us to consider included the magnitude of
stock size, the magnitude of yield, the variability in yield, and the probability of stock
collapse. The committee also noted that they were primarily interested in the
performance of these metrics over a 50 year time span. For this reason, performance
statistics considered included the percent of years the spawning biomass was below 20%
of unfished spawning biomass, hereafter termed “risk” for brevity, the average harvest
(per year), the median spawning biomass in the final 5 years (Final SB; as a % of

79

unfished level), and the absolute annual variation in yield (AAV), as defined in Punt et
al. (2008):

H H
AAV 
H
y

y 1

y

y

y

Where H y denotes harvest in a given year. These metrics were summarized in terms of
the medians, 25th and 75th percentiles of their distributions over simulations.
Many of the harvest control rules and performance metrics are defined in terms
of spawning biomass:

SBy   N y ,a ,s P( Fisha  250mm) wa ,s
s

a

where wa ,s is sex-specific average weight at age of a cisco estimated using a vonBertalanffy function (Chapter 1), and P( Fisha  250mm) is defined as the probability that
a cisco of a given age is greater than 250 mm (Chapter 1). We assume that fish greater
than 250 mm in length are mature. This definition of spawning biomass was chosen to
align with how the current control rule allocates TAC of cisco in Thunder Bay (biomass
of cisco > 250 mm). We defined the unfished spawning biomass as the average of the
median spawning biomass levels in the final 5 years after running the simulation model
with no harvest.
Model
We developed a stochastic projection model (SPM) based on an integrated
Statistical Catch-at-Age Assessment (SCAA) model developed in Chapter 1. For each

80

control rule, 1000 simulations of the SPM were run to obtain distributions of
performance metrics. The SPM is age- and sex-structured, beginning at age 2 and
forming a plus group at 15. The SCAA model ends in 2015 thus the SPM spans from
2016-2056 (50yr time horizon):

N y 1,a ,s

0.5Ry 1
if a  2

( M  F
)
  N y ,a 1,s e s y ,a 1,s
if 3  a  15 

( M  F
)
( M  F
)
 N y ,14,s e s y ,14,s  N y ,15 ,s e s y ,15 , s if a  15 

where N y ,a ,s is the number of cisco age a of sex s in year y , R y is recruitment in year
y , M s is the natural mortality for sex s (drawn from the SCAA posterior for each

simulation), and Fy ,a ,s refers to fishing mortality for a given year, age, and sex
combination. Each simulation began by drawing from the posterior distribution of sexspecific 2015 abundance at age from the SCAA.
Recruitment
Recruitment of cisco, at least over the past several decades in Lake Superior, has
been characterized by a highly variable, boom-or-bust pattern where a large year class is
produced, followed by many years of almost no recruitment, until the next big year class
is produced (Chapter 1 Figure 1.5; Stockwell et. al, 2009). In the SPM, we modeled this
process by drawing from a Bernoulli distribution each year that determined whether a
given year would be boom or bust. The parameter for this Bernoulli distribution was
drawn for each simulation from a uniform distribution with bounds l and u : U [l , u] . If
a given year within a simulation is characterized as a boom year, a stock-recruit (SR)

81

function is applied; if characterized as bust, the model draws from a lognormal
distribution derived using posterior estimates of bust years from the SCAA model. For
boom years, the SR function used was derived based on the Ricker functional form
(Ricker, 1975) using point estimates (medians) of the posterior distribution of
recruitment and stock size estimates in the SCAA as data. Projected recruitment is then:

Ry  S y 2e

 S y  2  y

e

 y ~ N (0,  r 2 )
Where  and  are parameters of the SR model, which are drawn at random for each
simulation of the SPM from the posterior distribution, and  y are multiplicative
deviations invoking stochastic recruitment over time within a simulation. We fixed  r
at a value of 0.711 based on Thorson et al.’s (2014) meta-analysis of recruitment
deviation for the family Salmonidae. This was done due to the large value of estimated

 r within the SR function (because of sparse data), which had the effect of producing
many unrealistically high projected recruitments when initially used in the SPM. In an
attempt to avoid using assessment output as data, we initially tried to estimate a SR
function within the SCAA however found that the model would not converge on a
solution. Specifics on the derivation of the SR function can be found in the appendix.
Given uncertainty in what level of recruitment constitutes a boom or a bust year,
and based on the fact that the SR function and bounds of the uniform distribution will
be defined by this, we specifically explored 2 different recruitment scenarios. These
scenarios are hereafter termed 7yr and 4yr (Figure 2.2), characterized by how we define

82

what constitutes a boom year. The 7yr scenario treats years in the SCAA that had a
median recruitment (age-2 abundance) over 200k as boom years (7/17 years in the
SCAA fit this criteria), while the 4yr scenario treats years that had a median recruitment
(age-2 abundance) over 1 million as boom years (4/17 years in the SCAA fit this criteria).
The bounds of the uniform distribution for each recruitment scenario were based on the
perceived frequency of “boom” year classes over the past 30 years using observations
from both the SCAA (Chapter 1) and Figure 15 in Yule et al., (2006). These bounds were
defined as U(0.25,0.40) for the 7yr scenario based on evidence 0f ~9-11 boom year
classes over the last 30 years and U(0.15,0.25) for the 4yr scenario, based on evidence of
~6 boom year classes over the last 30 years. Recruitment values in the SCAA that were
not characterized as “boom” recruitment years were placed in the “bust” category and
used to derive a lognormal distribution of “bust” recruitments.

83

Figure 2.2. SR curves for each recruitment scenario. “Data”, medians of the posterior
distribution of the SCAA, are plotted as points. The 7yr scenario SR curve uses all “data”
points while the 4yr scenario solely uses the points highlighted in green.

Fishing Mortality
Our approach to setting fishing mortality rates for each year of the simulation
was to set fishing rates so the resulting harvest matched a value obtained by applying

84

the control rule to the assessed spawning biomass (see Assessment Error below). Some
complexity is added because we are modeling dynamics as sex specific and although
cisco harvest is dominated by female fish (  80%, Chapter 1 - Figure 1.1), there is interannual variation. Our approach was to stochastically simulate the sex ratio of the
fishing intensities (fully selected fishing mortality) each year, and then solve for the
fishing intensity of females, (and given the ratio, the fishing intensity of males) that
produced the desired harvest. The sex ratio of fishing intensities is defined as:
f yr 

f y ,m
f y ,m  f y , f

Where f yr denotes the fishing intensity ratio in a given year, f y ,m is male fishing
intensity, and f y , f is female fishing intensity. Fishing intensity ratios for all 17 years of
the SCAA were drawn for each simulation in the SPM and used to define a beta
distribution. Each beta distribution was defined in terms of two shape parameters, here

  1    
 1 and
2
 


denoted as p and q which can be written as p   

  (1   ) 
q  1   
 1 , where  and  2 are the mean and variance of the ratio of
2
 

fishing intensities pulled from the posterior distribution of the SCAA for each
simulation. The corresponding beta distribution for each simulation was used to draw
fishing intensity ratios for each year within the SPM. Fishing intensity for a given
sex/year combination in each simulation was solved for using Newton-Raphson
iterations given a harvest control rule:

85



  M
s

a



Fa , y ,s
s

 Fa , y ,s



N a , y 1,s 1  e

( M s  Fa , y , s )

W

a ,s


  Hy


Fy ,a ,s  sa f y ,s

where sa refers to age-specific cisco fishery selectivity (parameters that define selectivity
function were drawn from the SCAA posterior distribution) and H y denotes harvest in a
given year and is defined based on a control rule. Female fishing intensity for a given
year was solved for and male fishing intensity was calculated using the fishing intensity
ratio and female fishing intensity:

f y ,m 

f yr * f y , f
1  f yr

We set a maximum fishing mortality rate of 3 to limit unrealistic scenarios that could
have fisherman catching every last fish in a given year.
Assessment Error
Assessment estimation error was simulated within the SPM through an
autoregressive process
y

SˆBy  SBy e

 e2
2

for y  1

 y
y  
2

 y 1  1    y for y  1

 y ~ N (0, e 2 )

Where ŜB y denotes the assessed spawning biomass and SB y is the true spawning
biomass.  and  e were specified as 0.7 and 0.16, assuming a lognormal assessment

86

 e2
error with a CV of about 0.22 (
). This was based on the CV of spawning biomass
1 2
in the final year of the SCAA (  0.22). Alternate values of rho and sigma (  =0.9,  e
=0.3) were explored to assess the sensitivity of results to levels of assessment error.
Similar procedures have been done in previous harvest policy projections (Deroba and
Bence 2012; Irwin et al. 2008, Punt et al. 2008). We did not model implementation error
within the SPM given license holders rarely, if ever, go over their TAC. Thus, assuming
fishers meet their TAC (unless fishing mortality limit is reached) is likely a conservative
assumption.
Sensitivity Analyses
Sensitivity to the bounds of the uniform distribution for the probability of a boom
year class was examined by shifting the distribution  0.05 for each recruitment
scenario. Several of the control rules we considered use unfished spawning stock
biomass, and this value, determined based on running the SPM with negligible harvest,
depends on the distribution for the probability of boom years. Therefore, we explored
two alternate scenarios for setting unfished spawning biomass when shifting the
distribution for boom years. First, we assumed unfished spawning biomass at values
estimated using the original uniform distribution bounds and second we assumed
unfished spawning biomass at new values calculated when the uniform distribution is
either shifted up or down by 0.05. The first of these scenarios explores the situation
where managers erroneously specify the unfished spawning biomass when the
frequency of boom years is shifted, i.e., the shifts represent a situation where system
productivity is both different and miss-specified in the control rule. The second

87

represents a case where the change in system productivity is accounted for in the control
rule.
Results
Equilibrium or unfished spawning biomass was calculated as 4,750,000 kg and
4,400,000 kg for the 4yr and 7yr recruitment scenarios, respectively.
Recruitment Scenario
Relative relationships between harvest control rules were largely robust to
recruitment scenarios. However, absolute values did differ, with results reflecting the
increased productivity for the 7yr scenario (i.e. higher yield, lower risk, higher Final SB,
and lower AAV). For this reason, hereafter in text we present the results solely for the
4yr recruitment scenario. Results for the 7yr recruitment scenario can be found in Table
1 and supplemental figures 2.1-2.4.
Average Yield
Constant exploitation rate and its derivatives (CU, CUT1, CUT2) outperformed
constant catch rules in terms of the maximum (over policy parameters) average yield
over the 50yr simulation period (Figure 2.3). Within CU control rules, as we would
expect, average yield was lowest for the 0.05 rate. As exploitation rate increased from
0.05 to 0.10-0.25 however, an asymptote was reached at around 255,000 kg of yield per
year (Table 1, Figure 2.3). While median yields for CU reached an asymptote, the spread
of the 25-75 quantile range slightly increased as exploitation rate increased from 0.100.25. Derivatives of the CU rule involving thresholds (CUT1 and CUT2) experienced
increases in yield (Figure 2.3) over their CU counterparts with similar exploitation rates,

88

with a CUT2 rule involving an exploitation rate of 0.20 to be linearly reduced at 50% of
unfished spawning stock biomass and made 0 at 30% of unfished spawning stock
biomass (Policy 1.3.10, Table 1) experiencing the largest average yield over the
simulation period at 357,804 kg per year. The constant catch control rules, even at their
highest catch levels (300,000 kg per year), were only able to produce average yields of
around 180,000 kg per year. In fact, when we increased catch levels above 300,000 kg
(up to 550,000 kg) within CC, an asymptote in average yield was reached at around
220,000 kg per year. When thresholds were included in constant catch control rules
(CCC1 and CCC2), yield for control rules with similar catch levels did not increase and in
fact slightly decreased in almost all cases (Exceptions are policies 2.1.3 vs 2.2.3, 2.1.4 vs
2.2.4, and 2.1.5 vs 2.4.9; Table 1, Figure 2.3).

89

Figure 2.3. Summary of the distributions of average harvest over the simulation period
for each respective control rule. Shown are medians (horizontal bar) and 25-75 quantiles
(box). Labels specify policy parameters that make up each control rule (CU = “U”; CUT1
= “U SBT”; CUT2 = “U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”).
Exploitation rates are presented as decimals and biomass thresholds as percentages. For
CUT2 control rules, a label of “0.10 50-20%” describes a control rule that has an
exploitation rate of 0.10 above 50% of the unfished spawning biomass, while that rate
linearly declines below that threshold to 0 at 20% of the unfished spawning biomass.
Catch levels are described in 100,000 kg (i.e. 100k = 100,000 kg).

90

Risk (% of years SB < 20% unfished level)
Where CU rules did not show much difference in yield at 0.10-0.25 exploitation
rates, they exhibited large differences in risk. As exploitation rate increased within the
CU control rule from 0.05-0.25, the amount of risk more than tripled from 18% of years
having a SB below 20% of the unfished level at an exploitation rate of 0.05 to 66% of
years under an exploitation rate of 0.25 (Table 1, Figure 2.4). The inclusion of thresholds
in constant exploitation rate control rules greatly decreased risk within a given
exploitation rate. For CUT1 rules, risk decreased both compared to the respective CU
rule with the same exploitation rate and within the CUT1 rule as the threshold was
increased from 20-50% of unfished SB for all exploitation rates. Risk was further
decreased with the inclusion of a lower threshold SB within the CUT2 rules. That is, for
exploitation rates of 0.10 and 0.20, risk was lower for the CUT2 rule than for its CUT1
and CU counterparts. For an exploitation rate of o.10, risk was 34% for CU, 26% at its
lowest in CUT1, and 20% at its lowest in CUT2 (Policies 1.1.2, 1.2.7, and 1.3.5; Table 1). A
similar result occurred for exploitation rates of 0.20, where under CU risk was 58%,
44% at its lowest under CUT1, and 30% at its lowest under CUT2 (Policies 1.1.4, 1.2.16,
and 1.3.10; Table 1).
Within CC rules, risk increased from 22% at a catch level of 100,000 kg a year to
56% at a catch level of 300,000 kg a year. Risk decreased with the inclusion of
exploitation rate thresholds within the CCC1 control rule. Within CCC1, risk increased as
the threshold exploitation rate increased (Figure 2.4). For each level of catch, the use of
biomass thresholds under the CCC2 rule decreased risk compared to CC control rules. In
addition, within CCC2 risk decreased as threshold SB levels increased. For example,

91

under a catch level of 200,000 kg a year (CC risk=44%), including a biomass threshold
at 20% of unfished SB decreased risk to 34% and including a biomass threshold at 30%
of unfished SB decreased risk to 32%. The lowest risk level over all control rules was
therefore under a CCC2 rule with the lowest catch level, 100,000 kg, and a threshold of
30% of the unfished spawning biomass at which point the catch level would be cut in
half (Policy 2.3.2, risk=16%).

92

Figure 2.4. Summary of the distributions of risk level for each respective control rule.
Risk is defined as the percentage of years in each simulation where SB is below 20% of
the unfished level. Shown are medians (horizontal bar) and 25-75 quantiles (box).
Labels specify policy parameters that make up each control rule (CU = “U”; CUT1 = “U
SBT”; CUT2 = “U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”). Exploitation
rates are presented as decimals and biomass thresholds as percentages. For CUT2
control rules, a label of “0.10 50-20%” describes a control rule that has an exploitation
rate of 0.10 above 50% of the unfished spawning biomass, while that rate linearly
declines below that threshold to 0 at 20% of the unfished spawning biomass. Catch
levels are described in 100,000 kg (i.e. 100k = 100,000 kg).

93

Absolute Annual Variation in Yield (AAV)
Over all control rules, AAV proved considerably smaller for the constant catch
control rules (Figure 2.5). For example, a CC rule with a catch level of 200,000 kg a year
(Policy 2.1.3) had an AAV of 0.06 while a CU rule with an exploitation rate of 0.15
(Policy 1.1.3) had an AAV of 0.32. Also, the inclusion of a threshold within any rule
(CUT1 & CUT2 as compared to CU and CCC1 & CCC2 as compared to CC) increased
AAV. Within constant exploitation rate control rules, AAV increased as exploitation rate
increased. The inclusion a threshold biomass levels for CUT1 rules increased AAV over
all exploitation rates, and the inclusion of a lower threshold biomass at which
exploitation rate would become zero (for CUT2) increased AAV further compared to
CUT1 and CU control rules with similar exploitation rates.
For constant catch control rules, AAV increased as catch level increased, from
0.01 at 100,000 kg a year (Policy 2.1.1) to 0.12 at 300,000 kg a year (Policy 2.1.5). The
inclusion of threshold exploitation rates for CCC1 increased AAV as well. For example, a
CC rule with a catch level of 200,000 kg a year (Policy 2.1.3) had an AAV of 0.06 while a
CCC1 rule with a catch level of 200,000 kg per year and a threshold exploitation rate of
0.15 (Policy 2.2.1) increased AAV to 0.10 . Within CCC1, AAV decreased as the threshold
exploitation rate increased for a given catch level. The inclusion of threshold SB levels
related to CCC2 also increased AAV compared to CC with similar catch levels, while
within CCC2, AAV increased as threshold SB increased for similar catch levels.

94

Figure 2.5. Summary of the distributions of absolute annual variation in yield for each
respective harvest control rule Shown are medians (horizontal bar) and 25-75 quantiles
(box). Labels specify policy parameters that make up each control rule (CU = “U”; CUT1
= “U SBT”; CUT2 = “U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”).
Exploitation rates are presented as decimals and biomass thresholds as percentages. For
CUT2 control rules, a label of “0.10 50-20%” describes a control rule that has an
exploitation rate of 0.10 above 50% of the unfished spawning biomass, while that rate
linearly declines below that threshold to 0 at 20% of the unfished spawning biomass.
Catch levels are described in 100,000 kg (i.e. 100k = 100,000 kg).

95

Spawning Biomass at the end of the SPM (Final SB)
Spawning biomass at the end of the simulation period, defined as the median SB
for the final 5 years of each simulation (Final SB, presented as a percentage of unfished
SB), was similar among baseline harvest control rules (CU & CC, Figure 2.6), however
the spread of the Final SB for constant catch control rules was much greater than that of
the constant exploitation rate control rules.
Within CU rules, Final SB decreased as exploitation rate increased, from 61% of
the unfished level at an exploitation rate of 0.05 (Policy 1.1.1) to 8% at an exploitation
rate of 0.25 (Policy 1.1.5). For any given exploitation rate, adding a SB threshold within
CUT1 ubiquitously increased Final SB, and CUT2 rules involving an additional lower
threshold further increased Final SB, performing best among constant exploitation rate
control rules. For example, a CU rule with an exploitation rate of 0.10 produced a Final
SB 33% of the unfished level (Policy 1.1.2) while a CUT2 rule with an exploitation rate of
0.10, an upper SB threshold of 50% of unfished SB, and a lower SB threshold of 30% of
unfished SB produced a Final SB of 51% of the unfished level (Policy 1.3.5, Table 1).
Within CUT1 rules of a given exploitation rate, Final SB increased as threshold biomass
increased. Similarly, within CUT2 rules given a level of exploitation rate, Final SB
increased as both upper and lower SB thresholds increased.
For constant catch, within the CC control rule Final SB declined as catch levels
increased, from 58% of the unfished level at 100,000 kg a year (Policy 2.1.1), to 13% at
300,000 kg a year (Policy 2.1.5). The inclusion of a threshold exploitation rate within
CCC1 mostly increased Final SB, and within CCC1 Final SB decreased as threshold
exploitation rate increased. For catch levels of 200,000 kg per year, the addition of

96

threshold exploitation rates increased Final SB for rates of 0.15 and 0.25 and decreased
Final SB at a rate of 0.20 (likely an anomalous result due to stochasticity). For catch
levels of 250,000 kg, including exploitation rate thresholds increased Final SB (from
18% to 33%, 25%, and 19%, Policy 2.1.4 compared to 2.2.4, 2.2.5, and 2.2.6). For all
catch levels, the inclusion of SB thresholds within CCC2 rules increased Final SB levels
compared to CC rules with similar catch levels. Final SB also increased as SB threshold
increased within CCC2 rules.

97

Figure 2.6. Summary of the distributions of final spawning biomass for each respective
control rule, with final spawning biomass defined as the median of the last 5 years
spawning biomass in each simulation, characterized as a percentage of the unfished
level. Shown are medians (horizontal bar) and 25-75 quantiles (box). Labels specify
policy parameters that make up each control rule (CU = “U”; CUT1 = “U SBT”; CUT2 =
“U SBUT-SBLT”; CC = “C”; CCC1 = “C UT”; CCC2 = “C SBT”). Exploitation rates are
presented as decimals and biomass thresholds as percentages. For CUT2 control rules, a
label of “0.10 50-20%” describes a control rule that has an exploitation rate of 0.10
above 50% of the unfished spawning biomass, while that rate linearly declines below
that threshold to 0 at 20% of the unfished spawning biomass. Catch levels are described
in 100,000 kg (i.e. 100k = 100,000 kg).

98

Performance metrics under similar levels of yield
For CU and CC control rules, an exploitation rate of 0.05 and a catch level of
250k produced similar average yields over the simulation period (Policy 1.1.1 = 172,608
kg/year, Policy 2.1.4 = 181,734 kg/year). For these two specific control rules with similar
yield, risk was greater (Figure 2.7 upper panel, Policy 1.1.1 = 18%, Policy 2.1.4 = 48%),
and spawning biomass at the end of the time series was lower (Figure 2.7 middle panel,
Policy 1.1.1 = 61%, Policy 2.1.4 = 18%) for CC with a catch level of 250k. However AAV
was lower for CC with a catch level of 250k compared to CU with an exploitation rate of
0.05 (Figure 2.7 lower panel, Policy 1.1.1 = 0.26, Policy 2.1.4 = 0.09). The inclusion of
thresholds in CCC2 did not alter this comparison. For example, while producing a
similar amount of yield, CU with an exploitation rate of 0.05 (Policy 1.1.1) produced
lower risk, greater Final SB, and greater AAV compared to CCC2 with a catch level of
250,000 kg per year and a threshold SB of 30% (Policy 2.3.8).

99

Figure 2.7. Comparison of performance metrics for CU rule with an exploitation rate of
0.05 (black points) and CC rule with a catch level of 250,000 kg (red points). Upper
panel plots risk (percentage of years SB < 20% unfished) versus average harvest (in
millions of kg) obtained in the same individual simulations. Middle panel compares
final spawning biomass (millions of kg) and average harvest. Lower panel compares
absolute annual variation in yield and average harvest. Only results from simulations
that produced an average harvest within the 25-75% quantiles (i.e., inside the
interquartile range) are plotted.

100

Table 2.1. Performance statistics are presented as results for the 4yr and 7yr recruitment
scenarios (4yr | 7yr). Values are presented as medians over simulations. Average yield
(kg) is mean yield over years. Risk is calculated as the percentage of years SB is below
20% of the unfished condition. AAV is defined in methods. Final spawning biomass is
the median SB of the last 5 years in a simulation. Catch levels for constant catch control
rules are presented in 100,000 kg (i.e. 100k=100,000 kg).
Harvest
Policy

Policy Parameters

Yield (kg)

Risk

AAV

Final SB

CU
1.1.1

U=0.05

172608 | 192664

0.18 | 0.06

0.26 | 0.23

0.61 | 0.84

1.1.2

U=0.10

250215 | 309911

0.34 | 0.12

0.29 | 0.25

0.33 | 0.61

1.1.3

U=0.15

258561 | 374726

0.48 | 0.22

0.32 | 0.29

0.19 | 0.47

1.1.4

U=0.20

254988 | 354457

0.58 | 0.36

0.34 | 0.31

0.11 | 0.25

1.1.5

U=0.25

250156 | 352739

0.66 | 0.48

0.36 | 0.32

0.08 | 0.15

CUT1
1.2.1

U=0.05, SBT=20%

191189 | 189855

0.16 | 0.06

0.26 | 0.23

0.68 | 0.89

1.2.2

U=0.05, SBT=30%

185240 | 196540

0.16 | 0.06

0.27 | 0.23

0.63 | 0.84

1.2.3

U=0.05, SBT=40%

179941 | 194120

0.16 | 0.04

0.28 | 0.24

0.63 | 0.85

1.2.4

U=0.05, SBT=50%

187503 | 195835

0.14 | 0.04

0.28 | 0.25

0.65 | 0.90

1.2.5

U=0.10, SBT=20%

248154 | 309159

0.30 | 0.12

0.31 | 0.27

0.38 | 0.65

1.2.6

U=0.10, SBT=30%

270291 | 322572

0.28 | 0.12

0.32 | 0.27

0.47 | 0.68

1.2.7

U=0.10, SBT=40%

278925 | 318348

0.26 | 0.10

0.33 | 0.28

0.45 | 0.69

1.2.8

U=0.10, SBT=50%

264902 | 325076

0.28 | 0.10

0.34 | 0.29

0.44 | 0.72

1.2.9

U=0.15, SBT=20%

282876 | 369314

0.44 | 0.22

0.36 | 0.30

0.24 | 0.46

1.2.10

U=0.15, SBT=30%

274787 | 376849

0.40 | 0.20

0.37 | 0.31

0.27 | 0.50

1.2.11

U=0.15, SBT=40%

294627 | 400517

0.36 | 0.16

0.38 | 0.32

0.27 | 0.52

1.2.12

U=0.15, SBT=50%

298119 | 403222

0.34 | 0.14

0.39 | 0.32

0.33 | 0.57

1.2.13

U=0.20, SBT=20%

284403 | 402488

0.52 | 0.28

0.39 | 0.32

0.16 | 0.37

1.2.14

U=0.20, SBT=30%

297737 | 393828

0.50 | 0.28

0.41 | 0.34

0.21 | 0.38

1.2.15

U=0.20, SBT=40%

309875 | 392536

0.44 | 0.26

0.43 | 0.36

0.21 | 0.40

1.2.16

U=0.20, SBT=50%

294041 | 421060

0.44 | 0.22

0.44 | 0.36

0.24 | 0.41

1.2.17

U=0.25, SBT=20%

306222 | 386964

0.56 | 0.40

0.41 | 0.35

0.14 | 0.28

1.2.18

U=0.25, SBT=30%

286845 | 397990

0.56 | 0.36

0.43 | 0.36

0.15 | 0.28

1.2.19

U=0.25, SBT=40%

310483 | 418427

0.52 | 0.32

0.45 | 0.38

0.17 | 0.33

101

Table 2.1. (cont’d)
1.2.20

295778 | 425424

0.50 | 0.30

0.47 | 0.39

0.20 | 0.35

265445 | 309333

0.26 | 0.10

0.34 | 0.28

0.45 | 0.70

281574 | 307611

0.22 | 0.10

0.35 | 0.30

0.49 | 0.69

280217 | 313174

0.22 | 0.08

0.36 | 0.31

0.49 | 0.75

288841 | 319980

0.21 | 0.08

0.36 | 0.30

0.53 | 0.70

276439 | 312400

0.20 | 0.06

0.37 | 0.31

0.51 | 0.77

315477 | 427299

0.42 | 0.24

0.46 | 0.37

0.26 | 0.44

305147 | 424042

0.40 | 0.20

0.47 | 0.39

0.28 | 0.48

347791 | 443634

0.36 | 0.16

0.48 | 0.39

0.30 | 0.52

336807 | 433390

0.34 | 0.16

0.49 | 0.40

0.31 | 0.51

357804 | 413688

0.30 | 0.16

0.50 | 0.42

0.34 |0.50

C=100k

99439 | 99999

0.22 | 0.06

0.01 | 0

0.58 | 0.91

2.1.2

C=150k

138031 |149996

0.32 | 0.08

0.04 | 0

0.43 | 0.74

2.1.3

C=200k

161574 | 197245

0.44 | 0.14

0.06 | 0.02

0.28 | 0.60

2.1.4

C=250k

181734 | 236464

0.48 | 0.18

0.09 | 0.04

0.18 | 0.49

2.1.5

C=300k

188024 | 260835

0.56 | 0.26

0.12 | 0.06

0.13 | 0.33

CCC1
2.2.1

C=200k, UT=0.15

155451 | 184377

0.34 | 0.10

0.10 | 0.05

0.38 | 0.75

2.2.2

C=200k, UT=0.20

155038 | 191569

0.40 | 0.10

0.10 | 0.03

0.25 | 0.68

2.2.3

C=200k, UT=0.25

163187 | 193259

0.38 | 0.10

0.08 | 0.03

0.30 | 0.65

2.2.4

C=250k, UT=0.15

182654 | 218018

0.34 | 0.12

0.13 | 0.07

0.33 | 0.65

2.2.5

C=250k, UT=0.20

173591 | 225785

0.44 | 0.14

0.12 | 0.06

0.25 | 0.65

2.2.6

C=250k, UT=0.25

180300 | 232787

0.44 | 0.14

0.10 | 0.04

0.19 | 0.60

90001 | 97996

0.20 | 0.04

0.04 | 0.02

0.71 | 0.86

85999 | 94000

0.16 | 0.04

0.05 | 0.03

0.69 | 0.92

131872 | 143999

0.24 | 0.08

0.06 | 0.02

0.47 | 0.82

123783 | 139498

0.22 | 0.06

0.07 | 0.04

0.54 | 0.85

CUT2
1.3.1
1.3.2
1.3.3
1.3.4
1.3.5
1.3.6
1.3.7
1.3.8
1.3.9
1.3.10
CC
2.1.1

CCC2
2.3.1
2.3.2
2.3.3
2.3.4

U=0.25, SBT=50%
U=0.10, SBUT=30%,
SBLT=20%
U=0.10, SBUT=40%,
SBLT=20%
U=0.10, SBUT=50%,
SBLT=20%
U=0.10, SBUT=40%,
SBLT=30%
U=0.10, SBUT=50%,
SBLT=30%
U=0.20, SBUT=30%,
SBLT=20%
U=0.20, SBUT=40%,
SBLT=20%
U=0.20, SBUT=50%,
SBLT=20%
U=0.20, SBUT=40%,
SBLT=30%
U=0.20, SBUT=50%,
SBLT=30%

C=100k, SBT=20%,
CL=50k
C=100k, SBT=30%,
CL=50k
C=150k, SBT=20%,
CL=75k
C=150k, SBT=30%,
CL=75k

102

Table 2.1. (cont’d)
2.3.5
2.3.6
2.3.7
2.3.8
2.4.9
2.3.10

C=200k, SBT=20%,
CL=100k
C=200k, SBT=30%,
CL=100k
C=250k, SBT=20%,
CL=125k
C=250k, SBT=30%,
CL=125k
C=300k, SBT=20%,
CL=150k
C=300k, SBT=30%,
CL=150k

159809 | 187998

0.34 | 0.12

0.08 | 0.04

0.31 | 0.70

149982 | 181999

0.32 |0.08

0.08 | 0.04

0.39 | 0.78

181532 | 232018

0.42 | 0.14

0.09 | 0.04

0.23 | 0.61

176569 | 219999

0.38 | 0.12

0.10 | 0.06

0.30 | 0.64

191462 | 256409

0.50 | 0.23

0.12 | 0.06

0.17 | 0.42

185590 | 252969

0.46 | 0.18

0.12 | 0.07

0.19 | 0.50

Sensitivity
Results were largely robust to higher levels of assessment error (  e =0.3) in
addition to increased levels of autocorrelation (  =0.9), as relative comparison of
harvest control rules based on performance statistics changed little compared to the
status quo results (Supplemental figures 2.5-2.12). For AAV, absolute values were higher
among all constant exploitation rate control rules for  e = 0.3 compared to status quo
scenario (Supplemental figure 2.7). Under scenarios where bounds of the uniform
distribution defining the probability of a boom year class are shifted by + or - 0.05 and
management correctly specifies unfished spawning biomass (according to runs of the
SPM with new uniform bounds), relative comparison of harvest control rules was
similar (Supplemental figures 2.13-2.20). However, absolute values of the performance
statistics did change significantly, where under a less productive scenario (bounds of the
uniform - 0.05), risk and AAV increased and Final SB and yield decreased, and for the
more productive counterpart (bounds of the uniform + 0.05), the opposite occurred.
Under a scenario where bounds of the uniform distribution defining the probability of a
boom year class are shifted by - 0.05 and management incorrectly specifies unfished
spawning biomass according to status quo scenario (4,750,000 kg), the only perceived
103

change from that where it correctly calculates unfished SB were increased AAV
(Supplemental figure 2.23), decreased risk (Supplemental figure 2.22), and a slight
increase in Final SB for control rules with biomass based thresholds (Supplemental
figure 2.24). Yield did not change appreciably across all harvest control rules. For a
scenario where bounds of the uniform distribution are shifted by + 0.05 and
management incorrectly specifies unfished spawning biomass according to status quo
scenario (4,750,000 kg), changes from that where it correctly calculates unfished SB
included a slight decrease in AAV (Supplemental figure 2.27) and an increase in risk
(Supplemental figure 2.26) for most control rules with biomass based thresholds. For
most harvest control rules changes between scenarios where bounds of the uniform
distribution are shifted by + 0.05 and unfished SB is either correctly or incorrectly
assumed were imperceptible in terms of yield and Final SB with the exception of
decreases in Final SB for most CUT2 control rules at exploitation rates of 0.10 and 0.20
(Supplemental Figure 2.28). It should be noted for both scenarios where unfished
biomass is assumed incorrectly that relative comparisons of all control rules are largely
unchanged.
Discussion
To address the first objective specified—to determine whether the current 10%
exploitation rate promotes sustainability of the Thunder Bay cisco fishery—we must
specify what constitutes “sustainability” of cisco in Thunder Bay. One simple way to look
at sustainability might be to observe the distribution of SB each year over the time series
and determine whether it is stable near the end, i.e. does the population distribution
crash or is it on a downward trajectory? In this case the 10% rate is “sustainable”, as the

104

trajectory over the 50yr time period for the 4yr recruitment scenario is seemingly stable
at around 1.5 million kg of SB (Figure 2.8).

Figure 2.8. Spawning biomass for the projection of the current harvest control rule, a
10% exploitation rate. Shown are medians (horizontal bar) and 25-75 quantiles (box).

A more robust way to explore the sustainability question may to examine it in
terms of maintaining SB above a threshold to ensure sufficient replenishment, as many
studies have presented arguments for maintaining SB above certain thresholds in fish

105

populations (Beddington and Cooke, 1983; Caddy and Mahon, 1995; Clark 1991; Francis
1993; Fujioka et al., 1997; Goodyear, 1993; Hollowed and Megrey, 1993; Leaman, 1993;
Quinn et al., 1990; Thompson, 1993). What is evident from these analyses is the
argument for maintenance of >20% of unfished spawning stock size. If we utilize this
criteria, the current 10% exploitation rate is usually “sustainable”, as the SPM projects a
median Final SB of 33% and 61% of the unfished level for the 4yr and 7yr scenarios
respectively. This “sustainability” designation is largely insensitive to reduced
productivity in terms of the probability of a boom year class. For example, when the
SPM is re-run with bounds of the uniform distribution defining the probability of a
boom year class shifted down by 0.05, Final SB is 27% and 52% of the new unfished
level (estimated using new bounds) under the 4yr and 7yr scenarios.
Alternatively, Mace (1994) and Myers et al. (1994) suggest maintaining spawning
biomass at a level that would produce 50% of the maximum recruitment level. If we
follow this convention and use the median estimates of our SR functions (Figure 2.2),
female SB should not decrease below ~24% of the unfished level for the 4yr scenario and
~15% for the 7yr scenario. Based on this criteria the current 10% rate is also seen as
generally sustainable. However, this designation may be risk-prone given the nature of
cisco year classes, as the SR function is only applied on average every few years
depending on the recruitment scenario. For our specific case study, we recommend
against this metric, given a large degree of uncertainty surrounding recruitment of cisco
in Thunder Bay.
In terms of our second objective: Can this control rule be improved upon to both
promote sustainability and meet fishery objectives? The answer is more complicated.

106

Within the framework of the CU control rule and levels of exploitation we considered,
the answer is no, as the current 10% rate effectively maximizes yield, maximizes Final
SB, and minimizes both risk and AAV compared to higher exploitation rates. However,
the adoption of a CUT1 or CUT2 rule will slightly increase yield, greatly decrease risk,
and increase Final SB. It is also possible that slight improvements could be obtained by
more fine evaluation of exploitation rates between 0.05 and 0.15. These results are
similar to those found by Deroba and Bence (2012) for lake whitefish, Coregonus
clupeaformis, in 1836 treaty waters of the Laurentian Great Lakes. The tradeoff lies in
the AAV, where adoption of a CUT2 rule will increase year-to-year variation in yield
most, followed by CUT1 rules compared to the current CU control rule. This is due to the
compensatory mechanism within these control rules that aims to change exploitation
rate below biomass thresholds. This difference averages around a ~4 unit increase in
AAV from CU to CUT1 and a ~7 unit increase from CU to CUT2 under an exploitation
rate of 0.10. If stakeholders are indifferent to this increase in AAV, and rather more
interested in magnitude of yield, decrease in risk, and increase in the Final SB, a CUT2
rule is likely most appropriate for cisco in Thunder Bay. Where we only ran the CUT2
rule for biomass thresholds ranging from 20-50% and exploitation rates of 0.10 and
0.20, one could surmise based on the relationship between the performance statistics,
biomass thresholds and exploitation rates, what effect alternate levels would produce.
For example had we chosen to include a lower biomass threshold (SBLT) at 10% of the
unfished level, this would likely have had the effect of decreasing AAV, Final SB and
yield (slightly) while increasing risk.

107

If low variation in yield is held in a high regard, or a higher regard than all other
performance metrics, a more appropriate control rule to adopt may be a constant catch
rule. These control rules did not outperform most exploitation rate rules (all but 0.05
rate) in terms of yield and exhibited more risk and smaller Final SB under rules
producing similar yield levels (for CC 250k vs CU 0.05, Figure 2.7), however they proved
far superior when comparing year-to-year variation in yield. When discussing the
efficacy of constant catch rules, it is important to note that if used indefinitely and
without a threshold, it has been argued constant catch will eventually lead to extinction
due to a stochastic environment (Punt 2010). Theoretically, this can be mitigated
through the inclusion of exploitation rate thresholds within CCC1 control rules. In this
study, the addition of thresholds within constant catch mostly decreased yield, increased
AAV, and usually increased Final SB and decreased risk. Out of the conditional constant
catch control rules, CCC2 rules were slightly more effective in decreasing risk, increasing
Final SB, while not costing much in yield and AAV compared to CC rules with similar
catch levels. Although it is important to note that where CCC2 performed better than
CCC1 over the 50yr time horizon in terms of AAV and both decreasing risk and
increasing Final SB compared to CC rules, over the long term the theory in Punt (2010)
predicts CCC2 will fail. The inclusion of a lower biomass threshold at which point catch
level is made zero within CCC2 could potentially mitigate this. It should also be noted
that CCC2 rules may perform better at different threshold levels (40-50%) and also
different lower catch levels (CL). We only simulated catch being reduced to half of its
original level at biomass thresholds of 20-30%. Similar to CUT2 rules, it is likely the
inclusion of higher biomass thresholds within CCC2 would greatly decrease risk and
increase Final SB, and the decrease in yield could potentially be mitigated by increasing

108

the lower catch level to 75% of its original level, although this would likely have the
effect of costing more in terms of AAV. If constancy in yield is held in a much higher
regard than other performance metrics, adoption of a constant catch control rule with a
threshold of the CCC2 type will most appropriately meet fishery objectives. If Punt
(2010)’s theory of extinction given a constant catch level and a stochastic environment is
considered, the addition of a lower biomass threshold at which point catch level is made
zero within CCC2 would be prudent. This would likely have the effect of decreasing risk,
increasing Final SB, and increasing AAV. Whether the addition of this threshold would
allow CCC2 to continue to outperform CCC1 in terms of AAV is unknown. In addition, it
is unclear what effect this lower threshold would have on yield, based on the fact that a
lower threshold for constant exploitation rate control rules actually proved to increase
yield over base CU rules.
Results comparing the four performance criteria were largely insensitive to
changes in the level and correlation of assessment error. This has been noted in similar
studies (Irwin et al., 2008; Deroba and Bence, 2012; Punt et al., 2008), where in others
it has proved consequential (Katsukawa 2004), largely in the direction of increased
assessment error decreasing the performance of control rules involving biomass
thresholds. It may be that the levels of assessment error we simulated (  e =0.3) are not
high enough to decrease the improvement of threshold-based control rules over those
without thresholds. One could imagine that as assessment error increases to infinity,
control rules based on changing exploitation or catch as a function of the assessed value
would diminish in performance. Our approach to simulating assessment error via
distributions instead of performing a full stock assessment simulation every year in the

109

SPM was primarily driven by time constraints for analysis. The lack of sensitivity to
assessment error suggests that results are likely robust to this simplifying assumption.
Similarly, our simulations assumed a stock assessment would be performed every year
for the stock. Our study could benefit from additional simulations including
management strategies where the control rule is applied to hydroacoustic estimates of
abundance (how TAC is currently set) to inform the utility of performing stock
assessments in each year.
Although relative comparison of the harvest control rules was largely unchanged
under different recruitment hypotheses/scenarios, it is important to note that if harvest
policy decisions are to be based in some part on the absolute values of certain metrics,
such as maintenance of a Final SB above 20% of the unfished level, more liberal
exploitation rates or catch levels may be employed under the 7yr scenario. For this
reason we recommend primarily comparing the relative performance of control rules
when making harvest policy decisions on cisco and when decisions necessitate
information on an absolute value, to err in a conservative fashion and use results from
the 4yr scenario. This subject is relevant once again when discussing sensitivity to
reduced productivity in terms of the probability of a boom year class. These sensitivity
runs which involved shifting the uniform distribution defining the probability of a boom
year class down by – 0.05 resulted in the same relative performance across all harvest
control rules. Although not surprisingly, absolute values differed, potentially resulting in
different conclusions as to which specific control rule meets sustainability criteria.
Although it should be noted that under reduced productivity, a CUT2 rule at an
exploitation rate of 0.10 can still achieve a Final SB > 20% of the unfished level.

110

In analyzing the sensitivity of our results to the probability of a boom year class,
the question of whether or not unfished biomass is correctly estimated pertaining to this
change in productivity was examined. What we may expect in a scenario where unfished
biomass is overestimated is altered performance of control rules utilizing biomass
thresholds, as the “true” threshold levels are actually higher than what is thought (i.e.
40% of the unfished level becomes 60% of the unfished level). This should result in an
increase in year-to-year variation in yield as catch and exploitation rates are changed
more frequently, and a decrease in yield as catch levels and exploitation rates decrease
at relatively high stock sizes, when the stock is in no apparent danger. This should also
result in decreased risk given more conservative thresholds. Our results indicated no
apparent decrease in yield for control rules with biomass-based thresholds however we
found increased AAV, decreased risk, and a slight increase in Final SB when unfished
biomass is overestimated compared to when it is correctly specified within the lower
productivity scenario. This may indicate that for the lower productivity scenario
(bounds on the probability of boom - 0.05), increasing biomass-based thresholds on
control rules would not cost in terms of reduced yield, however would result in lower
risk and higher Final SB at the cost of larger AAV. Alternatively what we might expect in
a scenario where unfished spawning biomass is underestimated is increased risk as
“true” biomass thresholds are more liberal. This should also result in decreased Final
SB, decreased yield, and decreased AAV. Our results corroborated all but decreased
Final SB and decreased yield (for most biomass-based control rules) as it pertains to
control rules with biomass-based thresholds when unfished biomass is underestimated.
This is likely due to the fact that the simulation was of the high-productivity scenario
(probability of a boom year class + 0.05), and exploitation rates may not have been high

111

enough to make a difference. Implications of over and under-estimation of unfished SB
on results when the productivity of the stock (bounds of the uniform) is shifted up or
down by 0.05 are largely trivial, as CUT1 and CUT2 rules continue to outperform CU
rules in terms of risk, Final SB, and yield, however at an increased cost of AAV (if
unfished SB is overestimated). In addition, CCC2 rules continued to outperform CC
rules in terms of risk and Final SB, at little to no additional cost to yield and AAV over
scenarios where unfished spawning biomass is correctly specified. Simulations of overand under-estimation of unfished spawning biomass on the status quo scenario (regular
binomial bounds) are needed to determine effects on the performance of harvest control
rules for the most probable recruitment scenario.
The reliability of estimated unfished biomass levels has been discussed in
previous studies, where life history characteristics of a species and temporal
autocorrelation in recruitment have been shown to alter estimation performance
(Haltuch et al., 2008, 2009). Our method of estimating unfished biomass most closely
relates to the “dynamic B0” method defined in (MacCall et al., 1985) and compared to
other methods (calculating the equilibrium point in the SR relationship and average
recruitment combined with spawning biomass per-recruit analyses) in Haltuch et al.,
(2008, 2009). Where Haltuch et al., (2008) and Haltuch et al., (2009) concluded that
calculating unfished biomass based on a fitted SR function was generally best (Haltuch
et al., 2009 – “if available catch and survey data do not span at least 50 years”), this
method was unavailable to us given the boom-or bust dynamics we specified in
projecting recruitment (i.e. SR function is only applied for select “boom” years). It
should also be noted that Haltuch et al., (2008) found that for all methods of estimating

112

unfished biomass examined, performance was generally poorer in the presence of high
recruitment variability, which cisco clearly exhibit. If the specification of unfished
biomass based on the SPM is of concern to managers, an alternative to setting biomass
thresholds based on an estimated unfished level is to set them based on some low
objective value, i.e., no harvest below 500,000kg of spawning biomass. This could have
the benefit of retaining some desirable characteristics of threshold policies (decreased
risk, increased Final SB) while not having to rely on correctly estimating the unfished
level of the stock.
As important as any findings of a study are the associated caveats based on its
assumptions and limitations. A critical assumption we make in our study is that the
probability of a boom year class, or the productivity of the stock, is static through time.
The dominant theory in the literature as it pertains to what is driving these sporadic
recruitment years for cisco is one of match-mismatch, where abiotic and biotic factors
are hypothesized to line up once every few years to allow for large cisco recruitment
events (Myers et al., 2015). Of these factors, high wind speeds have been shown to be
negatively correlated with cisco year class strength (Myers et al., 2015), the hypothesized
mechanism being an increased likelihood of advection of cisco larvae into colder and
less productive waters under high wind scenarios. It has also been hypothesized that
years with a large degree of ice cover and high temperatures directly after ice out benefit
cisco year classes through an increase prey availability and abundance (Myers et al.,
2015). Under anthropogenic climate change, surface temperatures are expected to
increase, potentially resulting in less ice cover in addition to increased wind speed
(Desai et al., 2009), each of which could result in decreasing the probability of a large

113

cisco year class. Simulations are necessary that take into account changing
environmental conditions through time in assessing the relative performance of harvest
control rules as they pertain to cisco.
In addition we must mention uncertainties about the process by which we
projected recruitment. First and foremost, we used assessment model output as data in
deriving the SR model. This has been called "doing statistics on statistics", or a “twostage analysis” (Link 1999) and potential issues resulting from this have been discussed
extensively in the fisheries literature (Brooks and Deroba, 2015; Maunder and Punt,
2013; Thorson et al., 2013). While we would have ideally avoided this, unfortunately
existing cisco SR functions that have been developed are inadequate based on known
bias in sampling gear (Stockwell et al., 2006), and we were unable to estimate a SR
function within the SCAA.
A second issue that applies to our stock-recruit analysis, as it does to many, is
that the data were sparse, or more specifically, few data were near the origin. This forced
us to rely on a heavily influential prior for the log of the maximum annual reproductive
rate parameter within the function, log( ~) . This prior was developed in Myers (1999)
based on stocks of the family Salmonidae, of which 100/106 were stocks of different
salmon species (others were 5 brook trout stocks and 1 lake trout). Where this prior
relates to the most taxonomically proximal group to cisco available in the meta-analyses,
it is safe to say the reproductive strategies of cisco and the species to which the prior was
derived are drastically different. Cisco generally aggregate in bays and nearshore areas
during the late-fall to spawn. They mate in the upper water column and then let their
eggs drop to the bottom showing no apparent preference of bottom substrate (Stockwell

114

et al. 2009). Cisco also spawn many times throughout their life span, whereas salmon
are generally anadromous and semelparous, meaning they travel up rivers to spawn
(usually over a preferred rocky substrate) and often die after doing so. What is
important in these descriptions of spawning habits is that while the species are relatively
close taxonomically, they actually have quite different reproductive habits. Our analyses
would greatly benefit from a more appropriate prior on the maximum annual
reproductive rate parameter, perhaps one based on species belonging to the subfamily
coregoninae, which includes bloater, kiyi, and other ciscoes that share similar
reproductive habits. Unfortunately stock-recruitment data for this group are sparse, and
a prior based on this subfamily is, to the best of our knowledge, nonexistent. The same
can be said for  r , the recruitment variability, which was chosen based on a metaanalysis of recruitment variation for the family Salmonidae (Thorson et al., 2014).
Depensation was also not considered in our analyses. The inclusion of
depensation at low stock sizes would have likely led to higher risk, lower Final SB, lower
yield, and larger AAV over all control rules as simulations that reached a depensatory
state may have hit a lower equilibrium point difficult to come out of. It is likely the
advantage of threshold control rules would have been increased over baseline CU and
CC with the inclusion of depensation.
Lastly, our treatment of “bust” recruitment years (i.e., drawing recruitment from
a lognormal distribution derived using recruitment estimates from the SCAA designated
“bust” years) may have presented a rescue for very low stock sizes. This treatment
presents the opposite effect that depensation would where recruitment under “bust”
designated years does not decrease as stock size decreases to very low levels. The effect

115

this had on our final results however is likely minimal, as bust recruitments are defined
at very low values and are highly unlikely to rescue a stock out of a state of very low
biomass.
It is important to note that where the current control rule in Thunder Bay is
defined as a function of the biomass of fish > 250mm. In Minnesota waters, the control
rule is defined in terms of the biomass of fish > 305mm. For this study we followed the
Thunder Bay convention in defining spawning biomass as cisco > 250mm given these
individuals are generally mature (Yule et al., 2006; Yule et al., 2008). If the results of
this comparison are to be used in determining harvest policies and control rules in other
cisco harvesting regions, the implication of different definitions for spawning biomass
should be considered.
In summary we have shown in this study that the current exploitation rate of 0.10
on Thunder Bay cisco is sustainable (given certain criteria). We have also simulated the
effects of a variety of alternate harvest control rules for managing cisco and found that,
compared to the current control rule, the inclusion of biomass thresholds within CUT1
or CUT2 control rules can greatly decrease risk and increase yield and spawning
biomass at the end of the time series, at a cost of increased year-to-year variation in
yield. Finally, if constancy in year-to-year yield is held in the highest regard, we have
shown that constant catch control rules greatly outperform constant exploitation rate
control rules in terms of this performance metric for cisco in Thunder Bay, and the
inclusion of biomass thresholds within CCC2 rules decreases risk and increases Final SB
at little cost to yield and AAV.

116

APPENDIX

117

The SR function used to project recruitment in the case of a “boom” year was
derived using spawning biomass (mature female kg) and recruitment data from median
point estimates of the posterior distribution of the SCAA (Chapter 1). Given spawning
biomass and recruitment are on a 2 year lag (i.e. SCAA has recruitment in 1999 and
2000) we calculated spawning biomass in 1997 and 1998 by hindcasting from the
estimated 1999 stock abundance using natural mortality and harvest in 1997-1998. Due
to the scarcity of stock-recruitment data (either 7 or 4 data points for each recruitment
scenario), we placed an informative prior on the log alpha parameter based on the
family Salmonidae in Myers et al. (1999): log( ~) ~ N (1.43,0.05) . The recruitment
estimates then had to be standardized

~
Ry  Ry SSBRF 0 (1  e M )
~
Where R
are the standardized recruitments, R y are the recruitment medians from the
y

SCAA, SSBRF 0 is spawning biomass (mature female kg) produced per recruit in the
unfished condition, and M is the SCAA median female natural mortality estimate. The
Ricker model is then fit as

~
 Ry 
  log( ~)   * SBy  2  
log
 SB 
 y 2 
Where SB denotes to spawning biomass, calculated as the weight of mature females.
This model was run for 10 million iterations saving every 500th and burning in 2500 of
the final iterations. When used in the SPM we must back transform ~

118

~

e

SSBRF 0 (1  e  M )
Then recruitments for boom years are projected by:

Ry   * SBy 2 * e

  *SBy  2

* e

Supplemental files containing figures related to the sensitivity to different levels
of assessment error, productivity, and incorrectly calculating unfished spawning
biomass (Figures 2.5-2.28) solely include the 4yr recruitment scenario.

119

REFERENCES

120

REFERENCES

Beddington, J.R., Cooke, J.G., 1983. The potential yield of fish stocks. FAO Fisheries
Technical Paper 242.
Brooks, E. N., & Deroba, J. J. 2015. When “data” are not data: the pitfalls of post hoc
analyses that use stock assessment model output. Canadian Journal of Fisheries
and Aquatic Sciences, 72(4), 634-641.
Caddy, J.F., Mahon, R., 1995. Reference points for fisheries management. FAO Fisheries
Technical Paper Number 347.
Clark, W.G., 1991. Groundfish exploitation rates based on life history parameters. Can.
J. Fish. Aquat. Sci. 48, 734–750.
Clark, W.G., Hare, S.R., 2004. A conditional constant catch policy for managing the
Pacific halibut fishery. N. Am. J. Fish. Manage. 24, 106–113.
Deroba, J.J., Bence, J.R., 2008. A review of harvest policies: understanding relative
performance of control rules. Fish. Res. 94, 210–223.
Deroba, J. J., & Bence, J. R. 2012. Evaluating harvest control rules for lake whitefish in
the Great Lakes: accounting for variable life-history traits. Fisheries Research,
121, 88-103.
Desai, A.R., Austin, J.A., Bennington, V., McKinley, G.A., 2009. Stronger winds over a
large lake in response to a weakening air-to-lake temperature gradient.
Nat.Geosci. 2, 855–858.
Ebener, M.P., Stockwell, J.D., Yule, D.L., Gorman, O.T., Hrabik, T.R., Kinnunen, R.E.,
Mattes, W.P., Oyadomari, J.K., Schreiner, D.R., and Geving, S. 2008. Status of
cisco (Coregonus artedi) in Lake Superior during 1970-2006 and management
and research considerations. Ann Arbor, Michigan: Great Lakes Fishery
Commission, Lake Superior Technical Report 1.
Francis, R.C., 1993. Monte Carlo evaluation of risks for biological reference points used
in New Zealand fishery assessments. In: Smith, S.J., Hunt, J.J., Rivard, D. (Eds.),
Risk Evaluation and Biological Reference Points for Fisheries Management, vol.
120. Canadian Special Publication of Fisheries and Aquatic Sciences, pp. 221–
230.
Fujioka, J.T., Heifetz, J., Sigler, M.F., 1997. Choosing a harvest strategy for sablefish
based on uncertain life-history parameters. In: NOAA Technical Report NMFS
130 Biology and Management of Sablefish; Papers from the International
Symposium on the Biology and Management of Sablefish, Seattle, pp. 247–251.
121

Goodyear, C.P., 1993. Spawning stock biomass per recruit in fisheries management:
foundation and current use. In: Smith, S.J., Hunt, J.J., Rivard, D. (Eds.), Risk
Evaluation and Biological Reference Points for Fisheries Management, vol. 120.
Canadian Special Publication of Fisheries and Aquatic Sciences, pp. 67–81.
Haltuch, M. A., Punt, A. E., & Dorn, M. W. 2008. Evaluating alternative estimators of
fishery management reference points. Fisheries Research, 94(3), 290-303.
Haltuch, M. A., Punt, A. E., & Dorn, M. W. 2009. Evaluating the estimation of fishery
management reference points in a variable environment. Fisheries Research,
100(1), 42-56.
Hollowed, A. B., & Megrey, B. A. 1993. Evaluation of risks associated with application of
alternative harvest strategies for Gulf of Alaska walleye pollock. In Proceedings of
the international symposium on management strategies for exploited fish
populations (pp. 291-320).
Irwin, B.J., Wilberg, M.J., Bence, J.R., Jones, M.L., 2008. Evaluating alternative harvest
policies for yellow perch in southern Lake Michigan. Fish. Res. 94, 267–281.
Katsukawa, T., 2004. Numerical investigation of the optimal control rule for decision
making in fisheries management. Fish. Sci. 70, 123–131.
Leaman, B.M., 1993. Reference points for fisheries management: the western Canadian
experience. In: Smith, S.J., Hunt, J.J., Rivard, D. (Eds.), Risk Evaluation and
Biological Reference Points for Fisheries Management, vol. 120. Canadian Special
Publication of Fisheries and Aquatic Sciences, pp. 15–30.
Link, W.A. 1999. Modeling pattern in collections of parameters. J. Wildl. Manage. 63:
1017–1027.
MacCall, A.D., Klingbeil, R.A., Methot, R.D., 1985. Recent increased abundance and
potential productivity of Pacific mackerel (Scomber japonicas). CalCOFI Rep. 26,
119–129.
Mace, P. M. 1994. Relationships between common biological reference points used as
thresholds and targets of fisheries management strategies. Canadian Journal of
Fisheries and Aquatic Sciences, 51(1), 110-122.
Maunder, M.N., and Punt, A.E. 2013. A review of integrated analysis in fisheries stock
assessment. Fish. Res.142: 61–74.
Myers, R.A., Rosenberg, A.A., Mace, P.M., Barrowman, N., Restrepo, V.R., 1994. In
search of thresholds for recruitment overfishing. ICES J. Mar. Sci. 51, 191–205.

122

Myers, R.A., Brown, K.G. and Barrowman, N.J. 1999. The maximum reproductive rate
of fish at low population sizes. Canadian Journal of Fisheries and Aquatic
Sciences 56, 2404–2419.
Myers, J. T., Yule, D. L., Jones, M. L., Ahrenstorff, T. D., Hrabik, T. R., Claramunt, R.
M., Ebener, M. P., Berglund, E. K. 2015. Spatial synchrony in cisco recruitment.
Fisheries Research, 165, 11-21.
Punt, A.E., Dorn, M.W., Haltuch, M.A., 2008. Evaluation of threshold management
strategies for groundfish off the US west coast. Fish. Res. 94, 251–266.
Punt, A. 2010. Harvest control rules and fisheries management. In: Grafton, RQ,
Hilborn R., Squires D., Tait M., & Williams A.(2010) Handbook of Marine
Fisheries Conservation and Management. Oxford University Press, Oxford,
England, 582-594.
Quinn II, T. J., Fagen, R., & Zheng, J. 1990. Threshold management policies for
exploited populations. Canadian Journal of Fisheries and Aquatic Sciences,
47(10), 2016-2029.
Quinn, T. J., & Deriso, R. B. 1999. Quantitative fish dynamics. Oxford University Press.
Ricker, W. E. 1975. Computation and interpretation of biological statistics of fish
populations. Bull. Fish. Res. Bd. Can., 191, 1-382.
Stockwell, J.D., Yule, D.L., Gorman, O.T., Isaac, E.J., Moore, S.A., 2006. Evaluation of
bottom trawls as compared to acoustics to assess adult lake herring (Coregonus
artedi) abundance in Lake Superior. J. Gt. Lakes Res. 32, 280–292.
Stockwell, J.D., Ebener, M.P., Black, J.A., Gorman, O.T., Hrabik, T.R., Kinnunen,
R.E.,Mattes, W.P., Oyadomari, J.K., Schram, S.T., Schreiner, D.R., Seider, M.J.,
Sitar,S.P., Yule, D.L., 2009. A synthesis of cisco recovery in Lake Superior:
implicationsfor native fish rehabilitation in the Laurentian Great Lakes. N. Am. J.
Fish. Manag.29, 626–652.
Thompson, G.G., 1993. A proposal for a threshold stock size and maximum fishing
mortality rate. In: Smith, S.J., Hunt, J.J., Rivard, D. (Eds.), Risk Evaluation and
Biological Reference Points for Fisheries Management, vol. 120. Canadian Special
Publication of Fisheries and Aquatic Sciences, pp. 303–320.
Thorson, J.T., Cope, J.M., Kleisner, K.M., Samhouri, J.F., Shelton, A.O., and Ward, E.J.
2013. Giants’ shoulders 15 years later: lessons, challenges, and guidelines in
fisheries meta-analysis. Fish Fish.
Thorson, J. T., Jensen, O. P., & Zipkin, E. F. 2014. How variable is recruitment for
exploited marine fishes? A hierarchical model for testing life history theory.
Canadian Journal of Fisheries and Aquatic Sciences, 71(7), 973-983.

123

Walters, C.J. and S.J. Martell, 2004. Fisheries ecology and management. Princeton
University Press.
Yule, D. L., Stockwell, J. D., Cholwek, G. A., Evrard, L. M., Schram, S., Seider, M., &
Symbal, M. 2006. Evaluation of methods to estimate lake herring spawner
abundance in Lake Superior. Transactions of the American Fisheries Society,
135(3), 680-694.
Yule, D. L., J. D. Stockwell, J. A. Black, K. I. Cullis, G. A. Cholwek, and J. T. Myers.
2008. How systematic age underestimation can impede understanding of fish
population dynamics: lessons learned from a Lake Superior cisco stock.
Transactions of the American Fisheries Society 137:481–495.

124