Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Variation in pink shrimp populations off the west coast of Vancouver Island : oceanographic and trophic.. Martell, Steven James Dean 2002

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
ubc_2003-792390.pdf [ 10.2MB ]
Metadata
JSON: 1.0074868.json
JSON-LD: 1.0074868+ld.json
RDF/XML (Pretty): 1.0074868.xml
RDF/JSON: 1.0074868+rdf.json
Turtle: 1.0074868+rdf-turtle.txt
N-Triples: 1.0074868+rdf-ntriples.txt
Original Record: 1.0074868 +original-record.json
Full Text
1.0074868.txt
Citation
1.0074868.ris

Full Text

V A R I A T I O N IN P I N K S H R I M P P O P U L A T I O N S O F F T H E W E S T C O A S T O F V A N C O U V E R ISLAND: O C E A N O G R A P H I C A N D TROPHIC INTERACTIONS by  STEVEN JAMES DEAN MARTELL  B . S c , University of British Columbia, 1997 M . S c , University of British Columbia, 1999  A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Zoology)  We accept this thesis as conforming to the required standard  U N I V E R S I T Y OF BRITISH C O L U M B I A December 2002 ® Steven James Dean Martell, 2002  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department  or  by  his  or  her  representatives.  It  is  understood  that  copying  or  publication of this thesis for financial gain shall not be allowed without my written permission.  Department The University of British Columbia Vancouver, Canada  DE-6 (2/88)  Abstract This thesis examines the historical dynamics of pink shrimp populations off the west coast of Vancouver Island. Fisheries for pink shrimp in this area started in 1973 and stocks collapsed on two occasions. Shrimp populations continued to decline in the absence of fishing. Using a single species stock assessment model, I reconstructed time series on shrimp abundance, recruitment, and natural mortality in two statistical areas to investigate environmental and trophic interactions that would explain historical variability in recruitment and natural mortality. In the southern statistical area, there was evidence of recruitment over-fishing, and in both statistical areas juvenile survival rate was density dependent. There was no evidence that northern flowing coastal current advected shrimp larvae from southern fishing grounds to northern fishing grounds. Incorporating oceanographic indices as environmental correlates did not explain much of the residual variation in the stock-recruitment relationships. Members of the family Gadidae primarily consume shrimp, aud I examined if variation in estimated shrimp mortality was correlated with Pacific cod and Pacific hake abundance. No significant negative correlations were found between shrimp abundance arid predator abundance; however, the Pacific hake time series was incomplete. Contrary to other studies, shrimp abundance and Pacific cod abundance were positively correlated, as well as shrimp recruitment and juvenile Pacific cod abundance. To examine pink shrimp dynamics from an ecosystem perspective I constructed a simplified ecosystem model that represented the main trophic interactions and commercial fisheries off the west coast of Vancouver Island. A time series of relative primary production off the west coast of Vancouver Island was estimated by fitting the ecosystem model to time series data,. Including trophic interactions and variation in primary production explained an additional 29% and 15% of the residual sum of squares for pink shrimp biomass and landings, respectively. Predation mortality by Pacific cod on pink shrimp declines in the mid 1980's and late 1990's due to depleted Pacific cod stocks, however, total mortality for pink shrimp increases due to increased shrimp fishing. The west coast Vancouver Island ecosystem is sensitive to variation in primary production, and estimated primary production shifts downward around 1977. ii  Contents Abstract  ii  Contents  iii  List of Tables  vi  List of Figures  ix  Acknowledgements 1  2  Introduction  1  1.1  Pink Shrimp Fisheries  2  1.2  Life history of Pavdalus jordanl  5  1.3  Research Survey Program  6  1.4  Thesis Goals  9  Shrimp Recruitment and Mortality  11  2.1  Introduction  11  2.2  Stock Reconstruction Model  12  2.2.1  Population dynamics model  13  2.2.2  Observation Models  14  2.2.3  Parameter Estimation aud Uncertainty  23  2.3  2.4 3  xvi  Results  26  2.3.1  Estimating Abundance  26  2.3.2  Parameter Estimates and Uncertainty  32  2.3.3  Mortality  40  2.3.4  Annual Recruitment and Stock-Recruitment Relationships . . . .  51  Discussion  55  Recruitment and Oceanography  60  3.1  60  Introduction  iii  iv  3.1.1  Summer and Winter Current. C i r c u l a t i o n Patterns  3.1.2  Distribution  of Shrimp Relative to Coastal Currents and Shelf  B r e a k Currents 3.2  4  Methods  64  3.2.1  67  E n v i r o n m e n t a l Indices  Results  69  3.4  Discussion  72  S h r i m p P r e d a t i o n by Gadids  77  4.1  Introduction  77  4.1.1  Stock Structure of Pacific C o d  78  4.1.2  Stock Structure of Pacific Hake  80  4.3  4.4  E s t i m a t i n g Predator Biomass  83  4.2.1  Pacific C o d Assessment  83  4.2.2  Pacific Hake Assessment  88  Results  90  4.3.1  Pacific Cod Biomass  90  4.3.2  C o r r e l a t i o n w i t h Predator Abundance  93  Discussion  95  Ecosystem Perspectives  99  5.1  Introduction  99  5.2  Methods  101  5.2.1  E c o p a t h w i t h Ecositn  101  5.2.2  F i t t i n g E c o s i m To T i m e Series Data  107  5.2.3  Description of M o d e l Groups  110  5.2.4  T i m e Series D a t a for West Coast Vancouver Island  11G  5.3  5.4 6  62  3.3  4.2  5  61  Results  119  5.3.1  F i t t i n g E c o s i m to T i m e Series Data  5.3.2  Changes i n Shrimp Mortality  Discussion  General Discussion and Summary  119 126 130 134  Bibliography  139  A  S i m u l a t i o n S t u d y for S h r i m p Assessment M o d e l  145  B  L i n g c o d Assessment  151  U s i n g changes i n diet c o m p o s i t i o n to estimate v u l n e r a b i l i t y  List of Tables 2.1  M o d e l parameter descriptions and symbols used i n population dynamics model  2.2  14  Parameter descriptions and symbols used i n length frequency observation model  2.3  17  Number of shrimp fishing vessels reporting fishing activity from 1987 to 2000. A vessel reporting > 100 shots per year are considered full time, and is roughly equal to 20 days of fishing activity per year  20  2.4  E s t i m a t e d growth parameters for area 124 and 125 pink shrimp populations. 32  2.5  E s t i m a t e d model parameters and standard deviations for statistical areas 124 and 125  2.6  32  E s t i m a t e d recruitment anomalies («>,) and standard deviations for statistical areas 124 and 125  2.7  33  E s t i m a t e d deviations i n mortality (;/ ) and standard deviations for statist  tical areas 124 and 125 2.8  34  Reported fishing effort i n commercial logbooks for statistical areas 124 and 125  44  2.9  N a t u r a l mortality rate statistics for areas 124 and 125  49  3.1  C o m p a r i n g the effects of environmental indices on stock recruitment relationships for areas 124-5. loijL — log-likelihood, \  2  = probability of the of  the environmental index relative to base model, and the % of the variance in the residuals explained by the enviromental correlate  71  3.2  Cross correlations between environmental indices shown in Figure 3.3 . .  71  3.3  Correlations between lo<i(H./S), fjog(R) and environmental indices shown in F i g 3.3  72  vi  vii 3.4  C o m p a r i n g the effects of environmental indices on stock recruitment relationships for area. 124 spa.wners producing area. 125 recruits and vice versa,. logL = log-likelihood, \ relative to base  2  — probability of the of the environmental index  model, and the % of the variance i n the residuals explained  by the enviroment.al correlate  4.1  72  P u b l i s h e d reports of Pacific hake (Merluccius producf;us) predation  on pink  shrimp (Pandalus jordnni) and reported % occurrence, or % by weight of pink s h r i m p i n hake stomachs 4.2  83  Parameter symbols and description of parameters for the delay-difference model  4.3  84  Biomass Estimates for Pacific  hake; for L a Peruse bank region (Ware and  McFarlane, 1995), and for hake i n the C a n a d i a n Zone, (Dorn et a l , 1999). 4.4  89  C o r r e l a t i o n coefficients (r) between pink shrimp abundance, mortality and recruitment, versus predator abundance from 1975 to 2000. N corresponds to the number of data  pairs, and significant F statistics indicated by * =  0.1, * < 0.05,** < 0.01. and * " < 0.001 5.1  E c o p a t h input  94  parameters, arid estimated parameters (in bold) for west  coast Vancouver Island model  If2  5.2  Diet m a t r i x for E c o p a t h model  114  5.3  Squared residuals for No M o d e l (or deviations from the mean of the data), a model w i t h no trophic' interactions (assume that predation mortality and consumption are constant), a, model w i t h  trophic interactions where pri-  mary production is held constant over time, and a, model that includes estimated p r i m a r y production anomalies (see Figure 5.7). Positive reductions in the square residuals (right hand  column) between trophic interactions  only model and trophic interactions w i t h variable p r i m a r y  productivity  indicate a better 'fit' to the observed data 5.4  C o r r e l a t i o n coefficients (/•) differences from  124  between hake abundance time series a n d Z  single species model and ecosystem model ( D a t a shown  in Figure 5.12) A.l  True a n d estimated  129  parameter values for catch-at-length model.  True  values were used to simulate time series data, w i t h no observation errors, and estimated values w i t h out and natural mortality rate  penalty functions on recruitment anomalies, 148  viii B.l  C o m m e r c i a l catch statistics for lingcod off southwest Vancouver Island (Source: K i n g and Surry, 2000)  B.2  152  E s t i m a t e d growth parameters from length and weight-at-age data from Cass et al,(1990)  153  List of Figures 1.1  Statistical area definitions used for smooth pink shrimp trawl fisheries.  1.2  P i n k shrimp Pandalus jordani  .  landings i n statistical areas 124 a n d 125  from 1975 to 2000 1.3  4  5  Survey biomass estimates for Pandalus jordani i n areas 124 & 125 from 1973 to 2001. Note no surveys  wore conducted i n 1974. 1984, and 1986 in  b o t h areas and no surveys were conducted i n area 125 i n 1989 and 1991. Biomass estimated as mean biomass per area, swept by survey times t o t a l area surveyed 1.4  8  L o c a t i o n of research trawl surveys from 1973 to 2001. A r e a 121 is southern most area (light grey), area 123 is darker grey, 124 = darkest grey, and area 125 is the most northern  area (black circles). Highest shrimp densities  on the shelf area occur between 100 and 150m depth contours 2.1  8  Reconstructed biomass and observed biomass indices for statistical areas 124 and 125. Biomass shown  here represents vulnerable biomass prior to  the start of the research survey, solid line is predicted biomass from stock assessment models, and points are the biomass indices from s h r i m p t r a w l surveys 2.2  27  Residuals for area 124 and 125 between observed biomass indices a n d model predicted indices (lo(j(Z)). Note that positive residuals indicate observed biomass index is greater then predicted biomass  2.3  P o p u l a t i o n age structure for pink shrimp i n areas 124 and area 125. W h i t e shaded polygon are age-1  shrimp and black shaded polygons axe age-4  shrimp 2.4  28  29  Predicted (solid line) and observed  (histograms) length frequencies for  area 124. E s t i m a t e d mean length-at-age is lines. Note there predicted  represented by vertical dotted  are no length frequency data for 1983-84, and 1986, a n d  distributions are shown as a, straight line  ix  30  X  2.5  Predicted (solid line) and observed (bars) length frequencies for area 125. mean length-at-age is represented by vertical dotted lines. No  Estimated  samples were taken in 1984, 1986, 1989, and 1991 in area  length frequency  2.6  P a i r plot of 500 random from an  2.7  sample estimates (/?„, L\. L . p, A , A ,'>, h) taken t  A  2  MCMC sample of length 250.000 for area 124  P a i r plot of 500 random  sample estimates (/?„, L j ,  1.4, p,  35  Ai, A , 7 , l ) taken 2  h  from an M C M C sample of length 250,000 for area 125 2.8  36  M a r g i n a l posterior distributions for A in statistical areas 124 and 125. 2  frequency data are relatively uninformative  Note that observed length  about how standard deviation in length-at-age changes w i t h age 2.9  37  Sensitivity of M to <r'f[ and a" ,. Left column represent the mean estimate t  of natural mortality as af, increases from 0.1 to 1, where the mean (line w i t h points) is calculated over all a' . Right column represents the mean 2  M  estimate of M as a  2 w  increases from 0.8 to 8, where the mean is calculated  over all cr' . I n the left columns, dotted lines represent estimates of M 2  when al, = 0.8 and dashed line represents estimates of M when a  2 r  = 8,  and in the right column dotted lines for a'j = 0.1 and dashed line for o  2  f  M  = 1 38  2.10 Sensitivity of B„ to a\ and <T ,. Left column represent the mean R as <j 2  2  t  D  M  increases from 0.1 to 1, where the mean (line w i t h points) is calculated over a l l a\ . f  Right column  represents the mean estimate of R as cr ,, 2  0  where the mean is calculated over all cr' . Dotted  increases from 0.8 to 8,  2  lines i n left column represent estimates of R„ when cr , = 0.8 and dashed 2  line represents estimates of R when cr , = 8, i n the right c o l u m n dotted 2  ()  lines represent a  2 M  = 0.1 and dashed  line o\, = 1  39  2.11 Contour plots for objective function values ( -'V"«'-'( )) versus a L  C  2  and a , 2  M  values. A r e a 125 clearly shows 3 local maxima 2.12 E s t i m a t e d area over 125 were distributed.  which shrimp populations in statistical areas 124Area, estimates were calculated using commercial  logbook data, from 1987 to 2000; these data are shown i n Figure 2.14  Shrimp abundance for  represents the  . .  41  abundance and stock area from 1987 to 2000.  2.13 Relationship between shrimp  abundance (D, =  40  each statistical area is represented by the total  N „ . j W ) . Note that abundance shown i n figure 2.1 a  vulnerable biomass  42  2.14 Spatial distribution of commercial fishing effort from 1987 to 2000 on the West Coast of Vancouver Island. Total number of tows from 1987 to 2000 was 59,9(52 2.15 E s t i m a t e d selectivity curves for areas 124 and 125, where dotted vertical line represents the length at 50% vulnerability 2.16 E s t i m a t e d fishing mortality rates from stock assessment model (lines in bottom panels) for areas 124-125. Circles represent area swept estimates of fishing mortality rates (F = ~-) from commercial log book data, and ,J  these data were used in the likelihood criterion for estimating model p a rameters. Top panels are the residuals between area swept estimates of F and estimated fishing mortality rates in the stock assessment model.  . .  2.17 M a x i m u m likelihood estimates of natural mortality rates (line w i t h circles) from stock assessment model for areas 124 (top) a n d 125 (bottom) from 1975 to 2000. G r e y points are a 1.000 randomly chosen samples from the posterior distribution generated by the Metropolis Hastings A l g o r i t h m , and the box plots represent the median, inner quartile, and sampling range 2.18 M a x i m u m likelihood estimate of age-f recruits (line) for statistical areas 124 and 125 from 1975 to 2000 and associated uncertainty (box plots). G r e y points are a 1000 random samples from the posterior distribution, and box plots represent median, inner quartiles and range of posterior samples. Average annual recruitment for area 124 and 125 is 623 and 704 m i l l i o n shrimp, respectively 2.19 Relationship between carapace length arid proportion female for pink shrimp in areas 124-5 combined. B l a c k dots represent the fraction of female shrimp in 1mm length intervals sampled between 1975 to 2000. A random sample of N = 10,000 from If2,732 shrimp measured a n d sexed was chosen for estimating parameters for a logistic curve. Parameters were estimated using a least squares criterion, where the maximum likelihood estimate for  7 = 0.89 and  l  h  = 19.71  2.20 Spawner biomass (t) at t — 2 years versus age-1 recruits for statistical areas 124 (top-left) a n d 125 (bottom-left) and corresponding Beverton-Holt (solid line) a n d R.icker model (dotted line) fits. Spawning stock biomass is calculated using the surviving age-2+ biomass from 2 years previous (see text). Figures to the right are corresponding log recruits per spawner versus spawncrs. where the regression line represents density-dependent juvenile survival rate  xii 3.1  3.2  Tow positions of research trawl surveys from 1973 to 2001. Statistical areas 125, 124, 123 and 121 are shown i n black, to lighter shades of grey, respectively. Depth contours are LOOm, 200m and lOOOin Deviations i n annual recruits per unit of spawniug biomass (log(R, jS\) log(R /S )) t  t  63  i n statistical areas 124 and 125. T h e dotted line is 3 standard  deviations for the  distribution and  points  beyond this  line were considered  outliers. E x t r e m e l y high juvenile survival rates occurred area 124 and 125 in 1999 and 2000, respectively 3.3  Standardized  65  oceanographie indices  used to compare as  an  environmental  correlate w i t h shrimp recruitment in area 124 o n west coast Vancouver Island. S S T = sea surface temperature, S L H = sea. level height at Bamfield M a r i n e S t a t i o n , S O I = southern oscillation index, P D O = Pacific decadal oscillation index. T i m e series spans 1977-2001 3.4  68  Beverton-Holt and Ricker models fit to stock-recruitment data from statistical area 124 using a. normal likelihood ( E q . 3.5) and robust likelihood ( E q . 3.7) E a c h model line represents the mean stock-recruitment relationship w i t h different environmental indices included as correlates  3.5  70  Stock-recruitment models for area 124 (top row) arid area 125 (bottom row) using spawning  biomass in adjacent  area  as  a correlate. Closed cir-  cles represent observed data, line is the model fit to the data, and open circles represent  predicted recruits  w i t h correlate included. Left panels are  Beverton-Holt models, right are Ricker models 3.6  S u r v i v a l index (log(R /S )) t  t  73  for age-1 recruits i n area 124 and area 125.  Dotted line represents mean log(R,/S) 4.1  75  G r o u n d f i s h trawl catch per unit effort of Pacific cod off the west coast of Vancouver Island (scale units are log(e.p.u.e). respond to groundfish  management areas;  Areas 3 C and 3 D cor-  shrimp management areas are  also shown. M a p courtesy of A l a n Sinclair, DFO Pacific Region 4.2  79  D i s t r i b u t i o n of Pacific hake i n 1992, 1995, a n d 1998 along survey transect lines.  measured using hydro acoustic supplied by M a r k Saunders, D F O )  Abundance is  techniques. (Figure  backscattering 82  xiii 4.3  P r o p o r t i o n of Pacific hake in the Canadian zone versus mean JanuaryFebruary sea level height anomalies (top panel). Positive anomalies indicate higher than average sea level heights in Bamfield, B r i t i s h C o l u m b i a 1  (data provided by M a r k Saunders, U F O ) . Estimated hake biomass (bottom panel) for total hake stock (solid line) and hake present i n the C a n a d i a n zone (dashed line) based on the sea level height relationship 4.4  90  E s t i m a t e d biomass time series of age 2+ Pacific cod and fishing mortality rate on the west coast of Vancouver Island from 1956 to 2000 (top row). Circles represent the shrimp trawl survey index scaled to cod biomass and dotted lines represent ± 2 standard deviations. B o t t o m row from left to right are recruitment anomalies, catch residuals, and mean weight residuals from the delay difference model  4.5  91  Likelihood profiles for W C V I Pacific cod population parameters where assumed variance's for natural mortality rates a  2 M  — 0.38 and recruitment  compensation <r| = 0.2 (solid line) versus doubling these variance terms (dashed lines) 4.6  93  Relative abundance time series of P i n k shrimp, Pacific C o d , and Pacific hake series I and II from 1975 to 2000. Lower panels are scatter plots of predator abundance versus pink shrimp abundance  5.1  95  A n example of the relationship between predation m o r t a l i t y rate a n d predator biomass for different values of k  lr  T h e dark circle represents  the base predation mortality rate m ,-,•,<>, calculated using E c o p a t h input  parameters. A s the; value of k  lt  Bj.o,niijs)  decreases the slope of the line through  decreases and the maximum predation rate (v ) of predator j tJ  on prey i decreases 5.2  107  F o o d web and fisheries for west coast Vancouver Island ecosystem model. T h e size of the box is proportional to biomass of each group, and the w i d t h of the connecting line is proportional to the How. Position on the y-axis is the trophic level of the group  •5.3  116  C o m m e r c i a l landings (t km"" ) for west coast Vancouver Island from 1950 2  to 2000 (left column). T o t a l mortality rate estimates for pink shrimp and adult Pacific herring.  P i n k shrimp total mortality rate estimates were  calculated from combined fishing and natural mortality rate estimates in chapter 2. Total mortality rates for Pacific herring were inferred from age composition data,. B o t t o m panel is mean body weights for Pacific cod, where mean body weight, is calculated using: \V = B /N f  t  t  117  xiv 5.4  Relative  abundance time series used as observational data for the ecosys-  tem model. Shrimp and eulachon biomass estimates are combined for two  statistical areas, herring abundance was constructed using a. catch-at-age model, lingcod abundance from stock assessment  model fit to commercial  C P U E data, and Pacific cod from assessment i n Chapter 4. E u p h a u s i i d abundance (from 1967 to 1998) estimated from a, trophodynamic model by R o b i n s o n arid Ware (1999) 5.5  Predicted (points) and  118  observed (vertical bars) landings for pink shrimp,  herring, lingcod and Pacific cod stocks off the west coast Vancouver Island shelf.  F i s h i n g mortality rate estimates for a l l four species were derived  from single species assessment models and used as input data into E c o s i m . 120 5.6  Predicted  biomass time series (line) from E c o s i m a n d observed relative  abundance time series (circles). Note that relative abundance time series are assumed proportional to true abundance, where Y = qB t  5.7  122  t  E s t i m a t e d p r i m a r y production anomalies from f 950 to 2000; where a spline function was fit to 30  equally spaced data points. T h e use of the spline  function helps filter out high frequency variation that occurs i n fitting the model for years 5.8  having only one or few times series w i t h large anomalies. 123  C o m p a r i s o n between estimated annual mean primary production anomalies from E c o s i m , and calculations of D i a t o m production (g C obtained by  2  x  Robinson and Ware (1999) using a, model based on upwelling  (nitrogen), water 5.9  m yr~ )  Comparison  temperature, arid solar radiation  between estimated annual mean primary production and the  negative Pacific Decadal Oscillation ( P D O ) 5.10 A l l mortality  125  126  components from E c o s i m compared to Z from single species  model ( S S M ) . T h e area between each line is the t o t a l mortality due to each  predator or fishing, and the area between t o t a l mortality and fishing  mortality is the unexplained mortality (1 — EE) 5.11 Standardized differences  127  between Z estimates from Ecosim versus Z es-  timates from single specie's model in  Chapter 2. Positive values indicate  years where E c o s i m predicted higher mortality rates than estimated using the S S M  128  5.12 T i m e series of standardized indices of hake abundance and residuals between ZEWE 1 was  and ZSSM- All indices have a mean 0 and a — 1. Hake series  derived from sea level height indices and hake series II and III were  calculated from area, swept, trawl surveys (see  Chapter 4 for details). . . . 129  XV  A.1  A.2  Example of an observed length frequency distribution and predicted frequencies from stock assessment model. Observed frequencies were drawn at random from the true multinomial distribution  1 46  Simulated time series data from stock assessment model (thin black lines), and reconstructed time series (thick grey lines) w i t h no  observation errors.  B t is biomass over time. F t is annual fishing rate, R t is annual recruitment, and Mt is animal natural mortality rate A.3  147  B o x plots of parameter estimates from 250 simulation-reconstructions using the catch-at-length model w i t h auxiliary information about trends in fishing mortality  (right, box-whisker  plots), and no information about  trends in fishing mortality (left box wisker plots). P l o t s titles correspond to parameters listed in Table 1. Hashed line is the true parameter value.  14!)  A . 4 Parameter correlations for 250 simulation-reconstructions using the catchat-length model described in Chapter II. Variables names defined in Table  B. l B.2  A.1  150  Retrospective projection of lingeod biomass Reconstructed biomass, observed landings, fishing  153 mortality rates,  residu-  als in catches and mean body weights, and estimates of age-2 recruits for southwest Vancouver Island lingeod  154  Acknowledgements No thanks would  be enough  for all the encouragement and  not only supported me financially  when times were  support of  my parents, w h o m  tough but morally. Just the  thought,  accomplishment enough.  alone of making my parents proud is  To my grandfather, Robert E . Robinson, a former trapper in B r i t i s h C o l u m b i a . Y o u r historical accounts of  fishing  and trapping, along w i t h the first rainbow trout you helped  me catch, has led me into a, new direction of assessing impacts on natural resources. Y o u r memories w i l l reach far beyond your lifetime. M a n y thanks to my committee members. Daniel P a u l y and helped me fishing  develop alternative ideas  about ecosystem modeling, as well as the occasional  west coast of Vancouver Island, and Also, I  would like to  Shelton Harley,  helping me  digest the Pacific cod history on the  great pasta dinners.  acknowledge the alternative  exposure provided  by R a n s o m Myers,  group is  Jeff Hutchings, arid Boris W o r m (Dalhousie University). T h i s  leaps and bounds ahead in computing power due to a  was Shelton Harley light years  who  ahead of  provided some excellent  (thanks  Ware, for exposing me to a complex area of oceanography  and climate affects. A l a n Sine lair for  that is  Christensen who  trip for harbor seals and sight seeing tours at F l o r i d a State University  V i l l y ) . J o h n Dower and D a n  It  Villy  pointed  out  Mircosoft".  feedback  0  general ban on  that, a thesis could Ana, Parma,  who  be  products.  written on free  software  shared many great stories and  on manuscripts from this thesis. xvi  Microsoft^  J o n Sclmute, who  xvii finally explained the black magic behind the Metropolis-Hastings algorithm, and kept me laughing during P S A R C meetings. R i k T h o m p s o n , who sent a l l of his reprints on the same flay I  requested 'one'. Dave Founder, for getting me up to speed on the latest  in non-linear parameter estimation software, and a bad hangover i n M a d i s o n W i s c o n s i n .  weekends surfing i n H o n o l u l u , poi, bass fishing i n  James K i t c h e l l , for the wonderful UNDER.C,  Superica in Santa Barbra and the next job. T h a n k s also to R o b Bison, for  some extra $$ and an opportunity to give back to the local steelhead community. M a n y thanks to m y supervisor, Carl Walters, the m a n who defines conservation as "kill something else in order to save it", aptly encourages graduate students w i t h sound advice such as "don't F 9 < ? K U P " , arid has his own icon for evolutionary ecology.  In addition to witnessing the  return of ' B e t s y ' , I've also witnessed many other great  mystery's w i t h C a r l , such as full moon's in have to thank C a r l for the many  broad day light and belly-button lint. I also  great, adventures, here in Vancouver, F l o r i d a , .Japan,  Honolulu and the G r a n d Canyon.  C a r l and his wife S a n d r a were kind in providing  accommodations for D a w n and myself, until W i l l y moved  back home and split ice tea on  Sandra's new carpet (sorry W i l l , I had to tell someone!). I a m forever grateful for a l l the love and support provided by m y girlfriend D a w n Cooper. T h a n k you so much for all your patience and understanding. Special thanks to my lab R i k Duckworth, Lisa  mates, Sean C o x , Robert A l l i e n s , B o b Lessard , P a u l Higgins,  Thompson, Leonardo Hua.to, N a t h a n Taylor, and D a v i d O ' B r i e n .  A l l of your discussions,  computer code, and feedback were very valuable. I especially  appreciate the 'shack' at Whistler provided by the Decks Family. Special thanks to the folks at D . F . O . ,  especially J i m Boutillicr who encouraged me to  xviii adopt  my own methods tor examining shrimp fisheries,  computing resources.  A l s o H a i Nguyen  and comprehending the data. Rob  and  Norm  and  provided access to  Olsen for  data and  assistance w i t h databases,  K r o n l u n d and M a r c Saunders for important discus-  sions about the problems with assessing Pacific hake abundance. F i n a l l y , the to crew of the C C V W . E . Ricker, the platform for  shrimp surveys  and the summer feeding grounds  for myself. F i n a c i a l support for  this project  was  administered through J i m Boutillier.  provided by  Department of Fisheries and Oceans,  Chapter 1 Introduction "Nothing is more common in the management of resources than surprise, usually unpleasant surprise." ... D o n a l d L u d w i g , 1996.  T h e fishery for pink shrimp (Pandalus jordani) is relatively new i n comparison to other trawl fisheries for fin fish off the west coast of Vancouver Island ( W C V I ) shelf. D u r i n g the rapid development of this fishery i n the early 1970's large catches were landed, and i n just a few short years stocks collapsed despite careful monitoring by the Department of Fisheries and Oceans ( D F O ) . T h e W C V I shrimp t r a w l fisheries is one of a few fisheries i n B r i t i s h C o l u m b i a where fishery independent surveys have been i n place since the i n ception of the fishery. Despite careful monitoring, routine stock assessments and fisheries management plans, shrimp populations off W C V I fail to respond to efforts to prevent over-fishing. O n two occasions populations continued to decline i n the absence of fishing. Historically, shrimp recruitment to W C V I has been highly variable, and unpredictable. T h i s thesis examines the history of shrimp p o p u l a t i o n dynamics off the W C V I and examines environmental and trophic interactions that lead to variability i n recruitment and natural mortality. In this chapter I give a brief overview of the history of shrimp t r a w l fisheries i n the W C V I region, a brief life-history description of pink shrimp, and a description of the types of data collected by fishery independent surveys. In the following chapter, I 1  2 reconstruct historical abundance for two statistical areas and estimate a time series of recruitment and natural mortality. In chapter 3, I then examine stock-recruitment data estimated from chapter 2 and determine if environmental variables could be used to better explain recruitment variation. N a t u r a l mortality rates i n pink shrimp populations tend to be highly variable, and I examine the variability i n natural mortality rates i n chapter 4. In chapter 4, I reconstruct the abundance of two shrimp predators and search for correlations between shrimp abundance, natural mortality, and recruitment w i t h predator abundance. T h e last part of this thesis consists of looking at shrimp population dynamics from an ecosystem perspective. In chapter 5, I use a dynamic model to reconstruct the W C V I ecosystem, which includes other major fisheries such as herring and ground fish, to examine how changes i n b o t h predator abundance and variation i n p r i m a r y productivity have contributed to variability i n per capita rates for pink shrimp.  T h e last chapter  summarizes all of the findings and discusses the relative importance of variability i n p r i m a r y production w i t h respect to shrimp production on the W C V I shelf.  1.1  Pink Shrimp  Fisheries  T r a w l fisheries for smooth pink shrimp i n B r i t i s h C o l u m b i a began i n Stuart C h a n n e l i n the Strait of Georgia during the 1950s. E x p a n s i o n occurred i n B a r k l e y Sound during the 1960s. S m a l l inshore beam trawlers were used at the time, but during the early 1970s, as markets for B C shrimp started to expand, larger offshore otter trawls were used. A joint provincial/federal survey was conducted i n 1972 to assess the potential of developing an offshore fishery for pink shrimp. C a t c h rates and preliminary assessments for pink shrimp i n areas 124 and 125 (see Figure 1.1) looked promising, and i n 1973 the first commercial landings were reported from an offshore fishery. Over the next 5 years, the fishery continued to develop and additional commercial vessels joined the unrestricted  3 boom.  F r o m 1974-1978 the U S A fishing fleet also participated i n shrimp fisheries off  W C V I , and i n 1978 w i t h the introduction of extended j u r i s d i c t i o n , U S A participation in the fishery discontinued. U p u n t i l 1977, the C a n a d i a n fishing fleet was composed of 12 m e d i u m sized otter trawlers, and there were no management restrictions i n place. Offshore fishing activities were targeted i n two areas, the Tofino grounds (area 124), and N o o t k a grounds (area 125). B y 1977, catch rates for p i n k shrimp i n area 124 had declined dramatically, much of the fleet had moved to more distant grounds i n area 125. In the following two years, catch rates i n area 125 had also declined. B y the early 1980's shrimp fisheries on the West Coast of Vancouver Island h a d essentially shut themselves down as shrimp had apparently vanished from the offshore area.  P r i o r to the collapse, the Department of  Fisheries and Oceans ( D F O ) had implemented quotas for area 124 i n 1978 i n response to record catches taken i n 1976 and indications of recruitment failure i n the 1977 research surveys. A t the time, officials experienced concern about how quickly the fishery was developing. A s a pseudo-experiment, area 125 was left open w i t h no restrictions i n place; they (the D F O ) intended to observe how effort would respond to open access.  After  failing to reach the quota i n area 124, most fishers shifted to the n o r t h where shrimp appeared more abundant than i n the previous year. It appeared the newly discovered resource off the West Coast of Vancouver Island had been seriously depleted i n area 124 and a wave of fishing effort was m o v i n g northward. D F O continued annual shrimp surveys throughout the early 1980's and a few fishermen, who dare to leave the protected waters of B a r k l e y sound, ventured out to test offshore waters. There was little sign of any significant recovery i n shrimp stocks i n areas 124 and 125 u n t i l 1985. Biomass estimates from D F O surveys continued to decline from 1981 to 1983. N o surveys were conducted i n 1984 and 1986. In 1985 a significant increase  4  —i 130  1  1  128  126  —i 124  1  r 122  Longitude  Figure 1.1: Statistical area definitions used for smooth pink shrimp trawl fisheries. i n abundance was noted for area 124. C o m m e r c i a l fishing activities i n area 124 resumed immediately and small, below average, catches were landed. Landings i n area 124 continued to increase over the next two years and it appeared the stocks were rebuilding. P r i o r to the collapse i n the late 1970's the majority of the fleet were using otter trawls to land shrimp, and during the rebuilding of the stock i n the m i d 1980's the majority of the fleet consisted of smaller, more profitable, beam trawlers. Landings i n area 124 through the late 1980's were relatively stable, and a few fishermen ventured further n o r t h to explore the grounds i n area 125. T h i s exploration was often done later i n the season, as catch rates declined i n area 124. Despite efforts, no significant concentrations of shrimp were detected i n area 125 u n t i l 1991. D F O had not surveyed shrimp populations i n area 125 i n 1991, but i n the previous year, assessment reports indicted an increase i n abundance over the 1988 assessments. B y 1992, record landings were being recorded for the N o o t k a grounds, and even larger catches were obtained i n 1993, mean while, landings i n area 124  5  Area 124 c/> a> c o 05 TO  Area 125  o o in  OJ CD  CM  c  o o o  <n u>  O £Z  CO  I.... I  "i 1975  r 1985  o o in  o o m  LZ  I.  CO  . , 11  i 1975  1995  I 11  1  1 1985  Year  r  "i 1995  i  Year  Figure 1.2: P i n k shrimp Pandalus jordani 1975 to 2000.  landings i n statistical areas 124 and 125 from  (or Torino grounds) were beginning to decline (Figure 1.2). In less t h a n two decades of fishing,  it appeared as if history were repeating itself, a b o o m fishery moving northward  followed by collapse. In addition to fishing mortality influencing shrimp abundance on the west coast of Vancouver Island, there is also evidence that ocean climates influence the b o t t o m of the food chain, which i n t u r n affects dynamics of zooplankton and pelagic fish communities ( M c G o w a n et al. 1998; Robinson and W a r e 1999). In the last two decades there have been increases i n sea surface temperature anomalies. Others have defined this change i n ocean conditions as a "regime shift" (Beamish et al. 1999), and there is some evidence of changes i n production and community reorganization i n response to ocean climate regime shifts (Anderson and P i a t t 1999). There have been no directed studies of the influence of climate variation on pink shrimp p o p u l a t i o n on the west coast of Vancouver Island.  1.2  L i f e h i s t o r y o f Pandalus  S m o o t h pink shrimp Pandalus jordani  jordani  are protandrous hermaphrodites, m a t u r i n g first as  males and then changing sex to female around 2 years of age. Some individuals develop  6 directly into females and this phenomenon has been observed more frequently at low population densities (Hannah et al.  1995).  Breeding occurs i n the fall (November).  Females lay eggs, which remain attached to abdominal appendages u n t i l hatching i n M a r c h - M a y . N a u p l i i or free-swimming larvae spend 2-3 months i n the water column and then settle by late June or early July. B y October, as males, they are 60-70 m m i n length. B y the following a u t u m n males reach 85 m m . M o s t adults become females i n the spring of their second year, and grow to 100 m m by October. M a x i m u m life span is 4 years and m a x i m u m sizes are 123 m m for males and 140 m m for females (Butler 1980). There are two different species of ocean pink shrimp present i n the waters off the coast of Vancouver Island: the smooth pink shrimp Pandalus jordani shrimp Pandalus  and the spiny pink  borealis. T h e west coast Vancouver Island ( W C V I ) shrimp trawl fishery  targets the more abundant smooth pink shrimp and B C is the northern limit of fishable populations (Butler 1980; Boutillier et al. 1998). S m o o t h pink shrimp fisheries i n B C were first reported i n 1952, i n Stuart C h a n n e l , and have also been operating i n Barkley Sound since the 1960s.  1.3  Research Survey  Program  A n n u a l research survey programs have been conducted off the West Coast of Vancouver Island since 1973, w i t h exception of a few years i n the 1980s.  Survey programs are  intended to provide an unbiased index of abundance, age-structure, total mortality rates (Boutillier et al.  1998).  and estimates of  In the past, assessing abundance  each  year was carried out using spatial interpolation of catch per unit area swept. T h u s trawl surveys were designed to map out the distribution of shrimp stocks on each ground. P r i o r to 1987, it was assumed that two m a i n aggregations existed, one i n area 124 (Tofino grounds) and one i n area 125 (Nootka grounds).  Recent survey d a t a suggest that the  7 offshore distribution of shrimp may be a continuous distribution that lies between the 100150 m depth contour (see Figure 1.4). T h i s survey program continues to monitor shrimp abundance i n areas 124-125, and since 1994, the program has expanded to southward to include areas 123 and 121. In each of the fishing grounds, fisheries research vessels ( G . B . Reed prior to 1987, and W . E . Ricker post 1997) use an otter t r a w l to sample to grounds. E a c h tow is roughly 30 minutes i n duration and sweeps roughly 1.0 nautical miles (1852 meters) at a speed of 2.0 knots. T h e length of the footrope on the otter t r a w l is 11 meters; therefore the average area swept per tow is roughly 20,372 m . T h e highest densities of shrimp are usually 2  found between the 100m and 150m depth contour ( F i g . 1.4). For each tow, the catch is sorted to species, weighed, and discarded. E a c h tow that contains > 1kg of shrimp is sampled to assess mean body size, proportion of males, transitional (stage between male and female), and females, as well as length frequency distribution. In most years surveys are conducted i n the m o n t h of M a y , and a systematic design has been employed to map the distribution of populations along shore. Inshore and offshore boundaries were determined by fishing along each transect u n t i l shrimp catches were negligible. Spatial catch per unit area swept data interpolated using a bicubic spline (Boutillier et al. 1998) was used to generate an estimate of absolute biomass for each statistical area ( F i g . 1.3). Age composition samples collected over time are useful for monitoring changes i n t o t a l m o r t a l i t y rates, assuming samples reflect the true p o p u l a t i o n age structure.  If  it is not possible to age the species, then length-frequency d a t a are generally used to infer population age structure. T h e standard length measurement for shrimp is carapace length, and because shrimp and other crustaceans grow by m o l t i n g their exoskeletons there are no hard body parts that can be used for aging. T h e general approach for using length frequency data to infer population age structure is to estimate  proportions-at-age  8  Area 124  w  CD  c o w  W CO  E o  Area 125 o o o  o o o o  w  CD  c o  -»—'  o o o  w  00  o o o  •3-  CO CC  00  E o  ll.li.  o 1975  I  i—r 1985 Year  T  1995  in  I LI I i I  I.I  I—I—I—I—I— 1975  1985  1995  Year  Figure 1.3: Survey biomass estimates for Pandalus jordani i n areas 124 k. 125 from 1973 to 2001. Note no surveys were conducted i n 1974, 1984, and 1986 i n b o t h areas and no surveys were conducted i n area 125 i n 1989 and 1991. Biomass estimated as mean biomass per area swept by survey times t o t a l area surveyed.  Figure 1.4: L o c a t i o n of research t r a w l surveys from 1973 to 2001. A r e a 121 is southern most area (light grey), area 123 is darker grey, 124 = darkest grey, and area 125 is the most northern area (black circles). Highest shrimp densities on the shelf area occur between 100 and 150m depth contours.  9 and growth parameters by comparing predicted length frequencies w i t h observed data (Schnute and Fournier 1980). G i v e n estimates of proportion-at-age, t o t a l mortality over time can be inferred from changes i n age composition (i.e. e~  Zi  = "^' N  l +l t  )•  L e n g t h frequency data collected from research trawl surveys are used to assess age structure i n the population. Past assessments have used the length frequency d a t a to assess relative recruitment and changes i n natural mortality rates over time by converting the length frequency observations into proportion-at-age using probability transition matrixes. In this thesis I develop a catch-at-length model i n w h i c h length frequency data are treated as observations for estimating model parameters.  I a m particularly inter-  ested i n estimating historical recruitment and mortality, and determining if variability i n recruitment and natural mortality rates are a result of variability i n climate and predation. I then investigate the effects of climate variability on shrimp recruitment, and if variability i n natural mortality is a result of changes i n predation mortality. Finally, I examine shrimp dynamics from an ecosystem perspective using a d y n a m i c ecosystem model that includes complex trophic interactions, other commercial fisheries, and variation i n p r i m a r y production.  1.4  Thesis Goals  T h i s thesis uses a combination of single species assessment methods and ecosystem modeling to reconstruct information about the historical dynamics of the W C V I shrimp stock. Specifically, the aims of the analysis were to: 1. Reconstruct the history of the W C V I shrimp trawl fishery, estimate shrimp abundance, recruitment, and changes i n natural mortality rates over time.  2. Determine if the observed variation i n recruitment over time is related to ocean conditions, and if there is evidence for advection of shrimp larvae.  10 3. Determine if variation i n natural mortality rates are related to predation by gadoid predators.  4. E x p l a i n how much of the observed variation can be attributed to trophic interactions; b o t h changes i n predation over time and increases i n production associated w i t h favorable climate conditions.  Chapter 2 Shrimp Recruitment and Mortality 2.1  Introduction  In this chapter I simultaneously estimate recruitment and total mortality for pink shrimp in 2 statistical areas off the West Coast of Vancouver Island using an age-structured stock assessment model. The assessment procedure is broken down into 3 parts: a population dynamics model, observation models, and a statistical fitting procedure. The population dynamics model is used to represent unobserved states, or quantities we are interested in such as numbers-at-age in the population. The observation model translates output from the dynamic model into predictions of data types that we have observed, such as catch rates, or length-frequency data. Population model parameters describe rates of change from one state to the next, and observation model parameters relate unobserved states to observed states. Examples of observation model parameters are catchability coefficient, or size-selectivity parameters. Finally, a statistical fitting criterion describes how likely the data are given model parameters, and parameters are estimated by finding the most likely combination. In the following sections I describe the models and parameter estimation procedures used to estimate annual recruitment and total mortality. Both recruitment and total mortality estimates will be used in subsequent chapters of the thesis. There are 3 important vital rates that determine how abundance changes over time:  11  12 recruitment rates, natural mortality and growth. For short lived, small bodied organisms it is important to determine how each of these v i t a l rates change over time.  2.2  Stock Reconstruction  Model  T h e goal of reconstructing shrimp populations on W C V I is to quantify recruitment and recruitment variation, total mortality and mortality due to fishing. To reconstruct shrimp populations i n different statistical areas, I use a method that does not assume an underlying stock-recruitment relationship. In other words, I treat annual recruitment as an arbitrary historical variable to be estimated from the data.  T h e advantage to this  approach is one does not have to assume an underlying stock recruitment relationship. A disadvantage, however, is that the method also incorporates additional parameters that may be confounded w i t h other model parameters. T h e three parts of this assessment model are: 1) a population dynamics model, 2) observation models and 3) the statistical methods used to estimate model parameters and parameter uncertainty. Here model parameters represent various population parameters including annual recruitment, mortality, and growth. T h e statistical objective function, also referred to as a likelihood function, is s i m p l y a function that equates a probability measure for observations (the data) given model parameters. Such a function is essential for comparing alternative models, allows for the use of numerical routines to search for better parameter combinations and exploring parameter uncertainty. Here I use a general non-linear search algorithm (Davidson, Fletcher, Powell algorithm, Press et al. 1996) to maximize the statistical objective function. Parameter uncertainty is explored by numerically i n tegrating a posterior distribution, w i t h uniform prior probabilities for the parameters, using M a r k o v C h a i n M o n t e C a r l o ( M C M C ) methods ( G e l m a n et al. 1995).  13  2.2.1  Population dynamics  model  I assume p a r t i a l recruitment to the fishery occurs as early as age 1, and that a l l shrimp die before age 5. A description of all model parameters and symbols used i n the following equations are found i n Table 2.1. Numbers of shrimp at age a + 1 and t + 1 are calculated using: N  = K  a + h t + 1  ^  -  ™  (2.1)  Recruitment each year is simply the number of age 1 shrimp, which is calculated using:  N  = Re~  (2.2)  wt  ht  where R represents mean recruitment. A n n u a l recruitment anomalies were assumed to be log-normal a n d constrained to sum to zero. I also assume that natural mortality rates are variable over time and have a mean M = 0.6 (Boutillier et al. 1998; M a r t e l l et al. 2000) and have a log-normal distribution. V a r i a t i o n i n natural mortality was implemented by estimating a sequence of deviations v such that M = 0.6e .  I assume recruitment to  Ut  t  t  the fishing gear is a function of length, and to represent mean vulnerability-at-age I use a simple logistic curve of the form:  where L  a  is the mean length-at-age a. I n general, annual fishing mortality F can be t  calculated i n one of two ways: either using fishing effort m u l t i p l i e d by a catchability coefficient (conditioned on effort), or total annual catch divided by total available biomass present i n that year (conditioned on catch). T h e time series of commercial fishing effort for W C V I shrimp t r a w l fisheries is incomplete, so I chose to use annual catches to calculate annual fishing mortalities:  F =-log t  where  C <VB t  t  (2.4)  14 Table 2.1: M o d e l parameter descriptions and symbols used i n population dynamics model  Description T i m e subscript Age subscript Numbers-at-age at time t N a t u r a l mortality rate i n year t F i s h i n g rate Average annual recruitment Recruitment anomaly M o r t a l i t y anomaly V u l n e r a b i l i t y at age L e n g t h at 50% vulnerability M e a n length at age a M e a n weight at age Shape parameter for logistic function A n n u a l landings Vulnerable biomass where vulnerable biomass (VB ) t  Symbol t a N = 0.6e"« at  M  t  Ft R  w Vt t  Va L  h  L  a  w  a  7 c  t  VB  t  i n each year is the sum of products of the numbers of  shrimp i n each age class, mean vulnerability-at-age, and mean weight-at-age: 4 VB  t  =  N  ^ V  a  W  a  (2.5)  a=l Average weight-at-age was calculated using a length-weight relationship (Butler 1980) which is described i n the following section. There are 55 population parameters, including variable parameters that describe annual recruitment and natural mortality rates (R,wi,... ,w _i,vi,.. n  2.2.2  .,i/ _i,L ,j). n  h  Observation Models  In the annual research t r a w l surveys there are two types of observations, an index of density (i.e. catch per unit area swept, or C P A ) and an index of population age structure (as measured by length frequency data). T h e purpose of an observation model is to generate a set of predicted observations based on predictions from the population dynamics model. These observed and predicted states can then be used to estimate model parameters using various statistical criteria. I use three separate observation models in this assessment: a model for predicting relative abundance, a model for predicting  15 length frequency distributions, and a model for predicting relative fishing rates from area swept calculations. Generating a set of predicted length frequencies also requires a growth curve, and here I estimate growth curve parameters simultaneously w i t h other observation model parameters.  Predicting Relative Abundance D F O estimates absolute abundance on the shrimp t r a w l grounds using spatially interpolated catch per area swept data from research trawls. A time series of absolute abundance estimates are available from 1975 to present day (estimates are shown i n F i g . 1.3 on page 8). I treat these estimates as a relative abundance index. A b u n d a n c e indices generated by D F O represent the fraction of the population that is vulnerable to fishing gear. Furthermore, it is reasonable to assume that the survey abundance information is proportional to stock size (no gear saturation, survey samples whole stock range) and that measurement errors have a lognormal distribution. T h e age structured p o p u l a t i o n dynamics model used i n this assessment expresses abundance i n numbers, while abundance indices are expressed i n biomass units; therefore, the assessment model requires a length-weight relationship to convert numbers-at-length to weight.  I used a length weight relationship developed by B u t l e r (1980) to convert  population numbers to biomass. T h e following relationship was used for the conversion.  W = aL  where a = 0.002257082 and b = 2.609395004  b  U s i n g the assumptions stated above, the observed relative abundance information is related to the true vulnerable stock biomass v i a : Y  = qVB e , Vt  t  t  (where v is n o r m a l l y t  distributed w i t h a mean 0 and an unknown variance). Here q represents the survey vessel catchability coefficient, Y is the observed index of relative abundance i n year t, and t  VB  t  is the vulnerable biomass i n year t (a predicted quantity). In this case q is an u n k n o w n  16 nuisance parameter, and is simply the slope of the zero intercept linear regression between Y and VB . t  t  A simple trick to eliminate q from the estimation procedure is always use  the conditional m a x i m u m likelihood estimate of q i n the calculation of the likelihood function. T h i s estimate is easily calculated using the "Z statistic":  (2.6) U s i n g this formulation, the m a x i m u m likelihood estimator for log(q) is simply the average Z of a l l Z statistics (see Walters and L u d w i g 1994). Note here that the vulnerable biomass term i n E q u a t i o n 2.6 is adjusted downward to include any commercial catches that may have been taken prior to the survey.  T o fit model predictions to observed abundance  indices I use a likelihood approach and minimize the following negative log-likelihood:  (2.7)  where N is the number of years of relative abundance observations.  Length Frequency Observations G r o w t h curves are generally estimated from length-at-age data.  Here I use length fre-  quency data to estimate growth parameters because no tagging or age data are available to independently estimate growth parameters. To estimate growth parameters I use a method that simultaneously analyzes multiple sets of length frequency data (Fournier et al. 1990). Here I fit predicted length-frequency distributions to the observed lengthfrequencies by tracking cohorts over time.  Note that symbols used i n the following  equations are described i n Table 2.2. M e a n length-at-age is described by a von Bertalanffy growth curve, and variation i n length-at-age is assumed normal. Therefore, the length of an i n d i v i d u a l at age a plus the  17 Table 2.2: Parameter descriptions and symbols used in length frequency observation model. Description Symbol Subscript for age index a Subscript for cohort index i Subscript for length interval X Time at which sample was taken t Mean length-at-age a for cohort i L Asymptotic length Loo Growth coefficient k Age at zero length for cohort i to Deviation from mean length-at-age 5a Proportion-at-age 4>a Coefficient of Variation in Length-at-age CV Survey selectivity at length x x Number of Observations in year t N Probability of length x at age a Px,a Observed Frequencies Qx,t Predicted Frequencies Qx,t a  s t  time t at which the sample was taken relative to its birthday t is: Q  L  = L^  A+T  (1 - - <»+«-'»>) + 5 fc  e  a  (2.8)  Schnute and Fournier (1980) offer an alternative parameterization of equation 2.8 that is more stable for parameter estimation and two of the growth parameters (li and L A ) can be initialized from modes in a length frequency distribution:  L = h + ( L - h)lZ A-i P  A  A  (-) 2  9  where li represents the mode of the first age class, L A represents the mode of the oldest age class, and p — e~ is the growth coefficient. Equation 2.9 is very convenient because k  growth parameters, especially li, can be estimated just by looking for modes in the length frequency data. For example, an estimate of first age class mode (l±) can be determined by eye just by looking at the length frequencies shown in Figures 2.4 and 2.5. Parameters  18 for equation 2.9 can also be back transformed into the parameters for equation 2.8 using: _  L  -  A  "  0 0  k  -  A xP  1  l - ^ - i  =  to  l  -log(p)  ~\ (it^)  =  ai  l09  where ai is the first age class in the length frequency data. The standard deviation for length-at-age is assumed to be a function of mean length: a  a  = \ exp { A x  2  -1 + 2  1-p , o - l 1p ~\  (2.10)  A  where X% and A are parameters to be estimated and the term inside the square brackets 2  represents the length dependency of the standard deviations independently of the length parameters l\ and LA- A constraint was imposed on A > 0 to ensure that standard 2  deviation in length-at-age increases with age. Note that if A = 0 then the standard 2  deviations are length independent (i.e. o~ — Ai for all ages). a  In this growth model, it is assumed that all cohorts have the same asymptotic length (Loo) and growth coefficient (p). However, to accommodate variability in the time of year at which length frequency samples were taken, a in equation 2.9 and 2.10 took non-integer values. Assuming shrimp lengths are normally distributed in each age group, then the probability of a shrimp of age a being in length interval x — | to x + | , where d is an interval width, is  Px,a  =  1  7= / a V2n J -d a  exp •  (x -  Lf a  dx  (2.11)  x  Equation 2.11 simply predicts the probability of age a shrimp in a length interval of width d around mean length x. The probability of sampling a shrimp of any age of length a; is a function of the selectivity and proportions-at-age in the population such that P  xt  — Y^ ipatS p a  x  Xta  Note that this term is normalized to ensure total probabil-  ity = 1, and proportions-at-age (ift ) are determined by the age-structured population a  19 dynamics model. T h e predicted number of observed shrimp i n each length interval is simply the total number of shrimp observed, multiplied by the probability of being i n length interval x summed over ages  (Q ,  x t  =  Yl  a  resulting i n a predicted length  ^tP ,a) x  frequency distribution. T h e contribution of length frequency data to the objective function was calculated using a robust likelihood method, which is described i n Fournier et a l . (1990). T h e log-likelihood equation for length frequency samples is:  N  L  Q  =  n  -l/2j2/Z ^ (^ lo  27T  N  -^nlogfa)  + 0.1/n)]  t-l x=l t=l x=l n N  t=l i = l  where P  xt  and P  t=l  H  2fe,  +  0.1/n)  (2.12) +  0  0  1  are the observed proportions and probability of sampling a shrimp i n  xt  length interval x i n year t, respectively. T o reduce the influence of large sample sizes , say > 400, the s u m of squares term for each length frequency sample i n year t is multiplied by r , where r 2  t  2 t  = mm(Q ,400). xt  T h e variance term £  x t  of being i n length interval x was estimated as (1 — P t)P t, x  x  for the expected probability therefore, if the probability  P t — 0 then the variance term = 0. A small number (0.1/n) is added to the variance x  term to ensure £  l t  does not approach zero, and this also reduces the influence of length  observations where the expected probability is near zero. One of the more difficult problems i n such analysis is determining the number of age-classes present i n the length frequency d a t a (see Schnute and Fournier 1980); since  Pandalus jordani are relatively short-lived and grow quickly, I assume a m a x i m u m of 4 age-classes present i n a l l years.  20  Relative Fishing Rates F i s h i n g rates based on area swept estimates were calculated using logbook records from commercial vessels reporting fishing activity i n areas 121 to 125. A commercial log book program for shrimp t r a w l fisheries became a mandatory license requirement i n 1987. D a t a for the number of vessels fishing i n areas 121-125 are presented i n Table 2.3 and estimates of total area swept were calculated using information on i n d i v i d u a l t r a w l durations and vessel speed (distance covered) from the logbook database maintained at the Pacific Biological Station, N a n a i m o , B C . Table 2.3: N u m b e r of shrimp fishing vessels reporting fishing activity from 1987 to 2000. A vessel reporting > 100 shots per year are considered full time, and is roughly equal to 20 days of fishing activity per year. N u m b e r of VesYear Total num- Area Area Area Area sels reporting > ber of ves- 121 123 124 125 100 Shots Per sels Year 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000  84 128 64 62 133 105 92 106 248 190 117 150 105 89  2 0 0 0 2 1 3 0 30 13 7 8 1 6  25 37 22 26 33 21 18 35 103 85 47 73 63 54  55 81 37 35 68 44 35 43 66 46 39 34 25 20  2 10 4 1 23 32 33 27 39 34 16 23 8 8  21 23 18 16 31 14 19 38 65 57 18 26 18 6  T o calculate fishing rate based on area swept estimates, first I define fishing rate as total catch divided by t o t a l vulnerable biomass  (F = C/VB).  A s s u m i n g catch is  proportional to stock size and fishing effort, then catch can be expressed as where E is fishing effort, and c is the catchability coefficient. equation into the fishing rate equation  C = cEVB  Substituting this catch  F is now expressed as F — cEB/VB,  i.e.  F = cE.  21 C a t c h a b i l i t y is a term that describes the effectiveness of one unit of fishing effort and is a function of the total area swept by one unit of fishing effort (a) divided by the total area over which the stock is distributed (A), or c = a/A.  T h i s is essentially the same  definition of catchability as i n Paloheimo and Dickie (1964), assuming the probability of capturing a shrimp w i t h i n the swept area a is equal to 1. Substituting this definition of c into the calculation of fishing rate, fishing rate is now expressed as F =  aE/A.  U s i n g information from commercial logbooks, fishing rate for any grid cell location i is calculated as F , = aEi/A,  where a is the total area swept per shot, A is the total area of  the cell and Ei is number of tows shot i n cell i. T o calculate fishing rates over the entire ground we simply sum over a l l grid cells fished and estimate the t o t a l area over which the stock is distributed each year: F = ^  (2.13)  t  T h e area swept by one unit of fishing effort is variable due to differences i n fishing gears between vessels and length of tow. Unfortunately vessel specific gear information is not available prior to 1998, and measuring tow distance accurately is also difficult, or not always reported i n the logbook information. Therefore, I assume an average area a swept by a l l vessels. Since 1998, the average foot rope length of all vessels reporting  fishing  activity was 65.48 feet (or 10.77*10~ nautical mile) and the average tow distance was 3  3.385 n a u t i c a l miles. E a c h tow roughly sweeps a — 3.65* 1 0  - 2  square nautical miles. T h e  t o t a l area swept is approximated by the number of tows (E )  multiplied by the average  t  area swept (numerator i n E q u a t i o n 2.13). It is also possible that only a fraction of the total shrimp w i t h i n an area swept are removed by a unit of effort, that is some may escape through the mesh, or slip over the head rope. To correct for this, the numerator i n E q u a t i o n 2.13 would be multiplied by a fraction of the total shrimp that escape the harvesting process. Unfortunately there are no d a t a available to calculate what fraction of  22 shrimp escape the trawl; therefore, I assume that the probability of capture has remained constant over time, a n d treat the F estimates as a relative index of the true fishing rate. T h e total area over which the stock is distributed (A ) is also estimated from comt  mercial logbook data using similar methods to those proposed by W i n t e r s a n d Wheeler (1985). T o calculate stock area each year and i n each statistical area, first I overlayed a 1 square nautical mile grid over the W C V I . T h e grid allows each tow location to be recorded i n a corresponding row and column, also known as rastering. For each year, a n d i n each statistical area, start locations for each tow were plotted o n the raster map. T h e total area was calculated by adding up the area of a l l cells that h a d non-zero catch rates. For example, if 80 cells were fished, and 75 cells h a d positive catch rates, then the total area over which the stock was distributed was 75 square nautical miles. T h i s method assumes that no area harboring commercial quantities of shrimp was left untouched, or access to any area was restricted. There have been no spatial closures for shrimp fisheries on the W C V I since 1987, a n d a l l areas accessible t o fishing have remained open for at least part of each fishing season. T h e estimated fishing rates based on area swept calculations are assumed to be proportional to the true fishing rate, and errors i n estimating F from area swept methods are assumed to be lognormal. T h e likelihood of fishing rates calculated from logbook records is calculated using similar methods as the relative abundance index. Here F from the logbook records is transformed into a n F statistic: F = Ln(F /F ), as  t  where F  as  is the  observed swept area F estimate from logbook data, a n d F is the calculated fishing rate t  using catch d a t a a n d E q u a t i o n 2.4. T h e likelihood of observed fishing rates is then:  (2.14) (this is a log-normal likelihood, integrated over the variance of observed F  as  around  23 Ft)-  2.2.3  Parameter Estimation and Uncertainty  M o d e l parameters were estimated using a Bayesian approach, where the distribution of errors for relative abundance data ( E q . 2.7) and area swept estimates of fishing mortality rates ( E q . 2.14) are assumed log normal. B o t h of these likelihood terms are weighted by the m a x i m u m likelihood estimate of the variance. length frequency data was assumed n o r m a l ( E q .  T h e error d i s t r i b u t i o n for the  2.14) and weighted by the relative  variance multiplied by a variance scale factor ( r ) . Large sample sizes were unweighted 2  t  to a m a x i m u m sample size of 400, allowing more emphasis on smaller samples sizes (r  2  determines the overall scale of the variance, see Fournier et al. 1990).  Normally,  likelihoods for length frequency data are based on the m u l t i n o m i a l d i s t r i b u t i o n because P  xt  are not independent, normally distributed random variables. Fournier et al. (1990)  go on to show that expression 2.12 approximates a m i n i m u m x  2  estimator, and is a  suitable "self-scaling" likelihood for analyzing length frequency data. In addition to the three criterion listed i n the previous two sections (equations 2.7, 2.12, and 2.14) prior distributions were used to constrain recruitment anomalies and natural m o r t a l i t y rates. A total of 60 model parameters were estimated for each statistical area by m i n i m i z i n g the following objective function:  L = (L  Z  + L  Q  + L  F  + P +P ) 1  2  (2.15)  where P ' s are prior distributions (equations 2.16a-2.16b). Priors were used for recruitment anomaly parameters and n a t u r a l mortality rate parameters and here I assume  24 recruitment and natural mortality for shrimp > age-1 are independent.  (2.16a) (2.16b)  T h e relative contribution of equations 2.16a and 2.16b to the overall posterior distribution (equation 2.15) is determined by prior assumptions about variation i n recruitment (a^j) and natural m o r t a l i t y ( a ) . 2  kf  Recruitment anomalies were constrained to s u m to 0 i n  order to estimate mean recruitment. T h e variance t e r m for recruitment anomaly prior distribution was set <7 = 4.0, which is equivalent to a coefficient of variation = 0.54 i n 2  area 124 and 1.66 i n area 125. T h e mean natural mortality rate was positively correlated w i t h mean recruitment, and i n preliminary assessments it was determined that very little information exists about natural mortality rates (i.e. the m a r g i n a l posterior distribution for M was flat). In order to resolve potential confounding of model parameters, I fixed mean natural m o r t a l i t y rates for each area at 0.6 and estimate annual deviations T h e variance specified for deviations i n natural m o r t a l i t y was a  2 M  {u ). t  = 0.4, which roughly  corresponds to a coefficient of variation equal to 0.2 and 0.35 i n areas 124 and 125, respectively. T h e variance weights used i n equations 2.16a and 2.16b reflect prior assumptions about how variable recruitment and natural m o r t a l i t y rates are i n shrimp populations. A s cr  2  and a  2 M  approach 0, the model becomes a observation error only model w i t h  constant recruitment and a fixed natural mortality rate equal to 0.6. <T and a 2  2 M  Conversely, as  increase i n magnitude, more variation is attributed to process error and  observation errors become increasingly smaller. Therefore, I also examine how sensitive  V"" 71 — 1 Q 1  estimates of mean recruitment and mean natural mortality rates ( range of o  2  and o . 2  w  M  t  ^_  1  ) are over a  M o d e l parameters were estimated over a 10x10 grid of <7 and a 2  2 M  25 values ranging from 0.8-8 i n cr^ axis and 0.1-1 i n o  2  of R  and M versus o  2  Q  w  and  M  axis. I then plot the mean estimates  a . 2  M  A t o t a l of 60 model parameters ( 9 = R , li, L A , p, A A , 7 , lh, w , u ) were estimated A  by m i n i m i z i n g equation 2.15.  1 (  2  t  t  M i n i m i z a t i o n was carried out using a phased approach  supplied i n the A D M o d e l Builder software ( © O t t e r Research L t d . ) , where some of the parameters are held fixed during the initial phase of the m i n i m i z a t i o n process. Parameter estimation was broken down into two phases; i n the first phase, a l l parameters w i t h the exception of 7 and v were estimated, then i n the second phase 7 and u were included. t  t  D u r i n g the first phase of m i n i m i z a t i o n additional priors were used on mean length at age one ((l\ —12.2) /0.001, 2  where 12.2 is the mode of the first age class i n the length frequency  data sets) and mean fishing mortality rate ((/og (/ /0.4)) /0.001). These additional priors ,  2  t  effectively reduce the number of iterations required for the non-linear search routine to converge while assuming that natural mortality rates was constant (M = 0.6) and the overall mean fishing mortality rate was 0.4.  In the second phase I relax these priors  and the objective function is no longer sensitive to deviations from the priors.  Using  a phased approach to parameter estimation i n n o n linear models usually provides more stable estimates, especially i n the case of over-parameterized models (Dave Fournier, pers comm., Otter Research L t d . , Sidney B C ) . T o examine parameter uncertainty i n a Bayes framework, assuming uniform prior distributions for all parameters, I sampled from the joint posterior d i s t r i b u t i o n using a Metropolis-Hastings A l g o r i t h m ( G e l m a n et al. 1995). T h e Metropolis-Hastings Algor i t h m is a M a r k o v C h a i n M o n t e C a r l o ( M C M C ) sampling routine that works well for high-dimensional problems. T h e algorithm works by taking a random walk through the parameter vector ( 0 ) , sampling parameter combinations at rates leading to a stationary distribution that is the joint posterior d i s t r i b u t i o n  P(Q\Y). G i v e n a target distribution  26 that is dependent on Y data, the algorithm creates a sequence of random points ©i, 82 ... whose distribution eventually converges to the target distribution. T h e algorithm starts by randomly choosing a candidate point 9*, and evaluates this point relative to the previous point t — 1, P ( 9 * | 9 _ i ) . If ©* is more probable t h a n 9 _ ! then it is autot  t  matically accepted as a point i n the posterior distribution. If 9 * is less probable than  9 _i  then it is accepted w i t h a probability of  (  P(Q*\Y)/P(Q -i\Y). t  Otherwise the point  is rejected, and a new candidate point 9 * is selected (note that if 9 * is rejected then the previous point 9 _ i is included i n the posterior sample again). Candidate points are t  selected using a j u m p i n g rule  J (Q*\@ _i). For the j u m p i n g d i s t r i b u t i o n I use a multivarit  t  ate normal distribution, using the "Cholesky" decomposition of the covariance m a t r i x generated by the non-linear search routine (Press et al. 1996). After the algorithm converges 10,000 random samples from the final posterior d i s t r i b u t i o n are used to generate marginal posterior probability distributions for each parameter. T h e assessment model described here was also tested using fake d a t a generated by the model to ensure parameter estimates were unbiased, and to ensure the model was working correctly. T h i s process involved using the population dynamics models, and observation models as a simulation model to generate time series data. T h e n , these fake data were analyzed w i t h the assessment model. Results of the simulation studies are shown i n A p p e n d i x A .  2.3 2.3.1  Results Estimating  A b u n d a n c e  B o t h predicted and observed indices of shrimp populations i n areas 124 and 125 declined precipitously from 1975 to 1978-79 ( F i g . 2.1). D u r i n g the early 1980's the fishery had essentially shut itself down because shrimp were v i r t u a l l y absent i n area 124 and low i n  27 area 125. B y 1985, stocks increased abruptly i n area 124. N o evidence of a significant increase i n area 125 was observed u n t i l 1991. In 1991, research surveys were not conducted in area 125, however, commercial landings indicated shrimp were present on N o o t k a grounds.  Research surveys i n area 125 resumed i n 1992, and a four-fold increase i n  abundance was observed over the 1990 survey.  Area 124  n 1975  i  i  1  1985  Area 125  1  \—  1995  Year  h 1975  1  1  1  1985  1  r  1995  Year  Figure 2.1: Reconstructed biomass and observed biomass indices for statistical areas 124 and 125. Biomass shown here represents vulnerable biomass prior to the start of the research survey, solid line is predicted biomass from stock assessment models, and points are the biomass indices from shrimp t r a w l surveys.  M o d e l predictions of exploitable biomass are i n general agreement w i t h biomass trends generated by D F O ' s spatial interpolation of area swept information from research t r a w l surveys. Residuals for the  log(Z) statistics for areas 124 and 125 are shown i n Figure  2.2. In area 124, the largest deviation was observed i n 1983. T h e largest area swept biomass estimate i n area 124 was 11505 tonnes observed i n 1976 ( F i g . 1.3), whereas predicted vulnerable biomass for that year was 3,653 tonnes. Recall that D F O biomass indices were treated as a relative abundance index, and only trend information was used for estimating model parameters. In 1982-3, the end of the 82-83 E l N i n o , shrimp stocks in b o t h areas had collapsed, and abundance indices were less than model predictions for area 124 i n 1983 and area 125 i n 1982.  28 Area 124  Area 125  CC a>  a  cc  1  1  1  1  1  1975  1980  1985  1990  1995  1— 2000  1  1  1  1  1  1975  1980  1985  1990  1995  Year  1— 2000  Year  Figure 2.2: Residuals for area 124 and 125 between observed biomass indices and model predicted indices (log(Z)). Note that positive residuals indicate observed biomass index is greater then predicted biomass. T h e largest deviation i n biomass estimates was i n area 125 i n 1998 ( F i g . 2.2), however, I suspect biomass was grossly under-estimated i n that area i n 1998.  Survey biomass  indices were estimated at 28 tonnes a month after the commercial fishery commenced. Unrestricted, commercial fishers removed a total of 207 tonnes from area 124, or 739% of the available biomass estimated using trawl survey data.  T h i s increase i n biomass  could not be accounted for by growth of biomass over the summer months.  Catches  exceeded estimated biomass i n area 125 on 3 occasions: 1978, 1995 and 1998. In 1978 and 1995 catches exceeded biomass estimates by 104% and 129%, respectively. In area 124, catches exceeded biomass estimates on 4 occasions: 1987, 1988, 1990 and 1994 where catch divided by estimated biomass was 116%, 305%, 123% and 102%, respectively. T h e estimated p o p u l a t i o n age structure over time i n b o t h offshore populations is shown i n F i g . 2.3. Since 1998, total biomass i n b o t h areas has been increasing, especially recently i n area 125. A l t h o u g h population numbers have been declining i n area 124 since 1999, biomass continues to increase due to reductions i n t o t a l m o r t a l i t y that have allowed for increased growth of the remaining biomass.  29  Area 124  1975  1985 Year  1995  Area 125  1975  1985  1995  Year  Figure 2.3: P o p u l a t i o n age structure for pink shrimp i n areas 124 and area 125. W h i t e shaded polygon are age-1 shrimp and black shaded polygons are age-4 shrimp. F i t t i n g the m o d e l to length frequency data over time largely influences model predictions of proportions at age. Figures 2.4 and 2.5 (see pages 30 and 31) show model fits to observed length frequency data sets for areas 124 and 125. Predicted length frequencies i n area 124 are i n general agreement w i t h the observed data w i t h exception of f998, 2000, and 2001. I n 1998, sample sizes were relatively small, and the youngest age-class is not apparent i n the observed data, whereas model predictions suggest a relatively strong age-1 cohort. Conversely, in 2000, a strong age-1 cohort was observed i n the data, but not apparent i n m o d e l predictions. Predicted length frequencies i n 2001 suggest that growth may be greater t h a n historical averages (this is also apparent i n 1985 and 1987 where the mode of age-1 animals is > > the mean length-at-age); nonetheless, 2- and 3-year old age classes are dominant i n both observed d a t a and model predictions. In area 125, the estimated coefficient of variation i n length-at-age was greater t h a n that observed i n area 124 (Table 2.4), resulting i n wider distributions for each cohort. A g a i n , observed and predicted length frequencies i n 1998 contradict each other and contradictions are probably attributed to the very small sample size taken i n 1998 ( N — 131 shrimp). Also,  30 i n 2001 shrimp appear to be growing at rates higher than the historical average. 1975  10  15  1985  20  25  30  10  15  1994  20  25  30  10  15  20  25  30  Length (mm) Figure 2.4: Predicted (solid line) and observed (histograms) length frequencies for area 124. E s t i m a t e d mean length-at-age is represented by vertical dotted lines. Note there are no length frequency d a t a for 1983-84, and 1986, and predicted distributions are shown as a straight line.  31  Figure 2.5: P r e d i c t e d (solid line) and observed (bars) length frequencies for area 125. E s t i m a t e d mean length-at-age is represented by vertical dotted lines. N o length frequency samples were taken i n 1984, 1986, 1989, and 1991 i n area 125.  32 Table 2.4: E s t i m a t e d growth parameters for area 124 and 125 pink shrimp populations. A r e a Loo (mm) k o\ to 124 125  29.97 23.25  0.26 0.52  -0.65 0.05  °\  1.295 1.362  1.297 1.362  1.298 1.362  °\  1.299 1.362  E s t i m a t e d growth parameters are shown i n Table 2.4. V a r i a t i o n i n length-at-age was slightly larger i n area 125. S h r i m p i n area 125 appear to grow at a faster rate, however attain a smaller mean asymptotic length  2.3.2 A  Parameter  Estimates  (Loo).  and  Uncertainty  t o t a l of 60 model parameters were estimated.  M a x i m u m likelihood estimates and  standard deviations for a l l model parameters i n both statistical areas are shown i n Tables 2.5-2.7.  Parameters listed i n Tables 2.6 and 2.7, are nuisance parameters and when  included i n the estimation procedure explain additional variation i n the data. O m i t t i n g these nuisance parameters from the estimation routine is equivalent to fitting the model assuming only observation errors, constant recruitment, and constant n a t u r a l mortality. In general standard deviations for estimated model parameters are smaller for area 124 t h a n 125, largely due to larger sample sizes and greater contrast i n the d a t a between these two areas. Table 2.5: E s t i m a t e d model parameters and standard deviations for and 125. A r e a 124 A r e a 125 Parameter MLE Estimate Std MLE Estimate 6.28 0.16 6.13 log(Ro) 10.60 0.21 9.04 h 0.13 20.27 21.18 IA 0.77 0.03 0.59 P 1.30 0.05 1.36 Ai 0.00 0.02 0.00 A 0.54 0.03 0.63 7 18.10 0.34 14.81 k 2  statistical areas 124  Std 0.32 0.35 0.11 0.03 0.05 0.00 0.08 0.60  M e a n recruitment (R ) is similar between the two areas, yet a l l other parameters listed 0  in table 2.5 differ markedly between the two areas. In general, shrimp i n area 125 appear to grow faster, attain a smaller asymptotic length and have more variation i n length-atage than shrimp i n area 124. It also appears that larger shrimp are caught i n area 124 versus area 125, based on estimated selectivity curve parameters Ih and 7. Uncertainty Table 2.6: E s t i m a t e d recruitment anomalies (w ) and standard deviations for statistical areas 124 and 125. A r e a 124 A r e a 125 Parameter N a m e M L E Estimate S t d M L E Estimate Std 0.91 0.23 0.39 0.38 -0.27 0.40 0.37 -0.09 ^1976 w1977 -0.37 0.41 -0.98 0.38 0.03 0.39 -0.11 0.40 ™1978 -0.15 0.39 0.02 0.40 ™1979 -0.41 ^1980 0.44 -0.15 0.46 -0.40 0.42 0.46 -0.70 ™1981 -0.07 0.40 -0.54 0.42 ^1982 0.18 0.33 -0.45 0.46 ™1983 0.88 0.26 0.29 0.37 ^1984 0.66 0.25 -0.49 0.42 ^1985 -0.17 0.73 0.23 0.33 ^1986 0.75 0.24 -0.30 0.34 ^1987 0.68 0.24 0.22 0.32 ™1988 0.21 0.26 ^1989 0.44 0.31 -0.04 0.28 0.24 0.35 ^1990 0.21 0.27 0.44 0.32 W1991 0.39 0.25 0.19 0.32 ^1992 0.26 0.26 0.46 0.31 ™1993 -0.68 0.31 -0.47 0.38 ^1994 -0.38 0.27 -0.66 0.50 ^1995 -1.12 0.33 -0.89 0.83 ^1996 -0.90 0.34 -1.09 0.32 W1997 0.10 0.34 0.37 0.59 ^1998 0.21 0.37 0.39 1.16 ^1999 -1.39 0.71 2.62 0.44 ^2000 t  and apparent correlations between parameters listed i n table 2.5 are illustrated i n Figures 2.6 and 2.7.  These figures are sub-samples from the posterior d i s t r i b u t i o n that was  constructed using the M a r k o v C h a i n M o n t e Carlo algorithm. In b o t h statistical areas, strong positive correlations exist between growth parameters l\,p,  and the length at  50% vulnerability (Ih)- T h e growth coefficient p determines how quickly shrimp grow,  therefore, the observed length frequency data could come from a population that grows quickly and are harvested at a smaller size, or vise versa (note higher p values equal slower growth rates). There was a strong negative correlation between the selectivity  Table 2.7: E s t i m a t e d deviations i n mortality (v ) and standard deviations for statistical areas 124 and 125 A r e a 124 A r e a 125 Parameter N a m e M L E Estimate S t d M L E Estimate S t d -0.06 0.30 0.15 0.24 ^1975 -0.20 0.40 -0.01 0.45 ^1976 0.19 0.43 0.23 0.45 ^1977 0.28 0.44 -0.01 0.46 ^1978 0.35 0.44 0.03 0.46 ^1979 0.24 0.43 0.12 0.48 ^1980 0.26 0.44 0.72 0.56 ^1981 0.22 0.44 0.00 0.48 ^1982 -0.24 0.38 0.08 0.49 ^1983 -0.24 0.37 0.12 0.49 ^1984 -0.28 0.36 0.01 0.47 ^1985 -0.28 0.36 0.47 0.06 ^1986 -0.26 0.35 -0.10 0.44 ^1987 -0.16 0.35 -0.25 0.41 ^1988 -0.18 0.35 -0.31 0.40 ^1989 -0.19 0.34 -0.43 0.39 ^1990 0.02 0.35 -0.42 0.40 ^1991 0.09 0.35 -0.22 0.44 ^1992 -0.06 0.33 -0.30 0.45 ^1993 0.10 0.31 -0.29 0.46 ^1994 0.24 0.32 0.10 0.50 ^1995 0.11 0.40 0.71 0.86 ^1996 0.19 0.41 0.87 0.52 ^1997 0.08 0.39 -0.17 0.45 ^1998 0.37 -0.13 -0.16 0.43 ^1999 -0.10 0.39 -0.19 0.43 ^2000 t  parameters 7 and lu- A s 7 increases i n value, the slope of the selectivity curve through lh increases (i.e. large values of 7 infer 'knife-edge' selection).  A s the length at 50%  vulnerability increases 7 decreases and a higher fraction of smaller shrimp are vulnerable to the fishing gear. T h e correlation was much stronger i n area 125 t h a n i n area 124 ( F i g . 2.6 and 2.7). Parameters for standard deviation i n length-at-age were also negatively correlated  Figure 2.6: P a i r plot of 500 r a n d o m sample estimates (R , Li, L A , from an M C M C sample of length 250,000 for area 124. Q  p,X\,^2,7,h)  8.0 I  6.0  7.5  .I I I . 9.0  19.8  0.50 ,  20.4  .  0.60 J_J  0.000 I  l_  1.25  .  ,  ,  0.020  13.5  I I. I I I 1  1.45  ,  0.5  ,  0.8  Figure 2.7: P a i r plot of 500 r a n d o m sample estimates (R , L\, L , p, A from an M C M C sample of length 250,000 for area 125. Q  15.5 I 1 1 I I I  A  1 ;  A , 7 , lh) 2  37  with growth parameters (specifically L\ and p versus A i ) . In other words, faster growth corresponds with more variation in length at age. There was very little information in the length frequency data about how standard deviation in length changes with age. In each statistical area, the marginal posterior distributions for A have several modes (Fig. 2  2.8). Recall that A was constrained to be greater than 0, such that standard deviation in 2  length-at-age increases with age, and if A = 0 then standard deviation in length-at-age 2  is the same for all ages.  Area 124  Area 125  1  r~—i 0.010  0.000  r 0.020  ^2  Figure 2.8: Marginal posterior distributions for A in statistical areas 124 and 125. Note that observed length frequency data are relatively uninformative about how standard deviation in length-at-age changes with age. 2  Sensitivity to a and 2  M  In principle it is not possible to estimate the ratio of observation errors to process errors, and the results presented thus far are based on a priori assumptions about the variability in nuisance parameters (w and v , process error terms). To explore sensitivity of M and t  t  R to prior assumptions about afj and er^, model parameters were estimated over a grid 0  of a  2 M  and  values. The results of the sensitivity analysis are shown in Figures 2.9  and 2.10. The mean value of M and R in each column or row in the matrix of possible 0  variance combinations were plotted against the vectors a  2  and a . For example, the line 2  M  w  with points in upper left panel of Figure 2.9 represents the mean estimate of M over all values of G , and the dotted or dashed lines represents mean estimates of M for a single 2  W  38 value of  a. 2  Area 124  Area 124  s  0.2  0.4  0.6  0.8  1  1.0  2  3  4  5  6  7  8  Figure 2.9: Sensitivity of M to cr ^ and c r . Left column represent the mean estimate of natural mortality as a increases from 0.1 to 1, where the mean (line w i t h points) is calculated over a l l a . Right c o l u m n represents the mean estimate of M as increases from 0.8 to 8, where the mean is calculated over a l l c r . In the left columns, dotted lines represent estimates of M when cr = 0.8 and dashed line represents estimates of M when cr^ — 8, and i n the right column dotted lines for a = 0.1 and dashed line for a = 1 2  2  2  M  2  M  2  2  2  2  M  As u  2 M  M  increases i n value, estimates of mean natural mortality rates tend to increase  for b o t h statistical areas (left c o l u m n of F i g . 2.9). For example, i n area 124 the mean natural mortality rate estimate when a  2 M  = 0.1 is 0.6, and as <r^ increases to 1, estimates  of natural mortality rates increase to 0.86. A s estimates of cr increase, mean natural 2  mortality estimates tend to remain relatively constant, however, as o  2 M  increases mean  estimates of M increase (as shown by dashed lines i n right column of F i g . 2.9). Estimates of mean recruitment R were less sensitive to changes i n variance weights for Q  the two prior distributions. Increasing  (variance i n recruitment anomalies) results i n  39 Area 124  0.2  0.4  0.6  Area 124  0.8  Area 125  CC  Area 125  o  —I— 0.2  i—  1.0  0.4  Figure 2.10: Sensitivity of R to a and c r . Left column represent the mean R as a increases from 0.1 to 1, where the mean (line w i t h points) is calculated over a l l a . Right column represents the mean estimate of R as a increases from 0.8 to 8, where the mean is calculated over a l l a^. D o t t e d lines i n left column represent estimates of R when cr^ = 0.8 and dashed line represents estimates of R when cr = 8, i n the right column dotted lines represent a = 0.1 and dashed line a = 1 . 2  Q  2  M  Q  2  2  M  M  2  a  Q  2  0  2  2  M  rather uniform estimates of R  M  i n area 124 (Figure 2.10). In area 125, R  D  Q  declines slightly  as variation i n recruitment anomalies increase. Estimates of R were fairly consistent over a  the range of variances i n natural mortality rate anomalies. V a r i a t i o n i n estimates of R  Q  and M over ranges of a  2  and a  2 M  shown i n Figure 2.9  and 2.10 was a result of a likelihood surface that has several local m a x i m a across the range of variances explored. For some combinations or a  2  converges to very different solutions for parameters.  M  and a^, the estimation routine  In other words, the shape of the  likelihood surface is not smooth and changes w i t h different values of a  2  2.11).  M  and <r^ (Figure  40 Area 124  0.2  0.4  0.6 «  Area 125  0.8  1.0  0.2  0.4  0.6 «  2  0.8  1.0  2  Figure 2.11: Contour plots for objective function values ( - / ( )) versus a and a values. Area 125 clearly shows 3 local maxima. L  e  2.3.3  max  L  2  2  M  w  Mortality  Total mortality for pink shrimp in each statistical area was estimated from combined fishing mortality and natural mortality. Fishing mortality rates for the single species model were calculated using equation 2.4, where F is conditioned on observed catches from 1975 to 2000. I also constructed a time series offishingmortality rates based on area swept calculations (equation 2.13) from commercial log book records from 1987 to 2000. These estimates of F were assumed to be proportional to the truefishingmortality rates and were treated as observations (equation 2.14) in estimating model parameters. For the results on mortality, Ifirstreport on the result of the log book analysis, where area swept calculations were used to construct an independent estimate offishingmortality. Then I report estimatedfishingmortality rates from the stock assessment model and compare these estimates with area swept estimates of F. Finally I present time varying natural mortality rates and associated uncertainty. Area swept estimates from logbook records  Estimated stock areas for statistical areas 124-125 are shown in Figure 2.12, and the spatial distribution of commercialfishingeffort is shown in Figure 2.14. Note that all  41 Area 124  1988  1994  Area 125  2000  Year  1988  1994  2000  Year  Figure 2.12: E s t i m a t e d area over w h i c h shrimp populations i n statistical areas 124-125 were distributed. A r e a estimates were calculated using commercial logbook data from 1987 to 2000; these d a t a are shown i n Figure 2.14 log-book records for the W C V I island are plotted on Figure 2.14, this also includes areas 121 a n d 123. Records i n areas 121 and 123 were not used for estimating stock areas or fishing mortality rates because shrimp populations i n these areas were not examined. O n average, shrimp are distributed over a much larger area i n statistical area 124 than i n statistical area 125. Stock area varies considerably from 1987 to 2000 (see Figure 2.12). D u r i n g periods of low abundance stock area is much reduced a n d during periods of high abundance stock area increases substantially, and as abundance declines so too does stock area. For example i n area 124, relative abundance (note, estimated abundance is shown on Figure 2.13) declined from 1,556 tonnes i n 1994 to 371 tonnes i n 1999, a n d estimated stock area declined from 348 square nautical miles to 30 square nautical miles (Figure 2.13). T h e tow start locations shown i n Figure 2.14 represent a l l of the tows logged since 1987 i n a l l statistical areas (121-125).  P r i o r to 1996, less than 50% of a l l tows were  spatially referenced, a n d therefore, I a m assuming that reported tow locations accurately represent t o t a l area  fished.  T h e two principal fishing areas between 1987-2000 are i n  statistical areas 123 a n d 124 ( F i g . 2.14, and Table 2.3). E s t i m a t e d fishing rates based o n  42  Area 124  2000 4000 6000 8000 Abundance (tonnes)  7)  Area 125  0  1000 2000 3000 4000 Abundance (tonnes)  5000  Figure 2.13: Relationship between shrimp abundance a n d stock area from 1987 to 2000. S h r i m p abundance for each statistical area is represented by the total abundance (B = 2~2 N ,tW ). Note that abundance shown i n figure 2.1 represents the vulnerable biomass. t  a  a  a  area swept calculations for statistical areas 124-125 are shown i n Figure 2.16 and were calculated using the total effort d a t a presented i n Table 2.5. In b o t h areas, fishing rates tend to increase w i t h stock area, as more fishing effort is attracted to areas where shrimp are more abundant. E s t i m a t e d fishing rates based o n log book d a t a i n areas 124 a n d 125 range from 0 to 1.06. F i s h i n g effort i n 1998 declined dramatically i n area 124, a n d increased slightly i n area 125. D u r i n g the late 1990s, eulachon (Thaleichthys pacificus) bycatch i n shrimp trawl fisheries was becoming a concern, a n d commercial fishermen began to voluntarily adopt the use of by catch reduction devices ( B R D s ) . B y year 2000, 100% of industry had voluntarily adopted the use of B R D s . A l t h o u g h industry had taken a preventative measure to reduce bycatch, commercial fisheries i n area 123 a n d 121 were suspended early i n the 2000-fishing season for eulachon conservation concerns, otter trawlers were also restricted i n areas 124-125, and Queen Charlotte Sound was closed to a l l shrimp fisheries. It should be noted that fishing rates shown i n Figure 2.16 were estimated assuming that no commercial quantities of shrimp were left undetected, and that the total area over  4:5  v<^^^^^  r  r " « ^ ^ ^ ^ ^  1989^-  ' "-  :  1997^-  .  V^S.  T-  .  r  V^S.  r  r  ^ W T ^J  r  r Area 125 1.  Area 124  <  <f  ^^Qfl^gS Area 123  r  *C^^  r  Figure 2.14: Spatial distribution of commercial fishing effort from 1987 to 2000 on the West Coast of Vancouver Island. T o t a l number of tows from 1987 to 2000 was 59,962.  44  Table 2.8: Reported fishing effort i n commercial logbooks for statistical areas 124 and 125 Year N u m b e r of Shots T o t a l Effort Average Tow T i m e S T D T o w T i m e A r e a 124 (Hours) (Hours) 1987 8702 1.67 5198 0.96 1988 6011 9013 0.76 1.5 1989 3911 5380 1.38 0.65 1990 3761 5464 1.45 0.76 1991 6520 10368 1.59 0.72 1992 2377 3162 1.33 0.62 1993 1422 1947 1.37 0.58 1994 5414 8279 1.53 0.61 1995 3109 4204 1.35 0.73 1996 3497 5932 0.59 1.7 1997 2109 3951 1.87 0.46 1998 279 484 1.73 0.54 1999 107 214 2 0.53 2000 304 523 1.72 0.59 A r e a 125 (Hours) (Hours) 1987 9 12 1.37 0.22 1988 39 41 1.05 0.68 1989 15 18 1.19 0.45 1990 1 2 N/A 1.5 1991 1249 1491 0.47 1.19 1992 1806 2244 1.24 0.49 1993 4852 6336 0.39 1.31 1994 2833 0.37 3536 1.25 1995 3834 4991 0.43 1.3 1996 1932 2573 1.33 0.7 1997 483 689 1.43 0.59 1998 940 1829 1.95 0.53 1999 263 432 1.64 0.46 2000 146 167 0.44 1.14  45 which the stock was distributed is well represented by non-zero commercial catch rates. T h i s last assumption is p a r t i c u l a r l y bothersome, i n that less t h a n half of the logbook records report precise start locations for each set prior to 1996. Therefore, estimates of stock area prior to 1996 are based on assuming that the small fraction of commercial vessels reporting tow locations represent the total area searched by commercial vessels. Furthermore, it also assumes that start locations were accurately recorded. A few tows reported i n the m i d to late 1990's were located on Vancouver Island itself (representing errors i n log book records, see F i g . 2.14), and tows were eliminated from area swept calculations. Finally, some areas simply cannot be trawled due to rough b o t t o m , depth limitations, or strong currents, and if these areas are large a n d harbor shrimp then the total area over which the stock was distributed is under-estimated.  Fishing Mortality Ftom Stock Assessment T h e instantaneous fishing m o r t a l i t y rates for areas 124-5 generated from the stock reconstruction model are shown i n Figure 2.16 on page 47. F i s h i n g mortality rates were calculated using observed catch a n d vulnerable biomass (i.e.  F — —log(l — C/VB)), not  total biomass. E s t i m a t e d selectivity curves for the two areas are shown i n Figure 2.15. Note that the estimated length at 50% vulnerability is much larger i n area 124 (If, = 18.1) versus area 125 (Ik — 14.8). A s a consequence, estimated fishing mortality rates i n area 124 tend to be much higher i n comparison to area 124. T h e highest estimated  fishing  rate occurred i n area 124 i n 1997 (F = 5). S h r i m p abundance i n 1997 was estimated to be at an a l l time low i n area 124, w i t h a high proportion of p a r t i a l l y vulnerable age-1 shrimp i n the population. F r o m 1976 to 1980, fishing m o r t a l i t y rates were declining i n area 124, corresponding w i t h the decline i n abundance i n area 124. In 1978 m u c h of the fishing effort had moved north into area 125, generating high fishing mortalities.  B y 1979, catches were also  46  Figure 2.15: E s t i m a t e d selectivity curves for areas 124 and 125, where dotted vertical line represents the length at 50% vulnerability. declining i n area 125, as off-shore stocks had declined to a point where it was no longer economical to fish (J. Boutillier, Pacific Biological Station, Nanaimo, B C , Pers C o m m ) . Throughout most of the 1980s, the remaining shrimp industry focused much of its efforts on inshore waters (including Strait of Georgia and B a r k l e y Sound).  B y 1985 offshore  shrimp populations had begun to rebuild ( F i g . 2.1) and unrestricted offshore  fishing  activity had resumed i n area 124. F r o m 1987 to 1990, v i r t u a l l y all of the offshore fishing effort occurred i n area 124 (Table 2.8), and fishing mortalities exceeded the recommended target exploitation rate of 33%. D u r i n g this same time period, biomass i n area 125 was beginning to rebuild; however, significant fishing effort d i d not resume i n this area u n t i l 1991. In 1991, fishing mortality rates increased to 0.82 and 0.31 i n areas 124 and 125, respectively. F r o m 1991 to 1994 the majority of fishing effort i n offshore areas alternated between area 124 and 125, as well as the fishing mortality rates.  Estimated  fishing  mortality rates i n area 124 and 125 during 1993-97 had averaged near 0.95, again much higher t h a n the target fishing rate of 0.4. B y 1998-99, fishing rates had dropped below 0.1 per year as a result of declining  47  Area 124  Area 125  LO  »•—*  "tf  o  CO  CD C  CM  -  At CO  co  CD  -  05  03  o  T  1975  T  1985 Year  1995  CD  i 1975  —i  r 1985  r  1995  Year  Figure 2.16: E s t i m a t e d fishing mortality rates from stock assessment model (lines i n bott o m panels) for areas 124-125. Circles represent area swept estimates of fishing mortality rates (F = j^) from commercial log book data, and these data were used i n the likelihood criterion for estimating model parameters. T o p panels are the residuals between area swept estimates of F and estimated fishing mortality rates i n the stock assessment model. 9  abundance, and declining prices for shrimp. A l s o , i n 2000 eulachon conservation concerns severely restricted fishing activities i n the offshore areas, and the fishery was shut down prematurely, even though substantial shrimp catches were being landed early i n the season. T h e trends i n area swept fishing mortality correspond well w i t h trends i n  fishing  mortality rates estimated from the stock assessment model (see F i g . 2.16). Including area swept estimates of F i n the likelihood calculations d i d influence the overall trends i n predicted F ' s . N o formal sensitivity analysis was carried out to examine the influence of area swept F ' s on parameter estimates; however, when equation 2.14 was omitted from the likelihood, the uncertainty i n n a t u r a l mortality rates increases substantially.  48 There are a few major differences i n the trends of F from the two calculations, especially mid-late 1990's. Typically, area swept estimates of fishing mortality are less than the predicted F from the stock assessment model during the 1993-97 time period. T h i s time period corresponds w i t h fairly intensive fishing activity on the offshore populations. It is likely that the t o t a l area A was over-estimated for years when fishing effort was high t  (i.e. fishermen were searching for new grounds as catch rates declined i n the p r i n c i p a l grounds). Over-estimating the total stock area resulted i n reduced estimates of fishing mortality rates because the denominator increased i n equation 2.13. Nonetheless, the commercial logbook program does provide some auxiliary information about trends i n exploitation rates over time. Relative abundance indices from research cruises indicate stocks i n areas 124 and 125 had reached a low i n T 9 9 8 . T h e stock assessment model and length frequency data indicate that this collapse was i n part due to high exploitation rates i n these two areas, just prior to 1998 i n area 124 and during 1998 i n area 125. Increased natural mortality rates also played a role i n this collapse (see next section). Whereas, area swept estimates of F indicate fishing rates were relatively moderate i n 1998 i n comparison to previous years. In fact, F's were estimated to be declining i n b o t h areas since 1995 (Figure 2.16) w i t h the exception of area 125 i n 1998. C o m m e r c i a l landings i n area 125 i n 1998 were estimated at 206 tonnes, a shocking 283% more than the 72.9 tonnes of available shrimp estimated by D F O . Clearly the biomass index i n area 125 i n 1998 under-estimated total biomass, and using this index of relative abundance to tune the stock assessment model resulted i n a high estimate of mortality. F i s h i n g rate estimates for 1996 and 1997 i n area 124 were substantially higher from the stock assessment model t h a n estimates calculated using commercial logbook data. In 1996, commercial catches also exceeded area swept biomass estimates for area 124.  49  Estimated Natural Mortality Rates Deviations i n annual natural mortality (y ) rates were treated as variable, and were estit  mated using a l l of the observed information from research surveys and relative estimates for fishing mortality based on area swept calculations. Estimates of natural mortality rates and uncertainty are shown i n Figure 2.14. T h e estimated error distributions around the mean estimate of M were calculated using the Metropolis-Hastings algorithm, randomly sampling from the posterior distribution 50,000 times. In b o t h statistical areas, natural mortality rates have been highly variable over time, but more so i n area 125 (Table 2.9). In area 125 the highest estimate of M was 1.23 i n 1981, and i n area 124 the highest estimate of M was 0.85 i n 1979. Qualitatively, trends i n natural m o r t a l i t y rates over time i n the two areas seem to be positively correlated, however, large changes i n M appear to occur a year earlier i n area 124 than i n area 125 (Fig.  2.17). T h i s pattern was also observed just prior to the 1997-98 E l N i n o , where  increases i n natural mortality occurred first i n area 124. Table 2.9: N a t u r a l mortality rate statistics for areas 124 and 125 A r e a 124 A r e a 125 Average 061 063 Coefficient of V a r i a t i o n 0.21 0.35 Minimum 0.45 0.39 Maximum 0.85 1.23  Uncertainty i n estimates of M is higher i n area 125 than i n area 124, largely due to sample size differences between the two areas. Between 1983 and 1990, n a t u r a l mortality rates i n area 124 were relatively stable, whereas i n area 125, natural mortality rates were declining from 1984 to 1990. Research surveys were not conducted i n either of these areas i n 1984 and 1986; therefore, no information about natural mortality rates from length frequency d a t a is available during this time period. Furthermore, relative estimates of fishing  mortality rates, the other component of t o t a l mortality, are not available u n t i l  50  A r e a 124 in co  o in  •=  CO  °  OJ  d o d  i  i  i  1975  i  ,  i ' "i <» K i w # y M V * M T $  i  i  1978  i  i  i  1981  i  i  i  i  1984  i  i  1987  • '  V  i  i  i  1990  i  i  i  1993  i  i  i  1996  i  r  1999  A r e a 125 in CO  o co' in  ~  CO  o  cj  in  d o d  J-  i 1975  i  J-  i  a.  i 1978  J.  i  -I. -j.  X  i  i 1981  i  i  -J  -J.  i 1984  j.  i  _i_  i  j_.  i  i  1987  i  "J.  i 1990  i  i  J-  j-'  - "  i  i  i  1993  1  i 1996  i  - Xi  J.  i  r  1999  Year  Figure 2.17: M a x i m u m likelihood estimates of natural mortality rates (line w i t h circles) from stock assessment model for areas 124 (top) and 125 (bottom) from 1975 to 2000. G r e y points are a 1000 randomly chosen samples from the posterior distribution generated by the Metropolis Hastings A l g o r i t h m , and the box plots represent the median, inner quartile, and sampling range.  51 1987.  F i s h i n g mortality rates for area 125 i n 1987 were low, and other auxiliary d a t a  suggest that populations were increasing w i t h relatively low total mortality rates.  2.3.4  A n n u a l Recruitment and Stock-Recruitment  Relationships  T h e m a x i m u m likelihood estimate for age-1 recruits and associated distributions from the posterior are shown i n F i g u r e 2.18. A s the spawning biomass declined from 1975 to 1977 i n area 124, recruitment declined almost proportionally and continued to decline for an additional 4 years ( F i g . 2.18). It appears that b o t h poor recruitment and increases in natural mortality led to declines i n shrimp abundance i n b o t h areas during the early 1980's and recently i n 1996 ( F i g . 2.17 and F i g . 2.18). Recently, natural mortality rates in b o t h areas have been below average and a combination of growth and recruitment has led  to increases i n biomass i n areas 124 and 125. A r e a 125 has apparently had a large  increase i n annual recruitment i n 2000 and 2001 and biomass i n area 125 has increased more than 10-fold since 1998. Stock-recruitment relationships for b o t h statistical areas were constructed using est i m a t e d age-1 recruits, and spawning stock biomass. Recall that Pandalus  jordani  are  protandrous hermaphrodites, and the age at transition to female begins around age-2 w i t h a small fraction of age-l's and a higher fraction of age-3's being female (Hannah et al. 1995). A relationship between length and proportion female was developed using  length frequency and sexing information from research surveys (see Figure 2.19).  To  estimate the proportion of female shrimp i n the population, numbers-at-length were  multiplied by the proportions-at-length j that are female (i.e. Nf i  ema e  — Nj * Pj).  T h e logistic curve estimated i n Figure 2.19 roughly approximates the transition from male to female.  length-at-  T h e fraction of the total stock that contributes to  52  Year  Figure 2.18: M a x i m u m likelihood estimate of age-1 recruits (line) for statistical areas 124 and 125 from 1975 to 2000 and associated uncertainty (box plots). G r e y points are a 1000 random samples from the posterior distribution, and box plots represent median, inner quartiles and range of posterior samples. Average annual recruitment for area 124 and 125 is 623 and 704 million shrimp, respectively.  53 P  _ 0) E  Q_  a> o P=  o  1+exp(-y(l-l )) h  d o o 5  10  15  20  25  Carapace Length (mm)  Figure 2.19: Relationship between carapace length and proportion female for pink shrimp in areas 124-5 combined. Black dots represent the fraction of female shrimp in 1mm length intervals sampled between 1975 to 2000. A random sample of N = 10,000 from 112,732 shrimp measured and sexed was chosen for estimating parameters for a logistic curve. Parameters were estimated using a least squares criterion, where the maximum likelihood estimate for 7 — 0.89 and lh = 19.71. larval production consists of fertilized females in the fall that survive to the following March. During this period female shrimp are exposed to bothfishingand natural mortalities, and the spawning biomass must be corrected downward to represent individuals that survive to spawn. The spawning biomass at time t was calculated using SB = t  ( y\ Nt-iSjPjWA  -Fi-i-Mt-i-OAM ^  e  t  w h e r e t h e t e r m  i n  Jackets is the sum of fe-  male biomass over length intervals j at the beginning of the previous year, multiplied by one year of survival plus survival to March when larvae are released. The Sj term is the probability of picking a shrimp of length j from the population, N ~i are the number of t  shrimp in the previous year, and Pj is the proportion of female shrimp in length interval j. The spawning biomass calculation assumes that all females are fertilized, and eggs only hatch on surviving females. Spawning biomass was used as the independent variable in the stock-recruitment plots shown in Figure 2.20. In area 124 one recruitment anomaly stands out in 1999, where the spawning biomass was calculated to be just less than a ton, resulting in a rather significant recruitment at just over 6.5 million age-1 recruits. This anomaly represents nearly a 4-fold increase in  54  Area 124  Area 124  1 -< CO  "i  1  1  200  400  600  1  r  800 1000  T  1  1  1  0  200  400  600  i  r  800 1000  Spawners  Spawners  A r e a 125  Area 125  o o o co  100  200  300  Spawners  400  500  0  100  200  300  400  500  Spawners  Figure 2.20: Spawner biomass (t) at t — 2 years versus age-1 recruits for statistical areas 124 (top-left) and 125 (bottom-left) and corresponding Beverton-Holt (solid line) and Ricker model (dotted line) fits. Spawning stock biomass is calculated using the surviving age-2+ biomass from 2 years previous (see text). Figures to the right are corresponding log recruits per spawner versus spawners, where the regression line represents densitydependent juvenile survival rate.  55 average juvenile survival rate (see bottom-right panel i n F i g . 2.20). Recently recruits per spawner have returned back to the historical average i n area 124. Similar improvements in juvenile survival rates have also been observed recently i n area 125 (see bottom-right i n F i g . 2.20). F r o m 1999 to 2001, survival of age-1 recruits has increased nearly 1-3 times over the historical average i n area 125. In b o t h statistical areas, the Beverton-Holt model suggests very high compensation i n juvenile survival rates when spawner abundance is low. However, the Ricker model produces more conservative estimates of m a x i m u m recruits per spawner i n area 124. O m i t t i n g the extreme outlier i n 1999 i n area 124 reduced the estimate of m a x i m u m recruits per spawner from 87 to 32 i n the Beverton-Holt model, whereas o m i t t i n g the outlier i n area 125 i n 2000 increases the estimate from 9 to 20. For the Ricker models, m a x i m u m recruits per spawner declined from 9.48 and 13.8 to 6.8 a n d 6.2 i n areas 124 and 125, respectively when o m i t t i n g extreme outliers.  2.4  Discussion  In stock assessment modeling, it is most often assumed that n a t u r a l m o r t a l i t y rates are constant over time (Walters 2000). T h i s simplifying assumption is usually made to reduce the number of estimated parameters. N a t u r a l mortality is confounded w i t h other model parameters m a k i n g it difficult to estimate from time series data, and is usually estimated using independent methods.  In this assessment, estimating annual n a t u r a l mortality  rates was made possible by the inclusion of length frequency data and relative indices of fishing  mortality based on area swept calculations. Changes i n length frequencies over  time provides information about the t o t a l mortality rate (both fishing and predation). If an estimate of fishing m o r t a l i t y is available then the other component of t o t a l mortality can be estimated from length frequency data, or proportions-at-age d a t a (Sinclair 1998).  56 The use of length-based methods greatly improves our ability to infer proportionsat-age i n populations that are difficult to age (Schnute and Fournier 1980; Q u i n n et al. 1998); however, the methods often assume constant growth rates over time. Also, tracking cohorts over time is easier for species that recruit once a year. T h i s is necessary to track modes i n the length frequency data (Smith et al.  1998).  There is no evidence that  pink shrimp populations spawn more t h a n once a year on the west Coast of Vancouver Island; however, it appears that variability i n the time of year at which they spawn can also bias length-based approaches. T h e problem shows up i n assessment methods that directly estimate variance i n length-at-age rather t h a n estimating a single coefficient of variation for a l l age classes.  Large variance i n the time of year at which eggs are  spawned, or young of the year emerge,,is apparent when variance i n length-at-age is larger for younger cohorts t h a n it is for older cohorts (if age-specific variances are estimated directly). In effect the wide range i n b i r t h dates artificially inflate variance i n length-atage calculations. In this assessment procedure, the function used to estimate variance i n length-at-age was constrained such that variance was either independent of age ( A = 0), 2  or increases w i t h age ( A > 0). Estimates of A were near zero for b o t h statistical areas 2  2  and length frequency d a t a were uninformative about variation i n length-at-age, probably due to large over-laps i n length-at-age and variability i n hatching times between cohorts. Similar findings for N o r t h e r n s h r i m p (Pandalus  borealis) were documented by Fournier  et al. (1991). E s t i m a t i n g fishing m o r t a l i t y rates based on area swept calculations are subject to errors i n estimating t o t a l area swept by commercial fisheries, and estimation of the t o t a l area over which the stock is distributed. For W C V I shrimp t r a w l fishers, estimates of t o t a l area swept were carried out using an average area swept m u l t i p l i e d by the total number of tows. T h i s was done because necessary information about tow distance and net w i d t h  57 was lacking for years prior to 1998. Unfortunately this averaging does not account for changes i n gear types over time (e.g. larger, more efficient, nets being developed over time, or the effect of bycatch reduction devices). T h e area swept estimates also assumes that 100% of shrimp w i t h i n the area swept were removed, and if this is not the case then the effective area swept is simply reduced by the fraction of shrimp that escape the net (Hannah 1995; W i n t e r s and Wheeler 1985; Paloheimo and Dickie 1964). A d j u s t i n g the numerator i n E q u a t i o n 2.13 to accommodate changes i n fishing technology over time is t r i v i a l compared to dealing w i t h errors i n estimating the t o t a l area over w h i c h the stock was distributed (Swain and Sinclair 1994). T h e total area over which the stock was distributed was estimated as the total area in which commercial fleets detected shrimp. T h i s was done by imposing an arbitrary grid of 1 N M over a map of reported tow locations and simply summing the number of 2  cells w i t h non-zero catch rates. A m o n g all of the problems w i t h the assumptions behind this approach for calculating catchability, lies the question of grid cell size.  Imposing  a finer grid on the map increases the chances of more cells containing a positive catch rate, possibly decreasing the total area estimate if the number of new empty cells is greater than the number of new positive catch rate cells. A s the grid becomes infinitely small, the total area estimated w i l l eventually reach an asymptote.  U s i n g today's high  tech G P S information, fishing grounds w i t h i n 1 N M may still contain a high fraction 2  of untrawlable ground. Close examination of spatial d a t a for ground fish t r a w l fisheries has revealed that fishermen actively target productive grounds on spatial scales much smaller t h a n 1 N M (see F i g . 3 i n Walters and Bonfil 1999). T h e total area calculations 2  assume that fishing effort is randomly distributed w i t h i n each block, and violations i n this assumption lead to over-estimating total area. There is another problem w i t h using commercial catch records to estimate stock area,  58 which arises through changes in fishing fleet size. W h e n there are good incentives to go fishing (i.e. high prices for shrimp, or high catch rates), fishing effort generally increases. Increased effort also leads to increased competition on the fishing grounds, and creates strong incentives for fishermen to explore new areas. It is generally thought that the effort distribution follows an ideal free d i s t r i b u t i o n (Gillis et al. 1993; Walters and Bonfil 1999). In years where fishing effort is high and shrimp prices are high, more and more fishermen  are willing to fish i n areas w i t h lower catch rates. T h e y may even go explore  new grounds i n hopes of discovering potentially unexploited areas. W h e n using spatial distribution of fishing effort to estimate total area over which the stock is distributed, high effort years w i l l probably lead to an over-estimate of stock area.  Conversely, low  effort years w i l l tend to under-estimate stock area, because remaining fishermen w i l l tend to concentrate their fishing effort i n core areas where history has shown the area to be productive.  A l s o , remaining fishermen also tend to have higher catch rates ( H i l b o r n  1985). I t r y to correct for this problem by only including records w i t h non-zero catch . rates to justify the total area over which the stock is distributed. B u t a non-zero catch rate can consist of catching only 1 k g of shrimp i n a 3-hour tow (hardly w o r t h staying i n the area). I n d i v i d u a l fishermen w i l l choose their own criterion for deciding to stay i n an area, or move to a new area. T h i s criterion is probably economically driven; therefore, variability i n annual shrimp prices, and variability i n the cost of fishing, w i l l influence decisions to stay or move, or simply quit fishing altogether. Including independent estimates of fishing mortality rates i n the estimation procedure greatly reduces the uncertainty associated w i t h estimating fishing mortality rates (even if the F estimates are used as relative estimates, see A p p e n d i x A ) . Changes i n length composition over time provide relatively good estimates of total mortality rates  (Z).  E s t i m a t i n g the F component of Z was based on total catches and the likelihood of area  59 swept estimates of F from commercial log book records. For the purposes of this thesis, I assumed any errors that arise through violations i n assumptions about estimating fishing mortality rates from area swept estimates are independent, and therefore not biased if used as a relative index. In Figure 2.16, the trends i n fishing mortalities are somewhat similar i n b o t h statistical areas. B u t is it really fair to make such a comparison? To determine if these additional observations were very informative about fishing and natural m o r t a l i t y rates, the same estimation procedure should be carried out w i t h little or no weight on the area swept data.  I d i d not carry out this analysis i n the thesis due to limitations i n computing  time,  (it takes 18 hours to estimate uncertainty i n model parameters on a P e n t i u m  I V computer).  In the simulation study ( A p p e n d i x A ) , including estimates of  fishing  mortality rates d i d not introduce any bias, provided errors i n area swept estimates of F were not auto correlated. Is there evidence of a stock-recruitment relationship i n these two areas? T o answer this question I fit Beverton-Holt and Ricker models to the data for each statistical area using a least squares fitting criterion. In b o t h areas the best Beverton-Holt fit to the d a t a shown on the left panels of F i g . 2.20 is a straight line, even when o m i t t i n g extreme outliers when spawning stock is small. T h e simplest interpretation of this result would be that recruitment is totally independent of stock size w i t h i n each statistical area. For example, there is little evidence that recruitment i n statistical area 125 is a result of spawner abundance i n area 125. A n alternative explanation is that observed recruitment i n each of the statistical areas is a function of spawner abundance elsewhere, and larvae or juveniles dispersed to new areas (a question I examine i n chapter 3). T h e Ricker model i n area 124 provided a more conservative estimate of m a x i m u m recruits per spawner, and provided evidence of recruitment over-fishing.  Chapter 3 Recruitment and Oceanography "We can't do it! We have nothing. There isn't a single predictive explanation in the scientific literature today that will reliably predict recruitment based on any environmental variable you want to pick. The environment undoubtedly affects recruitment, but we can't predict it." ... Jeff Hutchings i n L A M E N T FOR A N OCEAN (Harris 1999).  3.1  Introduction  In this chapter I examine oceanographic processes that may directly influence shrimp recruitment on the West Coast of Vancouver Island ( W C V I ) . T h e physical oceanography on the W C V I is strongly connected to seasonal primary p r o d u c t i v i t y regimes (Thomson et al. 1989), which would affect recruitment moderately through trophic interactions. T h e prominent oceanographic feature off the west coast of Vancouver Island is a persistent coastal current that flows northward. T h i s coastal current is confined landward of the 100 m depth contour (Thomson et al. 1989; L e B l o n d e et a l . 1986).  A l t h o u g h adult  shrimp i n the offshore areas are found between the 100m and 200m depth contours, the persistent coastal current still may provide a conduit for shrimp larvae to disperse northward.  Also, year classes tend to segregate on the fishing grounds; for example,  smaller northern pink shrimp are generally found at the edges of each fishing ground in shallower waters (Koeller 2000). It is possible that young of the year shrimp settle out on the banks then slowly move offshore to the shelf break as they grow. If currents 60  61 advect larvae, then I expect recruitment indices generated i n Chapter 2 to be correlated w i t h spawning population to the south.  O r if upwelling is responsible for advecting  recruits offshore, resulting i n recruitment failure, then I expect the recruitment indices to be negatively correlated w i t h upwelling indices. If, however, recruitment anomalies are positively correlated w i t h upwelling events then I would suspect that increases i n p r i m a r y p r o d u c t i v i t y led to increased recruitment. Effects of variation i n p r i m a r y p r o d u c t i v i t y w i l l be examined further i n Chapter 5, i n the context of an ecosystem model that explicitly represents trophic interactions i n the region.  3.1.1  S u m m e r  and Winter  Current Circulation  Patterns  T h e persistent coastal current is primarily driven by surface outflow from J u a n de Fuca Strait d u r i n g the summer, and southeasterly winds during the winter.  M u c h of the  current's energy originates from runoff of low salinity waters through J u a n de F u c a Strait. The coastal current off the W C V I during summer months (April-September) tends to spread offshore along the broad continental shelf N o r t h of the J u a n de F u c a Strait and south of B a r k l e y Sound. A s the current travels n o r t h it forms a narrower band north of B a r k l e y Sound. Northwesterly winds predominate on the W C V I d u r i n g a n o r m a l summer year, and these winds tend to reverse surface currents (10-20 m depth), therefore the coastal current is usually found at depth during summer. A t depths beyond the 100 m depth contour, the shelf break region, ocean current flows i n a southeasterly direction, parallel w i t h the continental shelf d u r i n g summer months. In winter months the shelf break current flows i n a northerly direction. T h e transition zone between these two currents varies around the 100-200 m depth contour (Thomson et al. 1989; Foreman and T h o m s o n 1997). A  second major oceanographic feature that dominates circulation patterns off the  62 west coast of Vancouver Island is summer upwelling and winter downwelling. Upwelling occurs v i a E k m a n transport, where northwesterly summer winds transport surface waters offshore.  A s the surface waters move offshore, deeper, cooler water is drawn up from  below. Upwelling water masses are rich i n nutrients and enhance p r i m a r y productivity, except i n E l N i n o years when the thermocline is much deeper ( M y s a k 1986). In years w i t h strong summer northwesterly winds, p r i m a r y p r o d u c t i v i t y can be significantly enhanced off the W C V I (Mackas et a l . 1987; R o b i n s o n and Ware 1994; R o b i n s o n and Ware 1999).  3.1.2  D i s t r i b u t i o n o f S h r i m p R e l a t i v e t oC o a s t a l C u r r e n t s Shelf Break  and  Currents  S h r i m p populations i n statistical areas 124 and 125 tend to be concentrated between the 100 m and 200 m depth contour (Figure 3.1). Other shrimp populations exist near the entrance to J u a n de F u c a strait and directly offshore of B a r k l e y Sound; these occurrences are not considered i n this study.  I n areas 124 and 125, however, shrimp tend to be  aggregated between the 100 m and 200 m depth contour, even i n years of low abundance (see maps shown i n Chapter 2). These aggregations appear to be right i n the transition zone between the coastal currents and shelf-break currents. The location of shrimp populations i n areas 124 and 125 relative to counter-flowing currents creates a potential mechanism for dispersal of larvae either north or south. For example, a floating particle at the surface could be dispersed i n a southerly direction if strong northwesterly winds persisted. A persistent strong northwesterly w i n d w i l l induce upwelling and drive surface currents offshore and if particles on or near the surface move far enough offshore they w i l l be advected southward by southern moving shelf-break currents (Parrish et a l . 1981). A l s o , strong northwesterly winds can cause inshore surface currents to change direction, and particles floating i n the upper 10 m surface waters  63  128  127  126  125  124  123  122  Longitude  Figure 3.1: Tow positions of research t r a w l surveys from 1973 to 2001. Statistical areas 125, 124, 123 and 121 are shown i n black, to lighter shades of grey, respectively. D e p t h contours are 100m, 200m and 1000m  w i l l also be transported to the south. In deeper inshore waters (> 20 m depth) coastal currents always displace particles to the north. T h u s , organisms that undergo diel vert i c a l migrations ( D V M ) have a mechanism to be transported i n opposite directions over large scales. A t smaller spatial scales, D V M specialists are also subject to displacement resulting from the ebb and flood of t i d a l currents ( S m i t h et al. 2001). In this C h a p t e r I examine the effects of various oceanographic indices as an environmental correlate for shrimp recruitment. Using the stock recruitment d a t a from C h a p t e r 2, I compare residuals from Beverton-Holt and Ricker models w i t h models that include an environmental variable, using methods similar to C a p u t i et al. (1998). I also examine if spawning stock biomass i n adjacent statistical areas is correlated w i t h recruitment, which would be indicative of larvae dispersing to different areas.  64  3.2  Methods  A s s u m i n g annual recruitment is a function of spawning stock biomass R = f{S ) t  t  then  deviations from an underlying mean function can be represented by additive process error terms: R = f(S )ei  (3.1)  Wt+Vt  t  t  where the vector W represents some standardized environmental index w i t h a mean 0 t  and a\  Vt  = 1 and v is the residual process error due to factors other than W . Note that t  t  I assume deviations from the underlying stock recruitment relationship are log-normally distributed (see R i n a l d i a n d G a t t o 1978). T h e scaling term 7 represents the magnitude of the effect the correlate has o n the resulting recruitment. If g a m m a is set to zero, then predicted recruitment (R ) is strictly a function of spawning stock biomass (S ). T w o t  t  reasonable stock recruitment models for the mean density-dependent relationships are: Beverton-Holt form:  =  " ^ - l e™** 1 + a/bS t  (3.2)  2  Ricker form:  Rt = aS e- - ^ bSt  2+  (3.3)  Wt+Vt  t  Residuals between observed a n d predicted recruits are:  vt = log(Rt) - (log(f(S )) + -yW ) t  (3.4)  t  where f(S ) is either the Beverton-Holt model or Ricker model and R , S are the observed t  t  t  recruits and spawners estimated from Chapter 2. Parameters are estimated by minimizing the following negative log-likelihood:  \log(a ) v  + \log^)  + ^  ^  (3.5)  65 T h e four model parameters a, b, a and 7 were estimated using a Newton-type algov  r i t h m (Dennis and Schnabel 1983; Schnabel and Weiss 1985) i n R © . Setting 7 = 0 reduces the model to process error only and the number of parameters from 4 to 3. These are nested models ( H i l b o r n and M a n g e l 1997), and twice the difference i n log-likelihoods has a x distribution w i t h the degrees of freedom equal to the difference 2  i n the number of parameters (df = 1 i n this case). T h e x  2  probability is used to assess  whether including the environmental index significantly improves our ability to predict recruitment.  Robust Estimation It is c o m m o n that ecological data contain large errors, or outliers, that make it difficult to estimate underlying ecological processes using likelihood functions (Tsou and R o y a l l 1995; H i l b o r n a n d M a n g e l 1997; C h e n a n d Fournier 1999).  I n the case of assessing  density dependent survival rates for pink shrimp, there were two extreme outliers i n the recruit per spawner estimates ( F i g . 3.2). Including these outliers tends to over-estimate  Area 124  Area 125  -o Lr  1 1980  1990 Year  2000  1980  1990  1  r 2000  Year  Figure 3.2: Deviations i n annual recruits per unit of spawning biomass (log(R /S ) — log(R /S )) i n statistical areas 124 and 125. T h e dotted line is 3 standard deviations for the distribution a n d points beyond this line were considered outliers. E x t r e m e l y high juvenile survival rates occurred area 124 and 125 i n 1999 and 2000, respectively. t  t  t  t  m a x i m u m juvenile survival rate (or the slope at the origin of the stock  recruitment  66 plot) when using equation 3.5 as the statistical criterion. I modify equation 3.5 w i t h a weighting term (ip ) that reduces the influence of extreme outliers o n the estimation t  stock recruitment parameters. E q u a t i o n 3.5 is sensitive to outliers because the normal d i s t r i b u t i o n has a small 'tail' and the probability of an extreme recruitment event at low spawner density drops off quickly after a few standard deviations from the mean. Large outliers tend to contaminate likelihood functions, and the estimation procedure "over-fits" the anomaly. There are two approaches to reducing the influence of outliers: 1) changing the probability distribution to have 'fat-tails'(see C h e n and Fournier (1999) for examples) or 2) reduce the influence of the outlier i n the likelihood by unweighting the residual. Residuals v were weighted using a modified Tukey's "biwieght": t  (3.6) = 0  if |&t| > c,  where the modification ensures that 0 < ip < 1 (the original equation is available i n t  Press et a l . 1996, page 702). T h e error term b for each t observation represents the t  residuals given initial estimates of a, b, a, 7 (b — v — log(R ) — (log(f(S )) t  t  t  t  + ~fW )) and t  c is a constant specifying how sensitive ip is to small departures from the mean. A s c increases i n value tp approaches 1, and a l l data are weighted equally. Note that for t  normally distributed errors the o p t i m a l value of c is 6.0 (Press et al. 1996). T h e "robust" likelihood equation is written as:  (3.7)  If deviations (b ) are greater that c then the observation is omitted from the likelihood t  calculation (tpt = 0). T h e value of c for each statistical area was set equal to 6.  67  3.2.1  Environmental  Indices  I examine six different environmental correlates: sea surface temperature ( S S T ) , B a k u n upwelling indices i n M a y - J u n e , Sea Level Height ( S L H ) i n M a r c h - A p r i l , and annual S L H anomalies, the Southern Oscillation Index (SOI), and the Pacific Decadal Oscillation (PDO)  i n March-November. A l l indices were standardized between 1977 and 2001 by  subtracting the mean index value and dividing by the standard deviation, such that each index has mean 0 and <x= 1. Standardized indices are shown i n Figure 3.3.  M o n t h l y mean sea surface temperatures are available from various lighthouse stations  1  since the early 1900's. S S T d a t a were extracted from A m p h i t r i t e point, which is located at 48°33'N L a t i t u d e , 1 2 5 ° 1 9 ' W Longitude (near Ucluelet, B r i t i s h C o l u m b i a ) . W a r m periods are represented by positive anomalies (e.g.  1982-3 and 1997-98 E l Nino's i n  Figure 3.3). Upwelling indices were calculated based on the theory of E k m a n , where 2  surface currents generated by w i n d stress shear 45° to the right of the w i n d direction in  the N o r t h e r n hemisphere.  T h e transport of the water c o l u m n is 90° to the w i n d  direction over the E k m a n layer. T y p i c a l summer northwesterly winds off the west coast of Vancouver Island drive surface waters offshore and are replaced by cooler upwelling water from depths of 50-100m. Upwelling waters are generally more productive; sustain larger biomasses of phytoplankton, zoolplankton, and fish communities. T h e upwelling index is a measure of net offshore transport of water, and units are usually expressed i n m s 3  _ 1  per 100m of shoreline. Positive anomalies indicate upwelling or offshore movement  of the surface layer, and negative anomalies indicate downwelling. T h e pacific decadal oscillation ( P D O ) is actually a sea surface temperature index for  the N o r t h Pacific Ocean ( M a n t u a et al. 1997). T h e term P D O actually refers to J  2  Data Available at http://www.ios.bc.ca/ios/osap/data/lighthouse/bcsop.htm. Data at available from http://www.pfeg.noaa.gov/products/las.html  SST (Amphitrite Point)  I, Upwelling (Mayjun)  '1  111 • • • -  •  -•  SLH (MarApr)  .1  SLH (Annual Mean)  — — I  • PDO (MarNov)  —I— 1980  I  1985  —I— 1990  I  1995  Year  Figure 3.3: Standardized oceanographic indices used to compare as an environmental correlate with shrimp recruitment in area 124 on west coast Vancouver Island. SST = sea surface temperature, SLH = sea level height at Bamfield Marine Station, SOI = southern oscillation index, PDO = Pacific decadal oscillation index. Time series spans 1977-2001. how productivity oscillates between northern areas, such as Alaska, and southern areas, such as the the Columbia River in Oregon/Washington. Several observations over the last century have noted that when salmon returns in the North are good, returns in the south are poor (Francis and Hare 1994). This corresponds to a positive PDO. When the opposite condition exist, the PDO index is negative. This index also correlates well with other indices that oscillate at interdecadal scales (Mantua et al. 1997; Ware and Thomson 2000). An alternative hypothesis for explaining variation in juvenile survival (represented by the jw term) is recruitment advection, where larvae are dispersed by currents from t  one area to a second area. To examine this hypothesis I first determined if there were positive correlations between shrimp abundance in one statistical area versus another statistical area. If larvae disperse north, then there should be a positive correlation between spawners in area 124 and corresponding age-1 recruits in area 125. If larvae  69 disperse south then the opposite pattern should exist.  It is not expected, however,  that these relationships are linear so I also examine spawner abundance i n the adjacent statistical area as a correlate of age-1 recruitment.  Spawning biomass i n the adjacent  area was standardized to an index w i t h mean = 0 and a = 1 using the same procedure that was done for environmental indices.  3.3  Results  Under the assumption that shrimp i n areas 124 and 125 are separate independent stocks, the coefficient of variation (a ) 2  ranged from 1.06 to 1.5 i n area 124, and 0.6 to 1.04  i n area 125 for Beverton-Holt and Ricker models, respectively. T h e robust likelihood greatly influenced estimates of of the slope at the origin i n the stock recruitment curves (see example shown Figure 3.4). In area 124, m a x i m u m improvements i n juvenile survival rates at low spawner density are roughly 4 and 5 fold for Ricker and Beverton-Holt models, respectively. In contrast, using the n o r m a l likelihood, improvements i n juvenile survival rates are 3 orders of magnitude higher. In other words, if a l l observations are treated equally then there is no information i n the data about compensation i n survival at low spawning densities, and the best statistical stock-recruitment fit is a flat line w i t h no S _2 effect. In area 125, the d a t a were uninformative about a using the Beverton-Holt t  model, and the slope for the Ricker model is roughly 5.65 using the robust likelihood and 10.6 using the normal likelihood. In area 124, the M a y - J u n upwelling indices significantly reduced the likelihood over the Ricker model w i t h no environmental indicies (Table 3.1).  However, including up-  welling as a correlate only explains a small fraction (1.36%) of the residual variance. For the Beverton H o l t model, the upwelling index increased variance i n the residuals (7.25%). In area 125, sea surface temperatures (SST) and the Pacific Decadal Oscillation  70 Normal Likelihood  0  200  600  1000  Spawners  Robust Likelihood  0  200  600  1000  Spawners  Figure 3.4: Beverton-Holt and Ricker models fit to stock-recruitment data from statistical area 124 using a normal likelihood (Eq. 3.5) and robust likelihood (Eq. 3.7) Each model line represents the mean stock-recruitment relationship with different environmental indices included as correlates. (PDO) as environmental correlates explained  12.89 and 28.73% of  the residual varia-  tion, respectively. SST and the P D O are positively correlated (Table 3.2).  In contrast,  the P D O index as a correlate in area 124 explains little (Ricker Model) or no variation in the residuals (Beverton-Holt Model). Upwelling indices did not explain any of the recruitment variation in area 125. There is no strong evidence that juvenile survival or absolute recruitment is tightly connected with environmental indices shown in Figure 3.3. Juvenile survival (log(R/S)) was negatively correlated with SST and P D O in both statistical areas (Table 3.3). No strong correlations existed between recruitment in area 124 and any of the environmental indices. Recruitment in area 125 was positively correlated with P D O index and negative with SST. Recruitment in area 125 was negatively correlated with spawner abundance in area  124 (r — —0.19), and 2  also negative in the reciprocal comparison  (r = —0.12). 2  These  weak correlations do not suggest a linear relationship between recruitment and spawner abundance in adjacent areas.  71 Table 3.1: C o m p a r i n g the effects of environmental indices on stock recruitment relationships for areas 124-5. logL — log-likelihood, x = probability of the of the environmental index relative to base model, and the % of the variance i n the residuals explained by the enviromental correlate. Beverton H o l t M o d e l Ricker M o d e l Model logL logL x x 2  2  Base SST Bak SLH mSLH SOI PDO  19.48 19.33 17.99 19.42 19.47 19.14 19.24  0.59 0.08 0.73 0.89 0.41 0.49  Base SST Bak SLH mSLH SOI PDO  20.72 19.12 20.70 20.63 20.03 20.49 18.34  0.07 0.88 0.67 0.24 0.50 0.03  —  —  Table 3.2: Cross correlations between SST Bak SST 1 -0.112 Bak -0.112 1 SLH 0.365 0.263 mSLH 0.624 0.192 SOI -0.646 -0.127 0.727 -0.0614 PDO  2  A r e a 124 — 21.96 3.14 21.68 0.46 -7.25 19.61 0.03 2.63 21.86 0.65 0.11 21.90 0.74 -2.88 21.92 0.77 -1.07 21.91 0.75 A r e a 125 — — 26.02 13.94 24.51 0.08 -0.34 26.01 0.91 -0.39 26.01 0.88 4.50 25.59 0.36 5.38 25.85 0.56 21.10 22.32 •0.01  —  5.00 1.36 -0.54 0.25 -1.63 2.50 —  12.89 -0.20 -0.32 4.74 3.92 28.73  environmental indices shown i n Figure 3.3 SLH mSLH SOI PDO 0.365 0.624 -0.646 0.727 0.263 0.192 -0.127 -0.061 0.720 1 -0.208 0.403 0.720 1 -0.475 0.505 -0.208 -0.475 1 -0.524 0.403 0.505 -0.524 1  U s i n g stock-recruitment models w i t h spawning biomass i n the adjacent statistical area as a correlate explained very little of the observed variation ( F i g . 3.5). T h e parameter 7 is the slope for a linear relationship between the environmental correlates and recruitment anomalies (7 represents the standard deviation i n w ). t  Estimates of 7 using spawning  biomass i n the adjacent area were a l l < 0.1 and negative, except for the Beverton-Holt model i n area 125. In contrast, 7 = —0.38 and -0.25 (Ricker model) for the P D O and S S T indices i n area 125, respectively. I also examined if recruitment i n area 125 was a result of spawner abundance i n area  72 Table 3.3: Correlations between log(R/S), Log(R) and environmental indices shown i n F i g 3.3 Area SST Bak SLH (Mar-Apr) SLH (.lan-Dec) SOI PDO (Mar-Nov)  Log(R/S) Log(R/S) Log(R) Log(R)i  -0.294 -0.338 -0.092 -0.395  12i  l25  124  25  0.095 -0.130 0.183 -0.063  0.097 -0.018 0.165 0.001  0.107 0.013 -0.025 -0.217  0.220 0.362 -0.135 0.305  -0.304 -0.393 0.163 0.486  124, and vice versa, by estimating stock recruitment parameters using spawner abundance i n the adjacent area as the independent variable. N o significant relationships were found, w i t h or without environmental variables (Table 3.4). Table 3.4: C o m p a r i n g the effects of environmental indices on stock recruitment relationships for area 124 spawners producing area 125 recruits and vice versa. logL = log-likelihood, % = probability of the of the environmental index relative to base model, and the % of the variance i n the residuals explained by the enviromental correlate. 2  Beverton Holt M o d e l Model  logL  Base SST Bak SLH mSLH SOI PDO  19.97 19.59 19.65 19.71 19.97 19.90 19.89  Base SST Bak SLH mSLH SOI PDO  Area 19.25 18.24 19.21 19.25 18.34 19.23 17.80  x  2  A r e a 124  3.4  —  0.38 0,42 0.47 0.96 0.71 0.68  Ricker M o d e l  logL  X % <r Sp awners A r e a 125 Recruits — — — 19.87 2.61 19.09 0.21 5.60 2.40 19.52 0.40 2.50 2.14 19.79 0.68 0.51 0.00 19.79 0.68 0.52 0.43 19.81 0.72 0.41 0.78 19.63 0.48 1.87 2  2  125 Spawners A r e a 124 Recruits —  0.15 0.78 0.95 0.18 0.83 0.09  —  12.32 -0.69 -0.01 4.62 2.02 18.87  28.82 28.32 28.14 28.60 28.21 28.78 27.04  —  0.32 0.24 0.50 0.27 0.77 0.06  —  8.16 -1.24 0.15 2.86 2.56 20.45  Discussion  T h e methods carried out i n this chapter to investigate the potential affect of various environmental indices on larval survival of pink shrimp assume that recruitment is a function of spawner abundance. A s mentioned i n chapter 2, the stock recruitment d a t a  73 Area 124  0  100  200  300  Spawners  Area 124  400  500  0  100  200  300  400  500  Spawners  Figure 3.5: Stock-recruitment models for area 124 (top row) and area 125 (bottom row) using spawning biomass i n adjacent area as a correlate. Closed circles represent observed data, line is the model fit to the data, and open circles represent predicted recruits w i t h correlate included. Left panels are Beverton-Holt models, right are Ricker models. from 1977 to 2001 have a lot of scatter and the best statistical fit to the d a t a is a flat line for the Beverton H o l t model, especially area 125. Recall that i n C h a p t e r 2, I d i d not assume a stock recruitment relationship and age-1 recruits were arbitrary historical values to be estimated from the data. In order to impose a realistic curve through the data, I adopted a robust likelihood approach that reduces the influence of unlikely data.  The  likelihood function (equation 3.7) performs equally well when the d a t a are informative and contain no outliers (Tsou and R o y a l l 1995; C h e n and Fournier 1999).  Using an  ordinary least squares approach resulted i n extremely high estimates of m a x i m u m juvenile survival rate (a) for b o t h statistical areas, fitting either model. T h e robust likelihood removed the effect of the two large anomalies i n 1999 and 2000 i n areas 124 and 125  74 (ip = 0 for these two observations). V a r i a b i l i t y i n annual surface temperatures d i d explain some recruitment variability off the west coast of Vancouver Island, at least i n area 125. Similar findings were also reported by H a n n a h (1993), where shrimp recruitment was negatively correlated w i t h early fall sea surface temperatures, roughly six months after larval release. Cooler water temperatures were also a significant abiotic factor i n the survival and d i s t r i b u t i o n of pink Pandalus  borealis off the Scotian shelf i n Nova Scotia (Koeller 2000). Poor survival i n  warm water has also been demonstrated for Pandalus jordani larvae reared i n laboratories (Rothlisberg and M i l l e r 1983). There are two possible outcomes resulting from strong upwelling events off the W C V I . First, larvae could potentially be transported offshore into shelf break currents and settle out i n unfavorable habitats resulting i n recruitment failure. T h e largest upwelling event observed was i n 1997 ( F i g . 3.3) and larval survival i n b o t h areas was slightly below average ( F i g . 3.6).  T h e strongest M a y - J u n e downwelling anomaly occurred i n 1982.  L a r v a l survival rates were near average for area 124, and well below average i n area 125 (see Chapter 2, F i g 2.20 on page 54). There was no consistent information i n the upwelling indices that would suggest that poor recruitment is a result of offshore advection of larvae. T h e second possible outcome resulting from strong upwelling is an increase i n survival due to increases i n p r i m a r y production and food availability. Including the upwelling index as a correlate for area 124 d i d improve the overall fit to the stock recruitment data, however, it explains very little of the residual variation. A g a i n there were no consistent relationship for both statistical areas that would suggest increased p r i m a r y production has resulted i n increased survival. However, shrimp are p r i m a r i l y benthic detritus feeders, and it is possible that good upwelling conditions may provide good feeding opportunities for benthic juveniles i n subsequent years. For example i n the years following the strong  75 upwelling event i n 1997, juvenile survival rates have been increasing i n b o t h statistical areas ( F i g . 3.6).  Area 124  1980  1990 Year  Area 125  2000  1980  1990  2000  Year  Figure 3.6: Survival index (log(Rt/St)) for age-1 recruits i n area 124 and area 125. D o t t e d line represents mean log(R/S).  There was no relationship between shrimp recruitment or survival rate and sea level height. Lower sea level heights tend to result i n higher t h a n average recruitment i n Oregon (Hannah 1993). Higher sea level heights generally mean stronger coastal currents and H a n n a h speculates that i n years when sea levels are extremely high (i.e. 82-83 & 97-98 E l N i n o years) larvae are advected northward resulting i n recruitment failure i n OregonWashington waters. There is some evidence of an age-1 recruitment event i n 1984 (as indicated by positive recruitment anomalies i n Table 2.6), and stocks began to rebuild i n area 124 following this recruitment event. Four years later, recruitment begins to improve i n area 125, which also led to an rebuilding of the stock i n that area. U p u n t i l 1997, the pattern appears to have been a progressive wave of strong cohorts propagating up the coast, growing to maturity, spawning and reseeding grounds further north. Since 1999, two strong year classes have shown up i n b o t h statistical areas, especially i n area 125 in more recent years. D u r i n g this period, sea level height anomalies have been negative and apparent juvenile survival rates have been high. T h i s is also the case for the above average juvenile survival rate i n area 124. T h u s it appears that large negative sea-level  76 height anomalies are a better indicator of juvenile survival rate t h a n they are as an index of advection from southern spawning grounds to more northern areas. T h e t i m i n g of the spring transition period from winter circulation patterns to summer circulation patters may also be very important for survival of larvae (Thomson et al. 1989). T h e spring transition occurs when there is a shift from southerly to northerly winds along the west coast, and off the west coast of Vancouver Island the shift occurs sometime between early February and late A p r i l .  A drop i n coastal sea level height  is also an indicator of the spring transition, where circulation conditions switch from downwelling to upwelling, and the northward coastal current increases i n velocity due to increased runoff from the Fraser river. If the spring transition is delayed relative to the t i m i n g of larval release, then ocean conditions for shrimp larvae w i l l be less productive, or w i l l transport shrimp larvae onshore into areas such as Barkely Sound. H a n n a h (1993) also suggested that the t i m i n g of the spring transition was important i n determining the strength and distribution of year classes i n Oregon populations. There is no clear evidence that shrimp i n area 124 are a seed source for recruitment in area 125. There was no significant correlation between recruitment i n area 125 and spawner abundance i n area 124, i n fact the sign of the correlation suggest that the opposite. In either case, using spawning abundance i n the adjacent area as a correlate for recruitment i n that area d i d not explain much of the observed variation. A l l of the results suggest that recruitment is probably occurring at much larger spatial scales t h a n those examined here.  Chapter 4 Shrimp Predation by Gadids 4.1  Introduction  Several recent studies have indicated strong predator-prey interactions between members of the family Gadidae and P a n d a l i d shrimps (Berenboim et al. 2000; H a n n a h 1995; L i l l y et al. 2000).  A s s u m i n g that gadids, such as Pacific cod (Gadus macrocephalus) and  Pacific hake (Merluccius productus), are p r i n c i p a l predators of pink shrimp species ( L i l l y et al.  1998; Westerheim 1996), here I examine correlations between Pacific cod and  Pacific hake abundance versus pink shrimp abundance, mortality and recruitment off the west coast of Vancouver Island. Previous work i n Oregon has demonstrated a positive correlation between hake abundance and pink shrimp mortality (Hannah 1995). If these two predators are strongly influencing the dynamics of shrimp populations (top-down control) the following predictions should hold: predator biomass and shrimp biomass should be negatively correlated, predator biomass and shrimp mortality should be positively correlated and shrimp recruitment and predator biomass should be negatively correlated. T h e later prediction represents predation on juvenile shrimp, or that juvenile survival rates for pink shrimp are correlated w i t h predator abundance. If predation plays a m i n i m a l role i n the dynamics then no negative correlation should exist, or if there is a positive correlation between predator and prey biomass, then dynamics may 77  78 be dominated, by bottom-up control processes or donor-controlled. I use time series estimates of natural mortality rates, recruitment and aggregated abundance from Chapter 2 to examine correlations between predator abundance and these estimates. In this chapter I first describe stock structure of the two predators investigated. I then reconstruct a biomass time series for Pacific cod on the west coast of Vancouver Island. This reconstruction also includes description of a stock assessment model that is used again in Chapter V for reconstructing the history of fishing mortality rates for other species in the W C V I ecosystem. Finally I calculate correlations between predators and pink shrimp, and discuss the importance of these findings.  4.1.1  Stock Structure of Pacific Cod  Pacific cod are generally captured in shallow shelf areas, in depths less than 150m, with bottom trawls. On the B C coast, Pacific cod has been separated into four stocks for management purposes (Westerheim 1996). Pacific cod off the west coast of Vancouver Island are treated as a single stock. There is no evidence that this stock is migratory, with the exception of seasonal spawning aggregations that form on small reef areas. Catch rates are generally higher during the spawning season. Much of the Pacific cod catch is taken north of the entrance to Barkley Sound, near Amphitrite Point, especially during the late winter spawning season (Sinclair et al. 2001). Spatial distribution of catch rates since 1996 for Pacific cod are shown in Figure 4.1 using data from the mandatory observer program for the B C trawl fishery (each shot in this fishery has been monitored by independent observers since 1996). Throughout the rest of the year, Pacific cod are generally dispersed over the shelf area. In the past, fishing closures from January through March (to protect spawner stock) have been used as a conservation tool during periods of low productivity (Stocker et al. 1995).  79  Figure 4.1: Groundfish t r a w l catch per unit effort of Pacific cod off the west coast of Vancouver Island (scale units are log(c.p.u.e). Areas 3 C and 3D correspond to groundfish management areas; shrimp management areas are also shown. M a p courtesy of A l a n Sinclair, D F O Pacific Region. Stock assessment for Pacific cod has relied on commercial C P U E indices for estimating abundance i n a l l areas. There have been no directed, fishery independent surveys for Pacific cod on the west coast of Vancouver Island. T h e most recent D F O assessment document (Sinclair et al. 2001) included an index of abundance generated from s h r i m p trawl surveys.  A l t h o u g h the index was not used for the D F O assessment, the trend  information from s h r i m p surveys seems to correlate well w i t h commercial indices since 1981. Because I a m interested i n calculating predation mortality on s h r i m p by Pacific cod, I choose to include the shrimp survey index of Pacific cod i n this assessment.  80  4.1.2 Two  Stock Structure of Pacific  H a k e  distinct types of Pacific hake (Merluccius productus) populations exist along the  West Coast of N o r t h A m e r i c a , a migratory coastal stock, and several isolated nonmigratory inshore stocks that are genetically distinct from coastal stocks (Utter 1971). Inshore stocks are generally confined to major inlets and fjords, including the Strait of Georgia and Puget Sound. These stocks are not considered potential predators of pink shrimp populations on the West Coast of Vancouver Island shelf and are subsequently ignored for this assessment. T h e coastal hake population is a predator on pink shrimp populations during the summer months off the West Coast of Vancouver Island (Rexstad and P i k i t c h 1986). After spawning during the winter months i n central California, adult coastal-stock hake move onshore i n the spring to feed and then migrate northward. B y m i d summer, adult hake form large midwater aggregations near the continental slope off Vancouver Island and feed on euphausiids, other pelagic schooling fishes, and P a n d a l i d shrimps ( D o r n et al. 1999; Mackas et al. 1997). Larger hake tend to migrate further,  become  more piscivorous, and consume large amounts of herring off Vancouver Island i n late summer (Ware and M c F a r l a n e 1995). T h e portion of the coastal stock that migrates into C a n a d i a n waters each summer consists of older (age 5+) hake, and the sex ratio is skewed to favor females ( D o r n et al. 1999); however, proportions of fish less than age-5 i n the C a n a d i a n zone have increased since 1996. D u r i n g E l N i n o events, a higher fraction of the stock migrates into C a n a d i a n waters, and ventures further n o r t h into A l a s k a n waters as shown i n Figure 4.2 ( D o r n et al. 1999; Ware and M c F a r l a n e 1995). D u r i n g the warmer years, the d i s t r i b u t i o n of juveniles also spreads north to Oregon and into C a n a d i a n waters. It is thought that this change in juvenile distribution has led to the break down of patterns i n d i c a t i n g among-cohort  81 juvenile cannibalism i n the stock-recruitment d a t a sets (i.e. strong cohorts appearing every 3-4 years, see D o r n et al. 1999). E s t i m a t i n g the fraction of mortality caused by Pacific hake on pink shrimp populations off Vancouver Island is difficult since the fraction of the hake stock i n C a n a d i a n waters must be estimated.  Previous research has documented a strong correlation be-  tween hake abundance i n the C a n a d i a n zone and E l N i n o events; more hake are found i n C a n a d i a n waters during E l N i n o years ( M a r k Saunders, pers comm., Pacific Biological Station, N a n a i m o , B C ) . Therefore, if hake are significant predators on pink shrimp, mortality rates for pink shrimp should be significantly higher i n E l N i n o years.  Pink  shrimp are not listed as one of the major diet items of hake feeding i n the C a n a d i a n zone (Tanasichuk et al. 1991). However, i n Washington and Oregon, predation mortality by Pacific hake on pink shrimp is thought to be significant (Rexstad and P i k i t c h 1986; H a n nah 1995) even though pink shrimp are a rare diet item for Pacific hake i n that region.  A l l observations of pink shrimp found i n the stomaches of Pacific hake report less than 10% of pink shrimp i n the diet (see Table 4.1), w i t h the exception of G o t s h a l l (1969). I n Gotshall's case, he was t r y i n g to develop a relative abundance index for pink shrimp and specifically sampling hake aggregations i n areas where shrimp are known to be present. P i n k shrimp are found more frequently i n the stomachs of larger hake (> 55 c m T L ) than i n smaller hake (Rexstad and P i k i t c h 1986), and larger hake are generally found closer to the b o t t o m ( M a r k Saunders, pers comm, Pacific Biological Station, N a n a i m o , B . C . ) whereas smaller hake generally form large pelagic schools. Rexstad and P i k i t c h (1986) estimate that 9.2 t / d of pink shrimp were consumed by hake off the Washington and Oregon coastlines during the summer of 1983. T h e hake migrations through this area takes roughly 40 days for larger hake and 80 days for smaller hake, and so the total  Figure 4.2: D i s t r i b u t i o n of Pacific hake i n 1992, 1995, and 1998 along survey transect lines. A b u n d a n c e is measured using hydro acoustic backscattering techniques. (Figure supplied by M a r k Saunders, D F O ) .  83 amount of shrimp consumed was estimated at 659.3 t. Table 4.1: Published reports of Pacific hake (Merluccius productus) predation on pink shrimp (Pandalus jordani) and reported % occurrence, or % by weight of pink shrimp in hake stomachs. Location  % Occurrence  California British Columbia Washington and Oregon Washington and Oregon Oregon to Vancouver Island Washington and Oregon  54% 3% 5.70%  4.2 4.2.1  0.30% 0.70%  1.70%  Estimating Predator Pacific C o d  Source  % by Weight  Gotshall (1969a, b) Outram and Haegele (1972) Alton and Nelson (1970) Livingston and Alton (1982) Livingston (1983) Rexstad and Pikitch (1986)  Biomass  Assessment  A time series of biomass for Pacific cod on the west coast of Vancouver Island was reconstructed using a delay-difference model (Deriso 1980; Schnute 1987). Two age groups are represented in this model, age-2 recruits and fish age 3 and older. The delaydifference model is developed by assuming that (1) growth can be adequately described by the Ford-Brody growth function  W i — a + pW , a+  a  (2) fish are fully recruited to the  fishery at age 3, and (3) natural andfishingmortality rates are independent of age for ages 3 and older (Hilborn and Walters 1992). Parameters and variables used in the delaydifference model are defined in Table 4.2. Model parameters were estimated by fitting the model to commercial catch-effort data, and to a cod abundance index derived from the shrimp trawl survey.  84  Table 4.2: Parameter symbols and description of parameters for the delay-difference model. Notation  Description  Estimated Parameters Unfished equilibrium age-2 recruits Ro B Ratio of initial biomass to unfished equilibrium biomass 5 Recruitment compensation parameter Catchability coefficient for commercial fishery Q M Instantaneous natural mortality rate Recruitment anomaly in year t Growth Parameters Value Asymptotic length 89.48 Growth coefficient 0.307 k to Time at length 0 -0.116 a Length-weight scaler 7.38E-06 Length-weight power 3.0963 b K Age of recruitment 2 Observed States E Annualfishingeffort Observed landings in tons c w Mean body weight of landed fish Unobservec States Unfished biomass, biomass in year t B, B Unfished numbers, numbers in year t No, Nt w Predicted mean body weight in year t Rt Age K recruits in year t Predicted annual catch in year t c Ft Fishing rate in year t Survival rate in year t s Derived Parameters Maximum recruits per unit of spawning stock biomass so Maximum recruitment scalar b Mean body weight at equilibrium W Weight-at-age of recruitment S Annual survival rate (e~ ) a Intercept of Ford-Brody growth equation P Slope of Ford-Brody growth equation r  t  t  t  t  0  t  t  t  K  M  85 Growth Parameters for the Ford-Brody growth equation can be calculated from parameters from the von Bertalanffy growth equation as:  oo  (1-/5)  Body weight at the age of recruitment is simply the weight at age K . Growth parameters and length-weight relationship parameters were taken from Westerheim (1996). I n i t i a l States Assuming unfished equilibrium conditions at the start of the model time period, the following derived initial states can be used to define stock recruitment parameters in terms of leading parameters to be estimated from the data. Given an initial estimate of unfished recruitment (R ), natural mortality (M), and recruitment compensation parameter (<5), 0  the unfished equilibrium numbers (N ), mean body weight (w ), biomass (B ) and stock 0  a  a  recruitment parameters (s ,(3) are initialized as follows: 0  N  0  =  Ro (1-5) {Sa + w (l-S)) 1-pS K  B  a  =  Nw 0  0  (4.1) (4.2) (4.3) (4.4) (4.5)  Note that Equation 4.5 describes the Ricker stock-recruitment model.  86  State D y n a m i c s P o p u l a t i o n numbers and biomass over time were calculated using equations 4.8 and 4.9, respectively.  A n n u a l survival (S ) is a function of fixed natural mortality and annual t  fishing mortality, where mortality associated w i t h fishing is conditioned on observed effort (E , E q u a t i o n 4.6). A n n u a l recruitment ( E q . 4.10) was calculated using a Ricker t  stock-recruitment curve, a n d annual recruitment anomalies (<f> ) were treated as model t  parameters to be estimated.  F  = qE  St  =  t  (4.6)  t  (- ~^  (4.7)  M  e  (4.8)  N  = NtS + Rt  B  = S (aN + B ) + w R  t+1  t  t+1  t  Rt  =  t  P  t  r  (4.9)  t  (4.10)  SoBt-KeWt-rtet*  P r e d i c t e d Observations Observed data consisted of annual catches and mean b o d y weights i n the catch. Predicted mean body weights were calculated using equation 4.11, and the Baranov catch equation (Eq. 4.12, Baranov 1918) was used for predicting annual catches.  Ct  = Tp^M)  ^~" e  (Ft+M)  )  t  B  (4.12)  Likelihoods Likelihoods for mean b o d y weights and observed catches are represented by equations 4.13 and 4.14. Here I assume that errors i n reported catches and observed mean weights  87 are lognormaliy distributed, and recruitment anomalies (process errors) also have a lognormal distribution and are penalized to have a mean 0 with variance cr ,. Equation 4.16 2  represents the likelihood of the shrimp trawl survey index, where Z — ln(Y /B ), t  t  t  which  also assumes that errors in relative abundance estimates are lognormaliy distributed (see Walters and Ludwig 1994).  L  w  =  —-—log  L  c  =  —-—log  (4.13)  5>(t)  (4.14)  !  (4.15) - 1 ) log  (Z - Z)'  (4.16)  t  Equations 4.13, 4.14 and 4.16 are self weighting by using the conditional maximum likelihood estimate of the variance (i.e. substituting ^ ^ " ^ likelihood results in ^^Log  (XX  1  _  for a in the negative log 2  ^) )> ignoring a constant scaling term). I further 2  use prior distributions on recruitment compensation (p (5)) and natural mortality rates 0  (po(M)y.  Po(M)  =  1 —7T- 2  ?rr z  Po(S)  =  c  r  M  , In  2  (4.17)  () M  V0.58 J  1 , / —-o In  2a  2  (4.18)  \ 2.74)  where the prior distributions of M and 5 are assumed to be lognormal. The best estimate of natural mortality rate available from the literature is 0.58, and the range of estimates is from 0.25-1.12 (Westerheim 1996). The parameters used for the prior distribution on S in equation 4.18 were taken from Table 1 in Myers et al. (1999), where mean compensation (S) for Gadiformes was used. Note that 5 « exp(loga) from Myers et al. (1999). The coefficient of variation in recruitment anomalies was set at a ^ — 0.6, for natural mortality 2  a  2 M  = 0.38 (Westerheim 1996), and for recruitment compensation erf — 0.2 (Myers et al.  88 1999). T h e objective function to be minimized is the sum of equations 4.13 through 4.18. L i k e l i h o o d profiles were used to explore uncertainty i n leading parameter (M,  estimates  R , and 5). Uncertainties i n nuisance parameters, such as recruitment anomalies, Q  were not explored i n this analysis. I also doubled the value of the variance terms i n the prior distributions to examine their influence on leading parameters.  4.2.2  Pacific H a k e  Assessment  C o a s t a l hake is treated as a single stock and is formally assessed by U . S . and C a n a d i a n agencies every 3 years using fishery independent trawl and acoustic surveys. A b u n d a n c e is assessed using age-structured models, and fishing m o r t a l i t y is partitioned into U . S . and C a n a d i a n components.  T h e relevant issue for my analysis is to determine what  fraction of the total estimated stock is present off southwest Vancouver Island i n the summer months, and how m u c h shrimp are consumed by hake. There are two sources of  information about hake abundance i n the C a n a d i a n zone: d a t a from the triennial  surveys, and volume swept estimates from Ware and M c F a r l a n e (1995). T h e area over which W a r e and M c F a r l a n e estimated hake biomass covers only the southern portion of statistical area 124, and no tow-by-tow information is available to estimate hake biomass in  area 124 and 125 independently.  T h e triennial surveys span the entire coast, from  California to A l a s k a but again there is no detailed information available that would allow for estimation of hake abundance i n areas 124 and 125 specifically. So  i n the absence of direct observation on the abundance of Pacific hake i n areas 124  and 125, I use the biomass estimates provided by Ware and M c F a r l a n e (1995) to assess possible predation impacts of hake on pink shrimp. A b u n d a n c e d a t a for Pacific hake from b o t h the triennial surveys and Ware and M c F a r l a n e (1995) are presented i n Table 4.3.  Note that the methods used to estimate abundance i n these two documents differ  89 Table 4.3: Biomass Estimates for Pacific hake McFarlane, 1995), and for hake i n the C a n a d i a n Year H a k e biomass in L a Peruse bank region (Kt) 1968 183.3 1969 226.42 1970 79 1974 53.25 1975 108.8 1980 1983 438.7 1985 65.4 1986 309.5 1987 256 1988 207.8 1989 138.3 1990 425.5 1991 251.8 1992 1995 1998  for L a Peruse bank region (Ware and Zone, ( D o r n et al, 1999). Hake biomass in Canadian Zone ( K t )  184 242 335  138  554 199 526  and cover different areas, so for years where abundance information is available from b o t h documents, estimates differ. There appears to be a significant relationship between the fraction of hake that enter the C a n a d i a n zone and winter sea level height ( M a r k Saunders pers comm.). Higher winter sea levels i n Bamfield are positively correlated w i t h the fraction of the t o t a l hake stock that enters the C a n a d i a n zone (Figure 4.3).  U s i n g this relationship I was also  able to construct a time series of hake abundances i n the C a n a d i a n zone, using sea-level height as a predictor for the proportion of hake present i n the C a n a d i a n zone from 1976 to 1998. Here I fit a logistic function to the observed time series d a t a shown i n the top panel i n Figure 4.3, where I assumed the asymptote of this function is 1 and estimated two parameters 7 and Xh assuming the deviations are n o r m a l l y distributed. Note that Xh corresponds to the sea level height anomaly at which 50% of the stock was present i n the C a n a d i a n zone, and 7 is a shape parameter that describes how steep the relationship  90 is. Proportion in Canadian Zone vs. Sea Level Height  i -100  i  i  i  i  i  -50  0  50  100  150  Sea  r~ 200  Level Anomaly (mm)  Hake Biomass  i  1  1980  1985  I 1990  1—  1995  Year  Figure 4.3: P r o p o r t i o n of Pacific hake i n the C a n a d i a n zone versus mean JanuaryFebruary sea level height anomalies (top panel). Positive anomalies indicate higher than average sea level heights i n Bamfield, B r i t i s h C o l u m b i a (data provided by M a r k Saunders, DFO). E s t i m a t e d hake biomass (bottom panel) for total hake stock (solid line) and hake present i n the C a n a d i a n zone (dashed line) based on the sea level height relationship.  4.3  Results  4.3.1 The  Pacific C o d  Biomass  reconstructed biomass and fishing m o r t a l i t y rates for Pacific cod are shown i n Figure  4.4 (top panel, residuals and recruitment anomalies shown i n b o t t o m panel). D u r i n g the mid  1950s to the 1970s, estimated cod biomass on the west coast of Vancouver Island  averaged around 13,986 tonnes, and increased to a m a x i m u m of 24,219 tonnes i n 1972. F r o m 1973 to 1986 biomass declined at a rate of 9% per year, reaching a low of 7,941 tonnes. D u r i n g the mid-late 1980's cod populations began to increase reaching another peak i n 1989, and reached an a l l time low of 1,808 tonnes i n 1996. T h e trend i n estimated biomass shown i n Figure 4.4 is consistent w i t h the results obtained by Haist and Fournier  91 (1995). Since 1996, the ground fish t r a w l fishery has had new quota regulations, as well as increased mesh sizes.  Pacific cod i n recent years is predicted to be increasing i n  abundance. Biomass  Fishing Mortality  o  1960  1980 Year  2000  1960  1980  2000  1960  Year  1980  2000  Year  Figure 4.4: E s t i m a t e d biomass time series of age 2+ Pacific cod and fishing mortality rate on the west coast of Vancouver Island from 1956 to 2000 (top row). Circles represent the shrimp t r a w l survey index scaled to cod biomass and dotted lines represent ± 2 standard deviations. B o t t o m row from left to right are recruitment anomalies, catch residuals, and mean weight residuals from the delay difference model.  F r o m the 1950's to 1990, fishing mortality rates (F ) t  dramatically increased to a m a x i m u m of 0 . 6 7 ~  yr  averaged near 0.l6~ , and yr  i n 1994. Since 1996, however, the new  trawl quota management plan has apparently reduced fishing mortality on Pacific cod to  om~ . yr  92 Management of the ground fish trawl fishery underwent a transformation i n 1996, where quota limits for each species were imposed on the fishery. T h i s change reduced targeting on spawning aggregations during the winter fishery, and one of the outcomes has been a decrease i n the mean body size landed. In the past Pacific cod assessments have been fit to catch-at-length data, and when I included this d a t a set i n the assessment, apparent recruitment i n recent years more than doubled. P r i o r to 1996 there was a low probability of capturing small cod because the fishery targeted larger cod concentrated i n spawning aggregations so I chose to omit the catch-at-length d a t a for this reconstruction. A l s o , recent commercial C P U E indices may not be informative about changes i n abundance if fishermen are effectively avoiding cod to prevent exceeding the species quota. G i v e n this behavior, C P U E indices would not reflect recent increases i n abundance that may  have occurred.  I tried to correct for this effect by including abundance indices  generated from the shrimp t r a w l survey, which indicate a small increase i n abundance since 1995. Furthermore, even without incorporating the catch-at-length data into the analysis, the reconstruction was fairly robust to changes i n prior distributions on natural mortality rates, and recruitment compensation (see F i g . 4.5). Catch-at-length data might be informative about changes i n age structure, and if used i n future assessments I would recommend including a second selectivity curve for the period post 1996. In general, the relative abundance d a t a and mean b o d y weight data were informative about key population parameters for W C V I cod. C h a n g i n g the mean of the prior distributions by 10% resulted i n less t h a n 1% changes i n key p o p u l a t i o n parameters. Results were also fairly robust to assumed variances on prior distributions for natural m o r t a l i t y rates and recruitment compensation parameters (Figure 4.5).  93 Recruitment Compensation  2  3  4  5  S Natural Mortality  0.45  1 0.50  0.55  0.60  0.65  0.70  0.75  o.ao  M Unfished Recruitment  i 0  2000  4000  6000  6000  10000  Figure 4.5: L i k e l i h o o d profiles for W C V I Pacific cod p o p u l a t i o n parameters where assumed variances for natural mortality rates a = 0.38 a n d recruitment compensation a — 0.2 (solid line) versus doubling these variance terms (dashed lines). 2  M  2  4.3.2  Correlation with Predator  A b u n d a n c e  P i n k shrimp a n d Pacific cod (adult a n d juvenile) abundances on the W C V I are positively correlated (Table 4.4). It appears that whenever conditions are good for shrimp populations to increase, they are also good for Pacific cod, at least i n statistical area 124.  S h r i m p a n d Pacific c o d abundances were high i n the m i d 1970s a n d reached lows  a decade later. D u r i n g the late 1980s a n d early 1990's b o t h shrimp a n d Pacific cod i n creased i n abundance.  Hake populations for the entire west Pacific h a d also increased  from the m i d 1970's and peaked i n abundance i n 1987 ( D o r n et al. 1999). N o significant correlation was found w i t h either Pacific hake time series presented i n Table 4.3 or w i t h the reconstructed abundance from sea-level height i n Figure 4.5. Correlations between predator abundance a n d natural mortality, a n d predator abundance and recruitment for shrimp i n both statistical areas are also presented i n Table 4.4. Evidence for interactions between predator abundance a n d shrimp m o r t a l i t y would  94 Table 4.4: Correlation coefficients (r) between pink shrimp abundance, mortality and recruitment, versus predator abundance from 1975 to 2000. N corresponds to the number of data pairs, and significant F statistics indicated by * = 0.1, * < 0.05,** < 0.01, and *** < 0.001 Area 124 Area 125 Areas Combined Predator  r  N  F  r  N  F  r  N  F  Abundance  Juv. Pacific cod Pacific cod Pacific Hake I Pacific Hake II Pacific Hake III a  6  C  0.66 25 18.07*** 0.24 25 1.41 0.55 25 10.54*** 0.59 25 12.95** 0.34 25 3.18* 0.56 25 10.89*** 0.04 22 0.03 -0.03 22 0.02 0.01 22 0.00 -0.39 8 1.23 -0.10 8 0.07 -0.31 8 0.73 -0.34 6 0.67 0.26 6 0.36 -0.05 6 0.01  Mortality  Juv. Pacific cod Pacific cod Pacific Hake I Pacific Hake II Pacific Hake III  -0.27 25 -0.04 25 -0.45 22 -0.17 8 0.06 6  1.90 0.05 5.47** 0.21 0.02  -0.05 25 0.07 -0.16 25 -0.04 25 0.04 -0.05 25 -0.16 22 0.55 -0.32 22 -0.17 8 0.21 -0.26 8 -0.41 6 0.99 -0.18 6  0.61 0.06 2.33 0.52 0.17  0.60 25 13.61** -0.25 25 1.59 -0.08 25 0.32 25 2.75 -0.21 25 1.10 -0.12 25 0.42 22 4.53** 0.09 22 0.18 0.35 22 -0.56 8 3.24 -0.22 8 0.37 -0.66 8 0.38 6 0.86 0.46 6 1.34 0.56 6  0.14 0.35 2.87 5.37* 2.23  Recruitment  Juv. Pacific cod Pacific cod Pacific Hake I Pacific Hake II Pacific Hake III  "Hake abundance determined from sea level heights Hake abundance data from Ware and McFarlane (1995) Hake abundance data from Dorn et al. (1999) 6  c  be indicated by positive correlations. No significant positive correlations existed between predator abundance and estimated natural mortality rates for pink shrimp. Similarly, evidence for interactions between predator abundance and shrimp recruitment would be indicated by negative correlations. Shrimp recruitment and juvenile Pacific cod were positively correlated, once again suggesting a shared bottom up affect for these two species (Figure 4.6). Although no significant correlations were found between Hake abundance and recruitment, the negative interaction in area 124 was fairly strong (Table 4.4).  95  0.5  1.0  1.5  Shrimp biomass  2.0  0.5  1.0  1.5  Shrimp biomass  2.0  0.5  1.0  1.5  2.0  Shrimp biomass  Figure 4.6: Relative abundance time series of P i n k shrimp, Pacific C o d , and Pacific hake series I and II from 1975 to 2000. Lower panels are scatter plots of predator abundance versus pink shrimp abundance.  4.4  Discussion  In contrast w i t h the east coast of C a n a d a , A l a s k a and the Barents Sea, shrimp abundance and Pacific cod abundance were positively correlated on the west coast of Vancouver Island.  D u r i n g the m i d 1980's b o t h Pacific cod and shrimp populations had declined  to low abundance, and by 1985, shrimp populations started to recover. W i t h roughly a 2-year lag Pacific cod populations were also recovering, and a similar lag was observed for the decline phase for b o t h species i n the early 1990's. T h e 2-year lag represents the age at which Pacific cod juveniles recruit to the adult population. O n the other hand, Pacific hake and shrimp abundance were not significantly correlated, but the sign of the correlation suggests that as hake abundance increases, shrimp abundance declines. T h e proximate agent for this decline could be predation by hake, causing increased mortality  96 rates for shrimp i n years when hake is abundant;  however there were no significant  positive correlations between hake abundance and shrimp mortality (see Table Tab4.5). D u r i n g w a r m years more hake move into the C a n a d i a n zone, apparently resulting i n increased predation m o r t a l i t y for many other species including herring and euphausiids (Ware and M c F a r l a n e 1995; Robinson and Ware 1999). T h i s is especially evident during the E l N i n o events of 1982-83 and 1997-98. I was unable to find direct evidence (i.e. stomach contents data) to suggest that either Hake or Pacific cod were eating juvenile shrimp or adult shrimp, although shrimp recruitment appears to decline w i t h increasing hake abundance. If Pacific cod were a key mortality agent on pink shrimp, I would have predicted substantial increases i n shrimp abundance when Pacific cod biomass had declined i n 1985 and 1995 (Figure 4.4), and significant positive correlations between shrimp n a t u r a l mortality rates and Pacific cod abundance. Such observations have occurred on the east coast of C a n a d a ( L i l l y et al. 2000, W o r m and Myers (Berenboim et al.  2000).  In Press), and i n the Barents Sea  O n the west coast of Vancouver Island, pink shrimp were  positively correlated w i t h Pacific cod, suggesting that b o t h species are responding to changes resulting from bottom-up forces (see Chapter 5). T h i s is further supported by the positive correlation between Pacific cod abundance and shrimp recruitment. Coastal Pacific Hake populations undergo a seasonal migration, and the range over which this migration occurs varies from year to year. T y p i c a l l y , i n w a r m E l N i n o years, the distribution of hake may extend as far north as A l a s k a (e.g. Figure 4.2, i n 1998), and there is no published information available about how long and how many hake are present i n statistical areas 124 and 125 during the summer. G i v e n the limited data set adjacent to area 124 (Ware and M c F a r l a n e 1995) there was an indication of a negative interaction between shrimp recruitment and hake abundance (Hake II series), but this  97 assumes the density of hake i n area 124 was the same as the adjacent area. T h e few stomach samples available for hake populations indicate that shrimp are a very small component of its diet (Table 4.1). If we were to apply the consumption of shrimp by hake estimated by R e x s t a d and P i k i t c h (1986) of 659 tonnes to the west coast of Vancouver Island population i n 1983, more shrimp would have been consumed than were available. A l o n g the west coast of Vancouver Island, the principal prey species for hake are euphausiids and herring (Robinson and Ware 1994). However, the ratio of mean hake biomass (using d a t a from Ware and M c F a r l a n e (1995)) to shrimp biomass is approximately 93:1, so if hake diets consisted of 1% shrimp this would translate into a mortality rate of 0.33 or approximately 655 tonnes y r  - 1  consumed (assuming hake were  present for 80 days, and its consumption per unit of biomass is equal to 1.6). Hake are thought to be present i n the C a n a d i a n zone for roughly 180 days d u r i n g w a r m years (Robinson and Ware 1994). B y comparison, hake on average consume 208 t of herring daily or 12,700 tonnes during the months of A u g u s t and September where the average diet fraction is 32% by weight (Ware and M c F a r l a n e 1995). R o b i n s o n and W a r e 1994 also noted that increases i n temperatures reduce p r o d u c t i v i t y of euphausiids, therefore i n warmer years hake tend to rely on other forage species, such as herring. G i v e n the available data, I was not able to conclude that migratory Pacific hake or Pacific cod play a significant predation role i n shrimp population dynamics. In the case of Pacific cod, it appears that abundance is i n phase w i t h pink shrimp p o p u l a t i o n along the west coast of Vancouver Island. There is no over-whelming evidence that Pacific cod rely on pink shrimp as a principal forage species, yet occurrence of p i n k shrimp i n the stomachs of Pacific cod is well documented i n this system as well as for gadoids i n other systems. Pacific hake have also been documented as foraging on pink shrimp and significant predation mortality is suspected i n this system as well as off the west  98 coast of Oregon (Hannah 1995). Unfortunately there is insufficient abundance and diet information about hake populations off the west coast of Vancouver Island to quantify the predation impacts of hake on pink shrimp populations; however, the relative differences i n abundance between hake and shrimp indicate that only a small fraction of shrimp i n hake diets would cause significant predation mortality. T h e variability i n hake abundance off the west coast of Vancouver Island appears to be driven by ocean temperature where warmer years leads to higher hake abundance, and shrimp mortality also appears to be related to temperature; it is possible that hake are responsible for at least some of the variability i n pink shrimp n a t u r a l mortality.  Chapter 5 Ecosystem Perspectives "We beiieve the food web modeling approach is hopeless as an aid to formulating management advice; the number of parameters and assumptions required are enormous." .. .Ray Hilborn and Carl Walters, 1992 (pg. 448).  5.1  Introduction  Much offisheriesscience in the past has focused on single species assessment, and ignores or assumes that the interactions that determine distribution and abundance remain constant over time. This view is rapidly changing. Recent observations in failures of stocks to recover (e.g. Northern cod), have raised questions about depensation effects by natural predators (Walters and Kitchell 2001). As we learn more about global climate regimes we are finding more and more correlations between productivity and climate indices (Hare et al. 1999). A central question is how do we sustain natural resources not knowing or understanding how the ecosystem will respond to changes in climate? The answer to this question comes through learning from experience and natural opportunities (Holling 2001). Attempts to understand variation in observed time series data usually consists of building more complex models that include interaction terms between predators or environmental indices (e.g. Walters et al. 1986; Hannah 1993; Caputi et al. 1998; Lilly 99  100 et al. 2000; B e r e n b o i m et al. 2000). Understanding the interactions between environmental variables and production is meaningless i n many management scenarios, because in order to make better predictions we would also have to be able to predict  future  environmental indices (Walters and Collie 1988). Furthermore, many of the environmental correlations break down over time (e.g. D r i n k w a t e r and Myers 1987). Alternatives to single species models driven by correlated information are more complex food web models or multi-species models (e.g. " B o r m i c o n " , Stefansson 1998; " E R S E M " , B a r e t t a et al. 1995, " E c o p a t h w i t h E c o s i m " Walters et al. 2000) where interactions (trophic and fishery) and variation i n p r i m a r y production can be explored. E x p l i c i t l y modeling the dynamic changes i n food webs over time allows for examination of direct interactions (such as predator-prey interactions), as well as indirect interactions (for example how herring fisheries might affect shrimp fisheries, even though there are no direct predator-prey linkages between herring and shrimp) and correlations created through shared responses to changes i n ecosystem p r o d u c t i v i t y and functions. In this Chapter I use a simplified ecosystem model to assess the relative impacts of mortality on west coast of Vancouver Island shrimp populations.  Here I use E c o p a t h  w i t h E c o s i m ( P a u l y et al. 2000) to p a r t i t i o n mortality among natural predators, and commercial fisheries and examine how total m o r t a l i t y changes over time. T h e first step in demonstrating the utility of an ecosystems model is to show that it can replay the history of disturbance i n the ecosystem.  T o do this I w i l l first fit time series d a t a to  E c o s i m predictions. T h i s accomplishes two goals, first it assures that E c o s i m can replay the general trends i n observation data (at least qualitatively), and second it helps i n parameter estimation for the model. Coastal upwelling i n the W C V I ecosystem is a very prominent feature and plays an important role i n p r i m a r y production. Historical annual  101 variation i n p r i m a r y production is estimated by fitting E c o s i m to observed abundance i n dices, total mortality, absolute catches and changes i n mean body weight over time. T h i s reconstructed variation as evidenced i n a n i m a l abundance is compared to other estimates of p r i m a r y production. B o t t o m up control i n E c o s i m is represented by specifying rates of flow from prey to predator i n a term k n o w n as flow control, or relative vulnerability of prey to its predators (see B u n d y 2001; Walters et al. 1997).  5.2 5.2.1  Methods Ecopath with  Ecosim  E c o p a t h is a trophic mass-balance modelling approach, which is based o n the assumption that t o t a l p r o d u c t i o n i n a particular reference year, or period of years, for any group i n an ecosystem can be expressed as:  n  (5.1)  where Bi is the biomass (t k m ) of group i, (P/B)i is the production to biomass ratio of - 2  group i, (Q/B)j is the consumption biomass ratio of group j (predators of group i), and  DCji is the average diet fraction of prey i for group j. EEi is the ecotrophic efficiency, or the fraction of total mortality that is accounted for w i t h i n the modelled system and Y", are fisheries yields. E q u a t i o n 5.1 presumes that the production of group i is equal to the consumption of group i by all predators j , plus non-predation losses (including fisheries landings) of group i over the time period i n question (usually one year). Here I define the P/B ratio as the total instantaneous mortality rate (Z = F + M) d u r i n g the 1950's. T h e defined area for the west coast Vancouver Island E c o p a t h model includes statistical  areas 121 through 126 out to the 200 m depth contour. T h e t o t a l area is approximately  11,000 k m . Biomass estimates and total landings were scaled to tonnes k m 2  - 2  for groups  102 where information from stock assessment reports could be obtained. For groups where no biomass information is available, I allowed E c o p a t h software to estimate  biomass  by specifying the fraction of total mortality being explained by the model (Ecotrophic Efficiency). E c o s i m is initialized using the E c o p a t h inputs, and biomass dynamics for most groups i n the model are then predicted using the differential equation:  SB  n  = f(B ) - M B ht  where % indexes group, other mortality),  M  ifi  ii0  iit  - FB iit  -  iit  d,j,t {B , B ) iit  itt  (5.2)  = (1 — EE/) * (P/B)i represents unaccounted mortality (e.g.,  F is the fishing m o r t a l i t y rate at time t and c i<t  i j ] t  is the consumption of  group i by predator j at time t. For p r i m a r y producers, the biomass production function  f{Bi ) is a simple saturating relationship of the form: it  f (B , ) = B /(1 z t  where  ri  + B hi)  i>t  (5.3)  ht  is the m a x i m u m intrinsic rate of growth, and hi is inversely proportional to  m a x i m u m density. Biomass production for consumer groups j is modelled as the product of t o t a l consumption and a growth efficiency g^.  n  f{B ) = gi^c {B B ) itt  where growth efficiency  = Bi(P/B)i/Qi  iti  itU  (5.4)  itt  is calculated from the balanced E c o p a t h esti-  mates. P r e d i c t i n g the dynamics of consumption, i.e. the "flow" rates Cj , i n E c o s i m is based ti  on the "Foraging A r e n a " concept (Walters and Juanes 1993). T h e basic premise of this theory is that each biomass pool is split between two behavioral states: a n "available state" (V ) representing prey that are vulnerable to consumption by predators, and an i<t  "unavailable state" (B  iit  — V ) representing individuals that are not vulnerable due to itt  103 prey behaviors such as schooling and h i d i n g i n refuge areas. Exchange of prey between available and unavailable states is a function of time that prey spend foraging, which may vary i n response to changing levels of food availability and abundance of competitors and predators.  E c o s i m assumes that prey exchange r a p i d l y between these two states, and  simply keeps track of the (time-varying) e q u i l i b r i u m biomass of available prey.  Thus,  consumption of prey i by predator j at any moment should vary as:  v  =  Vij + B  itt  2v^ + ciijBjt Cij (Bit,Bjt)  (5 5)  — ciijVijfBjf  For this model, a j is the rate of effective search for prey i by predator j and v^j is the it  rate of exchange of prey i between "available" and "unavailable" states.  T h e rate of  effective search a^j is estimated from re-arrangement of E q u a t i o n 5.5 using the E c o p a t h estimates of Qij o, -Bi,o t  a  n  d -Bj,o- Vulnerability exchange parameters v^j are calculated  from a user specified ratio of the m a x i m u m predation m o r t a l i t y to base predation mortality estimated i n E c o p a t h (see description on page 106). Large vulnerability exchange parameters (v^ »  1) represent top-down or L o t k a - V o l t e r r a type control; whereas, small  values represent bottom-up or donor controlled rates of flow from prey pools to predator pools. Note that E q u a t i o n 5.5 predicts a saturating relationship between predation mortality M y = Cij/Bi  and total predator biomass. T h i s form of relationship tends to  eliminate much of the unstable behavior that L o t k a - V o l t e r r a type predator-prey models typically produce. T h i s tends to conserve total biodiversity i n even the most complicated ecosystem models due to limitations on t o t a l predation m o r t a l i t y rates (Walters et al. 1997; Walters et al. 2000; Walters and K i t c h e l l 2001).  104  Age/Size structure in Ecosim Juvenile a n d adult fish typically feed on different prey resources and are vulnerable to different types of predators a n d fisheries. Simulating this trophic ontogeny cannot be accommodated using the simple biomass dynamics model i n E q u a t i o n 5.2. For species that are represented as juvenile a n d adult "pools", E c o s i m uses a Deriso-Schnute delaydifference model to represent adult dynamics a n d an explicit age-structured model that keeps track of monthly juvenile cohorts from b i r t h to the age at recruitment to the adult stock. T h e complete set of derivations and assumptions used i n this formulation can be found i n Walters et a l . (2000). Here I present the basic model structure for split-pool dynamics. Indexing adults as A a n d juveniles as J (group index i is implicit), the basic delaydifference model structure for all split pools is:  B ,t+i  =  e- *[a t{C t)N ,  N  Ait+1  =  N e- ^  Nj, i  =  R(B ,t,N ,t,C )  Nj, ,t i  =  e'^Nj^t  wj, ,t+i  =  wj -i t  A  llt+  a  +  a  + B , )+w Nj^  ZA  A  A  <a  A t  JAt  (5.6)  t  (5.7)  ikit  A  t  PA  + Nj  z  Ait  A  A t  (5.8)  A<t  l<a<k  (5.9)  + gCj /Nj }t  (5.10)  tt  where t is time (months); A^>,t is number of juveniles of age a (months); R(B , At  N ,  C )  A}t  Ait  predicts age-0 juvenile numbers as a function of adult numbers, biomass a n d food consumption; k is age (months) at graduation to adult pool; u>j, ,t is age a juvenile b o d y 0  weight at time t; Cj  is t o t a l food consumption by juveniles i n month t (2~2k k,J,t)'i' c  jt  Njj is total number of juveniles; g is juvenile growth efficiency; a ,t(C t) A  cept of F o r d - B r o d y growth model (w t+i  v s  Al  capita food consumption C  is the inter-  - A,t) which is assumed to depend on per w  (J2k k,A,t/N ); a n d p is the F o r d - B r o d y slope representing c  A}t  Aj  Ajt  105 metabolism. E c o s i m predicts t o t a l mortality rates of adults ZA,t and juveniles  Zj  >t  as the sum of un-  explained mortality M , t o t a l fishing mortality F (described below) and total predation 0  t  losses as:  Z ,t A  = M  + F  Afi  Cj,A,t/B  +  Ait  (5.11)  Ait  j  Zj, = M t  Jfi  + Fj, +  JM/Bj,t  C  t  (5.12)  3 where the s u m m a t i o n is over a l l predators j that eat group i adults and juveniles.  Density-dependent  mortality i n juvenile pools  E c o s i m produces density-dependent juvenile m o r t a l i t y from age-0 to age-fc i n two ways. F i r s t the user may specify a constant weight at recruitment w so that only those juvenile age classes that have body mass larger t h a n w recruit to the adult stock. If prey vulnerability rates to juveniles are set low, then per capita consumption rates decline rapidly w i t h increasing juvenile abundance. T h i s effectively means that under conditions of reduced prey availability, juveniles grow slower and spend more time exposed to high size-dependent predation mortality rates (i.e., longer time required to reach w). Second, simulated juveniles i n the model adjust the amount of time spent feeding i n response to changes i n prey availability and intraspecific competition. E c o s i m represents this effect by using a function that adjusts feeding time to t r y to m a i n t a i n near constant food intake rates per unit time, and feeding time is assumed to be time spent at risk to predation (i.e., adjusts t o t a l time required to consume fixed food ration). T h e model first computes the E c o p a t h base food intake rate copt w h i c h is assumed to optimize a feeding rate versus predation risk function.  A t each time step, E c o s i m computes the realized food intake  rate per juvenile biomass c and then adjusts time spent at risk to predation (i.e., time t  106 spent feeding) as: T = Tt-i(l-a  + ac t/ct-i)  t  (5.13)  op  where a is a user defined [0-1] value that defines the proportional increase i n time spent feeding per unit change i n realized food intake rate c ~\.  T h i s leads to a more complicated  t  model for predicting consumption of prey i by a l l its predators j (including juveniles as predators):  IT) -  C, JJ  { B l  V) \ — "  B j x )  ij i,jBi tBj tTi tTj t  a  }  ~ (v  hJ  / C I A \  v  +v T hJ  ht  i  i  j  + a /T^  (  5  1  4  )  hjBj  Initially, I used the feeding time adjustment method to represent density-dependent mortality for juvenile pools i n the W C V I model. Here I set feeding time adjusment rate a = 0.5, and a low v j for juvenile pools. I d i d not allow for foraging time adjustments it  (a = 0) i n adult pools and set Vij higher for the adults groups. T h e combined effects of low prey vulnerability and feeding time adjustment creates a compensatory relationship between juvenile abundance and survival w h i c h is supported by the substantial number of stock-recruitment d a t a sets now available (Myers et al. 1999). Further adjustments were made to exchange rate parameters (v^) to better explain observed time series data. T h e vulnerability exchange rate parameters take on a variety of positive values depending on E c o p a t h input parameters.  Input parameters specified i n E c o p a t h pro-  vide a base estimate of the predation m o r t a l i t y rate of predator i on prey j (My-o = Bjfi(Q/B)jDCij/Bifl  as represented by the dark circle i n Figure 5.1). T h e user specified  "flow control" parameters (kij) are used to specify a ratio of m a x i m u m predation mortality to the base mortality. T h e "vulnerability" parameters v j are specified as multiples t  of  M : ijfi  =  1 + k ••  t  Kij  Mi  ifi  kij[0,1]  where the user specifies values k^ between 0-1 (the "Flow control" parameter i n Ecosim). T h e Vij's represent the asymptote, or the m a x i m u m predation rate when either predator  107 (Bj) Vij  or prey biomass (Bi) becomes infinitely large.  A s the value of kij approaches 1,  tends to infinity and predation mortality increases linearly w i t h increased predator  biomass (slope of the line equals the encounter rate a^), or mass-action type interactions (Figure 5.1).  i  1  0  1  1  1  2  3  1 4  r 5  Predator Abundance  Figure 5.1: A n example of the relationship between predation mortality rate and predator biomass for different values of k^. T h e dark circle represents the base predation mortality rate m ^ o , calculated using E c o p a t h input parameters. A s the value of k^ decreases the slope of the line through S^n, i j , o decreases and the m a x i m u m predation rate (v j) of predator j on prey i decreases. m  t  5.2.2  Fitting E c o s i m T o T i m e Series  D a t a  Data Types Two  general classes of time series d a t a can be input for an E c o s i m scenario, these  are forcing data and observation data.  Forcing data can include abundance,  fishing  mortality, or relative fishing efforts, and are used i n calculating dynamic changes i n the model (e.g. F i s h i n g effort time series used to calculate fishing mortality rates over time). Observation data come i n two general forms, relative and absolute information, where relative information is assumed to be proportional to model predictions (Y = qX, where q is unknown). There are now 4 different observational time series d a t a types: abundance (biomass), total mortality, absolute catch, and mean body weight (catch and b o d y weight  108 data were added for this thesis, and are now available in the distribution version of EwE). Abundance may be treated as absolute or relative, and relative abundance indices are assumed proportional to the true abundance and are rescaled (using the conditional maximum likelihood estimate of q) to calculate a residual sum of squares. In this thesis, abundance indices and total mortality data are treated as relative indices. Absolute abundance and catch time series information are scaled to the Ecopath units (i.e. t km" ). 2  Changes in mean body weights over time provide information about changes in age composition or growth rates for any given stock. A rapid reduction in average body weight generally represents an increase in proportions of younger-age animals, i.e. in recruitment. Decreases in mean body weight may also occur due to increased mortality on older cohorts. A reduction in average body weight can also be indicative of reduced growth rates, which can occur with heavy competition for food resources, or reductions in time spent foraging due to increased predation risk. Mean body weights are only calculated for split-pool groups, where the needed age-structured accounting of both numbers and biomass is carried out. Mean body weight data collected from commercial catch samples also reflect changes in the proportion of age groups that are vulnerable to fishing. The accounting system in Ecosim assumes all individuals in the adult or juvenile group are fully vulnerable to fisheries, therefore, I choose to treat observed body weight data as a relative index. This may cause bias if changes in selectivity occur over time. For example, in Chapter 4 it was noted that changes in regulations for the groundfish trawl fishery have had dramatic effects on selectivity of juvenile cod. Because of quota regulations, spawning aggregations of Pacific cod are no longer targeted and as a result there has been a marked increase in landings of smaller cod. But if vulnerability is assumed constant over time, then the best interpretation of the data is an increase in  109  juvenile recruitment. Mean body weight data fromfisheryindependent surveys would be preferred over commercialfisherysamples where vulnerability-at-age should remain relatively constant over time and provide a better index of changes in age structure. At the heart of ecosystem modeling is specifying predator-prey interactions that determine the strength of trophic effects. Predator-prey interactions in EwE are initialized in Ecopath by specifying a diet matrix, then the strength of trophic effects are determined by Vij exchange rate parameters (calculated from user specified k^ parameters). Exchange rate parameters are specified by the user, and these (sometimes arbitrary) rates are often a point of debate. There are four alternative methods for estimating exchange rate parameters, including examining diet composition data from two separate time periods with two separate Ecopath models. In appendix C, I propose a method for fitting Ecosim to data on changes in diet compositions over time. Only a very limited set of diet composition data was available for the WCVI ecosystem, and I did not use methods outlined in appendix C to estimate kij parameters. Observation Models, Residuals and Likelihoods  Errors in all four data types are assumed to have a log normal distribution, and residuals are transformed using Z i = log ^ ° ^ ) • Note that the t subscript denotes time and t>  p  the i subscript denotes data set i. It is possible to have multiple data sets of the same type for each group in the model. The residual sum of squares between predicted and observed states for absolute information is then: (5.15)  and for relative abundance data we subtract the mean residual for rescaling: (5.16)  110 The  log-likelihood of a l l observed data given model predictions was then calculated using:  (5.17)  There is no l i m i t to the number of data sets that can be used i n the fitting criterion, and the contribution of each data set to the t o t a l likelihood can be modified by the user (wi). Hence, data that are known to have large errors relative to other d a t a sets should carry less weight.  Note, however, that equation 5.17 is evaluated at the conditional  m a x i m u m likelihood estimate of the variance for each d a t a set and setting a l l u>j = 1 assumes that a l l d a t a are weighted proportional to the inverse of the variance.  E s t i m a t i n g m o d e l parameters Two  types of model parameters can be varied by the search routine, which uses a m o d i -  fied M a r q u a r d t nonlinear search algorithm (see E c o p a t h Software Help M a n u a l ; " F i t t i n g E c o s i m to T i m e Series D a t a " for a description of the M a r q u a r d t A l g o r i t h m ) , to maximize the t o t a l log-likelihood function ( E q u a t i o n 5.17). These parameters are vulnerability exchange rate parameters (vij) and a time series of p r i m a r y production anomalies.  In  analogy w i t h single-species assessment models, vulnerability exchange rate parameters represent mean p r o d u c t i o n responses, and p r i m a r y production anomalies represent process errors.  5.2.3 The  Description of M o d e l  G r o u p s  ecosystem m o d e l was used to examine trophic interactions that contributed to vari-  ability i n pink shrimp abundance and mortality. Therefore, I chose to b u i l d a very simple model (fifteen groups including a detritus pool) that represents the m a i n trophic interactions and fisheries off the west coast of Vancouver Island. T h i s simplification was done to avoid creating a complex model that would include many groups for w h i c h data are  Ill lacking. M o d e l groups included known shrimp predators (Pacific cod, L i n g c o d , and D o g fish), p r i m a r y bycatch species i n shrimp t r a w l fisheries (eulachon and demersal fishes), principal forage species (euphausiids and herring) and the b o t t o m of the food web is represented by phytoplankton, herbivorous zoolplankton and copepods. A l l input model parameters are shown i n Table 5.1, and diet compositions for each group are shown i n Table 5.2.  Following is a brief documentation of each group i n the model and where  E c o p a t h input parameters were extracted from.  Lingcod T h e p r o d u c t i o n parameters (P/B and Q/B  or Z) for lingcod were taken from Cass et a l . (1990),  ratios were modified from Venier and K e l s o n (1996). L i n g c o d were split into j u -  venile and adult components, w i t h the juvenile phase being up to age 2. Biomass for each component was estimated using a delay-difference model fit to commercial catch-effort and mean b o d y weight data (results of stock assessment model are shown i n appendix B).  Dogfish A n estimated PjB  ratio for dogfish was taken from the natural m o r t a l i t y rate estimates  provided i n Walters and Bonfil (1999), assuming fishing mortality i n 1950 was negligible. A s s u m i n g individuals die at a constant instantaneous rate Z , then the productionbiomass ratio is equivalent to t o t a l mortality rate, regardless of metabolic rates ( A l l e n 1971). Q/B  T h e Q/B  ratio was estimated from a meta-analysis relationship that predicts  from aspect ratio (Froese and P a u l y 2002).  Pacific Cod Pacific cod were also separated into juvenile and adult groups, and the abundance data for each group was calculated using stock assessment results i n C h a p t e r 4. P/B  ratios  were taken from total mortality rate estimates i n (Haist and Fournier 1995), which also  112 Table 5.1: Ecopath input parameters, and estimated parameters (in bold) for west coast Vancouver Island model. Group  Lingeod Juv Lingeod Dogfish Adult Pacific Cod Juv. Pacific Cod Demersal Fish Adult Herring Juv. Herring Eulachon Euphausiids Pink Shrimp Copepods Herb. Zoolplankton Phytoplankton Detritus  Trophic level 4.51 4.24 3.89 4 3.68 3.67 3.61 3.54 3.34 2.61 2.04 2.4 2 1 1  Biomass (t/km ) 0.6 0.15 1.134 0.919 0.149 8.5 6.1 3.622 0.211 18.61 0.383 10.629 15.729 35.784 2  —  Prod/Biom (/year) 0.18 0.48 0.04 0.6 1.6 0.52 0.67 1.23 0.47 7.9 0.6 27 40 123  Cons/Biom (/year) 3.3 3.8 1.8 2.4 5.333 2.5 4 11.06 4.3 65 12 60 183.3  —-  —  —  Ecotrophic efficiency 0.331 0.636 0.595 0.243 0.29 0.985 0.794 0.43 0.8 0.6 0.95 0.95 0.95 0.8 0.121  correspond well with mortality estimates in chapter 4. I assumed juvenile P/B ratio was 2.6 time greater than adult total mortality. Pacific cod are fast growing species off the west coast of Vancouver Island, therefore I assumed conversion efficiencies (P/Q) are relatively high (0.25 and 0.3 for adults and juveniles respectively). Demersal fishes  This group is an aggregate group that represents benthicfishes,such has Plueronectids and Bothids (primarily Turbot and Pacific sanddabs), as well as a few shorter lived Sebastes species. Many of the demersalfishesare taken as bycatch in benthic trawl fisheries (shrimp trawls and groundfish trawls). It was assumed that during the 1950's bycatch was relatively high, as trawlfishingtechnology was relatively new at the time, and I assumed afishingmortality rate of roughly 20% for this group. Production and consumption parameters were adopted from Dalsgaard et al. (1998), and also data for turbot in Fargo and Starr (2001). Pacific Herring  Pacific herring is one of the principal forage species on the west coast of Vancouver island, and also supports an important commercialfishery.Prior to the collapse of herring in  113 the late 1960's, intensive harvesting took place prior to 1967 for a reduction fishery. Since the collapse, more valuable roe fisheries have developed and management has been much more conservative. A b u n d a n c e indices for juvenile and adult pools were taken from Schweigert (2000), and production parameters were estimated from t o t a l mortality rates i n a D F O stock assessment model (data provided i n Schweigert 2000). Eulachon There are no targeted commercial fisheries for eulachon on the west coast of Vancouver Island, but large commercial a n d first nations fisheries do exist i n large river estuaries. E u l a c h o n are anadromous, and spend part of their life cycle off the west coast of V a n couver island. Consistent observations of the 2- and 3-year age classes are made during research trawl surveys for shrimp. These d a t a were used to estimate abundance of eulachon i n the W C V I area using methods outlined i n H a y et al. (1997). E u l a c h o n diets were assumed to be similar to that of Pacific herring, however, I assumed a smaller fraction of euphausiids i n the diet, and that 10% of the diet is "imported" from outside the system (from prey species not included i n the model). Euphausiids T w o species of euphausiids (Euphausia  pacifica and Thysanoessa  spinifera)  are important  components of the zooplankton c o m m u n i t y off the west coast of Vancouver Island. B o t h species were aggregated into a single group, and were assumed to have similar production and consumption parameters. P r o d u c t i o n and consumption parameters were taken from Robinson and Ware (1994). N o euphausiid biomass estimates were available for 1950, therefore I allowed E c o p a t h to estimate biomass under the assumption that 60% of the mortality was being explained i n the ecosystem model (EE — 0.6). T h e biomass estimate from E c o p a t h is slightly less t h a n the 22 t k m during the m i d 1980's.  - 2  estimated by R o b i n s o n and Ware (1994)  114  uoi5(ireidioo2 'qraH  O  rH  d  d  oo  spodadoQ  °  d i u u q g >[ii!J  o d  spnsni3i]dng  o d  uoip-B|ng  O a  ex o o H  Atif  rt  a> d  CO  d d  rH  CO  LO rH  rH  LO  LO  d  d d o o o d d  i—i CO LO LO  o o d  CO H< CO  O LO o LO o o o o q q LO d d d d d d  n  LO  LO  o o q q o d o d  rH  d  rH LO O LO CN O rH  POO aHP^d 3 I P V  d LO  o d  a in  .  o  o d d  d  CO  Q  rn  CO  d  q S I J [13S.I9U.18Q  rH  in  O l CO H i—I 00 CO CN CN ^  d  Suujaj-f ^ n p v  0  CO  o o d d  in Sutjjajr  0  O  rH CO O CO O LO  r H f - L O C O O O C O C O t — O r H O J r H ^ f O O r H  q o r H q q q t - ^ q  IISIJSOQ  o o  pooSuiq Anr  pooSuiq  °  , I  LO °  S  o d °  o H d °  rH CN rH rH LO O O O O O O O O O O  CN rH CO  d> d d> czi  ci d ci d  -a o -o O o a cc o  -a o  'o  bO  S -  a,  00  ,.S  ta  J3  PH «  «H  '3 a3 C H  I—) nH  s  O T3  M ^a <  H  H  0  5  H  3 °  co H o d ° d  O l <M CO rH 00 CO O O  CD  O O  ci c  o bO  fa  0 ho C  ' H  a. S 3 rt 03 CO -3 cu  C rt _ O  a  a3  S  "a  o H : a. >> H CU O K -a a, Q  115 Pink Shrimp  Total mortality rates (P/B ratios) for pink shrimp were calculated in Chapter 2 for the 1975 to 2000 period. During the 1950s, there were no directed shrimp fisheries off the WCVI, so I used the natural mortality estimate for pink shrimp during the 1975-2000 period as the 1950 estimate of P/B. However, in Ecosim better statisticalfitsto the time series data were obtained by lowering the P/B estimate from 0.6 to 0.45. Consumption per unit biomass was calculated assuming a 5% conversion efficiency (P/Q — 0.05 is consistent with estimates provided in Dalsgaard et al. 1998). Biomass for shrimp during the 1950's was unknown, and I assume that 50% of the mortality for pink shrimp is being explained by predators listed in this ecosystem model. The estimated biomass in Ecopath is within the ranges of biomass estimates from the stock assessment model between 1975 and 2000. Copepods and Herbivorous Zoolplankton  Production and consumption parameters for copepods and herbivorous zoolplankton were taken from Robinson and Ware (1994), and Dalsgaard et al. (1998). Robinson and Ware reported biomass estimates for copepods and euphausiids, and the ratio between copepods and euphausiids during the mid 1980's was roughly 0.35. To try and achieve this relative difference in abundance for 1950's I adjusted EE's and diets to ensure that euphausiid biomass was at least greater than copepod biomass. Phytoplankton  Diatoms were the principal primary producer in the model. Historically, total production for diatoms in the La Peruse Bank region has ranged between 160-400 g C • m~ • year"" 2  depending on upwelling conditions. I assumed that most of this production is consumed by zoolplankton (EE — 0.95) and that the P/B ratio has been around 123 (Robinson and Ware 1994).  1  116  5.5 5.0  f 4.0  Herring Fishery!  iingcodl  4.5 ShiimpTraw)  _i  I Pacific Cod I  • iPogfishl  Demersal FishJ  nvmmsn  ._ 3.5 O  3.0 2.5  Copepods fPmkShrirtipl  2.0 1.5 Uettitusj  Phytoplankton  1.0 L  Figure 5.2: F o o d web a n d fisheries for west coast Vancouver Island ecosystem model. T h e size of the box is p r o p o r t i o n a l to biomass of each group, and the w i d t h of the connecting line is proportional to the flow. Position on the y-axis is the trophic level of the group.  5.2.4  T i m e Series D a t a for W e s t  Coast Vancouver  Island  T i m e series d a t a used i n the W C V I ecosystem model consisted of absolute catches, total mortality rates, relative abundance indices (either generated from stock assessment models or direct observations i n commercial catch rates) and fishing m o r t a l i t y rates used as time series forcing. Observed catch time series are shown i n F i g u r e 5.3. A l t h o u g h commercial fisheries for lingcod started i n the late 1800's, catch d a t a for southwest Vancouver Island was not separated u n t i l 1954, and I assumed catches were relatively stable between 1950 and 1954 (not shown i n the figure). C o m m e r c i a l fisheries for p i n k s h r i m p o n W C V I started i n 1973, and landings peaked i n 1975. Herring fisheries collapsed i n 1968 a n d the fishery was closed u n t i l 1971.  Recent catches of Pacific cod have been substantially reduced  because of new q u o t a regulations implemented i n 1994 and low stock size.  117 Shrimp Ct  1950  1960  1970  Shrimp Zt  1980  1990  2000  1950  1960  Lingeod Ct  -i 1950  —I 1960 1  I 1970  1960  1970  1980  1990  Herring Zt  I 1980  I 1990  2000  -r——"—i 1950 1960  Herring Ct  1950  1970  1 1970  1 1980  1 1990  2000  1990  2000  Pacific Cod Wt  1980  1990  2000  1950  1960  1970  1980  Year Pacific Cod Ct  oi o  i 1950  iii.,iinniiinninniiii.,ii ni....i t  1 1960  1—1970  1 1980  :  1 1990  r2000  1  Year  Figure 5.3: C o m m e r c i a l landings (t k m ) for west coast Vancouver Island from 1950 to 2000 (left column). T o t a l mortality rate estimates for pink shrimp a n d adult Pacific herring. P i n k shrimp total mortality rate estimates were calculated from combined fishing and natural mortality rate estimates i n chapter 2. T o t a l mortality rates for Pacific herring were inferred from age composition data. B o t t o m panel is mean b o d y weights for Pacific cod, where mean b o d y weight is calculated using: W = B /N . - 2  t  t  t  Relative abundance information for p i n k shrimp consists of the combined biomass estimates from area swept trawl surveys i n areas 124 a n d 125 (see Chapter 2). Pacific cod biomass estimates from Chapter 4, a n d lingeod biomass estimates were calculated using the same delay difference model (see A p p e n d i x B ) . Note that juvenile lingeod a n d Pacific cod are represented by the number of 2-year o l d fish, and adult indices are i n biomass units. A d u l t herring data were obtained from Jake Schweigert ( D F O , pers. comm.) a n d abundance was estimated using a catch-at-age model. T h e euphausiid index (1967-1998) was taken from Robinson a n d Ware (1999), w h i c h was generated using a trophodynamic simulation model driven by oceanographic d a t a (upwelling, sea-surface temperatures a n d  118 Shrimp Bt  Lingeod Bt  Juv Lincod Bt  Eulachon  Juv Herring Bt  Adult Herring Bt  Pacific Cod Bt  Juv Pacific cod Bt  Euphausiids  Year  Figure 5.4: Relative abundance time series used as observational d a t a for the ecosystem model. S h r i m p and eulachon biomass estimates are combined for two statistical areas, herring abundance was constructed using a catch-at-age model, lingeod abundance from stock assessment model fit to commercial C P U E data, and Pacific cod from assessment i n Chapter 4. E u p h a u s i i d abundance (from 1967 to 1998) estimated from a trophodynamic model by R o b i n s o n and Ware (1999). solar radiation). A n eulachon index was generated using the by-catch information from the shrimp t r a w l surveys, and is a combined index for statistical areas 123,124 and 125. A b u n d a n c e indices for eulachon were generated using spatial interpolation (bicubic spline) of area swept research t r a w l information. A l l indices were scaled to have a mean  =1.  119  5.3  Results  5.3.1  F i t t i n g E c o s i m to T i m e Series D a t a  Predicted catch time series from E c o s i m were consistent w i t h observed catch time series (Figure 5.5), despite not using the catch time series for lingcod and Pacific herring i n the fitting procedure. A total of 30 spline point parameters, representing variability i n p r i m a r y productivity, were estimated using the non-linear search routine provided i n E c o s i m . In a l l cases, the relative contribution of catch time series to the total likelihood was deliberately reduced to 10% of the biomass time series. After some initial  fitting  runs, it became apparent that the predicted lingcod biomass time series contradicted model predictions. W h e n the lingcod time series was included i n the t o t a l likelihood, the predicted changes i n herring and Pacific cod were inconsistent w i t h observed data given the estimated sequence of p r i m a r y production anomalies. A s a result of this apparent contradiction the lingcod biomass time series was omitted from the likelihood function (wi = 0.0). In addition to varying p r i m a r y productivity, better fits were obtained by increasing the c o l u m n of fe^ parameters for adult and juvenile herring to 0.85 and 0.75, respectively. Furthermore, the columns of % parameters for lingcod, juvenile lingcod and dogfish were reduced to 0.1, 0.05, and 0.1. Pacific cod (adult and juvenile) fc^'s were increased (0.5 and 0.45) to allow simulated abundance to r a p i d l y increase following the collapse i n the m i d 1980's. Demersal fish k^ were lowered to 0.1, and a l l remaining groups were held at the default value of 0.3. A general pattern is apparent for a l l of the relative abundance d a t a sets shown i n Fig.  5.4; abundance was high i n the 1970's followed by a collapse i n the early 1980's  (a period of low p r o d u c t i v i t y apparently associated w i t h the 1982-83 E l Nino) and a  120 Pink Shrimp Landings  Adult Herring Landings  lljllllllii 1950  1960  1970  1980  1990  2000  1950  Lingeod Landings  1950  1960  1970  1980  1990  1960  1970  1980  1990  2000  Adult Pacific C o d Landings  2000  i  1950  1 1960  1 1970  i  1980  1 1990  2000  Year  Figure 5.5: Predicted (points) and observed (vertical bars) landings for pink shrimp, herring, lingeod and Pacific cod stocks off the west coast Vancouver Island shelf. Fishing mortality rate estimates for all four species were derived from single species assessment models and used as input data into Ecosim. recovery period in the late 1980's early 1990's. Since the mid 1990's, biomass and landings from manyfisherieshave declined, and in some cases to record low abundances (e.g. Pacific cod, see Figs. 5.3 and 5.4). For most groups, trends in biomass time series predictions were consistent with observed data sets (Figure 5.6), with the exception of lingeod during the 1960-80 period. In the 1960's CPUE data indicate that lingeod stocks had increased in abundance, whereas Ecosim predictions indicate a decline in abundance. This "out of phase" pattern for lingeod continues until the mid 1980's, at which time fishing mortality rates increased sharply (see Appendix B) and lingeod abundance continued to decline. Adult lingeod biomass increased during the mid 1970's because of large, persistent recruitment events during this time period. Ecosim predicts favorable conditions for lingeod recruitment during this time period (high primary production, and abundant food resources for juvenile lingeod). Again, this contradicts the estimated numbers of age-2 recruits predicted from a single species assessment model (Appendix B).  121 M u c h of the t o t a l mortality estimated for Pacific cod was a result of fishing, and E c o s i m predictions correspond well w i t h estimates of abundance from the single species model shown i n C h a p t e r 4. Declines i n juvenile Pacific cod recruitment appear to be a result of declining spawner abundance (recruitment over-fishing), because predicted total mortality for juvenile Pacific cod is relatively stable during the decline phase. E c o s i m predictions of herring biomass over time show a similar pattern to the results of statistical catch-at-age models. A l t h o u g h E c o s i m fails to capture the large spike i n biomass during the 1970's following the collapse of the fishery, the predicted increase i n biomass was sufficient enough to fit well w i t h observed catch data (fishing mortalities from single species stock assessment models for Pacific herring were used as input d a t a for Ecosim). F i s h i n g mortalities during the 1960's were excessively high (0.7 down to 0.35) before the fishery was shut down. D u r i n g the closure, predicted predation mortalities (for b o t h juvenile and adult herring) remained relatively constant, preventing any rapid increases i n abundance. T h e relative abundance time series suggest abundance increased as much as 2 times i n comparison to predicted abundance during the m i d 1970's . E c o s i m predictions for eulachon and euphausiid do not fit the observed d a t a well (in the case of eulachon) or other model predictions (in the case of euphausiids). E c o s i m predicts that eulachon predation mortality declined i n the 1980's, due to declines i n predator abundances (lingcod and demersal fishes), however the eulachon stock continued to decline despite release from predation. Estimates of t o t a l mortality for eulachon actually increased due to bycatch i n the shrimp t r a w l fishery, and declining p r i m a r y production starting i n around the late 1970's also contributed to decreased eulachon production. There is no directed fishery for euphausiids off the west coast of Vancouver island, therefore no associated fishing mortality rates. Simulated euphausiid production was closely connected to variation i n p r i m a r y productivity. D u r i n g m i d 1970's simulated p r i m a r y  122 Pink Shrimp  ~t  1  1970  Eulachon  i  r  1  1  1950  1990  Adult Pacific Cod  1  1  1970  Lingeod  i  1990  Juv Herring  Juv Pacific Cod  E CO 03  E o  \  1950  1  1  1  1970  1  r  T  1990  Juv Lingeod  T  1970  1  1  1950  1  1  1  1970  r  1990  1950  1970  Adult Herring  r  i  1990  1950  i  i 1970  i  1990  Euphausiids  i 1990  r  \  1950  1  1  1970  1  1  r  1990  Year  Figure 5.6: P r e d i c t e d biomass time series (line) from E c o s i m and observed relative abundance time series (circles). Note that relative abundance time series are assumed proportional to true abundance, where Y = qB . t  t  production was at a n a l l time high, however E c o s i m predicted euphausiid abundance d i d not increase, due to large increases i n predation m o r t a l i t y by herring. D u r i n g the last decade, there was good agreement between predicted euphausiid biomass and the estimated euphausiid biomass from Robinson and Ware (1999), where b o t h models indicate a declining trend i n p r i m a r y productivity. A n n u a l p r i m a r y production anomalies were calculated i n E c o s i m by estimating 30 points for a spline function that spanned 1950-2000. These 30 parameters were estimated by fitting the E c o s i m model predictions to time series d a t a and the purpose of using a  123 spline function was to filter out high frequency variation that occurs when fitting the model to one, or few data sets that may have large anomalies. E s t i m a t e d mean primary production during the 1950-1977 period was 25% higher than mean p r i m a r y production from 1978-2000 period (Figure 5.7). E s t i m a t e d primary production reaches lows during the 1982-83 and 1997-98 E l N i n o years.  ~n 1950  1  1  1  1  1960  1970  1980  1990  r  -  2000  Year  Figure 5.7: E s t i m a t e d primary production anomalies from 1950 to 2000; where a spline function was fit to 30 equally spaced data points. T h e use of the spline function helps filter out high frequency variation that occurs i n fitting the model for years having only one or few times series w i t h large anomalies.  U s i n g the time series shown i n F i g . 5.7 to drive primary production results i n an overall 39% decline i n the total sum of squares between predicted and observed time series d a t a (Table 5.3). W i t h the exception of 5 d a t a sets, including simulated variation in primary p r o d u c t i v i t y resulted i n better model fits, especially for Pacific herring, Pacific cod and shrimp. For lingcod, adding variation i n primary production over time increased the residual variation between predicted and observed biomass (-189% represents an increase i n residual variation) catches and recruitment. Deviations i n mean b o d y weight for adult Pacific cod increased (4%) when including variation i n p r i m a r y production, and residual variation i n shrimp mortality increased 28%. P r i o r to 1973, there is no information about trends i n abundance for pink shrimp on  124 Table 5.3: Squared residuals for N o M o d e l (or deviations from the mean of the data), a model w i t h no trophic interactions (assume that predation mortality a n d consumption are constant), a model w i t h trophic interactions where p r i m a r y production is held constant over time, and a model that includes estimated primary production anomalies (see Figure 5.7). Positive reductions i n the square residuals (right hand column) between trophic interactions only model and trophic interactions w i t h variable primary productivity indicate a better 'fit' to the observed data. Data Group Shrimp Bt Shrimp Ct Shrimp Z Lingcod Bt Lingcod Ct Juv Lincod Eulachon Juv Herring Bt Adult Herring Bt Herring Ct Herring Z Pacific Cod Bt Pacific Cod Wt Pacific Cod Ct Juv. Pcod Bt Euphausiids Total  No model 10.44 71.04 1.84 0.67 15.36 6.76 24.07 34.07 19.77 69.52 22.89 20.27 1.23 49.11 23.53 7.06 377.64  No trophic interactions 9.97 16.43 0.60 0.82 10.61 6.77 21.04 87.35 62.85 55.29 18.13 12.50 1.22 12.53 18.36 7.06 341.55  the west coast of Vancouver Island shelf.  Trophic interactions with constant PP 9.32 9.01 0.84 1.17 11.94 12.23 22.73 39.12 26.47 30.49 17.89 8.97 1.27 8.97 14.09 6.35 220.85  'IVophic interactions and variable PP 6.57 7.64 1.07 3.38 12.99 16.12 18.28 11.06 8.39 13.89 17.48 2.33 1.32 2.39 6.12 6.16 135.19  Reduction in SS 29% 15% -28% -189% -9% -32% 20% 72% 68% 54% 2% 74% -4% 73% 57% 3% 39%  A s a result predicted biomass for shrimp is  relatively stable during this time period. Just prior to the start of the fishery, E c o s i m predicts an increase i n abundance associated w i t h increases i n p r i m a r y productivity. F r o m 1975 to 2000, simulated trends i n shrimp biomass are consistent w i t h combined abundance estimates from Chapter 2. I n comparison w i t h E c o s i m predictions, the single species model indicates more variability i n abundance, and the recent decline started later (1995) i n comparison to the decline predicted by E c o s i m (1989, see F i g . 5.6). Variations i n upwelling, sea surface temperatures, and solar radiation strongly influence production dynamics off the west coast of Vancouver Island (Robinson and Ware 1994), especially zoolplankton communities (Mackas et a l . 2001). E c o s i m was able to capture some of this variation by estimating relative trends i n p r i m a r y production using data only from higher trophic levels. N o oceanographic data or climate indices were used  125 to drive primary production i n the model. Instead, annual relative primary production was treated as an arbitrary time variable to be estimated from time series data on relative abundance, total mortality, mean body weights and absolute catches for higher trophic level species.  In effect, E c o s i m used d a t a from upper trophic levels to infer or back  calculate historical trends i n apparent p r i m a r y production. Recent downward trends i n calculated relative primary production correspond well w i t h estimated diatom production from Robinson and W a r e (1999) during the 1990's. D u r i n g the m i d 1970's Robinson and Ware estimate an increasing trend i n diatom production, whereas E c o s i m predicts a declining trend i n primary production (see F i g . 5.8). T h i s discrepancy occurs because abundance estimates for shrimp, herring, and Pacific cod were a l l declining during this time period, therefore the d a t a support a declining relative production scenario.  Primary Production Anomalies Diatom Production  1.5  o-  1.0  H  0.5  1970  1975  1980  1985  1990  1995  Year  Figure 5.8: Comparison between estimated annual mean p r i m a r y production anomalies 2 1 from Ecosim, and calculations of D i a t o m production (g C m yr~ ) obtained by Robinson and Ware (1999) using a model based on upwelling (nitrogen), water temperature, and solar radiation.  T h e relative primary production series estimated i n E c o s i m also resembles the negative pacific decadal oscillation (Figure 5.9). T h e P D O and estimated primary production are negatively correlated, thus during the negative phase of the P D O apparent primary  126 production off the west coast of Vancouver Island is higher t h a n average.  oo C  ir o  CM  Q CL  1950  1960  1970  1980  1990  2000  Year  Figure 5.9: Comparison between estimated annual mean p r i m a r y production and the negative Pacific Decadal Oscillation ( P D O ) .  5.3.2  Changes in S h r i m p  Mortality  D u r i n g the 1990's, total predation mortality on pink shrimp has declined ( F i g . 5.10) corresponding w i t h apparent declines i n predators. T h e p r i n c i p a l predators i n the ecosystem model, juvenile and adult Pacific cod, have been declining since the m i d 1970's. D u r i n g the period 1950-70, mean predation mortality rates induced by Pacific cod predation were roughly 0 . 5 1 ~ , or 78% of the t o t a l shrimp m o r t a l i t y rate. D u r i n g the 1990's the yr  calculated predation rate dropped to 0 . 2 2  _3/r  , or 30% of the t o t a l mortality rate. Other  pink shrimp predators (lingeod and dogfish) rarely consume pink shrimp, and less t h a n 20%  of the total mortality was attributed by E c o s i m to these predators.  Lingeod and  dogfish abundance have not declined nearly as much as Pacific cod abundance, and recently (1999-2000) dogfish predation has increased to 30% of the total shrimp mortality.  I assumed during the 1950's that 95% of the total mortality (Z) on pink shrimp populations was being explained by the predators included i n the model (this is done  127  1  1.4  H  1.2  H  — —  1.0  H  — •  0.8  J  0.4  -A  0.2  H  Total Mortality Dogfish Adult Pacific Cod Juv. Pacific Cod Lingcod Juv. Lingcod Fishing Mortality Z from SSM  3 c  o.o H 1950  1960  1970  1980  1990  2000  Year  Figure 5.10: A l l mortality components from E c o s i m compared to Z from single species model ( S S M ) . T h e area between each line is the t o t a l mortality due to each predator or fishing, and the area between total mortality and fishing mortality is the unexplained mortality (1 EE). by setting the EEi term equal to 0.95). There are a d d i t i o n a l predators that were not included i n the E c o p a t h model, and such additional predation is assumed to be constant over time i n E c o s i m . In contrast, the single species assessment model i n Chapter 2 treats natural mortality ( M ) as an arbitrary temporal variable to be estimated from the data. Overlaid on Figure 5.10 are the relative estimates of Z (F + M) from the single species model. There are two points to note about this comparison: 1) total mortality for b o t h E c o s i m and single species models increased 7% and 21%, respectively between 1975-85 and 1990-00, and 2) there are two negative residuals that correspond w i t h 1982-83 and 1997-98 E l N i n o events (see F i g 5.11). E c o s i m predicted an increase i n total mortality rates during the late 1980's and 1990's. T h i s increase i n t o t a l m o r t a l i t y rate was a result of increased fishing mortality. D u r i n g  128  1975  1980  i  1  1  1985  1990  1995  1— 2000  Year  Figure 5.11: Standardized differences between Z estimates from Ecosim versus Z estimates from single species model in Chapter 2. Positive values indicate years where Ecosim predicted higher mortality rates than estimated using the SSM. the same time period, predation mortality by Pacific cod was declining. In contrast, natural mortality rates estimated from the single species model were above average during the same time period (see Fig 2.17 on page 50). For example, in 1996, estimated natural mortality rates peaked in the single species model, whereas Ecosim predicted small predation mortality rates primarily associated with the large decline in Pacific cod. The apparent contradiction might be better explained if additional predators, such as migratory hake, were included in the model. The residuals in the Z estimates between the single species model (estimated from size composition data) and the ecosystem model (estimated from changes in trophic interactions) contain possible information about the impacts of predators that were not included in the ecosystem model. That is, some of the additional variation in Z might be explained if additional predators, such as hake, were included in Ecosim. Migratory Pacific hake were omitted from the ecosystem model because abundance information for the model area was incomplete. However, for years in which hake abundance data were available, I compared the residuals in the Z estimates from the two models with hake  129  •  Z differences  o c  o >  0  OC  -2  1975  1980  1985  1990  1995  2000  Year  Figure 5.12: Time series of standardized indices of hake abundance and residuals between ZEWE and ZSSM- All indices have a mean 0 and cr = 1. Hake series I was derived from sea level height indices and hake series II and III were calculated from area swept trawl surveys (see Chapter 4 for details). Table 5.4: Correlation coefficients (r) between hake abundance time series and Z differences from single species model and ecosystem model (Data shown in Figure 5.12). Pacific hake I Pacific hake II Pacific hake III Z difference 0.226 0.036 -0.229 abundance (Figure 5.12). During the 1982-83 El Nino event, Ecosim fails to capture a large increase in apparent shrimp mortality resulting in a large negative difference between the two estimates of Z. During this time period, two abundance indices suggest that hake were at an all time high. A negative deviation also occured during the most recent 1997-8 El Nino event as well. A negative correlation between the differences in Z (Fig. 5.12) and hake abundance would indicate that by including hake in the ecosystem model, some additional shrimp mortality could be explained by hake predation. In other words, are the differences in Z shown infigure5.12 a result of omitting hake from the ecosystem model? There were no strong correlations between the hake abundance indices from Chapter 4 and the differences in Z between the two models (Table 5.4). Including the hake III abundance in the ecosystem model might explain an additional 11% (r value from Table 5.4) of the difference in shrimp mortality between the two models.  2  130 5.4  Discussion  The goal of this chapter has been to examine shrimp population dynamics off the west coast of Vancouver Island (WCVI) in an ecosystem context. I used an ecosystem model to look at how changes in predation mortality and variation in primary production affect the population dynamics of pink shrimp. The ecosystem model suggests that reductions in predator abundance by commercial fishing has reduced predation mortality for pink shrimp. This reduction in predation mortality has apparently not resulted in an substantial increase in shrimp abundance, since the total mortality rate has slightly increased due to increased fishing activity for pink shrimp. Also, variation in primary production has apparently contributed substantially to variation in production of shrimp, as well as other species in the ecosystem model. Natural mortality rates for pink shrimp populations off the west coast are highly variable, and this is also true for other nearby systems (Hannah 1995), as well as for different species in northern waters (Fu and Quinn II 2000; Lilly et al. 2000). Calculated changes in predation mortality by Pacific cod can explain some of the estimated variability in shrimp mortality. Declining Pacific cod abundance during the 1980's and 90's has resulted in reduced overall predation mortality. In contrast to thefindingsof Worm and Myers (In Press) shrimp abundance on WCVI has been positively correlated with cod abundance, indicating that although shrimp mortality rates decline with declining cod abundance, there is no net increase in shrimp production. Both cod and shrimp have responded positively to increases in primary production at the bottom of the food chain, and the Pacific cod stocks have been insufficient, at least during the last 3 decades, to depress shrimp populations. Predation impact by lingcod and dogfish on shrimp populations has been relatively constant over the last 50 years. Shrimp are a minor component in the diet for both  131 of these predators, and combined affects of these predators account for roughly f 5% of the total mortality during the 1950-77 period and 20% during the 1978-2000 period. In contrast, during the same time periods, Pacific cod predation can account for roughly 78% and 30% of the predation mortality (see F i g . 5.10).  B o t h dogfish and lingeod  populations have been declining during the last 2 decades as a result of by-catch (in t r a w l fisheries) and directed fisheries, respectively.  Despite recent declines i n dogfish  and lingeod abundance predation rates on shrimp by these two predators has increased slightly due to declining Pacific cod and euphausiid abundance. T h e EE term i n E c o p a t h is the assumed proportion of t o t a l mortality that is explained by trophic interactions and fishing mortality.  T h e unexplained proportion (1 — EE)  i m p l i c i t l y represents other predation losses, disease, or senescence. If predation mortality is additive (i.e.  modelled predators are not consuming shrimp that are going to die  anyway), then unexplained loss rate (M  a  = 1 — EE) is assumed to be constant over time.  In the ecosystem model I assumed that losses due to predation by hake are represented i n the M  a  term, so these losses were treated as a constant times shrimp biomass over  time. Pacific hake were o m i t t e d from the Ecosystem model because time series data for the model area were incomplete and impossible to p a r t i t i o n from the coast wide hake stock. T h e migratory hake stock is large relative to prey groups specified i n the model. A s discussed i n Chapter 4 the spatial distribution of hake is highly variable, even w i t h i n the C a n a d i a n zone during summer months (Mackas et al.  1997), and it is extremely  difficult to estimate the mean biomass i n the model area. T h e influence of Pacific hake predation on the west coast of Vancouver Island ecosyst e m is significant for many forage species (Tanasichuk et al.  1991,Ware and M c F a r -  lane 1995,Robinson and Ware 1999). G i v e n the relative differences i n hake and shrimp biomass, and the potential for high predation mortality rates, hake should be included  132 i n the model. P r i m a r y and secondary prey species for hake are euphausiids and herring (Ware and M c F a r l a n e 1995), which are also p r i m a r y forage species for a l l large fish i n this ecosystem model. A s sea surface temperatures increase, euphausiid production declines (Robinson and Ware 1994), which may result i n increased predation on herring and perhaps shrimp by hake. A s noted i n Chapter 4, very small increases of shrimp i n the daily ration of Pacific hake could translate into very large changes i n predation mortality rate. A s suggested by the weak correlation between the residuals i n Z estimates between the two models and an estimate of hake abundance i n the C a n a d i a n zone, including hake i n the model may explain at least an additional 11% of the variability i n Z for shrimp. T h e results i n Table 5.3 are comparisons of 3 different models i n which (1) trophic i n teractions are assumed constant, (2) trophic interactions were included, and (3) b o t h trophic interactions and variability i n primary p r o d u c t i v i t y are included. A s s u m i n g trophic interactions are constant over time is equivalent to using a single species model w i t h constant natural m o r t a l i t y rates. B y comparing a m o d e l w i t h fixed natural mortality rate to a model w i t h variable natural mortality rates, we should be able to make inferences about how much of the dynamics are explained by trophic interactions. Less of the variation was explained for top predators (lingcod and Pacific cod) by including trophic interactions. In other words, much of the dynamics for lingcod and Pacific cod are determined by fishing mortality. In contrast, i n c l u d i n g trophic interactions explained a significant amount of variation i n Pacific herring abundance and mortality, and variations in natural mortality are likely, as was determined by Schweigert (2000). Including trophic interactions for pink shrimp d i d not result i n better explanations of the time series Z data; however, including b o t h trophic interactions and variation i n p r i m a r y productivity led to more of the variance being explained than could the fixed mortality model. In summary, the west coast Vancouver Island ecosystem model explains trends i n  133 catch and biomass time series for pink shrimp. These trends were consistent w i t h predictions generated i n the single species model. T h e recent decline i n Pacific cod explains some of the changes i n total mortality rates for shrimp; however, E c o s i m was unable to capture the two large natural mortality rate anomalies that occurred during the 198283 and 1997-98 E l N i n o events. Judging by the differences i n the estimates of Z from E c o s i m and the single species models, and evidence of more hake i n the C a n a d i a n zone during w a r m years, it appears that including Pacific hake i n the model would explain these large mortality anomalies.  T h e positive correlations i n shrimp and Pacific cod  abundance are best explained by variation i n p r i m a r y production. E s t i m a t e d variation i n primary production corresponds well w i t h other independent estimates of primary production, especially i n recent years.  Chapter 6 General Discussion and Summary This thesis has examined the dynamics of pink shrimp populations off the west coast of Vancouver Island from several perspectives. First, historical abundance, recruitment and natural mortality rates for pink shrimp were reconstructed from a single species perspective. From a stock-recruitment perspective, I examined if environmental correlates, such as upwelling or sea surface temperatures, explain additional variation in recruitment. I also examine if variation in estimated natural mortality are correlated with predator abundance. Finally, using a dynamic ecosystem model I combined all major interactions, including additionalfisheriesand variation in primary productivity, to explain variation in shrimp production and mortality. In each of the previous chapters, I have tried to develop more complex models in order to explain additional variation in observed data sets. The ultimate goal of explaining additional variation is to better manage the resource by constraining abundance predictions with known predator-prey or environmental relationships. Such models should recognize strong trophic interactions that directly affect natural mortality rates, examine tradeoffs between differentfisheries,and be able to generate alternative hypotheses about how variation in primary productivity cascades up through the ecosystem. In addition to building more complex models that explicitly represent environmental 134  135 variability and predator-prey interactions, a key aspect of sustainable fisheries management is estimating o p t i m a l harvest controls. In general, there are two key population parameters required for estimating o p t i m u m harvest objectives: a measure of how large the stock is (such as the unfished biomass, or carrying capacity), and rates of compensation, (such as the m a x i m u m intrinsic rate of growth, steepness of the stock-recruitment relationship, or parameter(s) that describe density dependent rate changes such as the E c o s i m Vij parameters).  Note that if the harvest control rule is a fixed exploitation  rate, then it is not always necessary to have an estimate of the carrying capacity, just a good estimate of present biomass or methods to directly monitor exploitation rates (Walters and M a r t e l l 2002). E v e n more complex age-structured models w i t h hundreds of nuisance parameters require these 2 key parameters for determining o p t i m a l harvest objectives. A key issue is to ensure that the assessment model and the type of data being used to challenge the model are informative, and unbiased, about rates of compensation and p o p u l a t i o n scale. If harvest rate objectives are met and p o p u l a t i o n parameters are well determined, then using fixed rules for managing the fishery are quite robust, over the long term, to environmental variability (Walters and P a r m a 1996). However, fixed harvest rules may also lead to pathological over-fishing because of delayed responses (i.e. Walters and K i t c h e l l 2001) or non-stationary effects, and here is where the u t i l i t y of ecosystem models can help identify these complex tradeoffs. In chapter 2, a stock assessment model was used to estimate average recruitment i n each statistical area assuming no underlying stock-recruitment relationship. T h i s was accomplished by directly estimating recruitment each year (an average recruitment level (R ) a  and a series of recruitment anomalies (w ) t  penalized to have a mean of 0).  In  other words, the model parametrization was meant to be informative about how large the stock has been but is uninformative about density-dependent changes i n juvenile  136  survival rates. To examine density-dependent changes in juvenile survival rates, I then compared recruits per spawner versus spawners. From a management perspective it is possible to proceed with making predictions about future recruitment using a derived stock-recruitment relationship from this assessment model. However, as noted in Chapter 2, past recruitment off the west coast of Vancouver Island have been highly variable, especially during periods of low spawner abundance. If the mechanism(s) that generates this variability can be identified, then additional information can be used to provide more informative predictions about future recruitment events (Caputi et al. 1998). For example, if it were known that spring sea-surface temperatures were negatively correlated with juvenile survival rates, then results from a multiple regression model could be used to provide additional information about juvenile survival rates given the recent sea surface temperatures. Using a nonlinear multiple regression framework in Chapter 3, I examined a few oceanographic indices as environmental correlates for shrimp recruitment. Working under the hypothesis that shrimp in each statistical area are separate stocks, I was unable tofinda single correlate that explained a significant amount of variation for both statistical areas. Although sea surface temperature, upwelling and PDO indices explain some variability in recruitment in each of the areas, the fraction of the total variance explained is small (< 20%). Even if a significant correlate was found to predict deviations in recruitment, it would be dangerous to rely on such relationship for predictive purposes because typical relationships between production and environmental indices break down over time (Drinkwater and Myers 1987). Furthermore, this method of improving predictions also assumes we can predict future environmental factors (Walters and Collie 1988).  137 V a r i a b i l i t y i n natural mortality rates over time also increases uncertainty about future abundance.  Using age-length composition information, I was able to estimate a  time series of natural m o r t a l i t y rates using the assessment model i n chapter 2. B y tracking cohorts over time, length frequency d a t a provide compositional information about population size structure. A s s u m i n g growth and selectivity have remained constant over time, t o t a l mortality rates (Z) can be inferred by tracking changes i n the length composition data. A l s o , the a d d i t i o n of independent estimates of relative fishing mortality rates from commercial log-book records allows for separation of Z into its components M and F. Estimates of changes i n natural mortality rates for pink shrimp off the west coast of Vancouver Island have been highly variable. N a t u r a l mortality rates for shrimp were estimated to increase during E l N i n o years and the mechanism for increases i n natural mortality for pink shrimp has never been explored directly, nor is there any direct information (i.e. stomach contents data) for elevated predation by predators. Indirect measures (correlations between predator biomass and prey abundance), as well as the small amount of information about hake distribution relative to sea surfaces height anomalies, does suggest that Pacific hake populations might play a significant role i n pink shrimp dynamics. These findings are also consistent w i t h findings about pink shrimp-hake interactions i n Oregon and W a s h i n g t o n (Hannah 1995). P r e d i c t i n g changes i n natural mortality rates for pink shrimp requires identification of the mechanisms responsible (e.g. hake predation) and establishing a functional relationship for this interaction. In chapter 5, I examine the shelf populations of pink shrimp from an ecosystem perspective. Using time series d a t a from several different fisheries and assessment programs, I was able to fit a d y n a m i c ecosystem model that included estimating a time series of relative p r i m a r y production anomalies. T h e goal of this chapter was to determine how much  138 of the variation i n production arid mortality could be explained by including trophic i n teractions and variability i n p r i m a r y production i n the assessment process. In summary, including trophic interactions while assuming constant production degraded statistical fits to the data; however, including variable p r i m a r y production explained an additional 29% and 15% of the variability i n relative shrimp abundance and observed shrimp catch time series, respectively (see Table 5.3 on page 124). Nearly all of the time series d a t a included i n the ecosystem assessment were consistent w i t h the hypothesis that p r i m a r y productivity has declined over the last 3 decades.  Coincidentally, long term trends i n  primary production (estimated i n E c o s i m by fitting model predictions to time series d a t a from fish populations) shifts downward around 1977, which corresponds well w i t h regimes shifts documented by B e a m i s h et al. (1999).  Bibliography A l l e n , K . R . (1971). R e l a t i o n between production and biomass. J. Fish. Canada 28, 1573-1581.  Res.  Bd.  Anderson, P. J . and J . F . P i a t t (1999). C o m m u n i t y reorganization i n the G u l f of A l a s k a following ocean climate regime shift. Mar. Ecol. Prog. Ser. 189, 117-123. Baranov, F . (1918). O n the question of of the biological basis of fisheries. Issled. Ikhtiologicheskii Inst. Izv. 1, 81-128.  Nauchn.  Baretta, J . , W . E b e n h o h , and P . R u a r d i j (1995). T h e E u r o p e a n regional seas ecosystem model, a complex marine ecosystem model. Netherlands Journal of Sea Research 33, 233-246. Beamish, R . J . , D . J . Noakes, G . A . M c F a r l a n e , L . K l y a s h t o r i n , V . V . Ivanov, and V . K u r a s h o v (1999). T h e regime concept and n a t u r a l trends i n the production of salmon. Can. J. Fish. Aquat. Sci. 56, 516-526. Berenboim, B . , A . Dolgov, V . K o r z h e v , and N . Y a r a g i n a (2000). T h e impact of cod on the dynamics of Barents Sea shrimp (Pandalus borealis) as determined by m u l t i species models. J. Northw. Atl. Fish. Sci. 27, 69-75. Boutillier, J . , I. Perry, B . W a d d e l l , and J . B o n d (1998). Assessment of the offshore Pandalus jordani fishery off the West Coast of Vancouver Island. PS ARC Working Paper no. 197-11 in M. Stocker [Ed.] Pacific Stock Assessment Review Committee (PSARC) Annual Report for 1998. Can. MS. Rep. Fish. Aquat. Sci., 231 p. Bundy, A . (2001). F i s h i n g on ecosystems: the interplay of fishing and predation i n Newfoundland-Labrador. Can. J. Fish. Aquat. Sci. 58, 1153-1167. Butler, T . (1980). Shirmps Fish. Aquat. Sci.  of the Pacific  Coast of Canada, Volume 202. C a n . B u l l .  C a p u t i , N . , J . W . Penn, L . M . Joll, and C . F . C h u b b (1998). Stock-recruitmentenvironment relationships for invertebrate species of Western A u s t r a l i a . Can. J. Fish. Aquat. Sci. 125, 247-255. Cass, A . J . , R . J . B e a m i s h , and G . A . M c F a r l a n e (1990). L i n g c o d (Opiodon elongatus). Can. Spec. Pub. Fish. Aquat. Sci., 40 p. Chen, Y . and D . Fournier (1999). Impacts of atyical data on Bayesian inference and robust Bayesian approach i n fisheries. Can. J. Fish. Aquat. Sci. 56, 1525-1533. Dalsgaard, J . , S. S. Wallace, S. Salas, and reconstructions of the Strait of Georgia: years ago. In D . Pauly, T . Pitcher, D . the future: reconstructing the Strait of Fisheries Centre Research Reports. 139  D . Preikshot (1998). Mass-balance mdoel the present, one hundred, and five hundred Preikshot, and J . Hearne (Eds.), Back to Georgia ecosystem., Volume 6, pp. 72-91.  140 Dennis, J . E . and R . B . Schnabel (1983). Numerical Methods for Unconstrained mization and Nonlinear Equations. P r e n t i c e - H a l l , Englewood Cliffs, N J .  Opti-  Deriso, R . B . (f980). Harvesting strategies and parameter estimation for an agestructured model. Can. J. Fish. Aquat. Sci. 37, 269-282. D o r n , M . W . , M . W . Saunders, C . D . W i l s o n , M . A . Guttormsen, K . Cooke, R . Kieser, and M . E . W i l k i n s (1999). Status of the coastal Pacific h a k e / w h i t i n g stock i n U . S . and C a n a d a i n 1998. Canadian Stock Assessment Secretariat Research Document 99, 101 p. Drinkwater, K . F . and R . A . Myers (1987). Testing predictions of marine fish and shellfish landings from environmental variables. Can. J. Fish. Aquat. Sci. 44, 15681573. Fargo, J . and P . J . Starr (2001). T u r b o t stock assessment for 2001 and recommendations for management i n 2002. PSARC Working Paper G2001-10. Foreman, M . G . G . and R . E . T h o m s o n (f997). Three-Dimensional M o d e l Simulations of Tides and B o u y a n c y Currents A l o n g the West Coast of Vancouver Island. Journal of Physical Oceanography 27, 1300-1325. Fournier, D . A . , J . R . Sibert, J . M a j k o w s k i , and J . H a m p t o n (1990). M u l t i f a n a likelihood-based method for estimating growth parameters and age composition from multiple length frequency d a t a sets illustrated using data for Southern bluefin t u n a (Thunnus maccoyii). Can. J. Fish. Aquat. Sci. 47, 301-317. Fournier, D . , J . Sibert, and M . Terceiro (1991). Analysis of L e n g t h Frequency Samples w i t h Relative Abundance D a t a for the G u l f of M a i n e N o r t h e r n S h r i m p (Pandalus borealis) by M U L T I F A N M e t h o d . Can. J. Fish. Aquat. Sci. 48, 591-598. Francis, R . C . and S. R . Hare (1994). Decadal-scale regime shifts i n the large marine ecosystems of the North-east Pacific: a case for historical science. Fish. Oceanogr. 3, 279-291. Froese, R . and D . P a u l y (2002). F i s h B a s e . World www.fishbase.org, M a y 22,2002.  Wide Web electronic  publication.  F u , C . and T . J . Q u i n n II (2000). E s t i m a b i l i t y of natural mortality and other population parameters i n a length-based model: Pandalus borealis i n K a c h e m a k Bay, Alaska. Can. J. Fish. Aquat. Sci 57, 2420-2432. G e l m a n , A . , J . C a r l i n , H . Stern, and D . R u b i n (1995). Bayesian man & Hall.  Data Analysis.  Chap-  G i l l i s , D . M . , R . M . Peterman, and A . V . T y l e r (1993). Movement dynamics i n a fishery: A p p l i c a t i o n of the ideal free d i s t r i b u t i o n to spatial allocation of effort. Can. J. Fish. Aquat. Sci. 50, 323-333. Gotshall, D . W . (1969). T h e use of predator food habits i n estimating relative abundance of the ocean shrimp, Pandalus jordani, R a t h b u n . FAO Fish. Rep. 57, 6 6 7 685. Haist, V . and D . Fournier (1995). Pacific cod stock assessment for f995 and recommended yield options for f 996. PSARC Working Paper G95-3. H a n n a h , R . W . , S. A . Jones, and M . R . L o n g (1995). Fecundity of the ocean shrimp (Pandalus jordani). Can. J. Fish. Aquat. Sci. 52, 2098-2107.  141 Hannah, R. W . (1993). Influence of environmental variation and spawning stock levels on recruitment of ocean shrimp (Pandalus jordani). Can. J. Fish. Aquat. Sci. 50, 612-622. Hannah, R. W . (1995). "Variation in geographic stock area, catchability, and natural mortality of ocean shrimp (Pandalus jordani): some new evidence for a trophic interaction with Pacific hake (Merluccius productus).". Can. J. Fish. Aquat. Sci. 52, 1018-1029. Hare, S. R., N . J. Mantua, and R. C. Francis (1999). Inverse production regimes: Alaskan and West Coast Salmon. Fisheries 24, 6-14. Harris, M . (1999). Lament for an Ocean. The Collapse of the Atlantic True Crime Story. McClelland & Stewart.  Cod Fishery:  A  Hay, D., J. Boutillier, M . Joyce, and G . Langford (1997). The eulachon (Thaleichthys pacificus) as an indicator species in the North Pacific. In F. Funk, T . J. Quinn, J. Heifetz, J. N . Ianelli, J. E . Powers, J. F. Schweigert, P. J. Sullivan, and C. I. Zhang (Eds.), Forage Fishes in Marine Ecosystems., pp. 509-530. Alaska Sea Grant College Program Report No. AK-SG-97-01, University of Alaska Fairbanks, 1997. Hilborn, R. and M . Mangel (1997). The Ecological Detective,  Confronting  Models with  Data. Princeton University Press. Hilborn, R. and C. Walters (1992)., Quantitative  Dynamics,  & Uncertainty.  Fisheries  Stock Assessment:  Choice,  Chapman & Hall.  Hilborn, R. (1985). Fleet Dynamics and Individual Variation: Why Some People Catch More Fish than Others. Can. J. Fish. Aquat. Sci. 42, 2-13. Holling, C. S. (2001). Understanding the complexity of economic, ecological, and social systems. Ecosystems 4, 390-405. Koeller, P. A . (2000). Relative importance of abiotic and biotic factors to the management of the Northern shrimp (Pandalus borealis) fishery on the Scotian shelf. J. Northw.  Atl. Fish. Sci. 27, 37-50.  LeBlonde, P. H . , B . M . Hickey, and R. E . Thomson (1986). Runoff driven coastal flow off British Columbia. NATO ASI Series G7, 309-317. Lilly, G . R., D. G . Parsons, and P. J. Vietch (1998). Spatial structure of northern shrimp (Pandalus borealis) off Labrador and eastern Newfoundland (northwest Atlantic). Can. J. Fish. Aquat. Sci. 125, 265-271. Lilly, G., D. Parsons, and D. Kulka (2000). Was the increase in shrimp biomass on the northeast Newfoundland shelf a consequence of a release in predation pressure from cod? J. Northw. Atl. Fish. Sci 27, 45-61. Mackas, D. L . , R. Kieser, M . Suanders, D. R. Yelland, R. M . Brown, and D. F . Moore (1997). Aggregation of euphausiids and Pacific Hake (Merluccius productus) along the outher continental shelf off Vancouver Island. Can. J. Fish. Aquat. Sci. 5J 2080-2096. h  Mackas, D., K . Denman, and A . Bennett (1987). Least square multiple tracer analysis of water mass composition. J. Geophys. Res. 92, 2907-2918. Mackas, D., R. Thomson, and M . Galbraith (2001). Changes in the zooplankton community of the British Columbia continental margin, 1985-1999, and their covariation with oceanographic conditions. Can. J. Fish. Aquat.Sci. 58, 685-702.  142 M a n t u a , N . J . , H . S. R . , Y . Zhang, W . J . M . , and R . C . Francis (1997). A Pacific interdecadal climate oscillation w i t h impacts of salmon production. Bull. Amer. Meteorol. Soc. 78, 1069-1079. M a r t e l l , S., J . Boutillier, H . Nguyen, and C . Walters (2000). Reconstructing the offshore Pandalus jordani trawl fishery off the west coast of Vancouver Island and simulating alternative management policies. Canadian Stock Assessment Secretariat 2000/149, 1-38. M c G o w a n , J . A . , D . R . C a y a n , and L . M . D o r m a n (1998). Climate-Ocean V a r i a b i l i t y and Ecosystem Response i n the Northeast Pacific. Science 10-Jul 281, 210-217. Myers, R . A . , K . G . Bowen, and N . J . B a r r o w m a n (1999). M a x i m u m reproductive rate of fish at low populations sizes. Can. J. Fish. Aquat. Sci. 56, 2404-2419. M y s a k , L . (1986). E l Nino, interannual variability and fisheries i n the northeast Pacific Ocean. Can. J. Fish. Aquat. Sci. 43, 464-497. Paloheimo, J . and L . Dickie (1964). A b u n d a n c e and fishing success. Rapp. P. -V. Reun Cons. Int. Explor. Mer. 155, 152-163. Parrish, R . , C . Neson, and A . B a k u n (1981). Transport mechanismsand reproductive successof fishes i n the California Current. Biol. Oceanogr. 1, 175-203. Pauly, D . , V . Christensen, and C . Walters (2000). " E c o p a t h , Ecosim, and Ecospace as tools for evaluating ecosystem impact of fisheries". ICES J. Mar. Sci. 57, 697-706. Press, W . , S. Teukolsky, W . Vetterling, and B . F l a n n e r y (1996). Numerical Recipies in Fortran 77, The Art of Scientific Computing (2nd ed.). Cambridge University Press. Q u i n n , T . , C . T u r n b u l l , and C . F u (1998). A Length-Based P o p u l a t i o n M o d e l for H a r d to-Age Invertebrate Populations. In F . Funk, T . J . Q u i n n , J . Heifetz, J . N . Ianelli, J . E . Powers, J . F . Schweigert, P. J . Sullivan, and C . I. Zhang (Eds.), Fishery Stock Assessment Models, pp. 531-556. A l a s k a Sea G r a n t College P r o g r a m Report N o . A K - S G - 9 8 - 0 1 , University of A l a s k a Fairbanks, 1998. Rexstad, E . A . and E . P i k i t c h (1986). Stomach contents and food consumption estimates of pacific hake, Merluccius productus. Fish. Bull. 84, 947-956. R i n a l d i , S. and M . G a t t o (1978). O n the determination of a commercial fishery production model. Ecological Modelling 8, 165-172. Robinson, C . L . K . and D . M . Ware (1994). M o d e l l i n g pelagic fish and plankton trophodynamics off Southwestern Vancouver Island, B r i t i s h C o l u m b i a . Can. J. Fish. Aquat. Sci. 51, 1737-1751. Robinson, C . L . K . and D . M . Ware (1999). Simulated and observed response of the southwest Vancouver Island pelagic ecosystem to oceanic conditions i n the 1990s. Can. J. Fish. Aquat. Sci. 56, 2433-2443. Rothlisberg, P. C . and C . B . M i l l e r (1983). Factors affecting the distribution, abundance, and survival of Pandalus jordani (Decapoda, Pandalidae) larvae off the Oregon coast. Fish. Bull. U.S. 81, 455-472. Schnabel, R . B . , K . J . E . and B . E . Weiss (1985). A modular system of algorithms for unconstrained m i n i m i z a t i o n . A CM Trans. Math. Software 11, 419-440. Schnute, J . and D . Fournier (1980). A new approach to length-frequency G r o w t h Structure. Can. J. Fish. Aquat. Sci. 37, 1337-1351.  analysis:  143 Schnute, J. (1987). A general fishery model for a size structured fish population. Can. J. Fish. Aquat. Sci. 44, 924-940.  Schweigert, J. (2000). Stock Assessment for British Columbia herring in 2000 and forcasts of the potential catch in 2001. tariat 2000/165, 74.  Candadian  Stock Assessment  Secre-  Sinclair, A . F. (1998). Estimating trends in fishing mortality at age and length directly from research survey and commercial catch data. Can. J. Fish. Aquat. Sci. 55, 1248-1263. Sinclair, A . , S. Martell, and J. Boutillier (2001). Assessment of Pacific cod off the west coast of Vancouver Island and in Hecate Strait, November 2001. Canadian Stock Assessment  Secretariat  2001/159,  60.  Smith, B . , L . Botsford, and S. Wing (1998). Estimation of growth and mortality parameters from size frequency distributions lacking age patterns: the red sea urchin (Strongylocentrotus franciscnus) as an example. Can. J. Fish. Aquat. Sci. 55, 1236— 1247. Smith, C , A . Hill, and M . Foreman (2001). Horizontal transport of marine ogranisms resulting from interactions between diel vertical migration and tidal currents off the west caost of Vancouver Island. Can. J. Fish. Aquat. Sci 58, 736-748. Stefansson, G . (1998). Comparing different information sources in a multispecies context. In Fishery Stock Assessment Models., pp. 741-758. Alaska Sea Grant College Program Report No. AK-SG-98-01, Fairbanks. Stocker, M . , C. Hand, and M . McAllister (1995). Pacific cod. In M . Stocker and J. Fargo (Eds.), Groundfish  stock assessments for the west coast of Canada  in 1994 o-nd  recommended yield options for 1995., Volume 2069, pp. 115-159. Can. Tech. Rep. Fish. Aquat. Sci. Swain, D. and A . Sinclair (1994). Fish distribution and catchability: what is the appropriate measure of distribution. Can. J. Fish. Aquat. Sci 51, 1046-1054. Tanasichuk, R., D. Ware, W . Shaw, and G . McFarlane (1991). Variations in diet, daily ration, and feeding periodicity of pacific hake (Merluccius productus) and spiny dogfish (Squalis acanthias) off the lower west coast of Vacouver Island. Can. J. Fish. Aquat. Sci. 48, 2118-2128.  Thomson, R., B . Hickey, and P. LeBlond (1989). The Vancouver Island coastal current: fisheries barrier and conduit. In R.J. Beamish and G.A. McFarlane [ed.] Effects of ocean variability on recruitment and an evaluation of parameters assessment models. Can. Spec. Publ. Fish. Aquat. Sci. 108.  used in stock  Tsou, T.-S. and R. M . Royall (1995). Robust Likelihoods. J. Amer. soc. 50(429), 316-320.  Stat. As-  Utter, F. (1971). Biochemical polymorphisms in Pacific hake (Merluccius productus). Cons. Perm.  Int. Explore.  Merr Rapp. P.-V.  Reun. 161, 87-89.  Venier, J. and J. Kelson (1996). The demersal fish box. In D. Pauly and V . Christensen (Eds.), Mass-balance  models of North-eastern  Pacific  ecosystems., Volume 4, pp.  67. Fisheries Centre Research Reports. Walters, C. J. and R. Bonfil (1999). Multispecies spatial assessment models for the British Columbia groundfish trawl fishery. Can. J. Fish. Aquat. Sci. 56, 601-628.  144 Walters, C. J. and J. S. Collie (1988). Is research on environmental factors useful to fisheries management? Can. J. Fish. Aquat. Sci. 45, 1848-1854. Walters, C. J. and F. Juanes (1993). Recruitment limitation as a consequence of natural selection for use of restricted feeding habitats and predation risk taking by juvenile fishes. Can. J. Fish. Aquat. Sci. 50, 2058-2070. Walters, C. J. and S. J. D. Martell (2002). Stock assessment needs for sustainable fisheries management. Bull. Mar. Sci. 70(2), 629-638. Walters, C , V . Christensen, and D. Pauly (1997). Structuring dynamic models of exploited ecosystems from trophic mass-balance assessments. Reviews in Fish Biology and Fisheries 7, 139-172. Walters, C. and J . F. Kitchell (2001). Cultivation/depensation effects on juvenile survival and recruitment: implications for the theory of fishing. Can. J. Fish. Aquat. Sci. 50, 1-12. Walters, C. and D . Ludwig (1994). Calculation of Bayes posterior probability distributions for key population parameters. Can. J. Fish. Aquat. Sci. 51, 713-722. Walters, C. and A . M . Parma (1996). Fixed exploitation rate strategies for coping with effects fo climate change. Can. J. Fish. Aquat. Sci. 53, 148-158. Walters, C , D . Pauly, V . Christensen, and J. F. Kitchell (2000). Representing Density Dependent Consequences of Life History Strategies in Aquatice Ecosystems: Ecosim II. Ecosystems 3, 70-83. Walters, C , M . Stocker, A. Tyler, and S. Westrheim (1986). Interaction between pacific cod (Gadus macrocephalus) and herring (Clupea harengus pallasi) in the Hecate Strait, British Columbia. Can. J. Fish. Aquat. Sci. 43. Walters, C. (2000). Natural selection for predation avoidance tactics: implications for marine population and community dynamics. Mar. Ecol. Prog. Ser. 208, 309-313. Ware, D. M . and G . A . McFarlane (1995). Climate-induced changes in Pacific hake (Merluccius productus) abundance and pelagic community interactions in the Vancouver Island upwelling system, p. 509-521. In R. J. Beamish [ed.] Climate change and norther fish populations. Can. Spec. Publ. Fish. Aquat. Sci. 121. Ware, D. M . and R. E . Thomson (2000). Interannual to multidecadal timescal climate variations in the Northeast Pacific. Journal of Climate 13, 3209-3220. Westerheim, S. (1996). On the Pacific cod (Gadus macrocephalus) in British Columbia waters, and a comparison with Pacific cod elsewhere, and Atlantic cod (Gadus morhua). Can. Tech. Rep. Fish. Aquat. Sci. 2092, 390. Winters, G . H . and J . P. Wheeler (1985). Interaction between stock area, stock abundance, and catchability coefficient. Can. J. Fish. Aquat. Sci. ^2(5), 989-998. Worm, B . and R. A . Myers (In Press). Top-down versus bottom-up control in oceanic food webs: a meta-analysis of cod-shrimp interactions in the North Atlantic. Ecology-  Appendix A Simulation Study for Shrimp Assessment M o d e l Traditional stock assessment models usually assume natural mortality rates remain relatively constant over time, and in Chapter 2, annual natural mortality rates were treated as arbitrary time variables. The stock assessment model was fit to 3 different types of observations: 1) relative abundance data from research trawl surveys, 2) relative fishing mortality rate data form commercial log-book data and area swept calculations, and 3) length-frequency data collected from research trawl surveys. To demonstrate that the assessment model used in Chapter 2 is appropriate for estimating annual natural mortality rates, I used the same model to simulate a set of time series data with known parameters, and then estimated population parameters using these data. Such a method is very important to ensure that the assessment model can estimate true parameter values, i.e. investigate potential bias (Hilborn and Mangel 1997). Length frequency data, much like catch-at-age data, provide information about relative proportions of each cohort (Hilborn and Walters 1992). Using a time series of length frequency data, an index of recruitment is apparent as new modes appearing in the data, and mortality is apparent by how quickly modes disappear over time. Using such data to investigate changes in natural mortality rate over time assumes that growth 145  146 over time is constant, that size selectivity remains constant, and that sampled length frequency distributions are representative of the entire p o p u l a t i o n . To determine if the assessment model was appropriate for estimating changes i n n a t u r a l mortality rates over time, I used the assessment model described i n Chapter 2 to generate 21-year time series of relative abundance data, relative fishing mortalities, and length frequency data, then estimated population parameters using the assessment m o d e l w i t h simulated data. Relative abundance data were assumed to have log-normal errors, and length frequency d a t a were taken by randomly sampling from the true m u l t i n o m i a l distribution i n each year. A n example of a length frequency sample over time is shown i n Figure A . l .  5  °  I i i—i—r 10 15 20 25 30 35 40 2  I—l^h—I—I I I 10 15 20 25 30 35 40  I T 1 1—^ 10 15 20 25 30 35 40 3  — i — i — i — i i 10 15 20 25 30 35 40 10  10 15 20 25 30 35 40 17  10 15 20 25 30 35 40 4  10 15 20 25 30 35 40 11  i—I—I—T I i 10 15 20 25 30 35 40 18  l—"PH—i—r *!—i  l i i — i — r ~ l — i 10 15 20 25 30 35 40 12  I I—I—I f I—I 10 15 20 25 30 35 40 19  \—r~]—I I I 10 15 20 25 30 35 40  r~n—I—I II i 10 15 20 25 30 35 40  J  5  10 15 20 25 30 35 40 5  U-  r"-"T—T—i—i I i 10 15 20 25 30 35 40 6 o  -I..  I i—i—r 10 15 20 25 30 35 40 7 I i — i — i — r 10 15 20 25 30 35 40  Length (mm)  Figure A . l : E x a m p l e of an observed length frequency d i s t r i b u t i o n and predicted frequencies from stock assessment model. Observed frequencies were d r a w n at random from the true m u l t i n o m i a l distribution.  If the model is working correctly, it should be able to accurately fit simulated datasets that have no observation errors, and correctly estimate recruitment anomalies and anomalies i n natural mortality rates (nuisance parameters) over time. T h i s was the case,  147 as shown in Figure A.2. Two penalty functions were used in Chapter 2 to help define the underlying mean estimates for annual recruitment and annual natural mortality rates. If these penalties are not included, time series anomalies can have a non-zero sum, which may result in biased parameters estimates.  For example, in Table A . 1 , estimates for mean annual  recruitment (R ) are over-estimated when the sum of recruitment anomalies can assume a  any value, whereas, estimates are slightly more accurate when penalties are included. Penalties on nuisance parameters are a convenient way of handling over-parametrization problems, but they may also mask the true underlying distribution for other parameters such as growth rates, or selectivity parameters. For example, when penalties are used to constrain nuisance parameters in the exercise carried out in Table A.1, other estimated parameters are biased.  Time (years)  <= - i  Time (years)  ,  5  I  I  I  5  10  15  Time (years)  -|  I—  1  1  1  20  5  10  15  1— 20  Time (years)  Figure A.2: Simulated time series data from stock assessment model (thin black lines), and reconstructed time series (thick grey lines) with no observation errors. Bt is biomass over time, Ft is annual fishing rate, Rt is annual recruitment, and M t is annual natural mortality rate.  When observation errors are included in the simulated data set, then it is not expected  148 Table A . 1 : True and estimated parameter values for catch-at-length model. True values were used to simulate time series d a t a w i t h no observation errors, and estimated values w i t h out penalty functions on recruitment anomalies, and natural m o r t a l i t y rate. Parameter True Value Estimate w / o penalty E s t i m a t e w penalty 26.00 R 33.12 38.36 M 0.7 0.699 0.698 15.03 15 14.99 u 30 29.99 30.01 0.685 0.68 0.679 P CV 0.08 0.08 0.079 0.775 0.8 0.803 7 18.114 18 17.998 L Q  H  that the exact parameter values w i l l be estimated correctly given a small time series. However, over an infinitely long time scale, or multiple simulated d a t a sets, on average the estimation procedure should reflect the true parameter values (i.e. unbiased). T o ensure the assessment model used i n Chapter 2 was unbiased, 250 independent data sets were simulated, generating 250 estimates of each model parameter.  I then plotted the  distribution of a l l estimated parameter to look for systematic bias (Figure A . 3 ) . W h e n length frequency samples are small and uninformative about changes i n mortality rates, auxiliary information about relative trends i n fishing m o r t a l i t y rates generally improves estimates of true population parameters, especially changes i n n a t u r a l mortality. For the simulations shown i n Figure A . 3 , annual samples of 10,000 were taken and were very informative about changes i n age structure over time. Smaller samples, or unrepresentative samples, led to biased parameter estimates. A disproportionate representation of smaller length classes w i l l over-estimate annual recruitment and under-estimate natural mortality, and vice versa when larger shrimp are over represented i n the sample. A s discussed i n Chapter 2, parameter correlation is of concern, especially between growth parameters and selectivity parameters.  For computational convenience I used  a variant of the von Bertalanffy growth equation to model length growth as L  A  L\ + (LN — L\) lZ N-\, P  where p = e~ , k  =  L\ and LN correspond to mean length of age  149 Y  Figure A.3: Box plots of parameter estimates from 250 simulation-reconstructions using the catch-at-length model with auxiliary information about trends in fishing mortality (right box-whisker plots), and no information about trends in fishing mortality (left box wisker plots). Plots titles correspond to parameters listed in Table 1. Hashed line is the true parameter value. 1 and age N individuals (Schnute and Fournier 1980). In this parameterization, the growth coefficient was positively correlated with the length at 50% vulnerability (Lh), and positively correlated with the shape parameter 7 (Figure A.4). The scaling parameter R appears to be positively correlated with mean natural mortality rates, indicating Q  that the observed data could come from a large population with a low natural mortality rate, or a small population with a high natural mortality rate. Therefore, in chapter 2 I fix the mean natural mortality rate to M — 0.6 and estimate a sequence of time varying anomalies v that were constrained to (VJ v — 0). t  t  t  150  0.68  Ro  m  0.72  29.95  30.10  0.0795  an « n in  17.8 18.6 I I I I,  fifl  i  M •  • — •  L1  X /  Ln  *1 I  '  rho  •  :  ML  H  •  I  S  CV1  • •  gamma  •  •  — _  • •  •r-  I I IT "II I 29 33  14.95  0.67 0.70  Ih  I I I I I 0.65 0.80  Figure A . 4 : Parameter correlations for 250 simulation-reconstructions using the catchat-length m o d e l described i n Chapter II. Variables names defined i n Table A . l .  Appendix B Lingeod Assessment Commercial catch statistics were used to reconstruct biomass andfishinghistory for lingcod in the southern Vancouver Island region (Table B.l). The assessment model described in Chapter 4 was used, and growth parameters for lingeod were estimated from mean length-at-age and weight-at-age presented in Cass et al. (1990). A prior distribution with a mean of 3.5 and a variance of 0.2 was used for the recruitment compensation parameter. Furthermore, parameters were estimated in two phases, in which thefirstphase included a penalty on the mean exploitation rate (0.25) to approximate correct scaling parameters (B ). In the second phase, this penalty was relaxed and nuisance parameters Q  (observation errors and process errors) were included in the estimated parameter set. I further assumed that natural mortality rates for lingeod remained constant over time, and natural mortality is equal to the von Bertalanffy growth coefficient k (Table B.2). Reconstructed biomass,fishingmortality rates, and age-2 recruits are shown in Figure B.2. Intensivefishingfor lingeod on the west coast of Vancouver Island began during the development of the groundfish trawlfisheryin early 1940's (Cass et al. 1990). Prior to the development of the trawlfisherylingeod werefishedquite intensively with handlines, and catches off the southwest coast of Vancouver island averaged near 2000 t/yr. By the 1950's average catches had declined and it is clear that the stock was well below its 151  152 Table B.l: Commercial catch statistics for lingeod off southwest Vancouver Island (Source: King and Surry, 2000). Year Catch (kg) CPUE (kg/hr) Effort (E ) Et/E Mean 1954 371469 123 3020 1.52 1955 539578 237 2277 1.14 1956 570520 310 1840 0.92 1957 431860 298 f449 0.73 1958 386548 356 1086 0.55 1959 284160 312 911 0.46 1960 353847 214 1653 0.83 1961 497840 312 1596 0.80 1962 180630 227 796 0.40 1963 171471 345 497 0.25 1964 309924 419 740 0.37 1965 512711 368 1393 0.70 1966 497307 308 1615 0.81 1967 716983 423 1695 0.85 1968 857458 603 1422 0.71 1969 472839 289 1636 0.82 1970 509759 282 1808 0.91 1971 649029 358 1813 0.91 1972 473373 286 1655 0.83 1973 607546 351 1731 0.87 1974 709254 303 2341 1.18 1975 1103985 370 2984 1.50 1976 582014 248 2347 1.18 5.09 1977 560852 288 1947 0.98 3.99 1978 294821 226 1305 0.66 5.16 1979 446212 260 1716 0.86 4.85 1980 363248 283 1284 0.65 3.95 1981 406054 262 1550 0.78 4.48 1982 1208108 412 2932 1.47 4.68 1983 558030 372 1500 0.75 4.52 1984 997055 445 2241 1.13 4.58 1985 1894282 531 3567 1.79 4.55 1986 455297 380 1198 0.60 1987 348270 206 1691 0.85 4.75 1988 475035 208 2284 1.15 4.26 1989 837611 316 2651 1.33 4.78 1990 1155190 278 4155 2.09 5.06 1991 1235360 260 4751 2.39 5.34 1992 961945 307 3133 1.57 4.55 1993 1413218 607 2328 1.17 1994 683642 219 3122 1.57 3.86 1995 78570 234 336 0.17 1996 747834 143 5230 2.63 6.11 1997 460515 226 2038 1.02 1998 523038 437 1197 0.60 5.09 1999 249687 234 1067 0.54 t  153 Table B . 2 : E s t i m a t e d growth parameters from length and weight-at-age data from Cass et al,(1990).  Parameter  Description A s y m p t o t i c length (cm) G r o w t h coefficient Length-weight scalar Length-weight power  Loo k a b  Value 106.3 0.212 1.12E-05 3.002  unfished state. T h e ratio of 1950's biomass relative to its unfished state was estimated at 29% i n this assessment. F r o m 1950 to 1969, biomass increased from as 9,500 tonnes to « 12,500 tonnes.  T h e most recent biomass estimate is roughly 9,100 tonnes.  The  retrospective projections show that more recent d a t a has resulted i n lower estimates of unfished biomass and p r o d u c t i v i t y (Table B . l ) . T h i s is likely due to the additional contrast that has been provided by recent declines i n stock abundance.  1 1960  1  1  1970  1980  1 1990  r 2000  Year  Figure B . l : Retrospective projection of lingcod biomass.  F i s h i n g m o r t a l i t y rates during the 1950-2000 period ranged from 0.01 i n 1995 to 0.18 i n 1996, and averaged roughly 0.07 during this p e r i o d (Figure B . l ) . Note that this is an adjustment downward after relaxing the penalty on mean fishing mortality rates. The early C P U E data indicate a strong rebuilding of the stock from 1950 to 1968, and here I made the very dangerous assumption that the qualified C P U E index has been proportional to abundance. A more realistic scenario might include an adjustment i n the  154  Figure B . 2 : Reconstructed biomass, observed landings, fishing mortality rates, residuals i n catches and mean body weights, and estimates of age-2 recruits for southwest Vancouver Island lingcod. catchability coefficient over time to reflect changes i n fishing technology that make the capture process more efficient. D u r i n g the m i d 1970's and early 1990's there was intensive fishing activity, and lingcod abundance declined markedly during these two periods.  Appendix C Using changes i n diet composition to estimate vulnerability One principal Ecosim parameter that is often debated is the vulnerability exchange rate parameter  (t>,j).  The parameter represents a behavioral process for which there are  almost no direct observations. In fact, it is probably impossible to directly measure this process because it requires that the observer knows when an individual is in a vulnerable and invulnerable state. There are methods however that allow for indirect measurement (4 to be exact), all of which are documented in the Ecopath Software manual. One method is to use Ecopath models from two different time periods, and determine what  Vij parameters were required to get from state 1 to state 2. If gross changes in prey biomass and predatory biomass, as well as changes in diet compositions, had occurred between the two states, then the data (as represented by two different mass-balance snapshots) were informative about v^.  The method outlined here works on the same  principle, however, rather than specify two Ecopath models, it is possible to actually fit Ecosim to changes in diet compositions over time. The predicted proportion of prey i in predator j's diet at each time step t is represented in Ecosim by: ^_  ij,t(Bj,t> Bj,t)  c  155  ^  c  156 Note that equation C l assumes any fixed imported diet has been scaled out in the denominator, and changes in consumption of i by j are within the defined system only. Here consumption of i by j is calculated using equation 5.5 at each time step, and is a function of Vij and a^. If  could be measured, or determined independently then we  could just solve equation 5.14 for v j and use this as an estimate. However, t  inseparable in consumption equations, and  and i»y are  actually represents a ratio (to be specified  by the user) of the base consumption of i by predator j specified in the Ecopath model, where ay is determined by this ratio (see Walters et al. 1997). Proportions of prey i in the diet of predator j, integrated over all size classes, is assumed to have a multinomial distribution with a mean an variance, and ignoring a scaling constant the log-likelihood of observed diet composition at time t is:  L(DCi \Bi , jit  it  B , ciij^ij) = rij * DC log (DC ) jit  ij>t  ijit  where  ]P DC = 1 (C.2) ij>t  i  where DCijj is an observed fraction of prey i in stomachs of predator j by weight and rij is the number of stomach samples. Information about changes in diet composition over time in this case should help resolve parameter confounding between predator-prey biomass and v^.  

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 20 1
China 16 20
Canada 12 0
Belgium 8 0
France 7 0
Netherlands 4 0
Russia 2 0
Germany 2 4
Poland 1 0
India 1 0
City Views Downloads
Unknown 19 5
Shenzhen 13 20
Ghent 8 0
Ashburn 7 0
Amsterdam 4 0
Wilkes Barre 3 0
Vancouver 3 0
Wilmington 3 0
Edmonton 2 0
Beijing 2 0
Mumbai 1 0
New York 1 0
Mountain View 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0074868/manifest

Comment

Related Items