VARIATION IN PINK SHRIMP POPULATIONS O F F T H E W E S T COAST OF 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 T R O P H I C INTERACTIONS by S T E V E N J A M E S D E A N M A R T E L L B.Sc , 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 R E Q U I R E M E N T S FOR T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y 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 inves- tigate 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 oceano- graphic indices as environmental correlates did not explain much of the residual variation in the stock-recruitment relationships. Members of the family Gadidae primarily con- sume 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 Pa- cific 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 pri- mary 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 xvi 1 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 2 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 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 2.4 Discussion 55 3 Recruitment and Oceanography 60 3.1 Introduction 60 iii iv 3.1.1 Summer and Winter Current. Ci rculat ion Patterns 61 3.1.2 Dis t r ibut ion of Shrimp Relative to Coastal Currents and Shelf Break Currents 62 3.2 Methods 64 3.2.1 Environmental Indices 67 3.3 Results 69 3.4 Discussion 72 4 Shr imp Preda t ion by Gadids 77 4.1 Introduction 77 4.1.1 Stock Structure of Pacific Cod 78 4.1.2 Stock Structure of Pacific Hake 80 4.2 Est imat ing Predator Biomass 83 4.2.1 Pacif ic C o d Assessment 83 4.2.2 Pacif ic Hake Assessment 88 4.3 Results 90 4.3.1 Pacific Cod Biomass 90 4.3.2 Corre lat ion wi th Predator Abundance 93 4.4 Discussion 95 5 Ecosystem Perspectives 99 5.1 Introduction 99 5.2 Methods 101 5.2.1 Ecopath wi th Ecositn 101 5.2.2 F i t t i ng Ecos im To T ime Series Data 107 5.2.3 Descr ipt ion of Mode l Groups 110 5.2.4 T i m e Series Da ta for West Coast Vancouver Island 11G 5.3 Results 119 5.3.1 F i t t i ng Ecos im to T ime Series Data 119 5.3.2 Changes in Shrimp Mortal i ty 126 5.4 Discussion 130 6 General Discussion and Summary 134 Bib l iography 139 A S imula t ion S tudy for Sh r imp Assessment M o d e l 145 B L ingcod Assessment 151 U s i n g changes i n d ie t c o m p o s i t i o n to e s t i m a t e v u l n e r a b i l i t y List of Tables 2.1 Mode l parameter descriptions and symbols used in populat ion dynamics model 14 2.2 Parameter descriptions and symbols used in length frequency observation model 17 2.3 Number of shr imp fishing vessels reporting f ishing act iv i ty from 1987 to 2000. A vessel reporting > 100 shots per year are considered ful l time, and is roughly equal to 20 days of fishing act iv i ty per year 20 2.4 Est imated growth parameters for area 124 and 125 pink shrimp populations. 32 2.5 Est imated model parameters and standard deviations for stat ist ical areas 124 and 125 32 2.6 Est imated recruitment anomalies («>,) and standard deviations for statis- t ical areas 124 and 125 33 2.7 Est imated deviations in mortality (;/t) and standard deviations for statis- t ical areas 124 and 125 34 2.8 Reported fishing effort in commercial logbooks for stat ist ical areas 124 and 125 44 2.9 Natura l mortal i ty rate statistics for areas 124 and 125 49 3.1 Compar ing the effects of environmental indices on stock recruitment rela- tionships for areas 124-5. loijL — log-l ikelihood, \ 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 ig 3.3 72 vi vii 3.4 Compar ing the effects of environmental indices on stock recruitment rela- tionships for area. 124 spa.wners producing area. 125 recruits and vice versa,. logL = log-l ikelihood, \ 2 — probabil i ty of the of the environmental index relative to base model, and the % of the variance in the residuals explained by the enviroment.al correlate 72 4.1 Publ ished reports of Pacific hake (Merluccius producf;us) predation on pink shrimp (Pandalus jordnni) and reported % occurrence, or % by weight of pink shr imp in hake stomachs 83 4.2 Parameter symbols and description of parameters for the delay-difference model 84 4.3 Biomass Est imates for Pacif ic hake; for L a Peruse bank region (Ware and McFarlane, 1995), and for hake in the Canadian Zone, (Dorn et a l , 1999). 89 4.4 Correlat ion coefficients (r) between pink shrimp abundance, mortal i ty 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 94 5.1 Ecopath input parameters, arid estimated parameters (in bold) for west coast Vancouver Island model I f2 5.2 Diet matr ix for Ecopath model 114 5.3 Squared residuals for No Mode l (or deviations from the mean of the data), a model w i th no trophic' interactions (assume that predation mortal i ty and consumption are constant), a, model wi th trophic interactions where pr i - mary product ion is held constant over t ime, and a, model that includes esti- mated pr imary production anomalies (see Figure 5.7). Posit ive reductions in the square residuals (right hand column) between trophic interactions only model and trophic interactions wi th variable pr imary product iv i ty indicate a better 'fit' to the observed data 124 5.4 Corre lat ion coefficients (/•) between hake abundance time series and Z differences from single species model and ecosystem model (Data shown in Figure 5.12) 129 A . l True and estimated parameter values for catch-at-length model. True values were used to simulate time series data, wi th no observation errors, and estimated values with out penalty functions on recruitment anomalies, and natural mortal i ty rate 148 viii B . l Commercia l catch statistics for lingcod off southwest Vancouver Island (Source: K i n g and Surry, 2000) 152 B.2 Est imated growth parameters from length and weight-at-age data from Cass et al,(1990) 153 List of Figures 1.1 Stat ist ical area definitions used for smooth pink shr imp trawl fisheries. . 4 1.2 P ink shr imp Pandalus jordani landings in statist ical areas 124 and 125 from 1975 to 2000 5 1.3 Survey biomass estimates for Pandalus jordani in areas 124 & 125 from 1973 to 2001. Note no surveys wore conducted in 1974. 1984, and 1986 in both areas and no surveys were conducted in area 125 in 1989 and 1991. Biomass estimated as mean biomass per area, swept by survey times total area surveyed 8 1.4 Locat ion of research trawl surveys from 1973 to 2001. Area 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 8 2.1 Reconstructed biomass and observed biomass indices for statist ical 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 shr imp t rawl surveys 27 2.2 Residuals for area 124 and 125 between observed biomass indices and model predicted indices (lo(j(Z)). Note that positive residuals indicate observed biomass index is greater then predicted biomass 28 2.3 Populat ion age structure for pink shrimp in areas 124 and area 125. Wh i te shaded polygon are age-1 shrimp and black shaded polygons axe age-4 shrimp 29 2.4 Predicted (solid line) and observed (histograms) length frequencies for area 124. Est imated mean length-at-age is represented by vert ical dotted lines. Note there are no length frequency data for 1983-84, and 1986, and predicted distributions are shown as a, straight line 30 ix X 2.5 Predicted (solid line) and observed (bars) length frequencies for area 125. Est imated mean length-at-age is represented by vertical dotted lines. No length frequency samples were taken in 1984, 1986, 1989, and 1991 in area 2.6 Pai r plot of 500 random sample estimates (/?„, L\. LA. p, A t , A2,'>, h) taken from an MCMC sample of length 250.000 for area 124 35 2.7 Pa i r plot of 500 random sample estimates (/?„, L j , 1.4, p, Ai, A 2 , 7 , lh) taken from an M C M C sample of length 250,000 for area 125 36 2.8 Marg ina l posterior distr ibutions for A 2 in statist ical areas 124 and 125. Note that observed length frequency data are relatively uninformative about how standard deviation in length-at-age changes w i th age 37 2.9 Sensit ivi ty of M to <r'f[ and a"t,. Left column represent the mean estimate of natural mortal i ty as af, increases from 0.1 to 1, where the mean (line wi th points) is calculated over all a'2M. Right column represents the mean estimate of M as a2w increases from 0.8 to 8, where the mean is calculated over al l cr'2. In the left columns, dotted lines represent estimates of M when al, = 0.8 and dashed line represents estimates of M when a2r = 8, and in the right column dotted lines for a'jf = 0.1 and dashed line for o2M = 1 38 2.10 Sensit ivi ty of B„ to a\t and <T2,. Left column represent the mean RD as <j2M increases from 0.1 to 1, where the mean (line wi th points) is calculated over al l a\f. Right column represents the mean estimate of R0 as cr2,, increases from 0.8 to 8, where the mean is calculated over al l cr'2. Dotted lines in left column represent estimates of R„ when cr2, = 0.8 and dashed line represents estimates of R() when cr2, = 8, in the right column dotted lines represent a2M = 0.1 and dashed line o\, = 1 39 2.11 Contour plots for objective function values ( C - 'V"«'- ' ( L )) versus a2M and a2, values. A rea 125 clearly shows 3 local maxima 40 2.12 Est imated area over which shrimp populations in stat ist ical areas 124- 125 were distr ibuted. Area, estimates were calculated using commercial logbook data, from 1987 to 2000; these data are shown in Figure 2.14 . . 41 2.13 Relat ionship between shrimp abundance and stock area from 1987 to 2000. Shr imp abundance for each statist ical area is represented by the total abundance (D, = N „ . jW a ) . Note that abundance shown in figure 2.1 represents the vulnerable biomass 42 2.14 Spat ial distr ibut ion 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 Est imated selectivity curves for areas 124 and 125, where dotted vertical line represents the length at 50% vulnerabil i ty 2.16 Est imated fishing mortality rates from stock assessment model (lines in bottom panels) for areas 124-125. Circles represent area swept estimates of fishing mortal i ty rates (F = ,J~-) from commercial log book data, and these data were used in the likelihood criterion for est imating model pa- 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 likel ihood estimates of natural mortal i ty rates (line wi th circles) from stock assessment model for areas 124 (top) and 125 (bottom) from 1975 to 2000. Grey points are a 1.000 randomly chosen samples from the posterior distr ibut ion generated by the Metropolis Hastings A lgor i thm, and the box plots represent the median, inner quartile, and sampling range 2.18 M a x i m u m l ikel ihood estimate of age-f recruits (line) for stat ist ical areas 124 and 125 from 1975 to 2000 and associated uncertainty (box plots). Grey points are a 1000 random samples from the posterior distr ibut ion, and box plots represent median, inner quartiles and range of posterior samples. Average annual recruitment for area 124 and 125 is 623 and 704 mil l ion shrimp, respectively 2.19 Relat ionship between carapace length arid proport ion 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 If2,732 shrimp measured and sexed was chosen for estimating parameters for a logistic curve. Parameters were estimated us- ing a least squares cri terion, where the maximum l ikel ihood estimate for 7 = 0.89 and lh = 19.71 2.20 Spawner biomass (t) at t — 2 years versus age-1 recruits for statist ical ar- eas 124 (top-left) and 125 (bottom-left) and corresponding Beverton-Holt (solid line) and 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 xi i 3.1 Tow positions of research trawl surveys from 1973 to 2001. Stat ist ical areas 125, 124, 123 and 121 are shown in black, to lighter shades of grey, respectively. Depth contours are LOOm, 200m and lOOOin 63 3.2 Deviat ions in annual recruits per unit of spawniug biomass (log(R, jS\) - log(Rt/St)) in statist ical areas 124 and 125. The dotted line is 3 standard deviations for the distribution and points beyond this line were considered outliers. Extremely high juvenile survival rates occurred area 124 and 125 in 1999 and 2000, respectively 65 3.3 Standardized oceanographie indices used to compare as an environmental correlate wi th shrimp recruitment in area 124 on west coast Vancouver Island. S S T = sea surface temperature, S L H = sea. level height at Bamfield Mar ine Stat ion, SOI = southern oscil lation index, P D O = Pacif ic decadal oscil lation index. T ime series spans 1977-2001 68 3.4 Beverton-Holt and Ricker models fit to stock-recruitment data from sta- t ist ical area 124 using a. normal l ikel ihood (Eq. 3.5) and robust l ikel ihood (Eq. 3.7) Each model line represents the mean stock-recruitment relation- ship wi th different environmental indices included as correlates 70 3.5 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 th correlate included. Left panels are Beverton-Holt models, right are Ricker models 73 3.6 Surv ival index (log(Rt/St)) for age-1 recruits in area 124 and area 125. Dotted line represents mean log(R,/S) 75 4.1 Groundf ish trawl catch per unit effort of Pacif ic cod off the west coast of Vancouver Island (scale units are log(e.p.u.e). Areas 3C and 3D cor- respond to groundfish management areas; shrimp management areas are also shown. M a p courtesy of A lan Sinclair, DFO Pacif ic Region 79 4.2 Dis t r ibut ion of Pacif ic hake in 1992, 1995, and 1998 along survey tran- sect lines. Abundance is measured using hydro acoustic backscattering techniques. (Figure supplied by Mark Saunders, D F O ) 82 xi i i 4.3 Proport ion of Pacific hake in the Canadian zone versus mean January- February sea level height anomalies (top panel). Posit ive anomalies indi- cate higher than average1 sea level heights in Bamfie ld, Br i t i sh Co lumbia (data provided by Mark Saunders, U F O ) . Est imated hake biomass (bottom panel) for total hake stock (solid line) and hake present in the Canadian zone (dashed line) based on the sea level height relationship 90 4.4 Est imated biomass time series of age 2+ Pacif ic cod and fishing mortal i ty 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. Bo t tom row from left to right are recruitment anomalies, catch residuals, and mean weight residuals from the delay difference model 91 4.5 L ikel ihood profiles for W C V I Pacific cod population parameters where assumed variance's for natural mortality rates a2M — 0.38 and recruitment compensation <r| = 0.2 (solid line) versus doubling these variance terms (dashed lines) 93 4.6 Relat ive abundance time series of P ink shrimp, Pacif ic C o d , and Pacif ic hake series I and II from 1975 to 2000. Lower panels are scatter plots of predator abundance versus pink shrimp abundance 95 5.1 A n example of the relationship between predation mortal i ty rate and predator biomass for different values of klr The dark circle represents the base predation mortal i ty rate m ,-,•,<>, calculated using Ecopath input parameters. As the; value of klt decreases the slope of the line through Bj.o,niijs) decreases and the maximum predation rate (vtJ) of predator j on prey i decreases 107 5.2 Food web and fisheries for west coast Vancouver Island ecosystem model. The size of the box is proport ional to biomass of each group, and the width of the connecting line is proportional to the How. Posi t ion on the y-axis is the trophic level of the group 116 •5.3 Commerc ia l landings (t km""2) for west coast Vancouver Island from 1950 to 2000 (left column). Tota l mortality rate estimates for pink shr imp and adult Pacif ic herring. P ink shrimp total mortal i ty rate estimates were calculated from combined fishing and natural mortal i ty rate estimates in chapter 2. Total mortal i ty rates for Pacif ic herring were inferred from age composit ion data,. Bo t tom panel is mean body weights for Pacif ic cod, where mean body weight, is calculated using: \Vf = Bt/Nt 117 xiv 5.4 Relative abundance time series used as observational data for the ecosys- tem model. Shr imp and eulachon biomass estimates are combined for two statistical areas, herring abundance was constructed using a. catch-at-age model, l ingcod abundance from stock assessment model fit to commercial C P U E data, and Pacific cod from assessment in Chapter 4. Euphausi id abundance (from 1967 to 1998) estimated from a, trophodynamic model by Robinson arid Ware (1999) 118 5.5 Predicted (points) and observed (vertical bars) landings for pink shr imp, herring, l ingcod and Pacific cod stocks off the west coast Vancouver Island shelf. F ish ing mortal i ty rate estimates for al l four species were derived from single species assessment models and used as input data into Ecosim. 120 5.6 Predicted biomass t ime series (line) from Ecos im and observed relative abundance t ime series (circles). Note that relative abundance t ime series are assumed proport ional to true abundance, where Yt = qBt 122 5.7 Est imated pr imary production anomalies from f 950 to 2000; where a spline function was fit to 30 equally spaced data points. The use of the spline function helps filter out high frequency variat ion that occurs in f i t t ing the model for years having only one or few times series wi th large anomalies. 123 5.8 Compar ison between estimated annual mean pr imary product ion anoma- lies from Ecos im, and calculations of Dia tom product ion (g C m2yr~x) obtained by Robinson and Ware (1999) using a, model based on upwelling (nitrogen), water temperature, arid solar radiat ion 125 5.9 Compar ison between estimated annual mean primary product ion and the negative Pacif ic Decadal Osci l lat ion ( P D O ) 126 5.10 A l l mortal i ty components from Ecosim compared to Z from single species model (SSM) . The area between each line is the total mortal i ty due to each predator or fishing, and the area between total mortal i ty and fishing mortal i ty is the unexplained mortal i ty (1 — EE) 127 5.11 Standardized differences between Z estimates from Ecosim versus Z es- timates from single specie's model in Chapter 2. Posit ive values indicate years where Ecosim predicted higher mortal i ty rates than estimated using the S S M 128 5.12 T ime series of standardized indices of hake abundance and residuals be- tween ZEWE and ZSSM- All indices have a mean 0 and a — 1. Hake series 1 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). . . . 129 XV A.1 Example of an observed length frequency distr ibut ion and predicted fre- quencies from stock assessment model. Observed frequencies were drawn at random from the true multinomial distr ibut ion 1 46 A.2 Simulated time series data from stock assessment model (thin black lines), and reconstructed time series (thick grey lines) wi th no observation errors. B t is biomass over time. F t is annual fishing rate, Rt is annual recruitment, and Mt is animal natural mortal i ty rate 147 A . 3 Box plots of parameter estimates from 250 simulation-reconstructions us- ing the catch-at-length model wi th auxi l iary information about trends in fishing mortal i ty (right, box-whisker plots), and no information about trends in fishing mortal i ty (left box wisker plots). Plots 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 catch- at-length model described in Chapter II. Variables names defined in Table A.1 150 B. l Retrospective projection of lingeod biomass 153 B.2 Reconstructed biomass, observed landings, fishing mortal i ty 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 support of my parents, whom not only supported me financially when times were tough but morally. Just the thought, alone of making my parents proud is accomplishment enough. To my grandfather, Robert E. Robinson, a former trapper in Br i t ish Columbia. Your historical accounts of fishing and trapping, along with the first rainbow trout you helped me catch, has led me into a, new direction of assessing impacts on natural resources. Your memories wi l l reach far beyond your lifetime. Many thanks to my committee members. Daniel Pau ly and Villy Christensen who helped me develop alternative ideas about ecosystem modeling, as well as the occasional fishing tr ip for harbor seals and sight seeing tours at F lo r ida State Universi ty (thanks Vi l l y ) . John Dower and Dan Ware, for exposing me to a complex area of oceanography and climate affects. A l a n Sine lair for helping me digest the Pacif ic cod history on the west coast of Vancouver Island, and great pasta dinners. A lso, I would like to acknowledge the alternative exposure provided by Ransom Myers, Shelton Harley, Jeff Hutchings, arid Boris Worm (Dalhousie Universi ty). Th is group is leaps and bounds ahead in computing power due to a general ban on Microsoft^ products. It was Shelton Harley who pointed out that, a thesis could be written on free software that is light years ahead of Mircosoft"0. Ana, Parma, who shared many great stories and provided some excellent feedback on manuscripts from this thesis. Jon Sclmute, who xv i xvii f inally explained the black magic behind the Metropol is-Hast ings algori thm, and kept me laughing during P S A R C meetings. R ik Thompson, who sent al l of his reprints on the same flay I requested 'one'. Dave Founder, for gett ing me up to speed on the latest in non-linear parameter estimation software, and a bad hangover in Madison Wisconsin. James Ki tche l l , for the wonderful weekends surfing in Honolu lu, poi, bass fishing in U N D E R . C , Superica in Santa Barbra and the next job. Thanks also to Rob Bison, for some extra $$ and an opportuni ty to give back to the local steelhead community. Many thanks to my supervisor, Carl Walters, the man who defines conservation as "ki l l something else in order to save it", aptly encourages graduate students wi th sound advice such as "don't F9<?KUP", arid has his own icon for evolutionary ecology. In addit ion to witnessing the return of 'Betsy ' , I've also witnessed many other great mystery's wi th Ca r l , such as full moon's in broad day light and bel ly-button lint. I also have to thank Ca r l for the many great, adventures, here in Vancouver, F lor ida, .Japan, Honolulu and the Grand Canyon. Car l and his wife Sandra were kind in providing accommodations for Dawn and myself, unti l 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 am forever grateful for al l the love and support provided by my girlfr iend Dawn Cooper. Thank you so much for all your patience and understanding. Special thanks to my lab mates, Sean Cox, Robert A l l iens , Bob Lessard ,Paul Higgins, R i k Duckworth, Lisa Thompson, Leonardo Hua.to, Nathan Taylor, and Dav id O 'Br ien . 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 im Boutillicr who encouraged me to xviii adopt my own methods tor examining shrimp fisheries, and provided access to data and computing resources. A lso Hai Nguyen and Norm Olsen for assistance w i th databases, and comprehending the data. Rob Kron lund and Marc Saunders for important discus- sions about the problems with assessing Pacif ic hake abundance. F inal ly , the to crew of the C C V W . E . Ricker, the platform for shrimp surveys and the summer feeding grounds for myself. F inac ia l support for this project was provided by Department of Fisheries and Oceans, administered through J im Bouti l l ier. Chapter 1 Introduction "Nothing is more common in the management of resources than surprise, usually unpleasant surprise." ... Donald Ludwig , 1996. The fishery for pink shrimp (Pandalus jordani) is relatively new in comparison to other trawl fisheries for fin fish off the west coast of Vancouver Island ( W C V I ) shelf. Dur ing the rapid development of this fishery in the early 1970's large catches were landed, and in just a few short years stocks collapsed despite careful monitoring by the Department of Fisheries and Oceans ( D F O ) . The W C V I shrimp trawl fisheries is one of a few fisheries in Br i t i sh Columbia where fishery independent surveys have been in place since the in - 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 in the absence of fishing. Historically, shrimp recruitment to W C V I has been highly variable, and unpredictable. This thesis examines the history of shrimp population dynamics off the W C V I and ex- amines environmental and trophic interactions that lead to variabili ty in recruitment and natural mortality. In this chapter I give a brief overview of the history of shrimp trawl fisheries in 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. Natura l mortali ty rates i n pink shrimp populations tend to be highly variable, and I examine the variabili ty in natural mortali ty rates in chapter 4. In chapter 4, I reconstruct the abundance of two shrimp predators and search for correlations between shrimp abundance, natural mortality, and recruitment wi th predator abundance. The 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 in both predator abundance and variation in pr imary productivity have contributed to variability in per capita rates for pink shrimp. The last chapter summarizes al l of the findings and discusses the relative importance of variabili ty i n primary production wi th respect to shrimp production on the W C V I shelf. 1.1 P i n k S h r i m p F i s h e r i e s Trawl fisheries for smooth pink shrimp i n Br i t i sh Columbia began in Stuart Channel in the Strait of Georgia during the 1950s. Expansion occurred in Barkley Sound during the 1960s. Smal 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 in 1972 to assess the potential of developing an offshore fishery for pink shrimp. Ca tch rates and preliminary assessments for pink shrimp in areas 124 and 125 (see Figure 1.1) looked promising, and in 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. From 1974-1978 the U S A fishing fleet also participated in shrimp fisheries off W C V I , and in 1978 wi th the introduction of extended jurisdict ion, U S A participation in the fishery discontinued. U p unt i l 1977, the Canadian fishing fleet was composed of 12 medium sized otter trawlers, and there were no management restrictions in place. Offshore fishing activities were targeted in two areas, the Tofino grounds (area 124), and Nootka grounds (area 125). B y 1977, catch rates for pink shrimp in area 124 had declined dramatically, much of the fleet had moved to more distant grounds in area 125. In the following two years, catch rates in area 125 had also declined. B y the early 1980's shrimp fisheries on the West Coast of Vancouver Island had essentially shut themselves down as shrimp had apparently vanished from the offshore area. Pr ior to the collapse, the Department of Fisheries and Oceans ( D F O ) had implemented quotas for area 124 in 1978 in response to record catches taken in 1976 and indications of recruitment failure in 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 wi th no restrictions in place; they (the D F O ) intended to observe how effort would respond to open access. After failing to reach the quota in area 124, most fishers shifted to the north where shrimp appeared more abundant than in the previous year. It appeared the newly discovered resource off the West Coast of Vancouver Island had been seriously depleted in area 124 and a wave of fishing effort was moving northward. D F O continued annual shrimp surveys throughout the early 1980's and a few fish- ermen, who dare to leave the protected waters of Barkley sound, ventured out to test offshore waters. There was litt le sign of any significant recovery in shrimp stocks in areas 124 and 125 unt i l 1985. Biomass estimates from D F O surveys continued to decline from 1981 to 1983. No surveys were conducted i n 1984 and 1986. In 1985 a significant increase 4 —i 1 1 —i 1 r 130 128 126 124 122 Longitude Figure 1.1: Statistical area definitions used for smooth pink shrimp trawl fisheries. in abundance was noted for area 124. Commercial fishing activities i n area 124 resumed immediately and small, below average, catches were landed. Landings in area 124 contin- ued to increase over the next two years and it appeared the stocks were rebuilding. Pr ior to the collapse in the late 1970's the majority of the fleet were using otter trawls to land shrimp, and during the rebuilding of the stock in the mid 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 north to explore the grounds in area 125. This exploration was often done later in the season, as catch rates declined in area 124. Despite efforts, no significant concentrations of shrimp were detected in area 125 unt i l 1991. D F O had not surveyed shrimp populations in area 125 in 1991, but in the previous year, assessment reports indicted an increase in abundance over the 1988 assessments. B y 1992, record landings were being recorded for the Nootka grounds, and even larger catches were obtained in 1993, mean while, landings in area 124 5 Area 124 Area 125 c/> a> c o 05 TO CO o o in CM o o o I . . . . I I. "i r 1975 1985 1995 Year OJ CD c O <n u> £Z LZ CO o o in o o m . , 11 I 11 i 1 1 r "i i 1975 1985 1995 Year Figure 1.2: P ink shrimp Pandalus jordani landings i n statistical areas 124 and 125 from 1975 to 2000. (or Torino grounds) were beginning to decline (Figure 1.2). In less than two decades of fishing, it appeared as if history were repeating itself, a boom fishery moving northward followed by collapse. In addit ion to fishing mortality influencing shrimp abundance on the west coast of Vancouver Island, there is also evidence that ocean climates influence the bot tom of the food chain, which in turn affects dynamics of zooplankton and pelagic fish communities (McGowan et al . 1998; Robinson and Ware 1999). In the last two decades there have been increases in sea surface temperature anomalies. Others have defined this change in ocean conditions as a "regime shift" (Beamish et al. 1999), and there is some evidence of changes in production and community reorganization in response to ocean climate regime shifts (Anderson and Pia t t 1999). There have been no directed studies of the influence of climate variation on pink shrimp populat ion on the west coast of Vancouver Island. 1.2 L i f e h i s t o r y o f Pandalus jordani Smooth pink shrimp Pandalus jordani are protandrous hermaphrodites, maturing 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 in the fall (November). Females lay eggs, which remain attached to abdominal appendages unt i l hatching in March-May. Naupl i i or free-swimming larvae spend 2-3 months in the water column and then settle by late June or early July. B y October, as males, they are 60-70 m m in length. B y the following autumn males reach 85 mm. Most adults become females in 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 maximum 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 in the waters off the coast of Vancouver Island: the smooth pink shrimp Pandalus jordani and the spiny pink shrimp Pandalus borealis. The west coast Vancouver Island ( W C V I ) shrimp trawl fishery targets the more abundant smooth pink shrimp and B C is the northern l imit of fishable populations (Butler 1980; Bouti l l ier et al. 1998). Smooth pink shrimp fisheries in B C were first reported in 1952, in Stuart Channel , and have also been operating in Barkley Sound since the 1960s. 1 .3 R e s e a r c h S u r v e y P r o g r a m Annua l research survey programs have been conducted off the West Coast of Vancouver Island since 1973, wi th exception of a few years in the 1980s. Survey programs are intended to provide an unbiased index of abundance, age-structure, and estimates of total mortali ty rates (Boutil l ier et al. 1998). In the past, assessing abundance each year was carried out using spatial interpolation of catch per unit area swept. Thus trawl surveys were designed to map out the distribution of shrimp stocks on each ground. Pr ior to 1987, it was assumed that two main aggregations existed, one in area 124 (Tofino grounds) and one in area 125 (Nootka grounds). Recent survey data suggest that the 7 offshore distr ibution of shrimp may be a continuous distribution that lies between the 100- 150 m depth contour (see Figure 1.4). Th is 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 rawl to sample to grounds. Each tow is roughly 30 minutes i n duration and sweeps roughly 1.0 nautical miles (1852 meters) at a speed of 2.0 knots. The length of the footrope on the otter t rawl is 11 meters; therefore the average area swept per tow is roughly 20,372 m 2 . The highest densities of shrimp are usually found between the 100m and 150m depth contour (Fig. 1.4). For each tow, the catch is sorted to species, weighed, and discarded. Each 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 distr ibution. In most years surveys are conducted in the month of May, 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 unt i l shrimp catches were negligible. Spatial catch per unit area swept data interpolated using a bicubic spline (Boutil l ier et al. 1998) was used to generate an estimate of absolute biomass for each statistical area (Fig. 1.3). Age composition samples collected over time are useful for monitoring changes in total mortal i ty rates, assuming samples reflect the true populat ion age structure. If it is not possible to age the species, then length-frequency data are generally used to infer population age structure. The standard length measurement for shrimp is carapace length, and because shrimp and other crustaceans grow by molt ing their exoskeletons there are no hard body parts that can be used for aging. The 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 0 0 o o o o o o o o l l . l i . I i—r T 1975 1985 1995 Year Area 125 w CD c o -»—' w CO CC E o in o o o 00 o o o •3- 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 in areas 124 k. 125 from 1973 to 2001. Note no surveys were conducted in 1974, 1984, and 1986 in both areas and no surveys were conducted in area 125 in 1989 and 1991. Biomass estimated as mean biomass per area swept by survey times total area surveyed. Figure 1.4: Locat ion of research t rawl 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 th observed data (Schnute and Fournier 1980). Given estimates of proportion-at-age, total mortali ty over time can be inferred from changes in age composition (i.e. e~Zi = N"^'lt+l )• Length frequency data collected from research trawl surveys are used to assess age structure in the population. Past assessments have used the length frequency data to assess relative recruitment and changes in natural mortali ty rates over time by convert- ing the length frequency observations into proportion-at-age using probabili ty transition matrixes. In this thesis I develop a catch-at-length model in which length frequency data are treated as observations for estimating model parameters. I am particularly inter- ested in estimating historical recruitment and mortality, and determining if variability in recruitment and natural mortal i ty rates are a result of variabil i ty in climate and pre- dation. I then investigate the effects of climate variability on shrimp recruitment, and if variabil i ty in natural mortal i ty is a result of changes in predation mortality. Final ly , I ex- amine shrimp dynamics from an ecosystem perspective using a dynamic ecosystem model that includes complex trophic interactions, other commercial fisheries, and variation in primary production. 1.4 Thesis Goals This thesis uses a combination of single species assessment methods and ecosystem mod- eling 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 abun- dance, recruitment, and changes i n natural mortali ty rates over time. 2. Determine if the observed variation in recruitment over time is related to ocean conditions, and if there is evidence for advection of shrimp larvae. 10 3. Determine if variation in natural mortali ty rates are related to predation by gadoid predators. 4. Exp la in how much of the observed variation can be attributed to trophic interac- tions; both changes in predation over time and increases i n production associated wi th favorable climate conditions. Chapter 2 Shrimp Recruitment and Mortality 2 . 1 I n t r o d u c t i o n 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 mortali ty and growth. For short lived, small bodied organisms it is important to determine how each of these v i ta l rates change over time. 2 . 2 S t o c k R e c o n s t r u c t i o n M o d e l The goal of reconstructing shrimp populations on W C V I is to quantify recruitment and recruitment variation, total mortal i ty and mortali ty due to fishing. To reconstruct shrimp populations in different statistical areas, I use a method that does not assume an un- derlying stock-recruitment relationship. In other words, I treat annual recruitment as an arbitrary historical variable to be estimated from the data. The 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 wi th other model parameters. The 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 re- cruitment, mortality, and growth. The statistical objective function, also referred to as a likelihood function, is s imply a function that equates a probabili ty measure for obser- vations (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 param- eter 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 in - tegrating a posterior distr ibution, wi th uniform prior probabilities for the parameters, using Markov Cha in Monte Car lo ( M C M C ) methods (Gelman et al . 1995). 13 2 . 2 . 1 P o p u l a t i o n d y n a m i c s m o d e l I assume part ial recruitment to the fishery occurs as early as age 1, and that al l shrimp die before age 5. A description of al l model parameters and symbols used in the following equations are found in Table 2.1. Numbers of shrimp at age a + 1 and t + 1 are calculated using: N a + h t + 1 = K ^ - ™ (2.1) Recruitment each year is simply the number of age 1 shrimp, which is calculated using: Nht = Re~wt (2.2) where R represents mean recruitment. Annua l recruitment anomalies were assumed to be log-normal and constrained to sum to zero. I also assume that natural mortali ty rates are variable over time and have a mean M = 0.6 (Boutil l ier et al. 1998; Mar te l l et al . 2000) and have a log-normal distribution. Var ia t ion i n natural mortality was implemented by estimating a sequence of deviations vt such that Mt = 0.6eUt. I assume recruitment to 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 La is the mean length-at-age a. In general, annual fishing mortali ty Ft can be calculated in one of two ways: either using fishing effort mult ipl ied by a catchability coefficient (conditioned on effort), or total annual catch divided by total available biomass present in that year (conditioned on catch). The time series of commercial fishing effort for W C V I shrimp trawl fisheries is incomplete, so I chose to use annual catches to calculate annual fishing mortalities: Ft =-log where Ct<VBt (2.4) 14 Table 2.1: Mode l parameter descriptions and symbols used in population dynamics model Descr ip t ion Symbo l Time subscript t Age subscript a Numbers-at-age at time t Nat Natura l mortal i ty rate in year t Mt = 0.6e"« Fishing rate Ft Average annual recruitment R Recruitment anomaly wt Morta l i ty anomaly Vt Vulnerabi l i ty at age Va Length at 50% vulnerabili ty Lh Mean length at age a La M e a n weight at age wa Shape parameter for logistic function 7 A n n u a l landings ct Vulnerable biomass VBt where vulnerable biomass (VBt) i n each year is the sum of products of the numbers of shrimp in each age class, mean vulnerability-at-age, and mean weight-at-age: 4 VBt = 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 in the following section. There are 55 population parameters, includ- ing variable parameters that describe annual recruitment and natural mortali ty rates (R,wi,... ,wn_i,vi,.. .,i/n_i,Lh,j). 2.2.2 O b s e r v a t i o n M o d e l s In the annual research t rawl 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 struc- ture (as measured by length frequency data). The purpose of an observation model is to generate a set of predicted observations based on predictions from the population dy- namics 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 wi th other observation model parameters. Predicting Relative Abundance D F O estimates absolute abundance on the shrimp trawl grounds using spatially interpo- lated 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 in F ig . 1.3 on page 8). I treat these estimates as a relative abundance index. Abundance indices generated by D F O represent the fraction of the population that is vulnerable to fishing gear. Further- more, 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. The age structured population dynamics model used in this assessment expresses abundance in numbers, while abundance indices are expressed in 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 Butler (1980) to convert population numbers to biomass. The following relationship was used for the conversion. W = aLb where a = 0.002257082 and b = 2.609395004 Using the assumptions stated above, the observed relative abundance information is related to the true vulnerable stock biomass via: Yt = qVBteVt, (where vt is normally distributed wi th a mean 0 and an unknown variance). Here q represents the survey vessel catchability coefficient, Yt is the observed index of relative abundance in year t, and VBt is the vulnerable biomass in year t (a predicted quantity). In this case q is an unknown 16 nuisance parameter, and is simply the slope of the zero intercept linear regression between Yt and VBt. A simple trick to eliminate q from the estimation procedure is always use the conditional maximum likelihood estimate of q in the calculation of the likelihood function. This estimate is easily calculated using the "Z statistic": Using this formulation, the maximum 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 in Equat ion 2.6 is adjusted downward to include any commercial catches that may have been taken prior to the survey. To fit model predictions to observed abundance indices I use a likelihood approach and minimize the following negative log-likelihood: where N is the number of years of relative abundance observations. Length Frequency Observations Growth 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 length- frequencies by tracking cohorts over time. Note that symbols used i n the following equations are described in Table 2.2. Mean length-at-age is described by a von Bertalanffy growth curve, and variation in length-at-age is assumed normal. Therefore, the length of an individual at age a plus the (2.6) (2.7) 17 Table 2.2: Parameter descriptions and symbols used in length frequency observation model. Descr ip t ion S y m b o l 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 La 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 sx Number of Observations in year t Nt Probability of length x at age a Px,a Observed Frequencies Qx,t Predicted Frequencies Qx,t time t at which the sample was taken relative to its birthday tQ is: LA+T = L^ (1 - e-fc<»+«-'»>) + 5a (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 A = h + ( L A - h)lZ PA-i (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~k is the growth coefficient. Equation 2.9 is very convenient because 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: _ LA - lxPA-1 0 0 " l - ^ - i k = -log(p) to = ai~\l09(it^) 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: aa = \xexp { A 2 - 1 + 2 1-p ,o - l (2.10) 1- pA~\ where X% and A 2 are parameters to be estimated and the term inside the square brackets represents the length dependency of the standard deviations independently of the length parameters l\ and LA- A constraint was imposed on A 2 > 0 to ensure that standard deviation in length-at-age increases with age. Note that if A 2 = 0 then the standard deviations are length independent (i.e. o~a — Ai for all ages). 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 1 Px,a = 7 = / exp • aaV2n Jx-d (x - Laf dx (2.11) 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 Pxt — Y^aipatSxpXta- Note that this term is normalized to ensure total probabil- ity = 1, and proportions-at-age (ifta) are determined by the age-structured population 19 dynamics model. The predicted number of observed shrimp in each length interval is simply the total number of shrimp observed, mult ipl ied by the probabil i ty of being in length interval x summed over ages (Qx,t = Yla ^tPx,a) resulting in a predicted length frequency distr ibution. The contr ibution of length frequency data to the objective function was calculated using a robust l ikelihood method, which is described in Fournier et a l . (1990). The log-l ikelihood equation for length frequency samples is: N n N LQ = -l/2j2/Zlo^27T(^ + 0.1/n)] -^nlogfa) t=l x=l N n -l l t=l t=l i = l H 2 f e , + 0.1/n) + 0 0 1 (2.12) where Pxt and Pxt are the observed proportions and probabi l i ty of sampling a shrimp in length interval x in year t, respectively. To reduce the influence of large sample sizes , say > 400, the sum of squares term for each length frequency sample in year t is mult ipl ied by r t 2, where r t 2 = mm(Qxt,400). The variance term £ x t for the expected probabil i ty of being in length interval x was estimated as (1 — Pxt)Pxt, therefore, if the probabil i ty Pxt — 0 then the variance term = 0. A small number (0.1/n) is added to the variance term to ensure £ l t does not approach zero, and this also reduces the influence of length observations where the expected probabi l i ty is near zero. One of the more difficult problems in such analysis is determining the number of age-classes present in the length frequency data (see Schnute and Fournier 1980); since Pandalus jordani are relatively short-l ived and grow quickly, I assume a max imum of 4 age-classes present in al l years. 20 Relative Fishing Rates Fishing rates based on area swept estimates were calculated using logbook records from commercial vessels reporting fishing activity in areas 121 to 125. A commercial log book program for shrimp trawl fisheries became a mandatory license requirement in 1987. D a t a for the number of vessels fishing in areas 121-125 are presented in Table 2.3 and estimates of total area swept were calculated using information on individual t rawl durations and vessel speed (distance covered) from the logbook database maintained at the Pacific Biological Station, Nanaimo, B C . Table 2.3: 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. Year Tota l num- A r e a Area A r e a A r e a Number of Ves- ber of ves- 121 123 124 125 sels reporting > sels 100 Shots Per Year 1987 84 2 25 55 2 21 1988 128 0 37 81 10 23 1989 64 0 22 37 4 18 1990 62 0 26 35 1 16 1991 133 2 33 68 23 31 1992 105 1 21 44 32 14 1993 92 3 18 35 33 19 1994 106 0 35 43 27 38 1995 248 30 103 66 39 65 1996 190 13 85 46 34 57 1997 117 7 47 39 16 18 1998 150 8 73 34 23 26 1999 105 1 63 25 8 18 2000 89 6 54 20 8 6 To calculate fishing rate based on area swept estimates, first I define fishing rate as total catch divided by total vulnerable biomass (F = C/VB). Assuming catch is proportional to stock size and fishing effort, then catch can be expressed as C = cEVB where E is fishing effort, and c is the catchability coefficient. Substi tuting this catch equation into the fishing rate equation F is now expressed as F — cEB/VB, i.e. F = cE. 21 Catchabi l i ty 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. Th is is essentially the same definition of catchability as in Paloheimo and Dickie (1964), assuming the probability of capturing a shrimp wi thin the swept area a is equal to 1. Substi tuting this definition of c into the calculation of fishing rate, fishing rate is now expressed as F = aE/A. Using 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 in cell i. To calculate fishing rates over the entire ground we simply sum over a l l grid cells fished and estimate the total area over which the stock is distributed each year: Ft = ^ (2.13) The area swept by one unit of fishing effort is variable due to differences in 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 in the logbook information. Therefore, I assume an average area a swept by a l l vessels. Since 1998, the average foot rope length of al l vessels reporting fishing activity was 65.48 feet (or 10.77*10~ 3 nautical mile) and the average tow distance was 3.385 nautical miles. Each tow roughly sweeps a — 3.65* 1 0 - 2 square nautical miles. The total area swept is approximated by the number of tows (Et) mult ipl ied by the average area swept (numerator in Equat ion 2.13). It is also possible that only a fraction of the total shrimp wi th in 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 in Equat ion 2.13 would be mult ipl ied by a fraction of the total shrimp that escape the harvesting process. Unfortunately there are no data available to calculate what fraction of 22 shrimp escape the trawl; therefore, I assume that the probabili ty of capture has remained constant over time, and treat the F estimates as a relative index of the true fishing rate. The total area over which the stock is distributed (At) is also estimated from com- mercial logbook data using similar methods to those proposed by Winters and Wheeler (1985). To calculate stock area each year and in each statistical area, first I overlayed a 1 square nautical mile grid over the W C V I . The grid allows each tow location to be recorded in a corresponding row and column, also known as rastering. For each year, and in each statistical area, start locations for each tow were plotted on the raster map. The total area was calculated by adding up the area of a l l cells that had non-zero catch rates. For example, if 80 cells were fished, and 75 cells had positive catch rates, then the total area over which the stock was distributed was 75 square nautical miles. This 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, and all areas accessible to fishing have remained open for at least part of each fishing season. The estimated fishing rates based on area swept calculations are assumed to be pro- portional to the true fishing rate, and errors in estimating F from area swept methods are assumed to be lognormal. The 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 an F statistic: F = Ln(Fas/Ft), where Fas is the observed swept area F estimate from logbook data, and Ft is the calculated fishing rate using catch data and Equat ion 2.4. The likelihood of observed fishing rates is then: (this is a log-normal likelihood, integrated over the variance of observed Fas around (2.14) 23 Ft)- 2.2.3 Parameter Estimation and Uncertainty M o d e l parameters were estimated using a Bayesian approach, where the distr ibution of errors for relative abundance data (Eq. 2.7) and area swept estimates of fishing mortality rates (Eq. 2.14) are assumed log normal. B o t h of these likelihood terms are weighted by the maximum likelihood estimate of the variance. The error distr ibution for the length frequency data was assumed normal (Eq. 2.14) and weighted by the relative variance mult ipl ied by a variance scale factor ( r t 2 ) . Large sample sizes were unweighted to a maximum 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 mult inomial distr ibution because Pxt are not independent, normally distributed random variables. Fournier et al . (1990) go on to show that expression 2.12 approximates a minimum x2 estimator, and is a suitable "self-scaling" likelihood for analyzing length frequency data. In addit ion 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 mortal i ty rates. A total of 60 model parameters were estimated for each statistical area by minimiz ing the following objective function: L = (LZ + LQ + LF + P1+P2) (2.15) where P ' s are prior distributions (equations 2.16a-2.16b). Priors were used for recruit- ment anomaly parameters and natural mortal i ty rate parameters and here I assume 24 recruitment and natural mortali ty for shrimp > age-1 are independent. (2.16b) (2.16a) The relative contribution of equations 2.16a and 2.16b to the overall posterior distribution (equation 2.15) is determined by prior assumptions about variation in recruitment (a^j) and natural mortal i ty(a 2 k f ) . Recruitment anomalies were constrained to sum to 0 i n order to estimate mean recruitment. The variance term for recruitment anomaly prior distribution was set <7 2 = 4.0, which is equivalent to a coefficient of variation = 0.54 i n area 124 and 1.66 in area 125. The mean natural mortali ty rate was positively correlated wi th mean recruitment, and in preliminary assessments it was determined that very little information exists about natural mortality rates (i.e. the marginal posterior distribution for M was flat). In order to resolve potential confounding of model parameters, I fixed mean natural mortal i ty rates for each area at 0.6 and estimate annual deviations {ut). The variance specified for deviations in natural mortal i ty was a2M = 0.4, which roughly corresponds to a coefficient of variation equal to 0.2 and 0.35 in areas 124 and 125, respectively. The variance weights used in equations 2.16a and 2.16b reflect prior assumptions about how variable recruitment and natural mortal i ty rates are i n shrimp populations. A s cr2 and a2M approach 0, the model becomes a observation error only model wi th constant recruitment and a fixed natural mortali ty rate equal to 0.6. Conversely, as <T2 and a2M increase in magnitude, more variation is attributed to process error and observation errors become increasingly smaller. Therefore, I also examine how sensitive V""1 71 — 1 Q estimates of mean recruitment and mean natural mortali ty rates ( t ^ _ 1 ) are over a range of o2w and o2M. Mode l parameters were estimated over a 10x10 grid of <7 2 and a2M 25 values ranging from 0.8-8 in cr^ axis and 0.1-1 in o2M axis. I then plot the mean estimates of RQ and M versus o2w and a2M. A total of 60 model parameters ( 9 = RA, li, L A , p, A 1 ( A 2 , 7 , lh, wt, ut) were estimated by minimiz ing equation 2.15. Min imiza t ion was carried out using a phased approach supplied in the A D Mode l Builder software ( © O t t e r Research Ltd . ) , where some of the parameters are held fixed during the ini t ia l phase of the minimizat ion process. Parameter estimation was broken down into two phases; i n the first phase, a l l parameters wi th the exception of 7 and vt were estimated, then in the second phase 7 and ut were included. Dur ing the first phase of minimizat ion addit ional priors were used on mean length at age one ((l\ —12.2)2/0.001, where 12.2 is the mode of the first age class in the length frequency data sets) and mean fishing mortality rate ((/og , (/ t /0.4)) 2 /0.001). These addit ional priors effectively reduce the number of iterations required for the non-linear search routine to converge while assuming that natural mortal i ty rates was constant (M = 0.6) and the overall mean fishing mortali ty 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 in non linear models usually provides more stable estimates, especially in the case of over-parameterized models (Dave Fournier, pers comm., Otter Research L t d . , Sidney B C ) . To examine parameter uncertainty in a Bayes framework, assuming uniform prior distributions for all parameters, I sampled from the joint posterior distr ibution using a Metropolis-Hastings Algor i thm (Gelman et al . 1995). The Metropolis-Hastings Algo- r i thm is a Markov Cha in Monte Carlo ( M C M C ) sampling routine that works well for high-dimensional problems. The algorithm works by taking a random walk through the parameter vector (0 ) , sampling parameter combinations at rates leading to a stationary distr ibution that is the joint posterior distr ibution P(Q\Y). Given a target distribution 26 that is dependent on Y data, the algorithm creates a sequence of random points ©i, 82 . . . whose distr ibution eventually converges to the target distribution. The algorithm starts by randomly choosing a candidate point 9*, and evaluates this point relative to the previous point t — 1, P(9*|9 t_i). If ©* is more probable than 9 t _ ! then it is auto- matically accepted as a point in the posterior distribution. If 9* is less probable than 9 ( _ i then it is accepted wi th a probability of P(Q*\Y)/P(Qt-i\Y). Otherwise the point is rejected, and a new candidate point 9* is selected (note that if 9* is rejected then the previous point 9 t _ i is included in the posterior sample again). Candidate points are selected using a jumping rule Jt(Q*\@t_i). For the jumping distr ibution I use a multivari- ate normal distribution, using the "Cholesky" decomposition of the covariance matr ix generated by the non-linear search routine (Press et al . 1996). After the algorithm con- verges 10,000 random samples from the final posterior distr ibution are used to generate marginal posterior probabili ty distributions for each parameter. The assessment model described here was also tested using fake data generated by the model to ensure parameter estimates were unbiased, and to ensure the model was working correctly. Th is process involved using the population dynamics models, and observation models as a simulation model to generate time series data. Then, these fake data were analyzed wi th the assessment model. Results of the simulation studies are shown in Appendix A . 2.3 Results 2 . 3 . 1 E s t i m a t i n g 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 (Fig . 2.1). Dur ing the early 1980's the fishery had essentially shut itself down because shrimp were vir tual ly absent in area 124 and low in 27 area 125. B y 1985, stocks increased abruptly in area 124. No evidence of a significant increase in area 125 was observed unt i l 1991. In 1991, research surveys were not conducted in area 125, however, commercial landings indicated shrimp were present on Nootka grounds. Research surveys in area 125 resumed in 1992, and a four-fold increase in abundance was observed over the 1990 survey. Area 124 Area 125 n i i 1 1 \— h 1 1 1 1 r 1975 1985 1995 1975 1985 1995 Year 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 trawl surveys. M o d e l predictions of exploitable biomass are in general agreement wi th biomass trends generated by D F O ' s spatial interpolation of area swept information from research trawl 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 in 1983. The largest area swept biomass estimate in area 124 was 11505 tonnes observed i n 1976 (Fig. 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 Nino, shrimp stocks in both areas had collapsed, and abundance indices were less than model predictions for area 124 i n 1983 and area 125 in 1982. 28 Area 124 Area 125 a> C a cc 1 1 1 1 1 1— 1975 1980 1985 1990 1995 2000 Year 1 1 1 1 1 1— 1975 1980 1985 1990 1995 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. The largest deviation in biomass estimates was in area 125 in 1998 (Fig . 2.2), however, I suspect biomass was grossly under-estimated in 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. Th i s increase in biomass could not be accounted for by growth of biomass over the summer months. Catches exceeded estimated biomass in 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. The estimated population age structure over time in both offshore populations is shown in F ig . 2.3. Since 1998, total biomass in both areas has been increasing, especially recently in area 125. Al though population numbers have been declining in area 124 since 1999, biomass continues to increase due to reductions in total mortal i ty that have allowed for increased growth of the remaining biomass. 29 Area 124 Area 125 1975 1985 1995 1975 1985 1995 Year Year Figure 2.3: Popula t ion age structure for pink shrimp in areas 124 and area 125. Whi t e shaded polygon are age-1 shrimp and black shaded polygons are age-4 shrimp. F i t t ing the model to length frequency data over time largely influences model predic- tions 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 in area 124 are in general agreement wi th the observed data w i th exception of f998, 2000, and 2001. In 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 in model predictions. Predicted length frequencies in 2001 suggest that growth may be greater than historical averages (this is also apparent in 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 in both observed data and model predictions. In area 125, the estimated coefficient of variation in length-at-age was greater than that observed in area 124 (Table 2.4), resulting in wider distributions for each cohort. Aga in , observed and predicted length frequencies in 1998 contradict each other and contradictions are probably attr ibuted to the very small sample size taken in 1998 (N — 131 shrimp). Also , 30 in 2001 shrimp appear to be growing at rates higher than the historical average. 1975 1985 1994 10 15 20 25 30 10 15 20 25 30 10 15 20 25 30 Length (mm) Figure 2.4: Predicted (solid line) and observed (histograms) length frequencies for area 124. Est imated mean length-at-age is represented by vertical dotted lines. Note there are no length frequency data for 1983-84, and 1986, and predicted distributions are shown as a straight line. 31 Figure 2.5: Predicted (solid line) and observed (bars) length frequencies for area 125. Est imated mean length-at-age is represented by vertical dotted lines. No length frequency samples were taken i n 1984, 1986, 1989, and 1991 in area 125. 32 Table 2.4: Est imated growth parameters for area 124 and 125 pink shrimp populations. A r e a Loo (mm) k to °\ o\ °\ 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 Est imated growth parameters are shown in Table 2.4. Varia t ion in length-at-age was slightly larger in area 125. Shrimp in area 125 appear to grow at a faster rate, however attain a smaller mean asymptotic length (Loo). 2 . 3 . 2 P a r a m e t e r E s t i m a t e s a n d U n c e r t a i n t y A total of 60 model parameters were estimated. M a x i m u m likelihood estimates and standard deviations for al l model parameters in both statistical areas are shown in Tables 2.5-2.7. Parameters listed in Tables 2.6 and 2.7, are nuisance parameters and when included in the estimation procedure explain additional variation i n the data. Omi t t ing these nuisance parameters from the estimation routine is equivalent to fitting the model assuming only observation errors, constant recruitment, and constant natural mortality. In general standard deviations for estimated model parameters are smaller for area 124 than 125, largely due to larger sample sizes and greater contrast i n the data between these two areas. Table 2.5: Est imated model parameters and standard deviations for statistical areas 124 and 125. A r e a 124 A r e a 125 Parameter M L E Estimate Std M L E Estimate Std log(Ro) 6.28 0.16 6.13 0.32 h 10.60 0.21 9.04 0.35 IA 21.18 0.13 20.27 0.11 P 0.77 0.03 0.59 0.03 Ai 1.30 0.05 1.36 0.05 A 2 0.00 0.02 0.00 0.00 7 0.54 0.03 0.63 0.08 k 18.10 0.34 14.81 0.60 Mean recruitment (R0) is similar between the two areas, yet al l other parameters listed in table 2.5 differ markedly between the two areas. In general, shrimp in area 125 appear to grow faster, at tain a smaller asymptotic length and have more variation in length-at- age than shrimp in area 124. It also appears that larger shrimp are caught in area 124 versus area 125, based on estimated selectivity curve parameters Ih and 7. Uncertainty Table 2.6: Estimated recruitment anomalies (wt) and standard deviations for statistical areas 124 and 125. A r e a 124 A r e a 125 Parameter Name M L E Estimate Std M L E Estimate Std 0.91 0.23 0.39 0.38 ^1976 -0.27 0.40 -0.09 0.37 w1977 -0.37 0.41 -0.98 0.38 ™1978 0.03 0.39 -0.11 0.40 ™1979 -0.15 0.39 0.02 0.40 ^1980 -0.41 0.44 -0.15 0.46 ™1981 -0.40 0.42 -0.70 0.46 ^1982 -0.07 0.40 -0.54 0.42 ™1983 0.18 0.33 -0.45 0.46 ^1984 0.88 0.26 0.29 0.37 ^1985 0.66 0.25 -0.49 0.42 ^1986 0.73 0.23 -0.17 0.33 ^1987 0.75 0.24 -0.30 0.34 ™1988 0.68 0.24 0.22 0.32 ^1989 0.21 0.26 0.44 0.31 ^1990 -0.04 0.28 0.24 0.35 W1991 0.21 0.27 0.44 0.32 ^1992 0.39 0.25 0.19 0.32 ™1993 0.26 0.26 0.46 0.31 ^1994 -0.68 0.31 -0.47 0.38 ^1995 -0.38 0.27 -0.66 0.50 ^1996 -1.12 0.33 -0.89 0.83 W 1 9 9 7 -0.90 0.34 -1.09 0.32 ^1998 0.10 0.34 0.59 0.37 ^1999 0.21 0.37 1.16 0.39 ^2000 -1.39 0.71 2.62 0.44 and apparent correlations between parameters listed in table 2.5 are illustrated in Figures 2.6 and 2.7. These figures are sub-samples from the posterior distr ibution that was constructed using the Markov C h a i n Monte Carlo algorithm. In both statistical areas, strong positive correlations exist between growth parameters l\,p, and the length at 50% vulnerability (Ih)- The 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: Est imated deviations i n mortali ty (vt) and standard deviations for statistical areas 124 and 125 A r e a 124 A r e a 125 Parameter Name M L E Estimate Std M L E Estimate Std ^1975 -0.06 0.30 0.15 0.24 ^1976 -0.20 0.40 -0.01 0.45 ^1977 0.19 0.43 0.23 0.45 ^1978 0.28 0.44 -0.01 0.46 ^1979 0.35 0.44 0.03 0.46 ^1980 0.24 0.43 0.12 0.48 ^1981 0.26 0.44 0.72 0.56 ^1982 0.22 0.44 0.00 0.48 ^1983 -0.24 0.38 0.08 0.49 ^1984 -0.24 0.37 0.12 0.49 ^1985 -0.28 0.36 0.01 0.47 ^1986 -0.28 0.36 0.06 0.47 ^1987 -0.26 0.35 -0.10 0.44 ^1988 -0.16 0.35 -0.25 0.41 ^1989 -0.18 0.35 -0.31 0.40 ^1990 -0.19 0.34 -0.43 0.39 ^1991 0.02 0.35 -0.42 0.40 ^1992 0.09 0.35 -0.22 0.44 ^1993 -0.06 0.33 -0.30 0.45 ^1994 0.10 0.31 -0.29 0.46 ^1995 0.24 0.32 0.10 0.50 ^1996 0.11 0.40 0.71 0.86 ^1997 0.19 0.41 0.52 0.87 ^1998 0.08 0.39 -0.17 0.45 ^1999 -0.13 0.37 -0.16 0.43 ^2000 -0.10 0.39 -0.19 0.43 parameters 7 and lu- A s 7 increases in 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% vulnerabili ty increases 7 decreases and a higher fraction of smaller shrimp are vulnerable to the fishing gear. The correlation was much stronger in area 125 than i n area 124 (Fig . 2.6 and 2.7). Parameters for standard deviation in length-at-age were also negatively correlated Figure 2.6: Pair plot of 500 random sample estimates (RQ, Li, L A , p,X\,^2,7,h) from an M C M C sample of length 250,000 for area 124. 8.0 9.0 0.50 0.60 0.000 0.020 13.5 15.5 I .I I I . , . J _ J I l _ . , , I I. I I I 1 , , I 1 1 I I I 6.0 7.5 19.8 20.4 1.25 1.45 0.5 0.8 Figure 2.7: Pair plot of 500 random sample estimates (RQ, L\, LA, p, A 1 ; A 2 , 7 , lh) from an M C M C sample of length 250,000 for area 125. 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 2 have several modes (Fig. 2.8). Recall that A 2 was constrained to be greater than 0, such that standard deviation in length-at-age increases with age, and if A 2 = 0 then standard deviation in length-at-age is the same for all ages. Area 124 Area 125 r ~ — i 1 r 0.000 0.010 0.020 ^2 Figure 2.8: Marginal posterior distributions for A 2 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. Sensitivity to a2M and 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 (wt and vt, process error terms). To explore sensitivity of M and R0 to prior assumptions about afj and er ,̂ model parameters were estimated over a grid of a2M and values. The results of the sensitivity analysis are shown in Figures 2.9 and 2.10. The mean value of M and R0 in each column or row in the matrix of possible variance combinations were plotted against the vectors a2M and a2w. For example, the line with points in upper left panel of Figure 2.9 represents the mean estimate of M over all values of G2W , and the dotted or dashed lines represents mean estimates of M for a single 38 value of a2. Area 124 Area 124 s 0.2 0.4 0.6 0.8 1.0 1 2 3 4 5 6 7 8 Figure 2.9: Sensitivity of M to cr 2^ and c r 2 . Left column represent the mean estimate of natural mortali ty as a2M increases from 0.1 to 1, where the mean (line wi th points) is calculated over al l a2M. Right column represents the mean estimate of M as increases from 0.8 to 8, where the mean is calculated over al l c r 2 . In the left columns, dotted lines represent estimates of M when cr 2 = 0.8 and dashed line represents estimates of M when cr^ — 8, and in the right column dotted lines for a2M = 0.1 and dashed line for a2M = 1 A s u2M increases in value, estimates of mean natural mortality rates tend to increase for both statistical areas (left column of F i g . 2.9). For example, in area 124 the mean natural mortali ty rate estimate when a2M = 0.1 is 0.6, and as <r̂ increases to 1, estimates of natural mortali ty rates increase to 0.86. A s estimates of cr2 increase, mean natural mortali ty estimates tend to remain relatively constant, however, as o2M increases mean estimates of M increase (as shown by dashed lines in right column of F ig . 2.9). Estimates of mean recruitment RQ were less sensitive to changes in variance weights for the two prior distributions. Increasing (variance in recruitment anomalies) results i n 39 Area 124 Area 124 0.2 0.4 0.6 0.8 CC o —I— 0.2 Area 125 i— 0.4 1.0 Area 125 Figure 2.10: Sensitivity of RQ to a2M and cr 2 . Left column represent the mean RQ as a2M increases from 0.1 to 1, where the mean (line wi th points) is calculated over al l a2M. Right column represents the mean estimate of Ra as a 2 increases from 0.8 to 8, where the mean is calculated over a l l a^. Dotted lines in left column represent estimates of RQ when cr^ = 0.8 and dashed line represents estimates of R0 when cr2 = 8, i n the right column dotted lines represent a2M = 0.1 and dashed line a2M = 1 . rather uniform estimates of RD in area 124 (Figure 2.10). In area 125, RQ declines slightly as variation in recruitment anomalies increase. Estimates of Ra were fairly consistent over the range of variances in natural mortali ty rate anomalies. Varia t ion in estimates of RQ and M over ranges of a 2 and a2M shown in Figure 2.9 and 2.10 was a result of a likelihood surface that has several local maxima across the range of variances explored. For some combinations or a2M and a^, the estimation routine converges to very different solutions for parameters. In other words, the shape of the likelihood surface is not smooth and changes wi th different values of a2M and <r̂ (Figure 2.11). 40 Area 124 Area 125 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 « 2 « 2 Figure 2.11: Contour plots for objective function values (e-L/max(L)) versus a2M and a2w values. Area 125 clearly shows 3 local maxima. 2 . 3 . 3 M o r t a l i t y 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 of fishing mortality 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 true fishing mortality rates and were treated as observations (equation 2.14) in estimating model parameters. For the results on mortality, I first report on the result of the log book analysis, where area swept calculations were used to construct an independent estimate of fishing mortality. Then I report estimated fishing mortality 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 commercial fishing effort is shown in Figure 2.14. Note that all 41 Area 124 Area 125 1988 1994 2000 1988 1994 2000 Year Year Figure 2.12: Est imated area over which shrimp populations in statistical areas 124-125 were distributed. A r e a estimates were calculated using commercial logbook data from 1987 to 2000; these data are shown in Figure 2.14 log-book records for the W C V I island are plotted on Figure 2.14, this also includes areas 121 and 123. Records in areas 121 and 123 were not used for estimating stock areas or fishing mortality rates because shrimp populations in these areas were not examined. O n average, shrimp are distributed over a much larger area in statistical area 124 than in statistical area 125. Stock area varies considerably from 1987 to 2000 (see Figure 2.12). Dur ing periods of low abundance stock area is much reduced and during periods of high abundance stock area increases substantially, and as abundance declines so too does stock area. For example in area 124, relative abundance (note, estimated abundance is shown on Figure 2.13) declined from 1,556 tonnes in 1994 to 371 tonnes i n 1999, and estimated stock area declined from 348 square nautical miles to 30 square nautical miles (Figure The tow start locations shown i n Figure 2.14 represent al l of the tows logged since 1987 in al l statistical areas (121-125). Pr ior to 1996, less than 50% of al l tows were spatially referenced, and therefore, I am assuming that reported tow locations accurately represent total area fished. The two principal fishing areas between 1987-2000 are i n statistical areas 123 and 124 (Fig. 2.14, and Table 2.3). Est imated fishing rates based on 2.13). 42 Area 124 Area 125 7) 2000 4000 6000 8000 0 1000 2000 3000 4000 5000 Abundance (tonnes) Abundance (tonnes) Figure 2.13: Relationship between shrimp abundance and stock area from 1987 to 2000. Shrimp abundance for each statistical area is represented by the total abundance (Bt = 2~2a Na,tWa). Note that abundance shown in figure 2.1 represents the vulnerable biomass. area swept calculations for statistical areas 124-125 are shown in Figure 2.16 and were calculated using the total effort data presented in Table 2.5. In both areas, fishing rates tend to increase wi th stock area, as more fishing effort is attracted to areas where shrimp are more abundant. Est imated fishing rates based on log book data in areas 124 and 125 range from 0 to 1.06. Fishing effort in 1998 declined dramatically in area 124, and increased slightly in area 125. Dur ing the late 1990s, eulachon (Thaleichthys pacificus) bycatch in shrimp trawl fisheries was becoming a concern, and 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 . Al though industry had taken a preventative measure to reduce bycatch, commercial fisheries in area 123 and 121 were suspended early in the 2000-fishing season for eulachon conservation concerns, otter trawlers were also restricted in areas 124-125, and Queen Charlotte Sound was closed to al l shrimp fisheries. It should be noted that fishing rates shown in Figure 2.16 were estimated assuming that no commercial quantities of shrimp were left undetected, and that the total area over 4:5 r v < ^ ^ ^ ^ ^ 1 9 9 7 ^ - V̂ S. ' "-: . . T- r r " « ^ ^ ^ ^ ^ 1 9 8 9 ^ - V̂ S. r r r r ^WT ^J < < f r Area 125 1. Area 124 ^^Qfl^gS Area 123 *Ĉ^ r Figure 2.14: Spat ial distribution of commercial fishing effort from 1987 to 2000 on the West Coast of Vancouver Island. Tota l number of tows from 1987 to 2000 was 59,962. 44 Table 2.8: Reported fishing effort in commercial logbooks for statistical areas 124 and 125 Year Number of Shots Tota l Effort Average Tow Time S T D Tow T ime A r e a 124 (Hours) (Hours) 1987 5198 8702 1.67 0.96 1988 6011 9013 1.5 0.76 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 1.7 0.59 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 1.5 N / A 1991 1249 1491 1.19 0.47 1992 1806 2244 1.24 0.49 1993 4852 6336 1.31 0.39 1994 2833 3536 1.25 0.37 1995 3834 4991 1.3 0.43 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 1.14 0.44 45 which the stock was distributed is well represented by non-zero commercial catch rates. This last assumption is part icularly bothersome, in that less than 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 in the mid to late 1990's were located on Vancouver Island itself (representing errors in log book records, see F i g . 2.14), and tows were eliminated from area swept calculations. Final ly, some areas simply cannot be trawled due to rough bottom, depth limitations, or strong currents, and if these areas are large and harbor shrimp then the total area over which the stock was distributed is under-estimated. Fishing Mortality Ftom Stock Assessment The instantaneous fishing mortal i ty rates for areas 124-5 generated from the stock re- construction model are shown in Figure 2.16 on page 47. Fishing mortali ty rates were calculated using observed catch and vulnerable biomass (i.e. F — —log(l — C/VB)), not total biomass. Est imated selectivity curves for the two areas are shown i n Figure 2.15. Note that the estimated length at 50% vulnerabili ty is much larger in area 124 (If, = 18.1) versus area 125 (Ik — 14.8). A s a consequence, estimated fishing mortali ty rates in area 124 tend to be much higher in comparison to area 124. The highest estimated fishing rate occurred in area 124 in 1997 (F = 5). Shrimp abundance in 1997 was estimated to be at an al l time low i n area 124, wi th a high proportion of part ial ly vulnerable age-1 shrimp in the population. From 1976 to 1980, fishing mortal i ty rates were declining in area 124, corresponding wi th the decline in abundance in area 124. In 1978 much of the fishing effort had moved north into area 125, generating high fishing mortalities. B y 1979, catches were also 46 Figure 2.15: Est imated selectivity curves for areas 124 and 125, where dotted vertical line represents the length at 50% vulnerability. declining in area 125, as off-shore stocks had declined to a point where it was no longer economical to fish (J. Bouti l l ier , Pacific Biological Station, Nanaimo, B C , Pers Comm). Throughout most of the 1980s, the remaining shrimp industry focused much of its efforts on inshore waters (including Strait of Georgia and Barkley Sound). B y 1985 offshore shrimp populations had begun to rebuild (Fig . 2.1) and unrestricted offshore fishing activity had resumed in area 124. From 1987 to 1990, vir tual ly al l of the offshore fishing effort occurred i n area 124 (Table 2.8), and fishing mortalities exceeded the recommended target exploitation rate of 33%. Dur ing this same time period, biomass in area 125 was beginning to rebuild; however, significant fishing effort d id not resume in this area unt i l 1991. In 1991, fishing mortali ty rates increased to 0.82 and 0.31 in areas 124 and 125, respectively. F rom 1991 to 1994 the majority of fishing effort in offshore areas alternated between area 124 and 125, as well as the fishing mortali ty rates. Est imated fishing mortali ty rates in area 124 and 125 during 1993-97 had averaged near 0.95, again much higher than 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 Area 124 •»—* o CD C 03 LO "tf - CO - CM 1975 T T 1985 1995 Year Area 125 47 A t co CO CD 05 o CD i r 1975 1985 — i r 1995 Year Figure 2.16: Est imated fishing mortali ty rates from stock assessment model (lines in bot- tom panels) for areas 124-125. Circles represent area swept estimates of fishing mortali ty rates (F = 9j^) from commercial log book data, and these data were used in the like- l ihood criterion for estimating model parameters. Top panels are the residuals between area swept estimates of F and estimated fishing mortali ty rates i n the stock assessment model. abundance, and declining prices for shrimp. Also , in 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. The trends in area swept fishing mortali ty correspond well wi th trends i n fishing mortali ty rates estimated from the stock assessment model (see F i g . 2.16). Including area swept estimates of F in the likelihood calculations did influence the overall trends in predicted F ' s . No 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 in natural mortal i ty rates increases substantially. 48 There are a few major differences i n the trends of F from the two calculations, espe- cially mid-late 1990's. Typical ly , area swept estimates of fishing mortali ty are less than the predicted F from the stock assessment model during the 1993-97 time period. This time period corresponds wi th fairly intensive fishing activity on the offshore populations. It is likely that the total area At was over-estimated for years when fishing effort was high (i.e. fishermen were searching for new grounds as catch rates declined in the principal grounds). Over-estimating the total stock area resulted in reduced estimates of fishing mortali ty rates because the denominator increased in equation 2.13. Nonetheless, the commercial logbook program does provide some auxiliary information about trends in exploitation rates over time. Relative abundance indices from research cruises indicate stocks in areas 124 and 125 had reached a low inT998 . The stock assessment model and length frequency data indicate that this collapse was i n part due to high exploitation rates in these two areas, just prior to 1998 in area 124 and during 1998 in area 125. Increased natural mortali ty rates also played a role in this collapse (see next section). Whereas, area swept estimates of F indicate fishing rates were relatively moderate i n 1998 in comparison to previous years. In fact, F's were estimated to be declining in both areas since 1995 (Figure 2.16) wi th the exception of area 125 in 1998. Commercial landings in area 125 in 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 in area 125 in 1998 under-estimated total biomass, and using this index of relative abundance to tune the stock assessment model resulted in a high estimate of mortality. Fishing rate estimates for 1996 and 1997 in area 124 were substantially higher from the stock assessment model than 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 in annual natural mortality (yt) rates were treated as variable, and were esti- mated using a l l of the observed information from research surveys and relative estimates for fishing mortali ty based on area swept calculations. Estimates of natural mortali ty rates and uncertainty are shown in Figure 2.14. The estimated error distributions around the mean estimate of M were calculated using the Metropolis-Hastings algorithm, ran- domly sampling from the posterior distr ibution 50,000 times. In both statistical areas, natural mortali ty rates have been highly variable over time, but more so in area 125 (Table 2.9). In area 125 the highest estimate of M was 1.23 in 1981, and in area 124 the highest estimate of M was 0.85 in 1979. Qualitatively, trends in natural mortal i ty rates over time in the two areas seem to be positively correlated, however, large changes in M appear to occur a year earlier i n area 124 than in area 125 (Fig. 2.17). Th i s pattern was also observed just prior to the 1997-98 E l Nino, where increases in natural mortali ty occurred first in area 124. Table 2.9: Natura l mortality rate statistics for areas 124 and 125 A r e a 124 A r e a 125 Average 061 0 6 3 Coefficient of Var ia t ion 0.21 0.35 M i n i m u m 0.45 0.39 M a x i m u m 0.85 1.23 Uncertainty in estimates of M is higher in area 125 than in area 124, largely due to sample size differences between the two areas. Between 1983 and 1990, natural mortali ty rates in area 124 were relatively stable, whereas in area 125, natural mortali ty rates were declining from 1984 to 1990. Research surveys were not conducted in either of these areas in 1984 and 1986; therefore, no information about natural mortali ty rates from length frequency data is available during this time period. Furthermore, relative estimates of fishing mortali ty rates, the other component of total mortality, are not available unt i l 50 Area 124 in co o in •= ° CO OJ d o d in CO o co' in ~ o CO cj in d o d , i ' "i <» K i w # y M V *M T $ V • ' i i i i i i i i i i i i i i i i i i i i i i i i i r 1975 1978 1981 1984 1987 1990 1993 1996 1999 Area 125 J- J- a. J. X -I. -j. -J. j . _i_ j_. "J. J- j-' - 1 " -J- X- J. i i i i i i i i i i i i i i i i i i i i i i i i i r 1975 1978 1981 1984 1987 1990 1993 1996 1999 Year Figure 2.17: M a x i m u m likelihood estimates of natural mortali ty rates (line wi th circles) from stock assessment model for areas 124 (top) and 125 (bottom) from 1975 to 2000. Grey points are a 1000 randomly chosen samples from the posterior distr ibution generated by the Metropolis Hastings Algor i thm, and the box plots represent the median, inner quartile, and sampling range. 51 1987. Fishing mortali ty rates for area 125 in 1987 were low, and other auxiliary data suggest that populations were increasing wi th relatively low total mortali ty rates. 2 . 3 . 4 A n n u a l R e c r u i t m e n t a n d S t o c k - R e c r u i t m e n t R e l a t i o n s h i p s The maximum likelihood estimate for age-1 recruits and associated distributions from the posterior are shown i n Figure 2.18. A s the spawning biomass declined from 1975 to 1977 in area 124, recruitment declined almost proportionally and continued to decline for an additional 4 years (Fig . 2.18). It appears that both poor recruitment and increases in natural mortali ty led to declines in shrimp abundance i n both areas during the early 1980's and recently in 1996 (Fig . 2.17 and F i g . 2.18). Recently, natural mortali ty rates in both areas have been below average and a combination of growth and recruitment has led to increases in biomass in areas 124 and 125. A r e a 125 has apparently had a large increase in annual recruitment in 2000 and 2001 and biomass in area 125 has increased more than 10-fold since 1998. Stock-recruitment relationships for both statistical areas were constructed using es- t imated age-1 recruits, and spawning stock biomass. Recal l that Pandalus jordani are protandrous hermaphrodites, and the age at transition to female begins around age-2 wi th 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 us- ing length frequency and sexing information from research surveys (see Figure 2.19). To estimate the proportion of female shrimp in the population, numbers-at-length were mult ipl ied by the proportions-at-length j that are female (i.e. Nfemaie — Nj * Pj). The logistic curve estimated in Figure 2.19 roughly approximates the length-at- transition from male to female. The 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). Grey 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 mil l ion shrimp, respectively. 53 P Q_ E _0) a> o o o o d P = 1+exp(-y(l-lh)) 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 follow- ing March. During this period female shrimp are exposed to both fishing and natural mortalities, and the spawning biomass must be corrected downward to represent indi- viduals that survive to spawn. The spawning biomass at time t was calculated using SBt = ( y\ Nt-iSjPjWA e-Fi-i-Mt-i-OAMt^ 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, Nt~i are the number of 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 1 r 200 400 600 800 1000 Spawners Area 125 T 1 1 1 i r 0 200 400 600 800 1000 Spawners 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 density- dependent juvenile survival rate. 55 average juvenile survival rate (see bottom-right panel in F ig . 2.20). Recently recruits per spawner have returned back to the historical average in area 124. Similar improvements in juvenile survival rates have also been observed recently in area 125 (see bottom-right in F i g . 2.20). From 1999 to 2001, survival of age-1 recruits has increased nearly 1-3 times over the historical average in area 125. In both statistical areas, the Beverton-Holt model suggests very high compensation in juvenile survival rates when spawner abundance is low. However, the Ricker model produces more conservative estimates of maximum recruits per spawner i n area 124. Omi t t i ng the extreme outlier in 1999 in area 124 reduced the estimate of maximum recruits per spawner from 87 to 32 in the Beverton-Holt model, whereas omit t ing the outlier in area 125 in 2000 increases the estimate from 9 to 20. For the Ricker models, maximum recruits per spawner declined from 9.48 and 13.8 to 6.8 and 6.2 in areas 124 and 125, respectively when omit t ing extreme outliers. 2.4 Discussion In stock assessment modeling, it is most often assumed that natural mortal i ty rates are constant over time (Walters 2000). Th is simplifying assumption is usually made to reduce the number of estimated parameters. Natura l mortality is confounded wi th other model parameters making it difficult to estimate from time series data, and is usually estimated using independent methods. In this assessment, estimating annual natural mortali ty rates was made possible by the inclusion of length frequency data and relative indices of fishing mortali ty based on area swept calculations. Changes in length frequencies over time provides information about the total mortali ty rate (both fishing and predation). If an estimate of fishing mortali ty is available then the other component of total mortali ty can be estimated from length frequency data, or proportions-at-age data (Sinclair 1998). 56 The use of length-based methods greatly improves our abili ty to infer proportions- at-age i n populations that are difficult to age (Schnute and Fournier 1980; Quinn 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. Th is is necessary to track modes in the length frequency data (Smith et al. 1998). There is no evidence that pink shrimp populations spawn more than once a year on the west Coast of Vancouver Island; however, it appears that variabil i ty in the time of year at which they spawn can also bias length-based approaches. The problem shows up i n assessment methods that directly estimate variance in length-at-age rather than estimating a single coefficient of variation for al l age classes. Large variance in the time of year at which eggs are spawned, or young of the year emerge,,is apparent when variance in length-at-age is larger for younger cohorts than it is for older cohorts (if age-specific variances are estimated directly). In effect the wide range in bi r th dates artificially inflate variance in length-at- age calculations. In this assessment procedure, the function used to estimate variance in length-at-age was constrained such that variance was either independent of age (A 2 = 0), or increases wi th age ( A 2 > 0). Estimates of A 2 were near zero for both statistical areas and length frequency data were uninformative about variation in length-at-age, probably due to large over-laps in length-at-age and variability in hatching times between cohorts. Similar findings for Northern shrimp (Pandalus borealis) were documented by Fournier et al. (1991). Est imat ing fishing mortali ty rates based on area swept calculations are subject to errors in estimating total area swept by commercial fisheries, and estimation of the total area over which the stock is distributed. For W C V I shrimp t rawl fishers, estimates of total area swept were carried out using an average area swept mult ipl ied by the total number of tows. This was done because necessary information about tow distance and net wid th 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). The area swept estimates also assumes that 100% of shrimp wi th in 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; Winters and Wheeler 1985; Paloheimo and Dickie 1964). Adjust ing the numerator in Equat ion 2.13 to accommodate changes i n fishing technology over time is t r iv ia l compared to dealing wi th errors in estimating the total area over which the stock was distributed (Swain and Sinclair 1994). The total area over which the stock was distributed was estimated as the total area in which commercial fleets detected shrimp. This was done by imposing an arbitrary grid of 1 N M 2 over a map of reported tow locations and simply summing the number of cells wi th non-zero catch rates. A m o n g all of the problems wi th 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 wi l l eventually reach an asymptote. Using today's high tech G P S information, fishing grounds wi th in 1 N M 2 may st i l l contain a high fraction of untrawlable ground. Close examination of spatial data for ground fish trawl fisheries has revealed that fishermen actively target productive grounds on spatial scales much smaller than 1 N M 2 (see F i g . 3 in Walters and Bonfi l 1999). The total area calculations assume that fishing effort is randomly distributed wi th in each block, and violations i n this assumption lead to over-estimating total area. There is another problem wi th 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 competit ion on the fishing grounds, and creates strong incentives for fishermen to explore new areas. It is generally thought that the effort distr ibution follows an ideal free distr ibution (Gil l is et al . 1993; Walters and Bonfil 1999). In years where fishing effort is high and shrimp prices are high, more and more fishermen are wil l ing to fish in areas wi th lower catch rates. They may even go explore new grounds in hopes of discovering potentially unexploited areas. W h e n using spatial distr ibution of fishing effort to estimate total area over which the stock is distributed, high effort years wi l l probably lead to an over-estimate of stock area. Conversely, low effort years wi l l tend to under-estimate stock area, because remaining fishermen w i l l tend to concentrate their fishing effort in core areas where history has shown the area to be productive. Also , remaining fishermen also tend to have higher catch rates (Hilborn 1985). I t ry to correct for this problem by only including records w i th 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 kg of shrimp in a 3-hour tow (hardly worth staying in the area). Individual fishermen wi l l choose their own criterion for deciding to stay in an area, or move to a new area. This criterion is probably economically driven; therefore, variability in annual shrimp prices, and variabil i ty in the cost of fishing, w i l l influence decisions to stay or move, or simply quit fishing altogether. Including independent estimates of fishing mortali ty rates in the estimation procedure greatly reduces the uncertainty associated wi th estimating fishing mortali ty rates (even if the F estimates are used as relative estimates, see Appendix A ) . Changes in length composition over time provide relatively good estimates of total mortal i ty rates (Z). Est imat ing the F component of Z was based on total catches and the l ikelihood 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 in assumptions about estimating fishing mortali ty rates from area swept estimates are independent, and therefore not biased if used as a relative index. In Figure 2.16, the trends in fishing mortalities are somewhat similar i n both 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 mortal i ty rates, the same estimation procedure should be carried out wi th litt le or no weight on the area swept data. I d id not carry out this analysis in the thesis due to l imitations in computing time, (it takes 18 hours to estimate uncertainty i n model parameters on a Pent ium I V computer). In the simulation study (Appendix A ) , including estimates of fishing mortali ty rates did not introduce any bias, provided errors in area swept estimates of F were not auto correlated. Is there evidence of a stock-recruitment relationship in these two areas? To answer this question I fit Beverton-Holt and Ricker models to the data for each statistical area using a least squares fitting criterion. In both areas the best Beverton-Holt fit to the data shown on the left panels of F i g . 2.20 is a straight line, even when omit t ing extreme outliers when spawning stock is small. The simplest interpretation of this result would be that recruitment is totally independent of stock size wi th in each statistical area. For example, there is lit t le evidence that recruitment in statistical area 125 is a result of spawner abundance in area 125. A n alternative explanation is that observed recruitment in each of the statistical areas is a function of spawner abundance elsewhere, and larvae or juveniles dispersed to new areas (a question I examine in chapter 3). The Ricker model in area 124 provided a more conservative estimate of maximum 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 in LAMENT 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 ) . The physical oceanography on the W C V I is strongly connected to seasonal primary productivi ty regimes (Thomson et al. 1989), which would affect recruitment moderately through trophic interactions. The prominent oceanographic feature off the west coast of Vancouver Island is a persistent coastal current that flows northward. This coastal current is confined landward of the 100 m depth contour (Thomson et al. 1989; LeBlonde et al . 1986). Al though adult shrimp in the offshore areas are found between the 100m and 200m depth contours, the persistent coastal current s t i l l 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 in Chapter 2 to be correlated wi th spawning population to the south. Or if upwelling is responsible for advecting recruits offshore, resulting in recruitment failure, then I expect the recruitment indices to be negatively correlated wi th upwelling indices. If, however, recruitment anomalies are positively correlated wi th upwelling events then I would suspect that increases in primary productivi ty led to increased recruitment. Effects of variation in primary productivity wi l l be examined further in Chapter 5, i n the context of an ecosystem model that explicit ly represents trophic interactions in the region. 3 . 1 . 1 S u m m e r a n d W i n t e r C u r r e n t C i r c u l a t i o n P a t t e r n s The persistent coastal current is pr imari ly driven by surface outflow from Juan de Fuca Strait during 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 Juan de Fuca Strait. The coastal current off the W C V I during summer months (April-September) tends to spread offshore along the broad continental shelf Nor th of the Juan de Fuca Strait and south of Barkley Sound. A s the current travels north it forms a narrower band north of Barkley Sound. Northwesterly winds predominate on the W C V I during a normal 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 in a southeasterly direction, parallel wi th the continental shelf during summer months. In winter months the shelf break current flows in a northerly direction. The transition zone between these two currents varies around the 100-200 m depth contour (Thomson et al. 1989; Foreman and Thomson 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 ia E k m a n transport, where northwesterly summer winds transport surface waters offshore. As the surface waters move offshore, deeper, cooler water is drawn up from below. Upwelling water masses are rich i n nutrients and enhance primary productivity, except in E l Nino years when the thermocline is much deeper (Mysak 1986). In years wi th strong summer northwesterly winds, primary productivi ty can be significantly enhanced off the W C V I (Mackas et al . 1987; Robinson and Ware 1994; Robinson 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 o C o a s t a l C u r r e n t s a n d S h e l f B r e a k C u r r e n t s Shrimp populations in 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 Juan de Fuca strait and directly offshore of Barkley Sound; these occurrences are not considered in this study. In 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 in Chapter 2). These aggregations appear to be right in the transition zone between the coastal currents and shelf-break currents. The location of shrimp populations in 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 in a southerly direction if strong northwesterly winds persisted. A persistent strong northwesterly wind wi l l induce upwelling and drive surface currents offshore and if particles on or near the surface move far enough offshore they wi l l be advected southward by southern moving shelf-break cur- rents (Parrish et al . 1981). Also , strong northwesterly winds can cause inshore surface currents to change direction, and particles floating in the upper 10 m surface waters 63 128 127 126 125 124 123 122 Longitude Figure 3.1: Tow positions of research trawl surveys from 1973 to 2001. Statist ical areas 125, 124, 123 and 121 are shown in black, to lighter shades of grey, respectively. Depth 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. Thus, organisms that undergo diel ver- t ical 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 idal currents (Smith et al. 2001). In this Chapter I examine the effects of various oceanographic indices as an environ- mental correlate for shrimp recruitment. Using the stock recruitment data from Chapter 2, I compare residuals from Beverton-Holt and Ricker models w i t h models that include an environmental variable, using methods similar to Capu t i et al. (1998). I also examine if spawning stock biomass in adjacent statistical areas is correlated wi th recruitment, which would be indicative of larvae dispersing to different areas. 64 3.2 Methods Assuming annual recruitment is a function of spawning stock biomass Rt = f{St) then deviations from an underlying mean function can be represented by additive process error terms: Rt = f(St)eiWt+Vt (3.1) where the vector Wt represents some standardized environmental index wi th a mean 0 and a\Vt = 1 and vt is the residual process error due to factors other than Wt. Note that I assume deviations from the underlying stock recruitment relationship are log-normally distributed (see Rina ld i and Gat to 1978). The scaling term 7 represents the magnitude of the effect the correlate has on the resulting recruitment. If gamma is set to zero, then predicted recruitment (Rt) is strictly a function of spawning stock biomass (St). Two reasonable stock recruitment models for the mean density-dependent relationships are: Beverton-Holt form: = " ^ - l e™** (3.2) 1 + a/bSt-2 Ricker form: Rt = aSte-bSt-2+^Wt+Vt (3.3) Residuals between observed and predicted recruits are: vt = log(Rt) - (log(f(St)) + -yWt) (3.4) where f(St) is either the Beverton-Holt model or Ricker model and Rt, St are the observed recruits and spawners estimated from Chapter 2. Parameters are estimated by minimizing the following negative log-likelihood: \log(av) + \log^) + ^ ^ (3.5) 65 The four model parameters a, b, av and 7 were estimated using a Newton-type algo- r i thm (Dennis and Schnabel 1983; Schnabel and Weiss 1985) in R © . Setting 7 = 0 reduces the model to process error only and the number of parameters from 4 to 3. These are nested models (Hi lborn and Mangel 1997), and twice the difference in log-likelihoods has a x2 distr ibution wi th the degrees of freedom equal to the difference in the number of parameters (df = 1 in this case). The x2 probability is used to assess whether including the environmental index significantly improves our abili ty to predict recruitment. Robust Estimation It is common that ecological data contain large errors, or outliers, that make it difficult to estimate underlying ecological processes using likelihood functions (Tsou and Royal l 1995; Hi lborn and Mangel 1997; Chen and Fournier 1999). In the case of assessing density dependent survival rates for pink shrimp, there were two extreme outliers in the recruit per spawner estimates (Fig. 3.2). Including these outliers tends to over-estimate Area 124 Area 125 -o Lr 1980 1990 Year 2000 1980 1 1 r 1990 2000 Year Figure 3.2: Deviations in annual recruits per unit of spawning biomass (log(Rt/St) — log(Rt/St)) in statistical areas 124 and 125. The dotted line is 3 standard deviations for the distr ibution and points beyond this line were considered outliers. Extremely high juvenile survival rates occurred area 124 and 125 in 1999 and 2000, respectively. maximum 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 wi th a weighting term (ipt) that reduces the influence of extreme outliers on the estimation stock recruitment parameters. Equat ion 3.5 is sensitive to outliers because the normal distr ibution has a small ' ta i l ' 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 probabili ty distribution to have 'fat-tails'(see Chen and Fournier (1999) for examples) or 2) reduce the influence of the outlier in the likelihood by unweighting the residual. Residuals vt were weighted using a modified Tukey's "biwieght": = 0 if |&t| > c, where the modification ensures that 0 < ipt < 1 (the original equation is available in Press et al. 1996, page 702). The error term bt for each t observation represents the c is a constant specifying how sensitive ip is to small departures from the mean. A s c increases in value tpt approaches 1, and al l data are weighted equally. Note that for normally distributed errors the optimal value of c is 6.0 (Press et al. 1996). The "robust" l ikelihood equation is writ ten as: If deviations (bt) are greater that c then the observation is omitted from the likelihood calculation (tpt = 0). The value of c for each statistical area was set equal to 6. (3.6) residuals given ini t ia l estimates of a, b, a, 7 (bt — vt — log(Rt) — (log(f(St)) + ~fWt)) and (3.7) 67 3.2.1 E n v i r o n m e n t a l I n d i c e s I examine six different environmental correlates: sea surface temperature (SST), Bakun upwelling indices in May-June, Sea Level Height (SLH) in M a r c h - A p r i l , and annual S L H anomalies, the Southern Oscil lat ion Index (SOI), and the Pacific Decadal Oscil lat ion ( P D O ) in 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 in Figure 3.3. Month ly mean sea surface temperatures are available from various lighthouse stations 1 since the early 1900's. S S T data were extracted from Amphi t r i t e point, which is located at 48°33'N Lati tude, 125°19 'W Longitude (near Ucluelet, Br i t i sh Columbia) . W a r m periods are represented by positive anomalies (e.g. 1982-3 and 1997-98 E l Nino's in Figure 3.3). Upwelling indices 2 were calculated based on the theory of Ekman , where surface currents generated by wind stress shear 45° to the right of the wind direction in the Northern hemisphere. The transport of the water column is 90° to the wind direction over the E k m a n layer. Typ ica 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. The upwelling index is a measure of net offshore transport of water, and units are usually expressed in m 3 s _ 1 per 100m of shoreline. Positive anomalies indicate upwelling or offshore movement of the surface layer, and negative anomalies indicate downwelling. The pacific decadal oscillation ( P D O ) is actually a sea surface temperature index for the Nor th Pacific Ocean (Mantua et al. 1997). The term P D O actually refers to JData Available at http://www.ios.bc.ca/ios/osap/data/lighthouse/bcsop.htm. 2Data 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 jwt term) is recruitment advection, where larvae are dispersed by currents from 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 in the adjacent statistical area as a correlate of age-1 recruitment. Spawning biomass in the adjacent area was standardized to an index wi th mean = 0 and a = 1 using the same procedure that was done for environmental indices. 3.3 Results Under the assumption that shrimp in areas 124 and 125 are separate independent stocks, the coefficient of variation (a2) ranged from 1.06 to 1.5 in area 124, and 0.6 to 1.04 in area 125 for Beverton-Holt and Ricker models, respectively. The robust likelihood greatly influenced estimates of of the slope at the origin in the stock recruitment curves (see example shown Figure 3.4). In area 124, maximum 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 normal likelihood, improvements in juvenile survival rates are 3 orders of magnitude higher. In other words, if al l observations are treated equally then there is no information in the data about compensation in survival at low spawning densities, and the best statistical stock-recruitment fit is a flat line wi th no St_2 effect. In area 125, the data were uninformative about a using the Beverton-Holt model, and the slope for the Ricker model is roughly 5.65 using the robust likelihood and 10.6 using the normal l ikelihood. In area 124, the M a y - J u n upwelling indices significantly reduced the likelihood over the Ricker model wi th 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 Holt model, the upwelling index increased variance in the residuals (- 7.25%). In area 125, sea surface temperatures (SST) and the Pacific Decadal Oscil lat ion 70 Normal Likelihood Robust Likelihood 0 200 600 1000 0 200 600 1000 Spawners Spawners Figure 3.4: Beverton-Holt and Ricker models fit to stock-recruitment data from statis- tical 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 environmen- tal 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 (r2 — —0.19), and also negative in the reciprocal comparison (r2 = —0.12). These weak correlations do not suggest a linear relationship between recruitment and spawner abundance in adjacent areas. 71 Table 3.1: Comparing the effects of environmental indices on stock recruitment relation- ships for areas 124-5. logL — log-likelihood, x2 = probability of the of the environmental index relative to base model, and the % of the variance in the residuals explained by the enviromental correlate. Beverton Holt M o d e l Ricker M o d e l Mode l logL x2 logL x 2 A r e a 124 Base 19.48 — 21.96 — — S S T 19.33 0.59 3.14 21.68 0.46 5.00 Bak 17.99 0.08 -7.25 19.61 0.03 1.36 S L H 19.42 0.73 2.63 21.86 0.65 -0.54 m S L H 19.47 0.89 0.11 21.90 0.74 0.25 SOI 19.14 0.41 -2.88 21.92 0.77 -1.63 P D O 19.24 0.49 -1.07 21.91 0.75 2.50 Area 125 Base 20.72 — — 26.02 — — S S T 19.12 0.07 13.94 24.51 0.08 12.89 Bak 20.70 0.88 -0.34 26.01 0.91 -0.20 S L H 20.63 0.67 -0.39 26.01 0.88 -0.32 m S L H 20.03 0.24 4.50 25.59 0.36 4.74 SOI 20.49 0.50 5.38 25.85 0.56 3.92 P D O 18.34 0.03 21.10 22.32 •0.01 28.73 Table 3.2: Cross correlations between environmental indices shown in Figure 3.3 S S T B a k S L H m S L H SOI P D O S S T 1 -0.112 0.365 0.624 -0.646 0.727 Bak -0.112 1 0.263 0.192 -0.127 -0.061 S L H 0.365 0.263 1 0.720 -0.208 0.403 m S L H 0.624 0.192 0.720 1 -0.475 0.505 SOI -0.646 -0.127 -0.208 -0.475 1 -0.524 P D O 0.727 -0.0614 0.403 0.505 -0.524 1 Using stock-recruitment models wi th spawning biomass in the adjacent statistical area as a correlate explained very li t t le of the observed variation (Fig . 3.5). The parameter 7 is the slope for a linear relationship between the environmental correlates and recruitment anomalies (7 represents the standard deviation in wt). Estimates of 7 using spawning biomass in the adjacent area were a l l < 0.1 and negative, except for the Beverton-Holt model in 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 in area 125 was a result of spawner abundance in area 72 Table 3.3: Correlations between log(R/S), Log(R) and environmental indices shown in F i g 3.3 Area S S T B a k S L H (Mar-Apr) S L H (.lan-Dec) SOI P D O (Mar-Nov) Log(R/S)12i -0.294 0.095 0.097 0.107 0.220 -0.304 Log(R/S)l25 -0.338 -0.130 -0.018 0.013 0.362 -0.393 Log(R)124 -0.092 0.183 0.165 -0.025 -0.135 0.163 Log(R)i25 -0.395 -0.063 0.001 -0.217 0.305 0.486 124, and vice versa, by estimating stock recruitment parameters using spawner abundance in the adjacent area as the independent variable. No significant relationships were found, wi th or without environmental variables (Table 3.4). Table 3.4: Comparing the effects of environmental indices on stock recruitment rela- tionships for area 124 spawners producing area 125 recruits and vice versa. logL = log-likelihood, % 2 = probabili ty of the of the environmental index relative to base model, and the % of the variance in the residuals explained by the enviromental correlate. Beverton Holt M o d e l Ricker M o d e l M o d e l logL x2 logL X 2 % <r2 A r e a 124 Sp awners A r e a 125 Recruits Base 19.97 — — 19.87 — — S S T 19.59 0.38 2.61 19.09 0.21 5.60 B a k 19.65 0,42 2.40 19.52 0.40 2.50 S L H 19.71 0.47 2.14 19.79 0.68 0.51 m S L H 19.97 0.96 0.00 19.79 0.68 0.52 SOI 19.90 0.71 0.43 19.81 0.72 0.41 P D O 19.89 0.68 0.78 19.63 0.48 1.87 A r e a 125 Spawners Area 124 Recruits Base 19.25 — — 28.82 — — S S T 18.24 0.15 12.32 28.32 0.32 8.16 Bak 19.21 0.78 -0.69 28.14 0.24 -1.24 S L H 19.25 0.95 -0.01 28.60 0.50 0.15 m S L H 18.34 0.18 4.62 28.21 0.27 2.86 SOI 19.23 0.83 2.02 28.78 0.77 2.56 P D O 17.80 0.09 18.87 27.04 0.06 20.45 3.4 D i s c u s s i o n The methods carried out in 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 data 73 Area 124 Area 124 0 100 200 300 400 500 0 100 200 300 400 500 Spawners Spawners Figure 3.5: Stock-recruitment models for area 124 (top row) and area 125 (bottom row) using spawning biomass in adjacent area as a correlate. Closed circles represent observed data, line is the model fit to the data, and open circles represent predicted recruits wi th 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 data is a flat line for the Beverton Holt model, especially area 125. Recall that in Chapter 2, I d id 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 data are informative and contain no outliers (Tsou and Royal l 1995; Chen and Fournier 1999). Using an ordinary least squares approach resulted in extremely high estimates of maximum juvenile survival rate (a) for both statistical areas, fitting either model. The robust likelihood removed the effect of the two large anomalies in 1999 and 2000 in areas 124 and 125 74 (ip = 0 for these two observations). Var iabi l i ty in annual surface temperatures did explain some recruitment variability off the west coast of Vancouver Island, at least in area 125. Similar findings were also reported by Hannah (1993), where shrimp recruitment was negatively correlated wi th early fall sea surface temperatures, roughly six months after larval release. Cooler water temperatures were also a significant abiotic factor in the survival and distr ibution of pink Pandalus borealis off the Scotian shelf in Nova Scotia (Koeller 2000). Poor survival in warm water has also been demonstrated for Pandalus jordani larvae reared in laboratories (Rothlisberg and Mi l l e r 1983). There are two possible outcomes resulting from strong upwelling events off the W C V I . Firs t , larvae could potentially be transported offshore into shelf break currents and settle out in unfavorable habitats resulting in recruitment failure. The largest upwelling event observed was in 1997 (Fig . 3.3) and larval survival in both areas was slightly below average (Fig . 3.6). The strongest May-June downwelling anomaly occurred in 1982. Larval survival rates were near average for area 124, and well below average in 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. The second possible outcome resulting from strong upwelling is an increase in survival due to increases in primary production and food availability. Including the upwelling index as a correlate for area 124 did improve the overall fit to the stock recruitment data, however, it explains very li t t le of the residual variation. A g a i n there were no consistent relationship for both statistical areas that would suggest increased primary production has resulted in increased survival. However, shrimp are pr imari ly 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 in the years following the strong 75 upwelling event in 1997, juvenile survival rates have been increasing in both statistical areas (Fig. 3.6). Area 124 Area 125 1980 1990 2000 1980 1990 2000 Year Year Figure 3.6: Survival index (log(Rt/St)) for age-1 recruits i n area 124 and area 125. Dot ted 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 in higher than average recruitment in Oregon (Hannah 1993). Higher sea level heights generally mean stronger coastal currents and Hannah speculates that in years when sea levels are extremely high (i.e. 82-83 & 97-98 E l Nino years) larvae are advected northward resulting in recruitment failure in Oregon- Washington 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 in area 124 following this recruitment event. Four years later, recruitment begins to improve in area 125, which also led to an rebuilding of the stock in that area. U p unt 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 in both statistical areas, especially i n area 125 in more recent years. Dur ing this period, sea level height anomalies have been negative and apparent juvenile survival rates have been high. This is also the case for the above average juvenile survival rate in area 124. Thus it appears that large negative sea-level 76 height anomalies are a better indicator of juvenile survival rate than they are as an index of advection from southern spawning grounds to more northern areas. The t iming 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). The 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 in 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 in velocity due to increased runoff from the Fraser river. If the spring transition is delayed relative to the t iming 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. Hannah (1993) also suggested that the t iming of the spring transition was important in determining the strength and distribution of year classes in Oregon populations. There is no clear evidence that shrimp in 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 in area 124, in fact the sign of the correlation suggest that the opposite. In either case, using spawning abundance in the adjacent area as a correlate for recruitment in that area did not explain much of the observed variation. A l l of the results suggest that recruitment is probably occurring at much larger spatial scales than those examined here. Chapter 4 Shrimp Predation by Gadids 4.1 I n t r o d u c t i o n Several recent studies have indicated strong predator-prey interactions between members of the family Gadidae and Pandal id shrimps (Berenboim et al. 2000; Hannah 1995; L i l l y et al . 2000). Assuming that gadids, such as Pacific cod (Gadus macrocephalus) and Pacific hake (Merluccius productus), are principal predators of pink shrimp species (L i l ly et al . 1998; Westerheim 1996), here I examine correlations between Pacific cod and Pacific hake abundance versus pink shrimp abundance, mortal i ty and recruitment off the west coast of Vancouver Island. Previous work in Oregon has demonstrated a positive correlation between hake abundance and pink shrimp mortali ty (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 mortali ty should be positively correlated and shrimp recruitment and predator biomass should be nega- tively correlated. The later prediction represents predation on juvenile shrimp, or that juvenile survival rates for pink shrimp are correlated wi th predator abundance. If preda- tion plays a minimal role in 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 es- timates 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 trawl catch per unit effort of Pacific cod off the west coast of Vancouver Island (scale units are log(c.p.u.e). Areas 3C 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 in a l l areas. There have been no directed, fishery independent surveys for Pacific cod on the west coast of Vancouver Island. The most recent D F O assessment document (Sinclair et al . 2001) included an index of abundance generated from shrimp trawl surveys. A l though the index was not used for the D F O assessment, the trend information from shrimp surveys seems to correlate well wi th commercial indices since 1981. Because I am interested in calculating predation mortali ty on shrimp by Pacific cod, I choose to include the shrimp survey index of Pacific cod i n this assessment. 80 4.1.2 S t o c k S t r u c t u r e o f P a c i f i c H a k e Two distinct types of Pacific hake (Merluccius productus) populations exist along the West Coast of Nor th Amer ica , a migratory coastal stock, and several isolated non- migratory 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. The 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 in central California, adult coastal-stock hake move onshore in the spring to feed and then migrate northward. B y mid summer, adult hake form large midwater aggregations near the continental slope off Vancouver Island and feed on euphausiids, other pelagic schooling fishes, and Pandal id shrimps (Dorn 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 in late summer (Ware and McFarlane 1995). The portion of the coastal stock that migrates into Canadian waters each summer consists of older (age 5+) hake, and the sex ratio is skewed to favor females (Dorn et al . 1999); however, proportions of fish less than age-5 in the Canadian zone have increased since 1996. Dur ing E l Nino events, a higher fraction of the stock migrates into Canadian waters, and ventures further north into Alaskan waters as shown in Figure 4.2 (Dorn et al. 1999; Ware and McFarlane 1995). Dur ing the warmer years, the distr ibution of juveniles also spreads north to Oregon and into Canadian waters. It is thought that this change in juvenile distribution has led to the break down of patterns indicating among-cohort 81 juvenile cannibalism in the stock-recruitment data sets (i.e. strong cohorts appearing every 3-4 years, see Dorn et al. 1999). Est imat ing the fraction of mortali ty caused by Pacific hake on pink shrimp popula- tions off Vancouver Island is difficult since the fraction of the hake stock in Canadian waters must be estimated. Previous research has documented a strong correlation be- tween hake abundance in the Canadian zone and E l Nino events; more hake are found in Canadian waters during E l Nino years (Mark Saunders, pers comm., Pacific Biologi- cal Station, Nanaimo, B C ) . Therefore, if hake are significant predators on pink shrimp, mortali ty rates for pink shrimp should be significantly higher in E l Nino years. P i n k shrimp are not listed as one of the major diet items of hake feeding in the Canadian zone (Tanasichuk et al . 1991). However, i n Washington and Oregon, predation mortali ty by Pacific hake on pink shrimp is thought to be significant (Rexstad and P ik i t ch 1986; Han- nah 1995) even though pink shrimp are a rare diet i tem for Pacific hake in that region. A l l observations of pink shrimp found in the stomaches of Pacific hake report less than 10% of pink shrimp in the diet (see Table 4.1), wi th the exception of Gotshal l (1969). In Gotshall 's case, he was t ry ing to develop a relative abundance index for pink shrimp and specifically sampling hake aggregations in areas where shrimp are known to be present. P ink shrimp are found more frequently in the stomachs of larger hake (> 55 cm T L ) than in smaller hake (Rexstad and P ik i t ch 1986), and larger hake are generally found closer to the bot tom (Mark Saunders, pers comm, Pacific Biological Station, Nanaimo, B . C . ) whereas smaller hake generally form large pelagic schools. Rexstad and P ik i t ch (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. The 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: Dis t r ibut ion of Pacific hake in 1992, 1995, and 1998 along survey transect lines. Abundance 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 % Occur- % by Source rence Weight California 54% Gotshall (1969a, b) British Columbia 3% Outram and Haegele (1972) Washington and Oregon 5.70% Alton and Nelson (1970) Washington and Oregon 0.30% Livingston and Alton (1982) Oregon to Vancouver Island 0.70% Livingston (1983) Washington and Oregon 1.70% Rexstad and Pikitch (1986) 4 . 2 E s t i m a t i n g P r e d a t o r B i o m a s s 4 . 2 . 1 P a c i f i c C o d A s s e s s m e n t 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 delay- difference model is developed by assuming that (1) growth can be adequately described by the Ford-Brody growth function Wa+i — a + pWa, (2) fish are fully recruited to the fishery at age 3, and (3) natural and fishing mortality rates are independent of age for ages 3 and older (Hilborn and Walters 1992). Parameters and variables used in the delay- difference 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 Ro Unfished equilibrium age-2 recruits Br Ratio of initial biomass to unfished equilibrium biomass 5 Recruitment compensation parameter Q Catchability coefficient for commercial fishery M Instantaneous natural mortality rate Recruitment anomaly in year t Growth Parameters Value Asymptotic length 89.48 k Growth coefficient 0.307 to Time at length 0 -0.116 a Length-weight scaler 7.38E-06 b Length-weight power 3.0963 K Age of recruitment 2 Observed States Et Annual fishing effort ct Observed landings in tons wt Mean body weight of landed fish Unobservec States B0, Bt Unfished biomass, biomass in year t No, Nt Unfished numbers, numbers in year t wt Predicted mean body weight in year t Rt Age K recruits in year t ct Predicted annual catch in year t Ft Fishing rate in year t st Survival rate in year t Derived Parameters so Maximum recruits per unit of spawning stock biomass b Maximum recruitment scalar Mean body weight at equilibrium WK Weight-at-age of recruitment S Annual survival rate (e~M) a Intercept of Ford-Brody growth equation P Slope of Ford-Brody growth equation 85 G r o w t h Parameters for the Ford-Brody growth equation can be calculated from parameters from the von Bertalanffy growth equation as: 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). In i t ia 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 (R0), natural mortality (M), and recruitment compensation parameter (<5), the unfished equilibrium numbers (N0), mean body weight (wa), biomass (Ba) and stock recruitment parameters (s0,(3) are initialized as follows: oo ( 1 - / 5 ) N0 = Ro (4.1) ( 1 - 5 ) {Sa + wK(l-S)) 1-pS (4.2) Ba = N0w0 (4.3) (4.4) (4.5) Note that Equation 4.5 describes the Ricker stock-recruitment model. 86 State Dynamics Populat ion numbers and biomass over t ime were calculated using equations 4.8 and 4.9, respectively. Annua l survival (St) is a function of fixed natural mortal i ty and annual fishing mortality, where mortal i ty associated wi th fishing is condit ioned on observed effort (Et, Equat ion 4.6). Annua l recruitment (Eq. 4.10) was calculated using a Ricker stock-recruitment curve, and annual recruitment anomalies (<f>t) were treated as model parameters to be estimated. Ft = qEt (4.6) St = e ( - M ~ ^ (4.7) Nt+1 = NtSt + Rt (4.8) Bt+1 = St(aNt + PBt) + wrRt (4.9) Rt = SoBt-KeWt-rtet* (4.10) Pred ic ted Observations Observed data consisted of annual catches and mean body weights in 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)) Bt (4.12) Likel ihoods Likel ihoods for mean body weights and observed catches are represented by equations 4.13 and 4.14. Here I assume that errors in reported catches and observed mean weights 87 are lognormaliy distributed, and recruitment anomalies (process errors) also have a log- normal distribution and are penalized to have a mean 0 with variance cr2,. Equation 4.16 represents the likelihood of the shrimp trawl survey index, where Zt — ln(Yt/Bt), which also assumes that errors in relative abundance estimates are lognormaliy distributed (see Walters and Ludwig 1994). Lw = —-—log Lc = —-—log 5>(t) ! - 1 ) log (Zt - Z)' (4.13) (4.14) (4.15) (4.16) Equations 4.13, 4.14 and 4.16 are self weighting by using the conditional maximum likelihood estimate of the variance (i.e. substituting ^ ^ " ^ for a2 in the negative log likelihood results in ^^Log (XX1 _ ^)2)> ignoring a constant scaling term). I further use prior distributions on recruitment compensation (p0(5)) and natural mortality rates (po(M)y. 2 Po(M) = Po(S) = 1 , (M) —7T- In ?rr2 z c r M V0.58 J 1 , / —-o In 2a2 \ 2.74) (4.17) (4.18) 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 a2^ — 0.6, for natural mortality a2M = 0.38 (Westerheim 1996), and for recruitment compensation erf — 0.2 (Myers et al. 88 1999). The objective function to be minimized is the sum of equations 4.13 through 4.18. Likel ihood profiles were used to explore uncertainty in leading parameter estimates ( M , RQ, and 5). Uncertainties i n nuisance parameters, such as recruitment anomalies, were not explored in this analysis. I also doubled the value of the variance terms in the prior distributions to examine their influence on leading parameters. 4 . 2 . 2 P a c i f i c H a k e A s s e s s m e n t Coastal hake is treated as a single stock and is formally assessed by U . S . and Canadian agencies every 3 years using fishery independent trawl and acoustic surveys. Abundance is assessed using age-structured models, and fishing mortali ty is partit ioned into U . S . and Canadian components. The relevant issue for my analysis is to determine what fraction of the total estimated stock is present off southwest Vancouver Island in the summer months, and how much shrimp are consumed by hake. There are two sources of information about hake abundance in the Canadian zone: data from the triennial surveys, and volume swept estimates from Ware and McFarlane (1995). The area over which Ware and McFarlane 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. The triennial surveys span the entire coast, from California to Alaska but again there is no detailed information available that would allow for estimation of hake abundance i n areas 124 and 125 specifically. So in the absence of direct observation on the abundance of Pacific hake in areas 124 and 125, I use the biomass estimates provided by Ware and McFarlane (1995) to assess possible predation impacts of hake on pink shrimp. Abundance data for Pacific hake from both the triennial surveys and Ware and McFarlane (1995) are presented in Table 4.3. Note that the methods used to estimate abundance in these two documents differ 89 Table 4.3: Biomass Estimates for Pacific hake for L a Peruse bank region (Ware and McFarlane, 1995), and for hake in the Canadian Zone, (Dorn et al, 1999). Y e a r H a k e b i o m a s s i n L a P e r u s e b a n k r e g i o n ( K t ) H a k e b i o m a s s i n C a n a d i a n Z o n e ( K t ) 1968 183.3 1969 226.42 1970 79 1974 53.25 1975 108.8 1980 184 1983 438.7 242 1985 65.4 1986 309.5 335 1987 256 1988 207.8 1989 138.3 138 1990 425.5 1991 251.8 1992 554 1995 199 1998 526 and cover different areas, so for years where abundance information is available from both documents, estimates differ. There appears to be a significant relationship between the fraction of hake that enter the Canadian zone and winter sea level height (Mark Saunders pers comm.). Higher winter sea levels i n Bamfield are positively correlated wi th the fraction of the total hake stock that enters the Canadian zone (Figure 4.3). Using this relationship I was also able to construct a time series of hake abundances i n the Canadian zone, using sea-level height as a predictor for the proportion of hake present in the Canadian zone from 1976 to 1998. Here I fit a logistic function to the observed time series data shown in the top panel in Figure 4.3, where I assumed the asymptote of this function is 1 and estimated two parameters 7 and Xh assuming the deviations are normally distributed. Note that Xh corresponds to the sea level height anomaly at which 50% of the stock was present in the Canadian zone, and 7 is a shape parameter that describes how steep the relationship 90 is. Proportion in Canadian Zone vs. Sea Level Height i i i i i i r~ -100 -50 0 50 100 150 200 Sea Level Anomaly (mm) Hake Biomass i 1 I 1— 1980 1985 1990 1995 Year Figure 4.3: Propor t ion of Pacific hake i n the Canadian zone versus mean January- February sea level height anomalies (top panel). Positive anomalies indicate higher than average sea level heights in Bamfield, Br i t i sh Columbia (data provided by M a r k Saunders, D F O ) . Est imated hake biomass (bottom panel) for total hake stock (solid line) and hake present in the Canadian zone (dashed line) based on the sea level height relationship. 4.3 Results 4.3.1 P a c i f i c C o d B i o m a s s The reconstructed biomass and fishing mortali ty rates for Pacific cod are shown in Figure 4.4 (top panel, residuals and recruitment anomalies shown in bottom panel). Dur ing 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 maximum of 24,219 tonnes in 1972. From 1973 to 1986 biomass declined at a rate of 9% per year, reaching a low of 7,941 tonnes. Dur ing the mid-late 1980's cod populations began to increase reaching another peak in 1989, and reached an all time low of 1,808 tonnes in 1996. The trend in estimated biomass shown in Figure 4.4 is consistent w i th the results obtained by Haist and Fournier 91 (1995). Since 1996, the ground fish t rawl fishery has had new quota regulations, as well as increased mesh sizes. Pacific cod in recent years is predicted to be increasing in abundance. Biomass Fishing Mortality o 1960 1980 2000 1960 1980 2000 1960 1980 2000 Year Year Year Figure 4.4: Est imated 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. Bot tom row from left to right are recruitment anomalies, catch residuals, and mean weight residuals from the delay difference model. From the 1950's to 1990, fishing mortali ty rates (Ft) averaged near 0.l6~yr, and dramatically increased to a maximum of 0 .67~ y r in 1994. Since 1996, however, the new trawl quota management plan has apparently reduced fishing mortali ty on Pacific cod to om~yr. 92 Management of the ground fish trawl fishery underwent a transformation in 1996, where quota limits for each species were imposed on the fishery. This change reduced targeting on spawning aggregations during the winter fishery, and one of the outcomes has been a decrease in the mean body size landed. In the past Pacific cod assessments have been fit to catch-at-length data, and when I included this data set in the assessment, apparent recruitment in recent years more than doubled. Pr ior to 1996 there was a low probability of capturing small cod because the fishery targeted larger cod concentrated in spawning aggregations so I chose to omit the catch-at-length data for this reconstruction. Also , recent commercial C P U E indices may not be informative about changes in abun- dance if fishermen are effectively avoiding cod to prevent exceeding the species quota. Given 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 rawl survey, which indicate a small increase in 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 natu- ral mortali ty rates, and recruitment compensation (see F i g . 4.5). Catch-at-length data might be informative about changes in age structure, and if used in future assessments I would recommend including a second selectivity curve for the period post 1996. In general, the relative abundance data and mean body weight data were informative about key population parameters for W C V I cod. Changing the mean of the prior distri- butions by 10% resulted in less than 1% changes in key populat ion parameters. Results were also fairly robust to assumed variances on prior distributions for natural mortal i ty rates and recruitment compensation parameters (Figure 4.5). 93 Recruitment Compensation 2 3 4 S Natural Mortality 5 0.45 1 0.50 0.55 0.60 0.65 M Unfished Recruitment 0.70 0.75 o.ao i 0 2000 4000 6000 6000 10000 Figure 4.5: Likel ihood profiles for W C V I Pacific cod population parameters where as- sumed variances for natural mortality rates a2M = 0.38 and recruitment compensation a2 — 0.2 (solid line) versus doubling these variance terms (dashed lines). 4 . 3 . 2 C o r r e l a t i o n w i t h P r e d a t o r A b u n d a n c e P i n k shrimp and Pacific cod (adult and juvenile) abundances on the W C V I are posi- tively 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 in statistical area 124. Shrimp and Pacific cod abundances were high i n the mid 1970s and reached lows a decade later. Dur ing the late 1980s and early 1990's both shrimp and Pacific cod in- creased in abundance. Hake populations for the entire west Pacific had also increased from the mid 1970's and peaked in abundance in 1987 (Dorn et al. 1999). No significant correlation was found wi th either Pacific hake time series presented in Table 4.3 or wi th the reconstructed abundance from sea-level height in Figure 4.5. Correlations between predator abundance and natural mortality, and predator abun- dance and recruitment for shrimp in both statistical areas are also presented in Table 4.4. Evidence for interactions between predator abundance and shrimp mortali ty 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 0.66 25 18.07*** 0.24 25 1.41 0.55 25 10.54*** Pacific cod 0.59 25 12.95** 0.34 25 3.18* 0.56 25 10.89*** Pacific Hake Ia 0.04 22 0.03 -0.03 22 0.02 0.01 22 0.00 Pacific Hake II6 -0.39 8 1.23 -0.10 8 0.07 -0.31 8 0.73 Pacific Hake IIIC -0.34 6 0.67 0.26 6 0.36 -0.05 6 0.01 Mortality Juv. Pacific cod -0.27 25 1.90 -0.05 25 0.07 -0.16 25 0.61 Pacific cod -0.04 25 0.05 -0.04 25 0.04 -0.05 25 0.06 Pacific Hake I -0.45 22 5.47** -0.16 22 0.55 -0.32 22 2.33 Pacific Hake II -0.17 8 0.21 -0.17 8 0.21 -0.26 8 0.52 Pacific Hake III 0.06 6 0.02 -0.41 6 0.99 -0.18 6 0.17 Recruitment Juv. Pacific cod 0.60 25 13.61** -0.25 25 1.59 -0.08 25 0.14 Pacific cod 0.32 25 2.75 -0.21 25 1.10 -0.12 25 0.35 Pacific Hake I 0.42 22 4.53** 0.09 22 0.18 0.35 22 2.87 Pacific Hake II -0.56 8 3.24 -0.22 8 0.37 -0.66 8 5.37* Pacific Hake III 0.38 6 0.86 0.46 6 1.34 0.56 6 2.23 "Hake abundance determined from sea level heights 6Hake abundance data from Ware and McFarlane (1995) cHake abundance data from Dorn et al. (1999) 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 pos- itively 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 2.0 0.5 1.0 1.5 2.0 0.5 1.0 1.5 2.0 Shrimp biomass Shrimp biomass Shrimp biomass Figure 4.6: Relative abundance time series of P ink 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 D i s c u s s i o n In contrast wi th the east coast of Canada, Alaska and the Barents Sea, shrimp abundance and Pacific cod abundance were positively correlated on the west coast of Vancouver Island. Dur ing the mid 1980's both 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 both species in the early 1990's. The 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. The proximate agent for this decline could be predation by hake, causing increased mortality 96 rates for shrimp in years when hake is abundant; however there were no significant positive correlations between hake abundance and shrimp mortali ty (see Table Tab4.5). Dur ing warm years more hake move into the Canadian zone, apparently resulting in increased predation mortali ty for many other species including herring and euphausiids (Ware and McFarlane 1995; Robinson and Ware 1999). Th is is especially evident during the E l Nino 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 wi th 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 in 1985 and 1995 (Figure 4.4), and significant positive correlations between shrimp natural mortali ty rates and Pacific cod abundance. Such observations have occurred on the east coast of Canada (L i l ly et al. 2000, W o r m and Myers In Press), and i n the Barents Sea (Berenboim et al. 2000). O n the west coast of Vancouver Island, pink shrimp were positively correlated wi th Pacific cod, suggesting that both species are responding to changes resulting from bottom-up forces (see Chapter 5). Th is 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. Typical ly , in warm E l Nino years, the distribution of hake may extend as far north as Alaska (e.g. Figure 4.2, in 1998), and there is no published information available about how long and how many hake are present in statistical areas 124 and 125 during the summer. Given the l imited data set adjacent to area 124 (Ware and McFarlane 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 in area 124 was the same as the adjacent area. The 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 Rexstad and P ik i t ch (1986) of 659 tonnes to the west coast of Vancouver Island population in 1983, more shrimp would have been consumed than were available. A long 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 data from Ware and McFarlane (1995)) to shrimp biomass is approximately 93:1, so if hake diets consisted of 1% shrimp this would translate into a mortali ty 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 in the Canadian zone for roughly 180 days during warm 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 August and September where the average diet fraction is 32% by weight (Ware and McFarlane 1995). Robinson and Ware 1994 also noted that increases in temperatures reduce productivi ty of euphausiids, therefore in warmer years hake tend to rely on other forage species, such as herring. Given the available data, I was not able to conclude that migratory Pacific hake or Pacific cod play a significant predation role in shrimp population dynamics. In the case of Pacific cod, it appears that abundance is in phase wi th pink shrimp populat ion 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 pink shrimp i n the stomachs of Pacific cod is well documented in this system as well as for gadoids in other systems. Pacific hake have also been documented as foraging on pink shrimp and significant predation mortali ty is suspected in 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 in abundance between hake and shrimp indicate that only a small fraction of shrimp in hake diets would cause significant predation mortality. The variabili ty in 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 mortali ty also appears to be related to temperature; it is possible that hake are responsible for at least some of the variabil i ty in pink shrimp natural mortality. Chapter 5 Ecosystem Perspectives "We beiieve the food web modeling approach is hopeless as an aid to for- mulating management advice; the number of parameters and assumptions required are enormous." .. .Ray Hilborn and Carl Walters, 1992 (pg. 448). 5.1 Introduction Much of fisheries science in the past has focused on single species assessment, and ignores or assumes that the interactions that determine distribution and abundance remain con- stant 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 en- vironmental indices (e.g. Walters et al. 1986; Hannah 1993; Caputi et al. 1998; Lilly 99 100 et al. 2000; Berenboim et al . 2000). Understanding the interactions between environ- mental variables and production is meaningless in 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 environmen- tal correlations break down over time (e.g. Drinkwater and Myers 1987). Alternatives to single species models driven by correlated information are more complex food web models or multi-species models (e.g. "Bormicon" , Stefansson 1998; " E R S E M " , Baret ta et al. 1995, "Ecopath wi th Ecosim" Walters et al . 2000) where interactions (trophic and fishery) and variation i n primary production can be explored. Exp l i c i t ly modeling the dynamic changes in food webs over time allows for examination of direct interactions (such as predator-prey interactions), as well as indirect interactions (for example how her- ring 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 in ecosystem productivi ty 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 Ecopath wi th Ecosim (Pauly et al . 2000) to part i t ion mortali ty among natural predators, and commercial fisheries and examine how total mortal i ty changes over time. The first step in demonstrating the ut i l i ty of an ecosystems model is to show that it can replay the history of disturbance in the ecosystem. To do this I w i l l first fit time series data to Ecosim predictions. Th i s accomplishes two goals, first it assures that Ecosim can replay the general trends in observation data (at least qualitatively), and second it helps in parameter estimation for the model. Coastal upwelling in the W C V I ecosystem is a very prominent feature and plays an important role in pr imary production. Historical annual 101 variation in primary production is estimated by fitting Ecosim to observed abundance in- dices, total mortality, absolute catches and changes in mean body weight over time. This reconstructed variation as evidenced in animal abundance is compared to other estimates of primary production. Bot tom up control in Ecosim is represented by specifying rates of flow from prey to predator in a term known as flow control, or relative vulnerabili ty of prey to its predators (see Bundy 2001; Walters et al. 1997). 5 . 2 . 1 E c o p a t h w i t h E c o s i m Ecopath is a trophic mass-balance modelling approach, which is based on the assumption that total production in a particular reference year, or period of years, for any group in an ecosystem can be expressed as: where Bi is the biomass (t k m - 2 ) of group i, (P/B)i is the production to biomass ratio of 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 mortali ty that is accounted for wi th in the modelled system and Y", are fisheries yields. Equat ion 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 in question (usually one year). Here I define the P/B ratio as the total instantaneous mortali ty rate (Z = F + M) during the 1950's. The defined area for the west coast Vancouver Island Ecopath model includes statisti- cal areas 121 through 126 out to the 200 m depth contour. The total area is approximately 11,000 k m 2 . Biomass estimates and total landings were scaled to tonnes k m - 2 for groups 5.2 Methods n (5.1) 102 where information from stock assessment reports could be obtained. For groups where no biomass information is available, I allowed Ecopath software to estimate biomass by specifying the fraction of total mortal i ty being explained by the model (Ecotrophic Efficiency). Ecosim is initialized using the Ecopath inputs, and biomass dynamics for most groups in the model are then predicted using the differential equation: SB n = f(Bht) - Mii0Biit - FiitBiit - d,j,t {Biit, Bitt) (5.2) where % indexes group, Mifi = (1 — EE/) * (P/B)i represents unaccounted mortality (e.g., other mortali ty), Fi<t is the fishing mortal i ty rate at time t and c i j ] t is the consumption of group i by predator j at time t. For primary producers, the biomass production function f{Biit) is a simple saturating relationship of the form: f (Bz,t) = riBi>t/(1 + Bhthi) (5.3) where is the maximum intrinsic rate of growth, and hi is inversely proportional to maximum density. Biomass production for consumer groups j is modelled as the product of total consumption and a growth efficiency g^. n f{Bitt) = gi^citi{BitUBitt) (5.4) where growth efficiency = Bi(P/B)i/Qi is calculated from the balanced Ecopath esti- mates. Predict ing the dynamics of consumption, i.e. the "flow" rates Cjti, in Ecosim is based on the "Foraging Arena" concept (Walters and Juanes 1993). The basic premise of this theory is that each biomass pool is split between two behavioral states: an "available state" (Vi<t) representing prey that are vulnerable to consumption by predators, and an "unavailable state" (Biit — Vitt) representing individuals that are not vulnerable due to 103 prey behaviors such as schooling and hiding in refuge areas. Exchange of prey between available and unavailable states is a function of time that prey spend foraging, which may vary in response to changing levels of food availability and abundance of competitors and predators. Ecosim assumes that prey exchange rapidly between these two states, and simply keeps track of the (time-varying) equil ibrium biomass of available prey. Thus, consumption of prey i by predator j at any moment should vary as: v = Vij + Bitt 2v^ + ciijBjt (5 5) Cij (Bit,Bjt) — ciijVijfBjf For this model, aitj is the rate of effective search for prey i by predator j and v^j is the rate of exchange of prey i between "available" and "unavailable" states. The rate of effective search a^j is estimated from re-arrangement of Equat ion 5.5 using the Ecopath estimates of Qijto, -Bi,o a n d -Bj,o- Vulnerabil i ty exchange parameters v^j are calculated from a user specified ratio of the maximum predation mortal i ty to base predation mor- tali ty estimated in Ecopath (see description on page 106). Large vulnerabili ty exchange parameters (v^ » 1) represent top-down or Lotka-Volterra type control; whereas, small values represent bottom-up or donor controlled rates of flow from prey pools to preda- tor pools. Note that Equat ion 5.5 predicts a saturating relationship between predation mortali ty M y = Cij/Bi and total predator biomass. Th is form of relationship tends to eliminate much of the unstable behavior that Lotka-Volterra type predator-prey models typically produce. Th is tends to conserve total biodiversity in even the most complicated ecosystem models due to limitations on total predation mortal i ty rates (Walters et al. 1997; Walters et al. 2000; Walters and Ki tche l l 2001). 104 Age/Size structure in Ecosim Juvenile and adult fish typically feed on different prey resources and are vulnerable to different types of predators and fisheries. Simulating this trophic ontogeny cannot be accommodated using the simple biomass dynamics model in Equat ion 5.2. For species that are represented as juvenile and adult "pools", Ecosim uses a Deriso-Schnute delay- difference model to represent adult dynamics and an explicit age-structured model that keeps track of monthly juvenile cohorts from bi r th to the age at recruitment to the adult stock. The complete set of derivations and assumptions used i n this formulation can be found in Walters et al . (2000). Here I present the basic model structure for split-pool dynamics. Indexing adults as A and juveniles as J (group index i is implic i t ) , the basic delay- difference model structure for al l split pools is: BA,t+i = e-ZA*[aAt{CAt)NA,t +PABA,t)+wJAtNj^t (5.6) NAit+1 = NAite-z^ + Njikit (5.7) Nj,llt+i = R(BA,t,NA,t,CA<t) (5.8) Nj,a,t+i = e ' ^ N j ^ t l<a<k (5.9) wj,a,t+i = wj<a-itt + gCj}t/Njtt (5.10) where t is time (months); A^>,t is number of juveniles of age a (months); R(BAt, NA}t, CAit) predicts age-0 juvenile numbers as a function of adult numbers, biomass and food con- sumption; k is age (months) at graduation to adult pool; u>j,0,t is age a juvenile body weight at time t; Cjjt is total food consumption by juveniles in month t (2~2kck,J,t)'i' Njj is total number of juveniles; g is juvenile growth efficiency; aA,t(CAjt) is the inter- cept of Ford-Brody growth model (wAlt+i v s - wA,t) which is assumed to depend on per capita food consumption CA}t (J2k ck,A,t/NAjt); and p is the Ford-Brody slope representing 105 metabolism. Ecosim predicts total mortality rates of adults ZA,t and juveniles Zj>t as the sum of un- explained mortali ty M 0 , total fishing mortali ty Ft (described below) and total predation losses as: ZA,t = MAfi + FAit + Cj,A,t/BAit (5.11) j Zj,t = MJfi + Fj,t + CJM/Bj,t (5.12) 3 where the summation is over al l predators j that eat group i adults and juveniles. D e n s i t y - d e p e n d e n t m o r t a l i t y i n j u v e n i l e p o o l s Ecosim produces density-dependent juvenile mortal i ty from age-0 to age-fc i n two ways. Firs t the user may specify a constant weight at recruitment w so that only those juve- nile age classes that have body mass larger than w recruit to the adult stock. If prey vulnerabili ty rates to juveniles are set low, then per capita consumption rates decline rapidly wi th increasing juvenile abundance. This effectively means that under conditions of reduced prey availability, juveniles grow slower and spend more time exposed to high size-dependent predation mortali ty rates (i.e., longer time required to reach w). Second, simulated juveniles in the model adjust the amount of time spent feeding i n response to changes in prey availability and intraspecific competition. Ecosim represents this effect by using a function that adjusts feeding time to t ry to maintain near constant food intake rates per unit time, and feeding time is assumed to be time spent at risk to predation (i.e., adjusts total time required to consume fixed food ration). The model first computes the Ecopath base food intake rate copt which is assumed to optimize a feeding rate versus predation risk function. A t each time step, Ecosim computes the realized food intake rate per juvenile biomass ct and then adjusts time spent at risk to predation (i.e., t ime 106 spent feeding) as: Tt = Tt-i(l-a + acopt/ct-i) (5.13) where a is a user defined [0-1] value that defines the proportional increase in time spent feeding per unit change in realized food intake rate ct~\. This leads to a more complicated model for predicting consumption of prey i by al l its predators j (including juveniles as predators): IT) V) \ — aijvi,jBi}tBjitTiitTjjt / C I A \ C,-JJ { B l " B j x ) ~ (vhJ + vhJTht + ahjBj/T^ ( 5 1 4 ) Initially, I used the feeding time adjustment method to represent density-dependent mortali ty for juvenile pools in the W C V I model. Here I set feeding time adjusment rate a = 0.5, and a low vitj for juvenile pools. I d id not allow for foraging time adjustments (a = 0) in adult pools and set Vij higher for the adults groups. The combined effects of low prey vulnerabili ty and feeding time adjustment creates a compensatory relationship between juvenile abundance and survival which is supported by the substantial number of stock-recruitment data sets now available (Myers et al. 1999). Further adjustments were made to exchange rate parameters (v^) to better explain observed time series data. The vulnerabili ty exchange rate parameters take on a variety of positive values de- pending on Ecopath input parameters. Input parameters specified in Ecopath pro- vide a base estimate of the predation mortal i ty rate of predator i on prey j (My-o = Bjfi(Q/B)jDCij/Bifl as represented by the dark circle in Figure 5.1). The user specified "flow control" parameters (kij) are used to specify a ratio of maximum predation mor- tality to the base mortality. The "vulnerability" parameters vtj are specified as multiples of Mijfi: 1 + k •• = Miifi kij[0,1] t Kij where the user specifies values k^ between 0-1 (the "Flow control" parameter in Ecosim). The Vij's represent the asymptote, or the maximum predation rate when either predator 107 (Bj) or prey biomass (Bi) becomes infinitely large. As the value of kij approaches 1, Vij tends to infinity and predation mortali ty increases linearly wi th increased predator biomass (slope of the line equals the encounter rate a^), or mass-action type interactions (Figure 5.1). i 1 1 1 1 r 0 1 2 3 4 5 Predator Abundance Figure 5.1: A n example of the relationship between predation mortali ty rate and predator biomass for different values of k^. The dark circle represents the base predation mortality rate m ^ o , calculated using Ecopath input parameters. A s the value of k^ decreases the slope of the line through S^n, m i j , o decreases and the maximum predation rate (vtj) of predator j on prey i decreases. 5 .2 .2 F i t t i n g E c o s i m T o T i m e S e r i e s D a t a Data Types T w o general classes of time series data can be input for an Ecosim scenario, these are forcing data and observation data. Forcing data can include abundance, fishing mortality, or relative fishing efforts, and are used in calculating dynamic changes in the model (e.g. Fishing effort time series used to calculate fishing mortali ty rates over time). Observation data come in 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 data types: abundance (biomass), total mortality, absolute catch, and mean body weight (catch and body 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 from fishery independent surveys would be preferred over commercial fishery samples 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 de- termine 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 deter- mined 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 pe- riods 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 Zt>i = log ̂ p ° ^ ) • Note that the t subscript denotes time and 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 al l observed data given model predictions was then calculated using: There is no l imi t to the number of data sets that can be used in the fitting criterion, and the contribution of each data set to the total l ikelihood can be modified by the user (wi). Hence, data that are known to have large errors relative to other data sets should carry less weight. Note, however, that equation 5.17 is evaluated at the conditional maximum likelihood estimate of the variance for each data set and setting all u>j = 1 assumes that al l data are weighted proportional to the inverse of the variance. Est imat ing model parameters Two types of model parameters can be varied by the search routine, which uses a modi- fied Marquardt nonlinear search algorithm (see Ecopath Software Help Manual ; "F i t t ing Ecosim to T ime Series Data" for a description of the Marquardt Algor i thm), to maximize the total log-likelihood function (Equation 5.17). These parameters are vulnerabili ty ex- change rate parameters (vij) and a time series of pr imary production anomalies. In analogy wi th single-species assessment models, vulnerabili ty exchange rate parameters represent mean production responses, and primary production anomalies represent pro- cess errors. 5 . 2 . 3 D e s c r i p t i o n o f M o d e l G r o u p s The ecosystem model was used to examine trophic interactions that contributed to vari- ability in pink shrimp abundance and mortality. Therefore, I chose to bui ld a very simple model (fifteen groups including a detritus pool) that represents the main trophic inter- actions and fisheries off the west coast of Vancouver Island. This simplification was done to avoid creating a complex model that would include many groups for which data are (5.17) I l l lacking. M o d e l groups included known shrimp predators (Pacific cod, Lingcod, and Dog- fish), pr imary bycatch species in shrimp trawl fisheries (eulachon and demersal fishes), principal forage species (euphausiids and herring) and the bottom of the food web is represented by phytoplankton, herbivorous zoolplankton and copepods. A l l input model parameters are shown in Table 5.1, and diet compositions for each group are shown in Table 5.2. Following is a brief documentation of each group in the model and where Ecopath input parameters were extracted from. Lingcod The production parameters (P/B or Z) for l ingcod were taken from Cass et al . (1990), and Q/B ratios were modified from Venier and Kelson (1996). Lingcod were split into j u - venile and adult components, wi th 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 body weight data (results of stock assessment model are shown in appendix B ) . Dogfish A n estimated PjB ratio for dogfish was taken from the natural mortal i ty rate estimates provided i n Walters and Bonfi l (1999), assuming fishing mortality i n 1950 was negligi- ble. Assuming individuals die at a constant instantaneous rate Z , then the production- biomass ratio is equivalent to total mortali ty rate, regardless of metabolic rates (Al len 1971). The Q/B ratio was estimated from a meta-analysis relationship that predicts Q/B from aspect ratio (Froese and Pauly 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 in Chapter 4. P/B ratios were taken from total mortali ty 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 Trophic Biomass Prod/Biom Cons/Biom Ecotrophic level (t/km2) (/year) (/year) efficiency Lingeod 4.51 0.6 0.18 3.3 0.331 Juv Lingeod 4.24 0.15 0.48 3.8 0.636 Dogfish 3.89 1.134 0.04 1.8 0.595 Adult Pacific Cod 4 0.919 0.6 2.4 0.243 Juv. Pacific Cod 3.68 0.149 1.6 5.333 0.29 Demersal Fish 3.67 8.5 0.52 2.5 0.985 Adult Herring 3.61 6.1 0.67 4 0.794 Juv. Herring 3.54 3.622 1.23 11.06 0.43 Eulachon 3.34 0.211 0.47 4.3 0.8 Euphausiids 2.61 18.61 7.9 65 0.6 Pink Shrimp 2.04 0.383 0.6 12 0.95 Copepods 2.4 10.629 27 60 0.95 Herb. Zoolplankton 2 15.729 40 183.3 0.95 Phytoplankton 1 35.784 123 — 0.8 Detritus 1 — —- — 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 benthic fishes, such has Plueronectids and Bothids (primarily Turbot and Pacific sanddabs), as well as a few shorter lived Sebastes species. Many of the demersal fishes are 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 trawl fishing technology was relatively new at the time, and I assumed a fishing mortality 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 commercial fishery. 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. Abundance indices for juvenile and adult pools were taken from Schweigert (2000), and production parameters were estimated from total mortality rates in 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 and first nations fisheries do exist i n large river estuaries. Eulachon are anadromous, and spend part of their life cycle off the west coast of Van- couver island. Consistent observations of the 2- and 3-year age classes are made during research trawl surveys for shrimp. These data were used to estimate abundance of eula- chon i n the W C V I area using methods outlined in Hay et al. (1997). Eulachon diets were assumed to be similar to that of Pacific herring, however, I assumed a smaller fraction of euphausiids in the diet, and that 10% of the diet is "imported" from outside the system (from prey species not included in the model). Euphausiids Two species of euphausiids (Euphausia pacifica and Thysanoessa spinifera) are important components of the zooplankton community 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. Product ion and consumption parameters were taken from Robinson and Ware (1994). No euphausiid biomass estimates were available for 1950, therefore I allowed Ecopath to estimate biomass under the assumption that 60% of the mortali ty was being explained i n the ecosystem model (EE — 0.6). The biomass estimate from Ecopath is slightly less than the 22 t k m - 2 estimated by Robinson and Ware (1994) during the mid 1980's. 114 O a ex o o H a Q in uoi5(ireidioo2 'qraH spodadoQ diuuqg >[ii!J spnsni3i]dng uoip-B|ng Sutjjajr Atif Suujaj-f ^ n p v q S I J [13S.I9U.18Q POO aHP^d 3I n PV IISIJSOQ pooSuiq Anr pooSuiq a, O rH d d oo ° 0 0 in O CO . o o d o d rH rn o o a> d d d Ol CO H i—I 00 CO CN CN ^ o d d CO CO CO rH d d d d in LO rH rH d d d LO LO o d o o d d rt i—i CO LO LO CO o o LO o O LO H< o o o o q q LO CO d d d d d d d d CO LO LO LO o o rH o q q d d o d o d rH LO O LO CN O rH rH CO O CO O LO o o 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 00 , I H H H ° LO ° S 0 5 3 ° o o H co H o d ° d ° d ° d rH CN rH rH LO O O O O O O O O O O d> d d> czi d CN Ol <M CO rH rH 00 O CO CO O O O ci ci d CD ci -a o bO S J3 - I—) nH ,.S s O T3 M ^ a < -a o -o O o o a cc 'o «H ta '3 P H a3 « C H fa bO 0 ho ' H C a. S c o rt a a3 C _ O 3 rt 03 CO -3 cu S "a o >> -a H CU O K a, Q H : a. 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 statistical fits to 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~2 • year""1 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). 116 5.5 5.0 4.5 f 4.0 _ i ._ 3.5 3.0 2.5 2.0 1.5 1.0 O iingcodl ShiimpTraw) • iPogfishl Demersal FishJ nvmmsn Herring Fishery! I Pacific Cod I fPmkShrirtipl Copepods Uettitusj Phytoplankton L Figure 5.2: Food web and fisheries for west coast Vancouver Island ecosystem model. The size of the box is proportional to biomass of each group, and the wid th of the connecting line is proportional to the flow. Posit ion on the y-axis is the trophic level of the group. 5 . 2 . 4 T i m e S e r i e s D a t a f o r W e s t C o a s t V a n c o u v e r I s l a n d Time series data used in the W C V I ecosystem model consisted of absolute catches, to- tal mortali ty rates, relative abundance indices (either generated from stock assessment models or direct observations in commercial catch rates) and fishing mortal i ty rates used as time series forcing. Observed catch time series are shown in Figure 5.3. Al though commercial fisheries for lingcod started in the late 1800's, catch data for southwest Vancouver Island was not separated unt i l 1954, and I assumed catches were relatively stable between 1950 and 1954 (not shown in the figure). Commercial fisheries for pink shrimp on W C V I started in 1973, and landings peaked in 1975. Herring fisheries collapsed in 1968 and the fishery was closed unt i l 1971. Recent catches of Pacific cod have been substantially reduced because of new quota regulations implemented in 1994 and low stock size. 117 Shrimp Ct Shrimp Zt 1950 1960 1970 1980 1990 2000 1950 1960 1970 1980 1990 Lingeod Ct Herring Zt - i 1—I I I I - r — — " — i 1 1 1 1950 1960 1970 1980 1990 2000 1950 1960 1970 1980 1990 2000 Herring Ct Pacific Cod Wt 1950 1960 1970 1980 1990 2000 1950 1960 1970 1980 1990 2000 Year Pacific Cod Ct o i iii.,iinniiinninniiii.t,ii ni....:i o i 1 1—- 1 1 r- 1 1950 1960 1970 1980 1990 2000 Year Figure 5.3: Commercial landings (t k m - 2 ) for west coast Vancouver Island from 1950 to 2000 (left column). Tota l mortal i ty rate estimates for pink shrimp and adult Pacific herring. P i n k shrimp total mortal i ty rate estimates were calculated from combined fishing and natural mortal i ty rate estimates in chapter 2. Tota l mortali ty rates for Pacific herring were inferred from age composition data. Bo t tom panel is mean body weights for Pacific cod, where mean body weight is calculated using: Wt = Bt/Nt. Relative abundance information for pink shrimp consists of the combined biomass estimates from area swept trawl surveys in areas 124 and 125 (see Chapter 2). Pacific cod biomass estimates from Chapter 4, and lingeod biomass estimates were calculated using the same delay difference model (see Appendix B ) . Note that juvenile lingeod and Pacific cod are represented by the number of 2-year old fish, and adult indices are in biomass units. A d u l t herring data were obtained from Jake Schweigert ( D F O , pers. comm.) and abundance was estimated using a catch-at-age model. The euphausiid index (1967-1998) was taken from Robinson and Ware (1999), which was generated using a trophodynamic simulation model driven by oceanographic data (upwelling, sea-surface temperatures and 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 data for the ecosystem model. Shrimp 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 in Chapter 4. Euphausi id abundance (from 1967 to 1998) estimated from a trophodynamic model by Robinson and Ware (1999). solar radiation). A n eulachon index was generated using the by-catch information from the shrimp trawl surveys, and is a combined index for statistical areas 123,124 and 125. Abundance indices for eulachon were generated using spatial interpolation (bicubic spline) of area swept research trawl information. A l l indices were scaled to have a mean =1. 119 5.3 R e s u l t s 5.3.1 F i t t i n g E c o s i m to T i m e Ser ies D a t a Predicted catch time series from Ecosim were consistent wi th observed catch time series (Figure 5.5), despite not using the catch time series for lingcod and Pacific herring in the fitting procedure. A total of 30 spline point parameters, representing variability in primary productivity, were estimated using the non-linear search routine provided i n Ecosim. In al 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 ini t ia l 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 in the total likelihood, the predicted changes in herring and Pacific cod were inconsistent w i th observed data given the estimated sequence of primary 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 primary productivity, better fits were obtained by increasing the column 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 rapidly increase following the collapse in the mid 1980's. Demersal fish k^ were lowered to 0.1, and all remaining groups were held at the default value of 0.3. A general pattern is apparent for al l of the relative abundance data sets shown in F i g . 5.4; abundance was high in the 1970's followed by a collapse i n the early 1980's (a period of low productivi ty apparently associated wi th the 1982-83 E l Nino) and a 120 Pink Shrimp Landings lljllllllii Adult Herring Landings 1950 1960 1970 1980 1990 2000 Lingeod Landings 1950 1960 1970 1980 1990 2000 Adult Pacific Cod Landings i 1 1 i 1 1950 1960 1970 1980 1990 2000 1950 1960 1970 1980 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 land- ings from many fisheries have 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 predic- tions 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 in- creased 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, per- sistent 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 total mortali ty estimated for Pacific cod was a result of fishing, and Ecosim predictions correspond well wi th estimates of abundance from the single species model shown in Chapter 4. Declines in juvenile Pacific cod recruitment appear to be a result of declining spawner abundance (recruitment over-fishing), because predicted total mortali ty for juvenile Pacific cod is relatively stable during the decline phase. Ecosim predictions of herring biomass over time show a similar pattern to the results of statistical catch-at-age models. Al though Ecos im fails to capture the large spike in biomass during the 1970's following the collapse of the fishery, the predicted increase in biomass was sufficient enough to fit well wi th observed catch data (fishing mortalities from single species stock assessment models for Pacific herring were used as input data for Ecosim). Fishing mortalities during the 1960's were excessively high (0.7 down to 0.35) before the fishery was shut down. Dur ing the closure, predicted predation mortalities (for both juvenile and adult herring) remained relatively constant, preventing any rapid increases in abundance. The relative abundance time series suggest abundance increased as much as 2 times in comparison to predicted abundance during the mid 1970's . Ecosim predictions for eulachon and euphausiid do not fit the observed data well (in the case of eulachon) or other model predictions (in the case of euphausiids). Ecos im pre- dicts that eulachon predation mortali ty declined i n the 1980's, due to declines in predator abundances (lingcod and demersal fishes), however the eulachon stock continued to de- cline despite release from predation. Estimates of total mortali ty for eulachon actually increased due to bycatch in the shrimp trawl fishery, and declining primary production starting in 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, there- fore no associated fishing mortali ty rates. Simulated euphausiid production was closely connected to variation in primary productivity. Dur ing m i d 1970's simulated primary 122 Pink Shrimp ~t 1 r 1970 1990 Lingeod E CO 03 E o \ 1 1 1 1 r 1950 1970 1990 Juv Lingeod T 1 r 1970 1990 Eulachon i 1 1 1 1 i 1950 1970 1990 Juv Herring T 1 1 1 1 r 1950 1970 1990 Adult Herring i i i i i r 1950 1970 1990 Adult Pacific Cod Juv Pacific Cod 1950 1970 1990 Euphausiids \ 1 1 1 1 r 1950 1970 1990 Year Figure 5.6: Predicted biomass time series (line) from Ecosim and observed relative abun- dance time series (circles). Note that relative abundance time series are assumed pro- portional to true abundance, where Yt = qBt. production was at an al l t ime high, however Ecos im predicted euphausiid abundance did not increase, due to large increases in predation mortali ty by herring. Dur ing the last decade, there was good agreement between predicted euphausiid biomass and the esti- mated euphausiid biomass from Robinson and Ware (1999), where both models indicate a declining trend in primary productivity. Annua l primary production anomalies were calculated i n Ecosim by estimating 30 points for a spline function that spanned 1950-2000. These 30 parameters were estimated by fitting the Ecosim model predictions to time series data 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. Est imated mean primary production during the 1950-1977 period was 25% higher than mean primary production from 1978-2000 period (Figure 5.7). Est imated primary production reaches lows during the 1982-83 and 1997-98 E l Nino years. ~n 1 1 1 1 r - 1950 1960 1970 1980 1990 2000 Year Figure 5.7: Est imated primary production anomalies from 1950 to 2000; where a spline function was fit to 30 equally spaced data points. The use of the spline function helps filter out high frequency variation that occurs in fitting the model for years having only one or few times series wi th large anomalies. Using the time series shown i n F i g . 5.7 to drive primary production results in an overall 39% decline in the total sum of squares between predicted and observed time series data (Table 5.3). W i t h the exception of 5 data sets, including simulated variation in primary productivi ty resulted in 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 body weight for adult Pacific cod increased (4%) when including variation in primary production, and residual variation in shrimp mortali ty increased 28%. Pr ior to 1973, there is no information about trends in abundance for pink shrimp on 124 Table 5.3: Squared residuals for No M o d e l (or deviations from the mean of the data), a model wi th no trophic interactions (assume that predation mortali ty and consump- tion are constant), a model wi th trophic interactions where primary production is held constant over time, and a model that includes estimated primary production anoma- lies (see Figure 5.7). Positive reductions in the square residuals (right hand column) between trophic interactions only model and trophic interactions wi th variable primary productivi ty indicate a better 'fit ' to the observed data. Data Group No model No trophic interactions Trophic interac- tions with con- stant PP 'IVophic interac- tions and vari- able PP Reduction in SS Shrimp Bt 10.44 9.97 9.32 6.57 29% Shrimp Ct 71.04 16.43 9.01 7.64 15% Shrimp Z 1.84 0.60 0.84 1.07 -28% Lingcod Bt 0.67 0.82 1.17 3.38 -189% Lingcod Ct 15.36 10.61 11.94 12.99 -9% Juv Lincod 6.76 6.77 12.23 16.12 -32% Eulachon 24.07 21.04 22.73 18.28 20% Juv Herring Bt 34.07 87.35 39.12 11.06 72% Adult Herring Bt 19.77 62.85 26.47 8.39 68% Herring Ct 69.52 55.29 30.49 13.89 54% Herring Z 22.89 18.13 17.89 17.48 2% Pacific Cod Bt 20.27 12.50 8.97 2.33 74% Pacific Cod Wt 1.23 1.22 1.27 1.32 -4% Pacific Cod Ct 49.11 12.53 8.97 2.39 73% Juv. Pcod Bt 23.53 18.36 14.09 6.12 57% Euphausiids 7.06 7.06 6.35 6.16 3% Total 377.64 341.55 220.85 135.19 39% the west coast of Vancouver Island shelf. A s a result predicted biomass for shrimp is relatively stable during this time period. Just prior to the start of the fishery, Ecosim predicts an increase in abundance associated wi th increases in pr imary productivity. From 1975 to 2000, simulated trends in shrimp biomass are consistent wi th combined abundance estimates from Chapter 2. In comparison wi th Ecosim predictions, the single species model indicates more variabil i ty i n abundance, and the recent decline started later (1995) in comparison to the decline predicted by Ecosim (1989, see F i g . 5.6). Variations in upwelling, sea surface temperatures, and solar radiation strongly influ- ence production dynamics off the west coast of Vancouver Island (Robinson and Ware 1994), especially zoolplankton communities (Mackas et al. 2001). Ecosim was able to capture some of this variation by estimating relative trends in primary production using data only from higher trophic levels. N o oceanographic data or climate indices were used 125 to drive primary production in 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, Ecosim used data from upper trophic levels to infer or back calculate historical trends in apparent primary production. Recent downward trends in calculated relative primary production correspond well wi th estimated diatom produc- t ion from Robinson and Ware (1999) during the 1990's. Dur ing the mid 1970's Robinson and Ware estimate an increasing trend in diatom production, whereas Ecosim predicts a declining trend in primary production (see F ig . 5.8). Th is discrepancy occurs because abundance estimates for shrimp, herring, and Pacific cod were a l l declining during this time period, therefore the data support a declining relative production scenario. 1.5 o- 1.0 H 0.5 Primary Production Anomalies Diatom Production 1970 1975 1980 1985 1990 1995 Year Figure 5.8: Comparison between estimated annual mean primary production anomalies from Ecosim, and calculations of Dia tom production (g C m2yr~1) obtained by Robinson and Ware (1999) using a model based on upwelling (nitrogen), water temperature, and solar radiation. The relative primary production series estimated i n Ecosim also resembles the nega- tive pacific decadal oscillation (Figure 5.9). The P D O and estimated primary production are negatively correlated, thus during the negative phase of the P D O apparent primary production off the west coast of Vancouver Island is higher than average. 126 C o o i r o Q CL CM 1950 1960 1970 1980 1990 2000 Year Figure 5.9: Comparison between estimated annual mean primary production and the negative Pacific Decadal Oscil lat ion ( P D O ) . 5 . 3 . 2 C h a n g e s i n S h r i m p M o r t a l i t y During the 1990's, total predation mortali ty on pink shrimp has declined (Fig. 5.10) cor- responding wi th apparent declines i n predators. The pr incipal predators in the ecosystem model, juvenile and adult Pacific cod, have been declining since the mid 1970's. Dur ing the period 1950-70, mean predation mortali ty rates induced by Pacific cod predation were roughly 0 .51~ y r , or 78% of the total shrimp mortal i ty rate. Dur ing the 1990's the calculated predation rate dropped to 0 .22 _ 3 / r , or 30% of the total mortality rate. Other pink shrimp predators (lingeod and dogfish) rarely consume pink shrimp, and less than 20% of the total mortali ty was attributed by Ecosim to these predators. Lingeod and dogfish abundance have not declined nearly as much as Pacific cod abundance, and re- cently (1999-2000) dogfish predation has increased to 30% of the total shrimp mortality. I assumed during the 1950's that 95% of the total mortal i ty (Z) on pink shrimp populations was being explained by the predators included in the model (this is done 127 1.4 H Total Mortality Dogfish Adult Pacific Cod 1.2 H — Juv. Pacific Cod — Lingcod Juv. Lingcod — Fishing Mortality • Z from SSM 1.0 H 1 0.8 J 3 c 0.4 -A 0.2 H o.o H 1950 1960 1970 1980 1990 2000 Year Figure 5.10: A l l mortal i ty components from Ecosim compared to Z from single species model (SSM). The area between each line is the total mortal i ty due to each predator or fishing, and the area between total mortali ty and fishing mortali ty is the unexplained mortali ty (1 - EE). by setting the EEi term equal to 0.95). There are addit ional predators that were not included in the Ecopath model, and such addit ional predation is assumed to be constant over time in Ecosim. In contrast, the single species assessment model in Chapter 2 treats natural mortali ty ( 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 mortali ty for both Ecosim and single species models increased 7% and 21%, respectively between 1975-85 and 1990-00, and 2) there are two negative residuals that correspond wi th 1982-83 and 1997-98 E l Nino events (see F i g 5.11). Ecosim predicted an increase in total mortal i ty rates during the late 1980's and 1990's. This increase in total mortal i ty rate was a result of increased fishing mortality. Dur ing 128 i 1 1 1— 1975 1980 1985 1990 1995 2000 Year Figure 5.11: Standardized differences between Z estimates from Ecosim versus Z es- timates 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 dur- ing 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 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 differ- ences 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 in figure 5.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% (r2 value from Table 5.4) of the difference in shrimp mortality between the two models. • Z differences 130 5 . 4 D i s c u s s i o n 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 substan- tial 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 the findings of 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 mortali ty (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 rawl 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. The EE term in Ecopath is the assumed proportion of total mortali ty that is explained by trophic interactions and fishing mortality. The unexplained proportion (1 — EE) impl ic i t ly represents other predation losses, disease, or senescence. If predation mortali ty is additive (i.e. modelled predators are not consuming shrimp that are going to die anyway), then unexplained loss rate (Ma = 1 — EE) is assumed to be constant over time. In the ecosystem model I assumed that losses due to predation by hake are represented in the Ma term, so these losses were treated as a constant times shrimp biomass over time. Pacific hake were omitted from the Ecosystem model because time series data for the model area were incomplete and impossible to part i t ion from the coast wide hake stock. The migratory hake stock is large relative to prey groups specified in the model. A s discussed in Chapter 4 the spatial distribution of hake is highly variable, even wi th in the Canadian zone during summer months (Mackas et al . 1997), and it is extremely difficult to estimate the mean biomass in the model area. The influence of Pacific hake predation on the west coast of Vancouver Island ecosys- tem is significant for many forage species (Tanasichuk et al . 1991,Ware and McFar - lane 1995,Robinson and Ware 1999). Given the relative differences in hake and shrimp biomass, and the potential for high predation mortali ty rates, hake should be included 132 in the model. P r imary and secondary prey species for hake are euphausiids and herring (Ware and McFarlane 1995), which are also primary forage species for a l l large fish i n this ecosystem model. As sea surface temperatures increase, euphausiid production declines (Robinson and Ware 1994), which may result in increased predation on herring and per- haps shrimp by hake. A s noted in Chapter 4, very small increases of shrimp in the daily ration of Pacific hake could translate into very large changes in predation mortali ty rate. A s suggested by the weak correlation between the residuals in Z estimates between the two models and an estimate of hake abundance in the Canadian zone, including hake in the model may explain at least an additional 11% of the variabili ty in Z for shrimp. The results in Table 5.3 are comparisons of 3 different models i n which (1) trophic in - teractions are assumed constant, (2) trophic interactions were included, and (3) both trophic interactions and variabili ty in primary productivi ty are included. Assuming trophic interactions are constant over time is equivalent to using a single species model wi th constant natural mortal i ty rates. B y comparing a model wi th fixed natural mor- tali ty rate to a model wi th variable natural mortali ty 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 l ingcod and Pacific cod are determined by fishing mortality. In contrast, including trophic interactions explained a significant amount of variation in Pacific herring abundance and mortality, and variations in natural mortali ty are likely, as was determined by Schweigert (2000). Including trophic interactions for pink shrimp d id not result in better explanations of the time series Z data; however, including both trophic interactions and variation in primary productivity led to more of the variance being explained than could the fixed mortali ty 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 wi th pre- dictions generated in the single species model. The recent decline in Pacific cod explains some of the changes i n total mortality rates for shrimp; however, Ecosim was unable to capture the two large natural mortali ty rate anomalies that occurred during the 1982- 83 and 1997-98 E l Nino events. Judging by the differences in the estimates of Z from Ecosim and the single species models, and evidence of more hake in the Canadian zone during warm years, it appears that including Pacific hake in the model would explain these large mortali ty anomalies. The positive correlations in shrimp and Pacific cod abundance are best explained by variation in primary production. Est imated variation in primary production corresponds well wi th 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 perspec- tive. 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 additional fisheries and 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 rela- tionships. Such models should recognize strong trophic interactions that directly affect natural mortality rates, examine tradeoffs between different fisheries, and be able to gen- erate 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 manage- ment is estimating optimal harvest controls. In general, there are two key population parameters required for estimating opt imum harvest objectives: a measure of how large the stock is (such as the unfished biomass, or carrying capacity), and rates of compensa- tion, (such as the maximum intrinsic rate of growth, steepness of the stock-recruitment relationship, or parameter(s) that describe density dependent rate changes such as the Ecosim 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 Mar te l l 2002). Even more complex age-structured models wi th hundreds of nuisance parameters require these 2 key parameters for determining opt imal 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 population scale. If harvest rate objectives are met and population parameters are well determined, then using fixed rules for managing the fishery are quite robust, over the long term, to environmental variabil i ty (Walters and Parma 1996). However, fixed harvest rules may also lead to pathological over-fishing because of delayed responses (i.e. Walters and Ki tche l l 2001) or non-stationary effects, and here is where the ut i l i ty of ecosystem models can help identify these complex tradeoffs. In chapter 2, a stock assessment model was used to estimate average recruitment in each statistical area assuming no underlying stock-recruitment relationship. This was accomplished by directly estimating recruitment each year (an average recruitment level (Ra) and a series of recruitment anomalies (wt) 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 in 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 un- der the hypothesis that shrimp in each statistical area are separate stocks, I was unable to find a single correlate that explained a significant amount of variation for both statis- tical 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 re- cruitment, 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 pre- dictions also assumes we can predict future environmental factors (Walters and Collie 1988). 137 Variabi l i ty in natural mortali ty rates over time also increases uncertainty about fu- ture abundance. Using age-length composition information, I was able to estimate a time series of natural mortal i ty rates using the assessment model in chapter 2. B y track- ing cohorts over time, length frequency data provide compositional information about population size structure. Assuming growth and selectivity have remained constant over time, total mortali ty rates (Z) can be inferred by tracking changes in the length com- position data. Also , the addit ion 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 in natural mortality rates for pink shrimp off the west coast of Vancouver Island have been highly variable. Natura l mortal i ty rates for shrimp were estimated to increase during E l Nino years and the mechanism for increases in 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. In- direct measures (correlations between predator biomass and prey abundance), as well as the small amount of information about hake distr ibution relative to sea surfaces height anomalies, does suggest that Pacific hake populations might play a significant role in pink shrimp dynamics. These findings are also consistent wi th findings about pink shrimp-hake interactions i n Oregon and Washington (Hannah 1995). Predict ing changes in natural mortal i ty 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 per- spective. Using time series data from several different fisheries and assessment programs, I was able to fit a dynamic ecosystem model that included estimating a time series of rela- tive primary production anomalies. The goal of this chapter was to determine how much 138 of the variation in production arid mortali ty could be explained by including trophic in- teractions and variability in primary production i n the assessment process. In summary, including trophic interactions while assuming constant production degraded statistical fits to the data; however, including variable pr imary production explained an additional 29% and 15% of the variabili ty in relative shrimp abundance and observed shrimp catch time series, respectively (see Table 5.3 on page 124). Nearly all of the time series data included i n the ecosystem assessment were consistent wi th the hypothesis that primary productivity has declined over the last 3 decades. Coincidentally, long term trends in primary production (estimated in Ecosim by fitting model predictions to time series data from fish populations) shifts downward around 1977, which corresponds well w i th regimes shifts documented by Beamish et al. (1999). Bibliography Al len , K . R. (1971). Relat ion between production and biomass. J. Fish. Res. Bd. Canada 28, 1573-1581. Anderson, P. J . and J . F . P ia t t (1999). Communi ty reorganization i n the Gu l f of Alaska 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. Nauchn. Issled. Ikhtiologicheskii Inst. Izv. 1, 81-128. Baretta, J . , W . Ebenhoh, and P. Ruardij (1995). The European 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 . McFarlane, L . Klyashtor in , V . V . Ivanov, and V . Kurashov (1999). The regime concept and natural trends in the production of salmon. Can. J. Fish. Aquat. Sci. 56, 516-526. Berenboim, B . , A . Dolgov, V . Korzhev, and N . Yaragina (2000). The impact of cod on the dynamics of Barents Sea shrimp (Pandalus borealis) as determined by mult i - species models. J. Northw. Atl. Fish. Sci. 27, 69-75. Bouti l l ier , J . , I. Perry, B . Waddel l , and J . Bond (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). Fishing on ecosystems: the interplay of fishing and predation in Newfoundland-Labrador. Can. J. Fish. Aquat. Sci. 58, 1153-1167. Butler , T . (1980). Shirmps of the Pacific Coast of Canada, Volume 202. Can . B u l l . F i sh . Aquat . Sci . Caput i , N . , J . W . Penn, L . M . Jol l , and C . F . Chubb (1998). Stock-recruitment- environment relationships for invertebrate species of Western Austra l ia . Can. J. Fish. Aquat. Sci. 125, 247-255. Cass, A . J . , R . J . Beamish, and G . A . McFarlane (1990). Lingcod (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 in fisheries. Can. J. Fish. Aquat. Sci. 56, 1525-1533. Dalsgaard, J . , S. S. Wallace, S. Salas, and D . Preikshot (1998). Mass-balance mdoel reconstructions of the Strait of Georgia: the present, one hundred, and five hundred years ago. In D . Pauly, T . Pitcher, D . Preikshot, and J . Hearne (Eds.), Back to the future: reconstructing the Strait of Georgia ecosystem., Volume 6, pp. 72-91. Fisheries Centre Research Reports. 139 140 Dennis, J . E . and R. B . Schnabel (1983). Numerical Methods for Unconstrained Opti- mization and Nonlinear Equations. Prentice-Hall , Englewood Cliffs, N J . Deriso, R . B . (f980). Harvesting strategies and parameter estimation for an age- structured model. Can. J. Fish. Aquat. Sci. 37, 269-282. Dorn , M . W . , M . W . Saunders, C . D . Wi l son , M . A . Guttormsen, K . Cooke, R . Kieser, and M . E . Wi lk in s (1999). Status of the coastal Pacific hake/whit ing stock i n U .S . and Canada in 1998. Canadian Stock Assessment Secretariat Research Docu- ment 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, 1568- 1573. Fargo, J . and P. J . Starr (2001). Turbot stock assessment for 2001 and recommenda- tions for management in 2002. PSARC Working Paper G2001-10. Foreman, M . G . G . and R. E . Thomson (f997). Three-Dimensional M o d e l Simulations of Tides and Bouyancy 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 . Majkowski , and J . Hampton (1990). Mul t i fan a likelihood-based method for estimating growth parameters and age composition from multiple length frequency data sets il lustrated using data for Southern bluefin tuna (Thunnus maccoyii). Can. J. Fish. Aquat. Sci. 47, 301-317. Fournier, D . , J . Sibert, and M . Terceiro (1991). Analysis of Length Frequency Samples wi th Relative Abundance Da ta for the G u l f of Maine Northern Shrimp (Pandalus borealis) by M U L T I F A N Method . Can. J. Fish. Aquat. Sci. 48, 591-598. Francis, R . C . and S. R . Hare (1994). Decadal-scale regime shifts in the large marine ecosystems of the North-east Pacific: a case for historical science. Fish. Oceanogr. 3, 279-291. Froese, R . and D . Pauly (2002). FishBase. World Wide Web electronic publication. www.fishbase.org, M a y 22,2002. Fu , C . and T . J . Qu inn II (2000). Es t imabi l i ty of natural mortali ty and other popu- lation parameters in a length-based model: Pandalus borealis in Kachemak Bay, Alaska. Can. J. Fish. Aquat. Sci 57, 2420-2432. Gelman, A . , J . Car l in , H . Stern, and D . R u b i n (1995). Bayesian Data Analysis. Chap- man & Ha l l . Gi l l i s , D . M . , R. M . Peterman, and A . V . Tyler (1993). Movement dynamics i n a fishery: Appl ica t ion of the ideal free distr ibution to spatial allocation of effort. Can. J. Fish. Aquat. Sci. 50, 323-333. Gotshal l , D . W . (1969). The use of predator food habits in estimating relative abun- dance of the ocean shrimp, Pandalus jordani, Rathbun. FAO Fish. Rep. 57, 667- 685. Haist, V . and D . Fournier (1995). Pacific cod stock assessment for f995 and recom- mended yield options for f 996. PSARC Working Paper G95-3. Hannah, R . W . , S. A . Jones, and M . R . Long (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 Cod Fishery: A True Crime Story. McClelland & Stewart. 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 Fisheries Stock Assessment: Choice, Dynamics, & Uncertainty. 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 manage- ment 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. 5Jh 2080-2096. 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 com- munity of the British Columbia continental margin, 1985-1999, and their covaria- tion with oceanographic conditions. Can. J. Fish. Aquat.Sci. 58, 685-702. 142 Mantua , N . J . , H . S. R. , Y . Zhang, W . J . M . , and R . C . Francis (1997). A Pacific interdecadal climate oscillation wi th impacts of salmon production. Bull. Amer. Meteorol. Soc. 78, 1069-1079. Mar te l l , S., J . Bouti l l ier , H . Nguyen, and C . Walters (2000). Reconstructing the off- shore Pandalus jordani trawl fishery off the west coast of Vancouver Island and simulating alternative management policies. Canadian Stock Assessment Secre- tariat 2000/149, 1-38. McGowan , J . A . , D . R . Cayan, and L . M . Dorman (1998). Climate-Ocean Variabi l i ty and Ecosystem Response in the Northeast Pacific. Science 10-Jul 281, 210-217. Myers, R . A . , K . G . Bowen, and N . J . Barrowman (1999). M a x i m u m reproductive rate of fish at low populations sizes. Can. J. Fish. Aquat. Sci. 56, 2404-2419. Mysak, L . (1986). E l Nino, interannual variabili ty and fisheries in the northeast Pacific Ocean. Can. J. Fish. Aquat. Sci. 43, 464-497. Paloheimo, J . and L . Dickie (1964). Abundance and fishing success. Rapp. P. -V. Reun Cons. Int. Explor. Mer. 155, 152-163. Parrish, R. , C . Neson, and A . Bakun (1981). Transport mechanismsand reproductive successof fishes in the California Current. Biol. Oceanogr. 1, 175-203. Pauly, D . , V . Christensen, and C . Walters (2000). "Ecopath , 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 . Flannery (1996). Numerical Recipies in Fortran 77, The Art of Scientific Computing (2nd ed.). Cambridge Universi ty Press. Quinn, T . , C . Turnbul l , and C . F u (1998). A Length-Based Populat ion Mode l for Hard- to-Age Invertebrate Populations. In F . Funk, T . J . Quinn , J . Heifetz, J . N . Ianelli, J . E . Powers, J . F . Schweigert, P. J . Sull ivan, and C . I. Zhang (Eds.), Fishery Stock Assessment Models, pp. 531-556. Alaska Sea Grant College Program Report No. A K - S G - 9 8 - 0 1 , University of Alaska Fairbanks, 1998. Rexstad, E . A . and E . P ik i t ch (1986). Stomach contents and food consumption esti- mates of pacific hake, Merluccius productus. Fish. Bull. 84, 947-956. Rina ld i , S. and M . Gat to (1978). O n the determination of a commercial fishery pro- duction model. Ecological Modelling 8, 165-172. Robinson, C . L . K . and D . M . Ware (1994). Model l ing pelagic fish and plankton trophodynamics off Southwestern Vancouver Island, Br i t i sh Columbia . 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 in the 1990s. Can. J. Fish. Aquat. Sci. 56, 2433-2443. Rothlisberg, P. C . and C . B . Mi l le r (1983). Factors affecting the distribution, abun- dance, 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 minimizat ion. A CM Trans. Math. Software 11, 419-440. Schnute, J . and D . Fournier (1980). A new approach to length-frequency analysis: Growth Structure. Can. J. Fish. Aquat. Sci. 37, 1337-1351. 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. Candadian Stock Assessment Secre- tariat 2000/165, 74. 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 pa- rameters 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 . Hil l , 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 con- text. 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 used in stock assessment models. Can. Spec. Publ. Fish. Aquat. Sci. 108. Tsou, T.-S. and R. M . Royall (1995). Robust Likelihoods. J. Amer. Stat. As- soc. 50(429), 316-320. 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 ex- ploited 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 sur- vival 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 distri- butions 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 Den- sity 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 Van- couver 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 abun- dance, 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. Ecol- ogy- Appendix A Simulation Study for Shrimp Assessment Mode l Traditional stock assessment models usually assume natural mortality rates remain rela- tively 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 as- sessment 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 parame- ters, 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 rel- ative 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 population. To determine if the assessment model was appropriate for estimating changes in natural mortali ty rates over time, I used the assessment model described in Chapter 2 to generate 21-year time se- ries of relative abundance data, relative fishing mortalities, and length frequency data, then estimated population parameters using the assessment model wi th simulated data. Relative abundance data were assumed to have log-normal errors, and length frequency data were taken by randomly sampling from the true mul t inomial distribution in each year. A n example of a length frequency sample over time is shown i n Figure A . l . I i i — i — r 10 15 20 25 30 35 40 2 I T 1 1—̂ 10 15 20 25 30 35 40 3 10 15 20 25 30 35 40 4 5 ° J l — " P H — i — r 5 * ! — i U- 10 15 20 25 30 35 40 5 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 I—l^h—I—I I I 10 15 20 25 30 35 40 — i — i — i — i i 10 15 20 25 30 35 40 10 10 15 20 25 30 35 40 11 l i i — i — r ~ l — i 10 15 20 25 30 35 40 12 10 15 20 25 30 35 40 17 i — I — I — T I i 10 15 20 25 30 35 40 18 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 Length (mm) r ~ n — I — I I I i 10 15 20 25 30 35 40 Figure A . l : Example of an observed length frequency dis tr ibut ion and predicted frequen- cies from stock assessment model. Observed frequencies were drawn at random from the true mult inomial distribution. If the model is working correctly, it should be able to accurately fit simulated data- sets that have no observation errors, and correctly estimate recruitment anomalies and anomalies in natural mortali ty rates (nuisance parameters) over time. This 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 (Ra) are over-estimated when the sum of recruitment anomalies can assume 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) Time (years) <= - i , 5 - | I I I I— 1 1 1 1— 5 10 15 20 5 10 15 20 Time (years) 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 Mt 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 data wi th no observation errors, and estimated values wi th out penalty functions on recruitment anomalies, and natural mortal i ty rate. Parameter True Value Estimate w / o penalty Estimate w penalty RQ 33.12 38.36 26.00 M 0.7 0.699 0.698 u 15 14.99 15.03 30 29.99 30.01 P 0.68 0.679 0.685 CV 0.08 0.079 0.08 7 0.8 0.803 0.775 LH 18 17.998 18.114 that the exact parameter values wi l l be estimated correctly given a small time series. However, over an infinitely long time scale, or multiple simulated data sets, on average the estimation procedure should reflect the true parameter values (i.e. unbiased). To ensure the assessment model used in Chapter 2 was unbiased, 250 independent data sets were simulated, generating 250 estimates of each model parameter. I then plotted the distr ibution of al l estimated parameter to look for systematic bias (Figure A .3 ) . W h e n length frequency samples are small and uninformative about changes in mortal- ity rates, auxil iary information about relative trends in fishing mortali ty rates generally improves estimates of true population parameters, especially changes i n natural mortality. For the simulations shown in Figure A . 3 , annual samples of 10,000 were taken and were very informative about changes in age structure over time. Smaller samples, or unrepre- sentative samples, led to biased parameter estimates. A disproportionate representation of smaller length classes wi l l over-estimate annual recruitment and under-estimate natu- ral mortality, and vice versa when larger shrimp are over represented in the sample. A s discussed in 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 LA = L\ + (LN — L\) lZPN-\, 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 parame- ter RQ appears to be positively correlated with mean natural mortality rates, indicating 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 vt that were constrained to (VJ t vt — 0). 150 0.68 0.72 29.95 30.10 0.0795 Ro m M •r- I I IT "II I 29 33 17.8 18.6 I I I I, an « n in fifl i L1 Ln *1 I ' rho • • X • — • / : ML H • S I CV1 • • gamma • — • _ • • Ih 14.95 0.67 0.70 I I I I I 0.65 0.80 Figure A .4 : Parameter correlations for 250 simulation-reconstructions using the catch- at-length model described in Chapter II. Variables names defined i n Table A . l . Appendix B Lingeod Assessment Commercial catch statistics were used to reconstruct biomass and fishing history for ling- cod 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 pa- rameter. Furthermore, parameters were estimated in two phases, in which the first phase included a penalty on the mean exploitation rate (0.25) to approximate correct scaling parameters (BQ). In the second phase, this penalty was relaxed and nuisance parameters (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, fishing mortality rates, and age-2 recruits are shown in Figure B.2. Intensive fishing for lingeod on the west coast of Vancouver Island began during the development of the groundfish trawl fishery in early 1940's (Cass et al. 1990). Prior to the development of the trawl fishery lingeod were fished quite 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 (Et) 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 153 Table B.2 : Est imated growth parameters from length and weight-at-age data from Cass et al,(1990). Descr ip t ion Parameter Value Asymptot ic length (cm) Loo 106.3 Growth coefficient k 0.212 Length-weight scalar a 1.12E-05 Length-weight power b 3.002 unfished state. The ratio of 1950's biomass relative to its unfished state was estimated at 29% in this assessment. From 1950 to 1969, biomass increased from as 9,500 tonnes to « 12,500 tonnes. The most recent biomass estimate is roughly 9,100 tonnes. The retrospective projections show that more recent data has resulted in lower estimates of unfished biomass and productivity (Table B . l ) . Th is is likely due to the additional contrast that has been provided by recent declines in stock abundance. 1 1 1 1 r 1960 1970 1980 1990 2000 Year Figure B . l : Retrospective projection of l ingcod biomass. Fishing mortali ty rates during the 1950-2000 period ranged from 0.01 in 1995 to 0.18 in 1996, and averaged roughly 0.07 during this period (Figure B . l ) . Note that this is an adjustment downward after relaxing the penalty on mean fishing mortali ty 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 in the 154 Figure B .2 : Reconstructed biomass, observed landings, fishing mortali ty rates, resid- uals in catches and mean body weights, and estimates of age-2 recruits for southwest Vancouver Island lingcod. catchability coefficient over time to reflect changes in fishing technology that make the capture process more efficient. Dur ing 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 in 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 repre- sented in Ecosim by: ^ _ cij,t(Bj,t> Bj,t) ^c 155 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 vtj and use this as an estimate. However, and i»y are inseparable in consumption equations, and 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(DCijit\Biit, Bjit, ciij^ij) = rij * DCij>tlog (DCijit) where ]P DCij>t = 1 (C.2) 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^.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Variation in pink shrimp populations off the west coast...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Variation in pink shrimp populations off the west coast of Vancouver Island : oceanographic and trophic.. Martell, Steven James Dean 2002
pdf
Page Metadata
Item Metadata
Title | Variation in pink shrimp populations off the west coast of Vancouver Island : oceanographic and trophic interactions |
Creator |
Martell, Steven James Dean |
Date | 2002 |
Date Created | 2009-11-12 |
Date Issued | 2009-11-12 |
Description | 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, and 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. |
Extent | 10692002 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-11-12 |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0074868 |
Degree |
Doctor of Philosophy - PhD |
Program |
Zoology |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 2002-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/14877 |
Aggregated Source Repository | DSpace |
Download
- Media
- ubc_2003-792390.pdf [ 10.2MB ]
- [if-you-see-this-DO-NOT-CLICK]
- 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
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 15 | 1 |
France | 6 | 0 |
Canada | 5 | 0 |
Netherlands | 4 | 0 |
Russia | 2 | 0 |
China | 2 | 20 |
Germany | 1 | 2 |
City | Views | Downloads |
---|---|---|
Unknown | 14 | 3 |
Ashburn | 5 | 0 |
Amsterdam | 4 | 0 |
Wilkes Barre | 3 | 0 |
Wilmington | 3 | 0 |
Beijing | 2 | 0 |
Vancouver | 2 | 0 |
Mountain View | 1 | 0 |
Atlanta | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
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>
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
Admin Tools
To re-ingest this item use button below, on average re-ingesting will take 5 minutes per item.
ReingestTo clear this item from the cache, please use the button below;
Clear Item cache