Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

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

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

Item Metadata

Download

Media
831-ubc_2003-792390.pdf [ 10.2MB ]
Metadata
JSON: 831-1.0074868.json
JSON-LD: 831-1.0074868-ld.json
RDF/XML (Pretty): 831-1.0074868-rdf.xml
RDF/JSON: 831-1.0074868-rdf.json
Turtle: 831-1.0074868-turtle.txt
N-Triples: 831-1.0074868-rdf-ntriples.txt
Original Record: 831-1.0074868-source.json
Full Text
831-1.0074868-fulltext.txt
Citation
831-1.0074868.ris

Full Text

VARIATION IN PINK SHRIMP POPULATIONS OFF THE WEST COAST OF VANCOUVER ISLAND: OCEANOGRAPHIC AND TROPHIC INTERACTIONS by STEVEN JAMES DEAN MARTELL B.Sc, University of British Columbia, 1997 M.Sc, University of British Columbia, 1999 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of Zoology) We accept this thesis as conforming to the required standard UNIVERSITY OF BRITISH COLUMBIA 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 iiList 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 12.2 Stock Reconstruction Model 2 2.2.1 Population dynamics model 13 2.2.2 Observation Models 4 2.2.3 Parameter Estimation aud Uncertainty 23 2.3 Results 26 2.3.1 Estimating Abundance 22.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 6iii iv 3.1.1 Summer and Winter Current. Circulation Patterns 61 3.1.2 Distribution of Shrimp Relative to Coastal Currents and Shelf Break Currents 62 3.2 Methods 64 3.2.1 Environmental Indices 7 3.3 Results 9 3.4 Discussion 72 4 Shrimp Predation by Gadids 77 4.1 Introduction4.1.1 Stock Structure of Pacific Cod 78 4.1.2 Stock Structure of Pacific Hake 80 4.2 Estimating Predator Biomass 83 4.2.1 Pacific Cod Assessment4.2.2 Pacific Hake Assessment 88 4.3 Results 90 4.3.1 Pacific Cod Biomass 94.3.2 Correlation with Predator Abundance 93 4.4 Discussion 95 5 Ecosystem Perspectives 99 5.1 Introduction5.2 Methods 101 5.2.1 Ecopath with Ecositn 105.2.2 Fitting Ecosim To Time Series Data 107 5.2.3 Description of Model Groups 110 5.2.4 Time Series Data for West Coast Vancouver Island 11G 5.3 Results 119 5.3.1 Fitting Ecosim to Time Series Data 115.3.2 Changes in Shrimp Mortality 126 5.4 Discussion 130 6 General Discussion and Summary 134 Bibliography 139 A Simulation Study for Shrimp Assessment Model 145 B Lingcod Assessment 151 Using changes in diet composition to estimate vulnerability List of Tables 2.1 Model parameter descriptions and symbols used in population dynamics model 14 2.2 Parameter descriptions and symbols used in length frequency observation model 7 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 20 2.4 Estimated growth parameters for area 124 and 125 pink shrimp populations. 32 2.5 Estimated model parameters and standard deviations for statistical areas 124 and 125 32.6 Estimated recruitment anomalies («>,) and standard deviations for statis tical areas 124 and 125 33 2.7 Estimated deviations in mortality (;/t) and standard deviations for statis tical areas 124 and 125 4 2.8 Reported fishing effort in commercial logbooks for statistical areas 124 and 125 42.9 Natural mortality rate statistics for areas 124 and 125 49 3.1 Comparing the effects of environmental indices on stock recruitment rela tionships for areas 124-5. loijL — log-likelihood, \2 = probability of the of the environmental index relative to base model, and the % of the variance in the residuals explained by the enviromental correlate 71 3.2 Cross correlations between environmental indices shown in Figure 3.3 . . 71 3.3 Correlations between lo<i(H./S), fjog(R) and environmental indices shown in Fig 3.3 72 vi vii 3.4 Comparing the effects of environmental indices on stock recruitment rela tionships for area. 124 spa.wners producing area. 125 recruits and vice versa,. logL = log-likelihood, \2 — probability of the of the environmental index relative to base model, and the % of the variance in the residuals explained by the enviroment.al correlate 72 4.1 Published reports of Pacific hake (Merluccius producf;us) predation on pink shrimp (Pandalus jordnni) and reported % occurrence, or % by weight of pink shrimp in hake stomachs 83 4.2 Parameter symbols and description of parameters for the delay-difference model 84 4.3 Biomass Estimates for Pacific hake; for La Peruse bank region (Ware and McFarlane, 1995), and for hake in the Canadian Zone, (Dorn et al, 1999). 89 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 94 5.1 Ecopath input parameters, arid estimated parameters (in bold) for west coast Vancouver Island model If2 5.2 Diet matrix for Ecopath model 114 5.3 Squared residuals for No Model (or deviations from the mean of the data), a model with no trophic' interactions (assume that predation mortality and consumption are constant), a, model with trophic interactions where pri mary production is held constant over time, and a, model that includes esti mated primary production anomalies (see Figure 5.7). Positive reductions in the square residuals (right hand column) between trophic interactions only model and trophic interactions with variable primary productivity indicate a better 'fit' to the observed data 124 5.4 Correlation 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, with no observation errors, and estimated values with out penalty functions on recruitment anomalies, and natural mortality rate 148 viii B.l Commercial catch statistics for lingcod off southwest Vancouver Island (Source: King and Surry, 2000) 152 B.2 Estimated growth parameters from length and weight-at-age data from Cass et al,(1990) 153 List of Figures 1.1 Statistical area definitions used for smooth pink shrimp trawl fisheries. . 4 1.2 Pink shrimp Pandalus jordani landings in statistical 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 Location 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 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 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 Population age structure for pink shrimp in areas 124 and area 125. White 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. Estimated 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 30 ix X 2.5 Predicted (solid line) and observed (bars) length frequencies for area 125. Estimated 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 Pair plot of 500 random sample estimates (/?„, L\. LA. p, At, A2,'>, h) taken from an MCMC sample of length 250.000 for area 124 35 2.7 Pair plot of 500 random sample estimates (/?„, Lj, 1.4, p, Ai, A2,7, lh) taken from an MCMC sample of length 250,000 for area 125 36 2.8 Marginal posterior distributions for A2 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 37 2.9 Sensitivity of M to <r'f[ and a"t,. Left column represent the mean estimate of natural mortality as af, increases from 0.1 to 1, where the mean (line with 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 all 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 Sensitivity 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 with points) is calculated over all a\f. Right column represents the mean estimate of R0 as cr2,, increases from 0.8 to 8, where the mean is calculated over all 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. Area 125 clearly shows 3 local maxima 40 2.12 Estimated area over which shrimp populations in statistical areas 124-125 were distributed. Area, estimates were calculated using commercial logbook data, from 1987 to 2000; these data are shown in Figure 2.14 . . 41 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 (D, = N„.jWa). Note that abundance shown in figure 2.1 represents the vulnerable biomass 42 2.14 Spatial distribution of commercial fishing effort from 1987 to 2000 on the West Coast of Vancouver Island. Total number of tows from 1987 to 2000 was 59,9(52 2.15 Estimated selectivity curves for areas 124 and 125, where dotted vertical line represents the length at 50% vulnerability 2.16 Estimated fishing mortality rates from stock assessment model (lines in bottom panels) for areas 124-125. Circles represent area swept estimates of fishing mortality rates (F = ,J~-) from commercial log book data, and these data were used in the likelihood criterion for estimating 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 Maximum likelihood estimates of natural mortality rates (line with 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 distribution generated by the Metropolis Hastings Algorithm, and the box plots represent the median, inner quartile, and sampling range 2.18 Maximum likelihood estimate of age-f 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 million shrimp, respectively 2.19 Relationship between carapace length arid 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 If2,732 shrimp measured and sexed was chosen for estimating parameters for a logistic curve. Parameters were estimated us ing a least squares criterion, where the maximum likelihood estimate for 7 = 0.89 and lh = 19.71 2.20 Spawner biomass (t) at t — 2 years versus age-1 recruits for statistical 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 xii 3.1 Tow positions of research trawl surveys from 1973 to 2001. Statistical 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 Deviations in annual recruits per unit of spawniug biomass (log(R, jS\) -log(Rt/St)) in statistical 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 63.3 Standardized oceanographie 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 68 3.4 Beverton-Holt and Ricker models fit to stock-recruitment data from sta tistical area 124 using a. normal likelihood (Eq. 3.5) and robust likelihood (Eq. 3.7) Each model line represents the mean stock-recruitment relation ship with 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 with correlate included. Left panels are Beverton-Holt models, right are Ricker models 73 3.6 Survival index (log(Rt/St)) for age-1 recruits in area 124 and area 125. Dotted line represents mean log(R,/S) 75 4.1 Groundfish trawl catch per unit effort of Pacific 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. Map courtesy of Alan Sinclair, DFO Pacific Region 79 4.2 Distribution of Pacific hake in 1992, 1995, and 1998 along survey tran sect lines. Abundance is measured using hydro acoustic backscattering techniques. (Figure supplied by Mark Saunders, DFO) 82 xiii 4.3 Proportion of Pacific hake in the Canadian zone versus mean January-February sea level height anomalies (top panel). Positive anomalies indi cate higher than average1 sea level heights in Bamfield, British Columbia (data provided by Mark Saunders, UFO). Estimated 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 Estimated 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. Bottom row from left to right are recruitment anomalies, catch residuals, and mean weight residuals from the delay difference model 91 4.5 Likelihood profiles for WCVI 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 Relative abundance time series of Pink shrimp, Pacific Cod, and Pacific hake series I and II from 1975 to 2000. Lower panels are scatter plots of predator abundance versus pink shrimp abundance 95 5.1 An example of the relationship between predation mortality rate and predator biomass for different values of klr The dark circle represents the base predation mortality 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 proportional to biomass of each group, and the width of the connecting line is proportional to the How. Position on the y-axis is the trophic level of the group 116 •5.3 Commercial landings (t km""2) for west coast Vancouver Island from 1950 to 2000 (left column). Total mortality rate estimates for pink shrimp and adult Pacific herring. Pink shrimp total mortality rate estimates were calculated from combined fishing and natural mortality rate estimates in chapter 2. Total mortality rates for Pacific herring were inferred from age composition data,. Bottom panel is mean body weights for Pacific 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. Shrimp and eulachon biomass estimates are combined for two statistical areas, herring abundance was constructed using a. catch-at-age model, lingcod abundance from stock assessment model fit to commercial CPUE data, and Pacific cod from assessment in Chapter 4. Euphausiid 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 shrimp, herring, lingcod 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. 120 5.6 Predicted biomass time series (line) from Ecosim and observed relative abundance time series (circles). Note that relative abundance time series are assumed proportional to true abundance, where Yt = qBt 122 5.7 Estimated primary 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 variation that occurs in fitting the model for years having only one or few times series with large anomalies. 123 5.8 Comparison between estimated annual mean primary production anoma lies from Ecosim, and calculations of Diatom production (g C m2yr~x) obtained by Robinson and Ware (1999) using a, model based on upwelling (nitrogen), water temperature, arid solar radiation 125 5.9 Comparison between estimated annual mean primary production and the negative Pacific Decadal Oscillation (PDO) 126 5.10 All mortality components from Ecosim compared to Z from single species model (SSM). The area between each line is the total mortality due to each predator or fishing, and the area between total mortality and fishing mortality is the unexplained mortality (1 — EE) 127 5.11 Standardized differences between Z estimates from Ecosim versus Z es timates from single specie's model in Chapter 2. Positive values indicate years where Ecosim predicted higher mortality rates than estimated using the SSM 128 5.12 Time 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 distribution and predicted fre quencies from stock assessment model. Observed frequencies were drawn at random from the true multinomial distribution 1 46 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 animal natural mortality rate 147 A.3 Box plots of parameter estimates from 250 simulation-reconstructions us ing 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. 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 mortality rates, residu als in catches and mean body weights, and estimates of age-2 recruits for southwest Vancouver Island lingeod 154 Acknowledgements No thanks would be enough for all the encouragement and 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 British 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 will reach far beyond your lifetime. Many thanks to my committee members. Daniel Pauly and Villy Christensen who helped me develop alternative ideas about ecosystem modeling, as well as the occasional fishing trip for harbor seals and sight seeing tours at Florida State University (thanks Villy). John Dower and Dan Ware, for exposing me to a complex area of oceanography and climate affects. Alan Sine lair for helping me digest the Pacific cod history on the west coast of Vancouver Island, and great pasta dinners. Also, I would like to acknowledge the alternative exposure provided by Ransom Myers, Shelton Harley, Jeff Hutchings, arid Boris Worm (Dalhousie University). This 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 xvi xvii finally explained the black magic behind the Metropolis-Hastings algorithm, and kept me laughing during PS ARC meetings. Rik Thompson, who sent all of his reprints on the same flay I requested 'one'. Dave Founder, for getting me up to speed on the latest in non-linear parameter estimation software, and a bad hangover in Madison Wisconsin. James Kitchell, for the wonderful weekends surfing in Honolulu, poi, bass fishing in UNDER.C, Superica in Santa Barbra and the next job. Thanks also to Rob Bison, for some extra $$ and an opportunity to give back to the local steelhead community. Many thanks to my supervisor, Carl Walters, the man who defines conservation as "kill something else in order to save it", aptly encourages graduate students with sound advice such as "don't F9<?KUP", arid has his own icon for evolutionary ecology. In addition to witnessing the return of 'Betsy', I've also witnessed many other great mystery's with Carl, such as full moon's in broad day light and belly-button lint. I also have to thank Carl for the many great, adventures, here in Vancouver, Florida, .Japan, Honolulu and the Grand Canyon. Carl and his wife Sandra were kind in providing accommodations for Dawn and myself, until Willy moved back home and split ice tea on Sandra's new carpet (sorry Will, I had to tell someone!). I am forever grateful for all the love and support provided by my girlfriend Dawn Cooper. Thank you so much for all your patience and understanding. Special thanks to my lab mates, Sean Cox, Robert Alliens, Bob Lessard ,Paul Higgins, Rik Duckworth, Lisa Thompson, Leonardo Hua.to, Nathan Taylor, and David O'Brien. All 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 Jim Boutillicr who encouraged me to xviii adopt my own methods tor examining shrimp fisheries, and provided access to data and computing resources. Also Hai Nguyen and Norm Olsen for assistance with databases, and comprehending the data. Rob Kronlund and Marc Saunders for important discus sions about the problems with assessing Pacific hake abundance. Finally, the to crew of the CCV W.E. Ricker, the platform for shrimp surveys and the summer feeding grounds for myself. Finacial support for this project was provided by Department of Fisheries and Oceans, administered through Jim Boutillier. 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 (WCVI) shelf. During 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 (DFO). The WCVI shrimp trawl fisheries is one of a few fisheries in British 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 WCVI fail to respond to efforts to prevent over-fishing. On two occasions populations continued to decline in the absence of fishing. Historically, shrimp recruitment to WCVI has been highly variable, and unpredictable. This thesis examines the history of shrimp population dynamics off the WCVI and ex amines environmental and trophic interactions that lead to variability in recruitment and natural mortality. In this chapter I give a brief overview of the history of shrimp trawl fisheries in the WCVI 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. Natural mortality rates in pink shrimp populations tend to be highly variable, and I examine the variability in natural mortality 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 with 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 WCVI ecosystem, which includes other major fisheries such as herring and ground fish, to examine how changes in both predator abundance and variation in primary productivity have contributed to variability in per capita rates for pink shrimp. The last chapter summarizes all of the findings and discusses the relative importance of variability in primary production with respect to shrimp production on the WCVI shelf. 1.1 Pink Shrimp Fisheries Trawl fisheries for smooth pink shrimp in British Columbia began in Stuart Channel in the Strait of Georgia during the 1950s. Expansion occurred in Barkley Sound during the 1960s. Small inshore beam trawlers were used at the time, but during the early 1970s, as markets for BC 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. Catch 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 USA fishing fleet also participated in shrimp fisheries off WCVI, and in 1978 with the introduction of extended jurisdiction, USA participation in the fishery discontinued. Up until 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). By 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. By 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. Prior to the collapse, the Department of Fisheries and Oceans (DFO) 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. At the time, officials experienced concern about how quickly the fishery was developing. As a pseudo-experiment, area 125 was left open with no restrictions in place; they (the DFO) 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. DFO 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 little sign of any significant recovery in shrimp stocks in areas 124 and 125 until 1985. Biomass estimates from DFO surveys continued to decline from 1981 to 1983. No surveys were conducted in 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 in 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. Prior 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 in 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 until 1991. DFO 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. By 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: Pink shrimp Pandalus jordani landings in 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 addition to fishing mortality influencing shrimp abundance on the west coast of Vancouver Island, there is also evidence that ocean climates influence the bottom 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 Piatt 1999). There have been no directed studies of the influence of climate variation on pink shrimp population on the west coast of Vancouver Island. 1.2 Life history of 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 until hatching in March-May. Nauplii or free-swimming larvae spend 2-3 months in the water column and then settle by late June or early July. By October, as males, they are 60-70 mm in length. By the following autumn males reach 85 mm. Most adults become females in the spring of their second year, and grow to 100 mm by October. Maximum life span is 4 years and maximum sizes are 123 mm for males and 140 mm 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 (WCVI) shrimp trawl fishery targets the more abundant smooth pink shrimp and BC is the northern limit of fishable populations (Butler 1980; Boutillier et al. 1998). Smooth pink shrimp fisheries in BC were first reported in 1952, in Stuart Channel, and have also been operating in Barkley Sound since the 1960s. 1.3 Research Survey Program Annual research survey programs have been conducted off the West Coast of Vancouver Island since 1973, with 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 mortality rates (Boutillier 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. Prior 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 distribution of shrimp may be a continuous distribution that lies between the 100-150 m depth contour (see Figure 1.4). This survey program continues to monitor shrimp abundance in 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 trawl to sample to grounds. Each tow is roughly 30 minutes in 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 trawl is 11 meters; therefore the average area swept per tow is roughly 20,372 m2. 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 distribution. 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 until shrimp catches were negligible. Spatial catch per unit area swept data interpolated using a bicubic spline (Boutillier et al. 1998) was used to generate an estimate of absolute biomass for each statistical area (Fig. 1.3). Age composition samples collected over time are useful for monitoring changes in total mortality rates, assuming samples reflect the true population 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 molting 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 00 o o o o o o o o ll.li. 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: Location 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. 9 and growth parameters by comparing predicted length frequencies with observed data (Schnute and Fournier 1980). Given estimates of proportion-at-age, total mortality 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 mortality rates over time by convert ing the length frequency observations into proportion-at-age using probability 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 mortality rates are a result of variability in climate and pre dation. I then investigate the effects of climate variability on shrimp recruitment, and if variability in natural mortality is a result of changes in predation mortality. Finally, 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 WCVI shrimp stock. Specifically, the aims of the analysis were to: 1. Reconstruct the history of the WCVI shrimp trawl fishery, estimate shrimp abun dance, recruitment, and changes in natural mortality 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 mortality rates are related to predation by gadoid predators. 4. Explain how much of the observed variation can be attributed to trophic interac tions; both changes in predation over time and increases in production associated with favorable climate conditions. Chapter 2 Shrimp Recruitment and Mortality 2.1 Introduction In this chapter I simultaneously estimate recruitment and total mortality for pink shrimp in 2 statistical areas off the West Coast of Vancouver Island using an age-structured stock assessment model. The assessment procedure is broken down into 3 parts: a population dynamics model, observation models, and a statistical fitting procedure. The population dynamics model is used to represent unobserved states, or quantities we are interested in such as numbers-at-age in the population. The observation model translates output from the dynamic model into predictions of data types that we have observed, such as catch rates, or length-frequency data. Population model parameters describe rates of change from one state to the next, and observation model parameters relate unobserved states to observed states. Examples of observation model parameters are catchability coefficient, or size-selectivity parameters. Finally, a statistical fitting criterion describes how likely the data are given model parameters, and parameters are estimated by finding the most likely combination. In the following sections I describe the models and parameter estimation procedures used to estimate annual recruitment and total mortality. Both recruitment and total mortality estimates will be used in subsequent chapters of the thesis. There are 3 important vital rates that determine how abundance changes over time: 11 12 recruitment rates, natural mortality and growth. For short lived, small bodied organisms it is important to determine how each of these vital rates change over time. 2.2 Stock Reconstruction Model The goal of reconstructing shrimp populations on WCVI is to quantify recruitment and recruitment variation, total mortality and mortality 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 with 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 simply a function that equates a probability 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 distribution, with uniform prior probabilities for the parameters, using Markov Chain Monte Carlo (MCMC) methods (Gelman et al. 1995). 13 2.2.1 Population dynamics model I assume partial recruitment to the fishery occurs as early as age 1, and that all shrimp die before age 5. A description of all 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: Na+ht+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. Annual recruitment anomalies were assumed to be log-normal and constrained to sum to zero. I also assume that natural mortality rates are variable over time and have a mean M = 0.6 (Boutillier et al. 1998; Martell et al. 2000) and have a log-normal distribution. Variation in 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 mortality Ft can be calculated in one of two ways: either using fishing effort multiplied 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 WCVI 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: Model parameter descriptions and symbols used in population dynamics model Description Symbol Time subscript t Age subscript a Numbers-at-age at time t Nat Natural mortality rate in year t Mt = 0.6e"« Fishing rate Ft Average annual recruitment R Recruitment anomaly wt Mortality anomaly Vt Vulnerability at age Va Length at 50% vulnerability Lh Mean length at age a La Mean weight at age wa Shape parameter for logistic function 7 Annual landings ct Vulnerable biomass VBt where vulnerable biomass (VBt) in 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^VaWa (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 mortality rates (R,wi,... ,wn_i,vi,.. .,i/n_i,Lh,j). 2.2.2 Observation Models In the annual research trawl surveys there are two types of observations, an index of density (i.e. catch per unit area swept, or CPA) 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 with other observation model parameters. Predicting Relative Abundance DFO 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 Fig. 1.3 on page 8). I treat these estimates as a relative abundance index. Abundance indices generated by DFO 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 with 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 all Z statistics (see Walters and Ludwig 1994). Note here that the vulnerable biomass term in Equation 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 in 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. Description Symbol Subscript for age index a Subscript for cohort index i Subscript for length interval X Time at which sample was taken t Mean length-at-age a for cohort i 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 LA) can be initialized from modes in a length frequency distribution: LA = h + (LA - h)lZ PA-i (2-9) where li represents the mode of the first age class, LA 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 00 " 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 { A2 -1 + 2 1-p ,o-l (2.10) 1- pA~\ where X% and A2 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 A2 > 0 to ensure that standard deviation in length-at-age increases with age. Note that if A2 = 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, multiplied by the probability of being in length interval x summed over ages (Qx,t = Yla ^tPx,a) resulting in a predicted length frequency distribution. The contribution of length frequency data to the objective function was calculated using a robust likelihood method, which is described in Fournier et al. (1990). The log-likelihood equation for length frequency samples is: N n N LQ = -l/2j2/Zlo^27T(^ + 0.1/n)] -^nlogfa) t=l x=l N n t-l x=l t=l t=l i=l H 2fe, + 0.1/n) +001 (2.12) where Pxt and Pxt are the observed proportions and probability 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 multiplied by rt2, where rt2 = mm(Qxt,400). The variance term £xt for the expected probability of being in length interval x was estimated as (1 — Pxt)Pxt, therefore, if the probability Pxt — 0 then the variance term = 0. A small number (0.1/n) is added to the variance term to ensure £lt does not approach zero, and this also reduces the influence of length observations where the expected probability is near zero. One of the more difficult problems 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-lived and grow quickly, I assume a maximum of 4 age-classes present in all 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. Data 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 trawl durations and vessel speed (distance covered) from the logbook database maintained at the Pacific Biological Station, Nanaimo, BC. 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 Total num Area Area Area Area 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. Substituting this catch equation into the fishing rate equation F is now expressed as F — cEB/VB, i.e. F = cE. 21 Catchability 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. This is essentially the same definition of catchability as in Paloheimo and Dickie (1964), assuming the probability of capturing a shrimp within the swept area a is equal to 1. Substituting this definition of c into the calculation of fishing rate, fishing rate is now expressed as F = aE/A. 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 all 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 all vessels. Since 1998, the average foot rope length of all 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* 10-2 square nautical miles. The total area swept is approximated by the number of tows (Et) multiplied by the average area swept (numerator in Equation 2.13). It is also possible that only a fraction of the total shrimp within 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 Equation 2.13 would be multiplied 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 probability 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 WCVI. 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 all 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 WCVI 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 Equation 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 Model parameters were estimated using a Bayesian approach, where the distribution of errors for relative abundance data (Eq. 2.7) and area swept estimates of fishing mortality rates (Eq. 2.14) are assumed log normal. Both of these likelihood terms are weighted by the maximum likelihood estimate of the variance. The error distribution for the length frequency data was assumed normal (Eq. 2.14) and weighted by the relative variance multiplied by a variance scale factor (rt2). Large sample sizes were unweighted to a maximum sample size of 400, allowing more emphasis on smaller samples sizes (r2 determines the overall scale of the variance, see Fournier et al. 1990). Normally, likelihoods for length frequency data are based on the multinomial distribution 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 addition to the three criterion listed in the previous two sections (equations 2.7, 2.12, and 2.14) prior distributions were used to constrain recruitment anomalies and natural mortality rates. A total of 60 model parameters were estimated for each statistical area by minimizing 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 mortality rate parameters and here I assume 24 recruitment and natural mortality 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 mortality(a2kf). Recruitment anomalies were constrained to sum to 0 in order to estimate mean recruitment. The variance term for recruitment anomaly prior distribution was set <72 = 4.0, which is equivalent to a coefficient of variation = 0.54 in area 124 and 1.66 in area 125. The mean natural mortality rate was positively correlated with 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 mortality rates for each area at 0.6 and estimate annual deviations {ut). The variance specified for deviations in natural mortality 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 mortality rates are in shrimp populations. As cr2 and a2M approach 0, the model becomes a observation error only model with constant recruitment and a fixed natural mortality 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 mortality rates ( t^_1 ) are over a range of o2w and o2M. Model parameters were estimated over a 10x10 grid of <72 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, LA, p, A1( A2,7, lh, wt, ut) were estimated by minimizing equation 2.15. Minimization was carried out using a phased approach supplied in the AD Model Builder software (©Otter Research Ltd.), where some of the parameters are held fixed during the initial phase of the minimization process. Parameter estimation was broken down into two phases; in the first phase, all parameters with the exception of 7 and vt were estimated, then in the second phase 7 and ut were included. During the first phase of minimization additional 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 additional priors effectively reduce the number of iterations required for the non-linear search routine to converge while assuming that natural mortality rates was constant (M = 0.6) and the overall mean fishing mortality rate was 0.4. In the second phase I relax these priors and the objective function is no longer sensitive to deviations from the priors. Using a phased approach to parameter estimation in non linear models usually provides more stable estimates, especially in the case of over-parameterized models (Dave Fournier, pers comm., Otter Research Ltd., Sidney BC). To examine parameter uncertainty in a Bayes framework, assuming uniform prior distributions for all parameters, I sampled from the joint posterior distribution using a Metropolis-Hastings Algorithm (Gelman et al. 1995). The Metropolis-Hastings Algo rithm is a Markov Chain Monte Carlo (MCMC) 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 distribution that is the joint posterior distribution 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 distribution 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*|9t_i). If ©* is more probable than 9t_! 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 with 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 9t_i is included in the posterior sample again). Candidate points are selected using a jumping rule Jt(Q*\@t_i). For the jumping distribution I use a multivari ate normal distribution, using the "Cholesky" decomposition of the covariance matrix generated by the non-linear search routine (Press et al. 1996). After the algorithm con verges 10,000 random samples from the final posterior distribution are used to generate marginal posterior probability 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. This 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 with the assessment model. Results of the simulation studies are shown in Appendix A. 2.3 Results 2.3.1 Estimating Abundance Both predicted and observed indices of shrimp populations in areas 124 and 125 declined precipitously from 1975 to 1978-79 (Fig. 2.1). During the early 1980's the fishery had essentially shut itself down because shrimp were virtually absent in area 124 and low in 27 area 125. By 1985, stocks increased abruptly in area 124. No evidence of a significant increase in area 125 was observed until 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. Model predictions of exploitable biomass are in general agreement with biomass trends generated by DFO's spatial interpolation of area swept information from research trawl surveys. Residuals for the log(Z) statistics for areas 124 and 125 are shown in 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 in 1976 (Fig. 1.3), whereas predicted vulnerable biomass for that year was 3,653 tonnes. Recall that DFO 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 El Nino, shrimp stocks in both areas had collapsed, and abundance indices were less than model predictions for area 124 in 1983 and area 125 in 1982. 28 Area 124 Area 125 a> CC 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 in 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. This 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 Fig. 2.3. Since 1998, total biomass in both areas has been increasing, especially recently in area 125. Although population numbers have been declining in area 124 since 1999, biomass continues to increase due to reductions in total mortality 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: Population age structure for pink shrimp in areas 124 and area 125. White shaded polygon are age-1 shrimp and black shaded polygons are age-4 shrimp. Fitting 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 with the observed data with exception of f998, 2000, and 2001. In 1998, sample sizes were relatively small, and the youngest age-class is not apparent in the observed data, whereas model predictions suggest a relatively strong age-1 cohort. Conversely, in 2000, a strong age-1 cohort was observed in 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. Again, observed and predicted length frequencies in 1998 contradict each other and contradictions are probably attributed 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. Estimated 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. Estimated mean length-at-age is represented by vertical dotted lines. No length frequency samples were taken in 1984, 1986, 1989, and 1991 in area 125. 32 Table 2.4: Estimated growth parameters for area 124 and 125 pink shrimp populations. Area 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 Estimated growth parameters are shown in Table 2.4. Variation 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 Parameter Estimates and Uncertainty A total of 60 model parameters were estimated. Maximum likelihood estimates and standard deviations for all 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 in the data. Omitting 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 in the data between these two areas. Table 2.5: Estimated model parameters and standard deviations for statistical areas 124 and 125. Area 124 Area 125 Parameter MLE Estimate Std MLE 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 A2 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 all other parameters listed in table 2.5 differ markedly between the two areas. In general, shrimp in area 125 appear to grow faster, attain 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. Area 124 Area 125 Parameter Name MLE Estimate Std MLE 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 W1997 -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 distribution that was constructed using the Markov Chain 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: Estimated deviations in mortality (vt) and standard deviations for statistical areas 124 and 125 Area 124 Area 125 Parameter Name MLE Estimate Std MLE 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- As 7 increases in value, the slope of the selectivity curve through lh increases (i.e. large values of 7 infer 'knife-edge' selection). As the length at 50% vulnerability 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 in 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, LA, p,X\,^2,7,h) from an MCMC 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, A1; A2,7, lh) from an MCMC sample of length 250,000 for area 125. 37 with growth parameters (specifically L\ and p versus Ai). 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 A2 have several modes (Fig. 2.8). Recall that A2 was constrained to be greater than 0, such that standard deviation in length-at-age increases with age, and if A2 = 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 A2 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 cr2^ and cr2. Left column represent the mean estimate of natural mortality as a2M increases from 0.1 to 1, where the mean (line with points) is calculated over all a2M. Right column represents the mean estimate of M as increases from 0.8 to 8, where the mean is calculated over all cr2. In the left columns, dotted lines represent estimates of M when cr2 = 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 As u2M increases in value, estimates of mean natural mortality rates tend to increase for both statistical areas (left column of Fig. 2.9). For example, in area 124 the mean natural mortality rate estimate when a2M = 0.1 is 0.6, and as <r^ increases to 1, estimates of natural mortality rates increase to 0.86. As estimates of cr2 increase, mean natural mortality estimates tend to remain relatively constant, however, as o2M increases mean estimates of M increase (as shown by dashed lines in right column of Fig. 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 in 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 cr2. Left column represent the mean RQ as a2M increases from 0.1 to 1, where the mean (line with points) is calculated over all a2M. Right column represents the mean estimate of Ra as a2 increases from 0.8 to 8, where the mean is calculated over all a^. Dotted lines in left column represent estimates of RQ when cr^ = 0.8 and dashed line represents estimates of R0 when cr2 = 8, in 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 mortality rate anomalies. Variation in estimates of RQ and M over ranges of a2 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 with 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 Mortality Total mortality for pink shrimp in each statistical area was estimated from combined fishing mortality and natural mortality. Fishing mortality rates for the single species model were calculated using equation 2.4, where F is conditioned on observed catches from 1975 to 2000. I also constructed a time series 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: Estimated area over which shrimp populations in statistical areas 124-125 were distributed. Area estimates were calculated using commercial logbook data from 1987 to 2000; these data are shown in Figure 2.14 log-book records for the WCVI 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. On 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). During 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 in 1999, and estimated stock area declined from 348 square nautical miles to 30 square nautical miles (Figure The tow start locations shown in Figure 2.14 represent all of the tows logged since 1987 in all statistical areas (121-125). Prior to 1996, less than 50% of all 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 in statistical areas 123 and 124 (Fig. 2.14, and Table 2.3). Estimated 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 with stock area, as more fishing effort is attracted to areas where shrimp are more abundant. Estimated 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. During 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 (BRDs). By year 2000, 100% of industry had voluntarily adopted the use of BRDs. Although 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 all 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<^^^^^ 1997^- V^S. ' "-: . . T-r r "«^^^^^ 1989^- V^S. r r r r ^WT ^J < <f r Area 125 1. Area 124 ^^Qfl^gS Area 123 *C^^ r Figure 2.14: Spatial distribution of commercial fishing effort from 1987 to 2000 on the West Coast of Vancouver Island. Total number of tows from 1987 to 2000 was 59,962. 44 Table 2.8: Reported fishing effort in commercial logbooks for statistical areas 124 and 125 Year Number of Shots Total Effort Average Tow Time STD Tow Time Area 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 Area 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 particularly 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 Fig. 2.14), and tows were eliminated from area swept calculations. Finally, 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 mortality rates for areas 124-5 generated from the stock re construction model are shown in Figure 2.16 on page 47. Fishing mortality rates were calculated using observed catch and vulnerable biomass (i.e. F — —log(l — C/VB)), not total biomass. Estimated selectivity curves for the two areas are shown in Figure 2.15. Note that the estimated length at 50% vulnerability is much larger in area 124 (If, = 18.1) versus area 125 (Ik — 14.8). As a consequence, estimated fishing mortality 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 all time low in area 124, with a high proportion of partially vulnerable age-1 shrimp in the population. From 1976 to 1980, fishing mortality rates were declining in area 124, corresponding with the decline in abundance in area 124. In 1978 much of the fishing effort had moved north into area 125, generating high fishing mortalities. By 1979, catches were also 46 Figure 2.15: Estimated 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. Boutillier, Pacific Biological Station, Nanaimo, BC, 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). By 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, virtually all of the offshore fishing effort occurred in area 124 (Table 2.8), and fishing mortalities exceeded the recommended target exploitation rate of 33%. During this same time period, biomass in area 125 was beginning to rebuild; however, significant fishing effort did not resume in this area until 1991. In 1991, fishing mortality rates increased to 0.82 and 0.31 in areas 124 and 125, respectively. From 1991 to 1994 the majority of fishing effort in offshore areas alternated between area 124 and 125, as well as the fishing mortality rates. Estimated fishing mortality 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. By 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 At co CO CD 05 o CD i r 1975 1985 —i r 1995 Year Figure 2.16: Estimated fishing mortality rates from stock assessment model (lines in bot tom panels) for areas 124-125. Circles represent area swept estimates of fishing mortality rates (F = 9j^) from commercial log book data, and these data were used in the like lihood criterion for estimating model parameters. Top panels are the residuals between area swept estimates of F and estimated fishing mortality rates in the stock assessment model. abundance, and declining prices for shrimp. Also, in 2000 eulachon conservation concerns severely restricted fishing activities in the offshore areas, and the fishery was shut down prematurely, even though substantial shrimp catches were being landed early in the season. The trends in area swept fishing mortality correspond well with trends in fishing mortality rates estimated from the stock assessment model (see Fig. 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 mortality rates increases substantially. 48 There are a few major differences in the trends of F from the two calculations, espe cially mid-late 1990's. Typically, area swept estimates of fishing mortality are less than the predicted F from the stock assessment model during the 1993-97 time period. This time period corresponds with 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 mortality 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 in 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 mortality rates also played a role in this collapse (see next section). Whereas, area swept estimates of F indicate fishing rates were relatively moderate in 1998 in comparison to previous years. In fact, F's were estimated to be declining in both areas since 1995 (Figure 2.16) with 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 DFO. 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 all of the observed information from research surveys and relative estimates for fishing mortality based on area swept calculations. Estimates of natural mortality rates and uncertainty are shown 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 distribution 50,000 times. In both statistical areas, natural mortality 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 mortality rates over time in the two areas seem to be positively correlated, however, large changes in M appear to occur a year earlier in area 124 than in area 125 (Fig. 2.17). This pattern was also observed just prior to the 1997-98 El Nino, where increases in natural mortality occurred first in area 124. Table 2.9: Natural mortality rate statistics for areas 124 and 125 Area 124 Area 125 Average 061 063 Coefficient of Variation 0.21 0.35 Minimum 0.45 0.39 Maximum 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 mortality rates in area 124 were relatively stable, whereas in area 125, natural mortality 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 mortality rates from length frequency data is available during this time period. Furthermore, relative estimates of fishing mortality rates, the other component of total mortality, are not available until 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 <» Kiw#yMV*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: Maximum likelihood estimates of natural mortality rates (line with 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 distribution generated by the Metropolis Hastings Algorithm, and the box plots represent the median, inner quartile, and sampling range. 51 1987. Fishing mortality rates for area 125 in 1987 were low, and other auxiliary data suggest that populations were increasing with relatively low total mortality rates. 2.3.4 Annual Recruitment and Stock-Recruitment Relationships The maximum likelihood estimate for age-1 recruits and associated distributions from the posterior are shown in Figure 2.18. As 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 mortality led to declines in shrimp abundance in both areas during the early 1980's and recently in 1996 (Fig. 2.17 and Fig. 2.18). Recently, natural mortality 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. Area 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 timated age-1 recruits, and spawning stock biomass. Recall that Pandalus jordani are protandrous hermaphrodites, and the age at transition to female begins around age-2 with 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 multiplied 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: Maximum 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 million 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^ where the term in 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 Fig. 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 Fig. 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 in area 124. Omitting 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 omitting 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 omitting extreme outliers. 2.4 Discussion In stock assessment modeling, it is most often assumed that natural mortality rates are constant over time (Walters 2000). This simplifying assumption is usually made to reduce the number of estimated parameters. Natural mortality is confounded with 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 mortality rates was made possible by the inclusion of length frequency data and relative indices of fishing mortality based on area swept calculations. Changes in length frequencies over time provides information about the total mortality rate (both fishing and predation). If an estimate of fishing mortality is available then the other component of total mortality can be estimated from length frequency data, or proportions-at-age data (Sinclair 1998). 56 The use of length-based methods greatly improves our ability to infer proportions-at-age in 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. This 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 variability in the time of year at which they spawn can also bias length-based approaches. The problem shows up in assessment methods that directly estimate variance in length-at-age rather than estimating a single coefficient of variation for all 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 birth 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 (A2 = 0), or increases with age (A2 > 0). Estimates of A2 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). Estimating fishing mortality 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 WCVI shrimp trawl fishers, estimates of total area swept were carried out using an average area swept multiplied by the total number of tows. This was done because necessary information about tow distance and net width 57 was lacking for years prior to 1998. Unfortunately this averaging does not account for changes in 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 within 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). Adjusting the numerator in Equation 2.13 to accommodate changes in fishing technology over time is trivial compared to dealing with 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 NM2 over a map of reported tow locations and simply summing the number of cells with non-zero catch rates. Among all of the problems with 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. As the grid becomes infinitely small, the total area estimated will eventually reach an asymptote. Using today's high tech GPS information, fishing grounds within 1 NM2 may still 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 NM2 (see Fig. 3 in Walters and Bonfil 1999). The total area calculations assume that fishing effort is randomly distributed within each block, and violations in this assumption lead to over-estimating total area. There is another problem with using commercial catch records to estimate stock area, 58 which arises through changes in fishing fleet size. When there are good incentives to go fishing (i.e. high prices for shrimp, or high catch rates), fishing effort generally increases. Increased effort also leads to increased competition on the fishing grounds, and creates strong incentives for fishermen to explore new areas. It is generally thought that the effort distribution follows an ideal free distribution (Gillis et al. 1993; Walters and Bonfil 1999). In years where fishing effort is high and shrimp prices are high, more and more fishermen are willing to fish in areas with lower catch rates. They may even go explore new grounds in hopes of discovering potentially unexploited areas. When using spatial distribution of fishing effort to estimate total area over which the stock is distributed, high effort years will probably lead to an over-estimate of stock area. Conversely, low effort years will tend to under-estimate stock area, because remaining fishermen will 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 try to correct for this problem by only including records with non-zero catch . rates to justify the total area over which the stock is distributed. But 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 will 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 variability in the cost of fishing, will influence decisions to stay or move, or simply quit fishing altogether. Including independent estimates of fishing mortality rates in the estimation procedure greatly reduces the uncertainty associated with estimating fishing mortality 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 mortality rates (Z). Estimating the F component of Z was based on total catches and the likelihood of area 59 swept estimates of F from commercial log book records. For the purposes of this thesis, I assumed any errors that arise through violations in assumptions about estimating fishing mortality rates from area swept estimates are independent, and therefore not biased if used as a relative index. In Figure 2.16, the trends in fishing mortalities are somewhat similar in both statistical areas. But is it really fair to make such a comparison? To determine if these additional observations were very informative about fishing and natural mortality rates, the same estimation procedure should be carried out with little or no weight on the area swept data. I did not carry out this analysis in the thesis due to limitations in computing time, (it takes 18 hours to estimate uncertainty in model parameters on a Pentium IV computer). In the simulation study (Appendix A), including estimates of fishing mortality 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 Fig. 2.20 is a straight line, even when omitting extreme outliers when spawning stock is small. The simplest interpretation of this result would be that recruitment is totally independent of stock size within each statistical area. For example, there is little evidence that recruitment in statistical area 125 is a result of spawner abundance in area 125. An 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 AN 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 (WCVI). The physical oceanography on the WCVI is strongly connected to seasonal primary productivity 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). Although adult shrimp in the offshore areas are found between the 100m and 200m depth contours, the persistent coastal current still may provide a conduit for shrimp larvae to disperse northward. Also, year classes tend to segregate on the fishing grounds; for example, smaller northern pink shrimp are generally found at the edges of each fishing ground in shallower waters (Koeller 2000). It is possible that young of the year shrimp settle out on the banks then slowly move offshore to the shelf break as they grow. If currents 60 61 advect larvae, then I expect recruitment indices generated in Chapter 2 to be correlated with 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 with upwelling indices. If, however, recruitment anomalies are positively correlated with upwelling events then I would suspect that increases in primary productivity led to increased recruitment. Effects of variation in primary productivity will be examined further in Chapter 5, in the context of an ecosystem model that explicitly represents trophic interactions in the region. 3.1.1 Summer and Winter Current Circulation Patterns The persistent coastal current is primarily driven by surface outflow from Juan de Fuca Strait during the summer, and southeasterly winds during the winter. Much of the current's energy originates from runoff of low salinity waters through Juan de Fuca Strait. The coastal current off the WCVI during summer months (April-September) tends to spread offshore along the broad continental shelf North of the Juan de Fuca Strait and south of Barkley Sound. As the current travels north it forms a narrower band north of Barkley Sound. Northwesterly winds predominate on the WCVI 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. At depths beyond the 100 m depth contour, the shelf break region, ocean current flows in a southeasterly direction, parallel with 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 via Ekman 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 in nutrients and enhance primary productivity, except in El Nino years when the thermocline is much deeper (Mysak 1986). In years with strong summer northwesterly winds, primary productivity can be significantly enhanced off the WCVI (Mackas et al. 1987; Robinson and Ware 1994; Robinson and Ware 1999). 3.1.2 Distribution of Shrimp Relative to Coastal Currents and Shelf Break Currents 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 in 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 will induce upwelling and drive surface currents offshore and if particles on or near the surface move far enough offshore they will 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. Statistical areas 125, 124, 123 and 121 are shown in black, to lighter shades of grey, respectively. Depth contours are 100m, 200m and 1000m will 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 tical migrations (DVM) have a mechanism to be transported in opposite directions over large scales. At smaller spatial scales, DVM specialists are also subject to displacement resulting from the ebb and flood of tidal 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 with models that include an environmental variable, using methods similar to Caputi et al. (1998). I also examine if spawning stock biomass in adjacent statistical areas is correlated with 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 with 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 Rinaldi and Gatto 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 rithm (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 (Hilborn and Mangel 1997), and twice the difference in log-likelihoods has a x2 distribution with 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 ability 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 Royall 1995; Hilborn 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 distribution 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 with a weighting term (ipt) that reduces the influence of extreme outliers on the estimation stock recruitment parameters. Equation 3.5 is sensitive to outliers because the normal distribution has a small 'tail' and the probability of an extreme recruitment event at low spawner density drops off quickly after a few standard deviations from the mean. Large outliers tend to contaminate likelihood functions, and the estimation procedure "over-fits" the anomaly. There are two approaches to reducing the influence of outliers: 1) changing the probability distribution to have 'fat-tails'(see 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. As c increases in value tpt approaches 1, and all data are weighted equally. Note that for normally distributed errors the optimal value of c is 6.0 (Press et al. 1996). The "robust" likelihood equation is written 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 initial estimates of a, b, a, 7 (bt — vt — log(Rt) — (log(f(St)) + ~fWt)) and (3.7) 67 3.2.1 Environmental Indices I examine six different environmental correlates: sea surface temperature (SST), Bakun upwelling indices in May-June, Sea Level Height (SLH) in March-April, and annual SLH anomalies, the Southern Oscillation Index (SOI), and the Pacific Decadal Oscillation (PDO) in March-November. All 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. Monthly mean sea surface temperatures are available from various lighthouse stations1 since the early 1900's. SST data were extracted from Amphitrite point, which is located at 48°33'N Latitude, 125°19'W Longitude (near Ucluelet, British Columbia). Warm periods are represented by positive anomalies (e.g. 1982-3 and 1997-98 El Nino's in Figure 3.3). Upwelling indices2 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 Ekman layer. Typical 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 m3s_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 (PDO) is actually a sea surface temperature index for the North Pacific Ocean (Mantua et al. 1997). The term PDO 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 with 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 in 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 all 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 with 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 likelihood. In area 124, the May-Jun upwelling indices significantly reduced the likelihood over the Ricker model with 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 Oscillation 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 PDO are positively correlated (Table 3.2). In contrast, the PDO 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 PDO 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 PDO 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 Model Ricker Model Model logL x2 logL x2 Area 124 Base 19.48 — 21.96 — — SST 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 SLH 19.42 0.73 2.63 21.86 0.65 -0.54 mSLH 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 PDO 19.24 0.49 -1.07 21.91 0.75 2.50 Area 125 Base 20.72 — — 26.02 — — SST 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 SLH 20.63 0.67 -0.39 26.01 0.88 -0.32 mSLH 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 PDO 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 SST Bak SLH mSLH SOI PDO SST 1 -0.112 0.365 0.624 -0.646 0.727 Bak -0.112 1 0.263 0.192 -0.127 -0.061 SLH 0.365 0.263 1 0.720 -0.208 0.403 mSLH 0.624 0.192 0.720 1 -0.475 0.505 SOI -0.646 -0.127 -0.208 -0.475 1 -0.524 PDO 0.727 -0.0614 0.403 0.505 -0.524 1 Using stock-recruitment models with spawning biomass in the adjacent statistical area as a correlate explained very little 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 all < 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 PDO and SST indices in 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 Fig 3.3 Area SST Bak SLH (Mar-Apr) SLH (.lan-Dec) SOI PDO (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, with 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 = 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 Model Ricker Model Model logL x2 logL X2 % <r2 Area 124 Sp awners Area 125 Recruits Base 19.97 — — 19.87 — — SST 19.59 0.38 2.61 19.09 0.21 5.60 Bak 19.65 0,42 2.40 19.52 0.40 2.50 SLH 19.71 0.47 2.14 19.79 0.68 0.51 mSLH 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 PDO 19.89 0.68 0.78 19.63 0.48 1.87 Area 125 Spawners Area 124 Recruits Base 19.25 — — 28.82 — — SST 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 SLH 19.25 0.95 -0.01 28.60 0.50 0.15 mSLH 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 PDO 17.80 0.09 18.87 27.04 0.06 20.45 3.4 Discussion 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. As mentioned in 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 with 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 did 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 Royall 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). Variability 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 with early fall sea surface temperatures, roughly six months after larval release. Cooler water temperatures were also a significant abiotic factor in the survival and distribution 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 Miller 1983). There are two possible outcomes resulting from strong upwelling events off the WCVI. First, 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, Fig 2.20 on page 54). There was no consistent information in 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 little of the residual variation. Again there were no consistent relationship for both statistical areas that would suggest increased primary production has resulted in increased survival. However, shrimp are primarily benthic detritus feeders, and it is possible that good upwelling conditions may provide good feeding opportunities for benthic juveniles in 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 in area 124 and area 125. Dotted 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 El Nino years) larvae are advected northward resulting in recruitment failure in Oregon-Washington waters. There is some evidence of an age-1 recruitment event in 1984 (as indicated by positive recruitment anomalies in 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. Up until 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 in area 125 in more recent years. During 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 timing 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 April. 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 timing of larval release, then ocean conditions for shrimp larvae will be less productive, or will transport shrimp larvae onshore into areas such as Barkely Sound. Hannah (1993) also suggested that the timing 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 in 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. All 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 Introduction Several recent studies have indicated strong predator-prey interactions between members of the family Gadidae and Pandalid shrimps (Berenboim et al. 2000; Hannah 1995; Lilly et al. 2000). Assuming that gadids, such as Pacific cod (Gadus macrocephalus) and Pacific hake (Merluccius productus), are principal predators of pink shrimp species (Lilly et al. 1998; Westerheim 1996), here I examine correlations between Pacific cod and Pacific hake abundance versus pink shrimp abundance, mortality and recruitment off the west coast of Vancouver Island. Previous work in Oregon has demonstrated a positive correlation between hake abundance and pink shrimp mortality (Hannah 1995). If these two predators are strongly influencing the dynamics of shrimp populations (top-down control) the following predictions should hold: predator biomass and shrimp biomass should be negatively correlated, predator biomass and shrimp mortality should be positively correlated and shrimp recruitment and predator biomass should be nega tively correlated. The later prediction represents predation on juvenile shrimp, or that juvenile survival rates for pink shrimp are correlated with 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 WCVI 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 BC 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 BC 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. Map courtesy of Alan Sinclair, DFO Pacific Region. Stock assessment for Pacific cod has relied on commercial CPUE indices for estimating abundance in all areas. There have been no directed, fishery independent surveys for Pacific cod on the west coast of Vancouver Island. The most recent DFO assessment document (Sinclair et al. 2001) included an index of abundance generated from shrimp trawl surveys. Although the index was not used for the DFO assessment, the trend information from shrimp surveys seems to correlate well with commercial indices since 1981. Because I am interested in calculating predation mortality on shrimp by Pacific cod, I choose to include the shrimp survey index of Pacific cod in this assessment. 80 4.1.2 Stock Structure of Pacific Hake Two distinct types of Pacific hake (Merluccius productus) populations exist along the West Coast of North America, 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 Pikitch 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. By mid summer, adult hake form large midwater aggregations near the continental slope off Vancouver Island and feed on euphausiids, other pelagic schooling fishes, and Pandalid 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. During El 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). During the warmer years, the distribution 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). Estimating the fraction of mortality 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 El Nino events; more hake are found in Canadian waters during El Nino years (Mark Saunders, pers comm., Pacific Biologi cal Station, Nanaimo, BC). Therefore, if hake are significant predators on pink shrimp, mortality rates for pink shrimp should be significantly higher in El Nino years. Pink shrimp are not listed as one of the major diet items of hake feeding in the Canadian zone (Tanasichuk et al. 1991). However, in Washington and Oregon, predation mortality by Pacific hake on pink shrimp is thought to be significant (Rexstad and Pikitch 1986; Han nah 1995) even though pink shrimp are a rare diet item for Pacific hake in that region. All 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), with the exception of Gotshall (1969). In Gotshall's case, he was trying to develop a relative abundance index for pink shrimp and specifically sampling hake aggregations in areas where shrimp are known to be present. Pink shrimp are found more frequently in the stomachs of larger hake (> 55 cm TL) than in smaller hake (Rexstad and Pikitch 1986), and larger hake are generally found closer to the bottom (Mark Saunders, pers comm, Pacific Biological Station, Nanaimo, B.C.) whereas smaller hake generally form large pelagic schools. Rexstad and Pikitch (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: Distribution of Pacific hake in 1992, 1995, and 1998 along survey transect lines. Abundance is measured using hydro acoustic backscattering techniques. (Figure supplied by Mark Saunders, DFO). 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 Estimating Predator Biomass 4.2.1 Pacific Cod Assessment A time series of biomass for Pacific cod on the west coast of Vancouver Island was reconstructed using a delay-difference model (Deriso 1980; Schnute 1987). Two age groups are represented in this model, age-2 recruits and fish age 3 and older. The 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 Growth 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). Initial 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 Population numbers and biomass over time were calculated using equations 4.8 and 4.9, respectively. Annual survival (St) is a function of fixed natural mortality and annual fishing mortality, where mortality associated with fishing is conditioned on observed effort (Et, Equation 4.6). Annual 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.7Nt+1 = NtSt + Rt (4.8) Bt+1 = St(aNt + PBt) + wrRt (4.9) Rt = SoBt-KeWt-rtet* (4.10) Predicted 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) Likelihoods Likelihoods 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 zcrM 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. Likelihood profiles were used to explore uncertainty in leading parameter estimates (M, RQ, and 5). Uncertainties in 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 Pacific Hake Assessment 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 mortality is partitioned 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 in 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 La Peruse bank region (Ware and McFarlane, 1995), and for hake in the Canadian Zone, (Dorn et al, 1999). Year Hake biomass in La Peruse bank region (Kt) Hake biomass in Canadian Zone (Kt) 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 in Bamfield are positively correlated with 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 in 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: Proportion of Pacific hake in the Canadian zone versus mean January-February sea level height anomalies (top panel). Positive anomalies indicate higher than average sea level heights in Bamfield, British Columbia (data provided by Mark Saunders, DFO). Estimated 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 Pacific Cod Biomass The reconstructed biomass and fishing mortality rates for Pacific cod are shown in Figure 4.4 (top panel, residuals and recruitment anomalies shown in bottom panel). During 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. During 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 with the results obtained by Haist and Fournier 91 (1995). Since 1996, the ground fish trawl 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: Estimated 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. Bottom 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 mortality rates (Ft) averaged near 0.l6~yr, and dramatically increased to a maximum of 0.67~yr in 1994. Since 1996, however, the new trawl quota management plan has apparently reduced fishing mortality on Pacific cod to om~yr. 92 Management of the ground fish trawl fishery underwent a transformation 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. Prior 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 CPUE 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, CPUE indices would not reflect recent increases in abundance that may have occurred. I tried to correct for this effect by including abundance indices generated from the shrimp trawl 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 in prior distributions on natu ral mortality rates, and recruitment compensation (see Fig. 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 WCVI cod. Changing the mean of the prior distri butions by 10% resulted in less than 1% changes in key population parameters. Results were also fairly robust to assumed variances on prior distributions for natural mortality 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: Likelihood profiles for WCVI 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 Correlation with Predator Abundance Pink shrimp and Pacific cod (adult and juvenile) abundances on the WCVI 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 in the mid 1970s and reached lows a decade later. During 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 with either Pacific hake time series presented in Table 4.3 or with 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 mortality 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 Pink shrimp, Pacific Cod, and Pacific hake series I and II from 1975 to 2000. Lower panels are scatter plots of predator abundance versus pink shrimp abundance. 4.4 Discussion In contrast with 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. During the mid 1980's both Pacific cod and shrimp populations had declined to low abundance, and by 1985, shrimp populations started to recover. With 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. On 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 mortality (see Table Tab4.5). During warm years more hake move into the Canadian zone, apparently resulting in increased predation mortality for many other species including herring and euphausiids (Ware and McFarlane 1995; Robinson and Ware 1999). This is especially evident during the El 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 with increasing hake abundance. If Pacific cod were a key mortality agent on pink shrimp, I would have predicted substantial increases in shrimp abundance when Pacific cod biomass had declined in 1985 and 1995 (Figure 4.4), and significant positive correlations between shrimp natural mortality rates and Pacific cod abundance. Such observations have occurred on the east coast of Canada (Lilly et al. 2000, Worm and Myers In Press), and in the Barents Sea (Berenboim et al. 2000). On the west coast of Vancouver Island, pink shrimp were positively correlated with Pacific cod, suggesting that both species are responding to changes resulting from bottom-up forces (see Chapter 5). This 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. Typically, in warm El 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 limited 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 Pikitch (1986) of 659 tonnes to the west coast of Vancouver Island population in 1983, more shrimp would have been consumed than were available. Along 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 mortality rate of 0.33 or approximately 655 tonnes yr-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). By 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 productivity 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 with pink shrimp population 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 in 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 mortality 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 variability 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 mortality also appears to be related to temperature; it is possible that hake are responsible for at least some of the variability 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; "ERSEM", Baretta et al. 1995, "Ecopath with Ecosim" Walters et al. 2000) where interactions (trophic and fishery) and variation in primary production can be explored. Explicitly 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 productivity 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 with Ecosim (Pauly et al. 2000) to partition mortality among natural predators, and commercial fisheries and examine how total mortality changes over time. The first step in demonstrating the utility of an ecosystems model is to show that it can replay the history of disturbance in the ecosystem. To do this I will first fit time series data to Ecosim predictions. This 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 WCVI ecosystem is a very prominent feature and plays an important role in primary 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. Bottom up control in Ecosim is represented by specifying rates of flow from prey to predator in a term known as flow control, or relative vulnerability of prey to its predators (see Bundy 2001; Walters et al. 1997). 5.2.1 Ecopath with Ecosim 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 km-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 mortality that is accounted for within the modelled system and Y", are fisheries yields. Equation 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 mortality 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 km2. Biomass estimates and total landings were scaled to tonnes km-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 mortality 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 mortality), Fi<t is the fishing mortality rate at time t and cij]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. Predicting 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) equilibrium 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 Equation 5.5 using the Ecopath estimates of Qijto, -Bi,o and -Bj,o- Vulnerability exchange parameters v^j are calculated from a user specified ratio of the maximum predation mortality to base predation mor tality estimated in Ecopath (see description on page 106). Large vulnerability 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 Equation 5.5 predicts a saturating relationship between predation mortality My = Cij/Bi and total predator biomass. This form of relationship tends to eliminate much of the unstable behavior that Lotka-Volterra type predator-prey models typically produce. This tends to conserve total biodiversity in even the most complicated ecosystem models due to limitations on total predation mortality rates (Walters et al. 1997; Walters et al. 2000; Walters and Kitchell 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 Equation 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 birth to the age at recruitment to the adult stock. The complete set of derivations and assumptions used in 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 implicit), the basic delay-difference model structure for all split pools is: BA,t+i = e-ZA*[aAt{CAt)NA,t +PABA,t)+wJAtNj^t (5.6) NAit+1 = NAite-z^ + Njikit (5.7Nj,llt+i = R(BA,t,NA,t,CA<t) (5.8) Nj,a,t+i = e'^Nj^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 vs- 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 mortality M0, total fishing mortality 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 all predators j that eat group i adults and juveniles. Density-dependent mortality in juvenile pools Ecosim produces density-dependent juvenile mortality from age-0 to age-fc in two ways. First 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 vulnerability rates to juveniles are set low, then per capita consumption rates decline rapidly with 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 mortality rates (i.e., longer time required to reach w). Second, simulated juveniles in the model adjust the amount of time spent feeding in response to changes in prey availability and intraspecific competition. Ecosim represents this effect by using a function that adjusts feeding time to try 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. At 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., time 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 all its predators j (including juveniles as predators): IT) V) \ — aijvi,jBi}tBjitTiitTjjt /CIA\ C,-JJ {Bl" Bjx) ~ (vhJ + vhJTht + ahjBj/T^ (514) Initially, I used the feeding time adjustment method to represent density-dependent mortality for juvenile pools in the WCVI model. Here I set feeding time adjusment rate a = 0.5, and a low vitj for juvenile pools. I did 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 vulnerability 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 vulnerability 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 mortality 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 mortality increases linearly with 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: An example of the relationship between predation mortality 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. As the value of k^ decreases the slope of the line through S^n, mij,o decreases and the maximum predation rate (vtj) of predator j on prey i decreases. 5.2.2 Fitting Ecosim To Time Series Data Data Types Two 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 mortality 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 all observed data given model predictions was then calculated using: There is no limit to the number of data sets that can be used in the fitting criterion, and the contribution of each data set to the total likelihood 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 all data are weighted proportional to the inverse of the variance. Estimating 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; "Fitting Ecosim to Time Series Data" for a description of the Marquardt Algorithm), to maximize the total log-likelihood function (Equation 5.17). These parameters are vulnerability ex change rate parameters (vij) and a time series of primary production anomalies. In analogy with single-species assessment models, vulnerability exchange rate parameters represent mean production responses, and primary production anomalies represent pro cess errors. 5.2.3 Description of Model Groups The ecosystem model was used to examine trophic interactions that contributed to vari ability in pink shrimp abundance and mortality. Therefore, I chose to build 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) Ill lacking. Model groups included known shrimp predators (Pacific cod, Lingcod, and Dog fish), primary 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. All 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 lingcod were taken from Cass et al. (1990), and Q/B ratios were modified from Venier and Kelson (1996). Lingcod were split into ju venile and adult components, with 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 An estimated PjB ratio for dogfish was taken from the natural mortality rate estimates provided in Walters and Bonfil (1999), assuming fishing mortality in 1950 was negligi ble. Assuming individuals die at a constant instantaneous rate Z, then the production-biomass ratio is equivalent to total mortality rate, regardless of metabolic rates (Allen 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 mortality rate estimates in (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 DFO stock assessment model (data provided in 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 in 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 in the WCVI 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. Both species were aggregated into a single group, and were assumed to have similar production and consumption parameters. Production 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 mortality was being explained in the ecosystem model (EE — 0.6). The biomass estimate from Ecopath is slightly less than the 22 t km-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 ^npv qSIJ [13S.I9U.18Q POO aHP^d 3InPV IISIJSOQ pooSuiq Anr pooSuiq a, O rH d d oo ° 00 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 ^ odd 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 rHf-LOCOOOCOCOt— OrHOJrH^fOOrH qorHqqqt-^q 00 , I H H H ° LO ° S 05 3 ° o o H co H o d ° d ° d ° d rH CN rH rH LO OOOOO OOOOO 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 PH a3 « CH 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 width of the connecting line is proportional to the flow. Position on the y-axis is the trophic level of the group. 5.2.4 Time Series Data for West Coast Vancouver Island Time series data used in the WCVI ecosystem model consisted of absolute catches, to tal mortality rates, relative abundance indices (either generated from stock assessment models or direct observations in commercial catch rates) and fishing mortality rates used as time series forcing. Observed catch time series are shown in Figure 5.3. Although commercial fisheries for lingcod started in the late 1800's, catch data for southwest Vancouver Island was not separated until 1954, and I assumed catches were relatively stable between 1950 and 1954 (not shown in the figure). Commercial fisheries for pink shrimp on WCVI started in 1973, and landings peaked in 1975. Herring fisheries collapsed in 1968 and the fishery was closed until 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 oi 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 km-2) for west coast Vancouver Island from 1950 to 2000 (left column). Total mortality rate estimates for pink shrimp and adult Pacific herring. Pink shrimp total mortality rate estimates were calculated from combined fishing and natural mortality rate estimates in chapter 2. Total mortality rates for Pacific herring were inferred from age composition data. Bottom 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. Adult herring data were obtained from Jake Schweigert (DFO, 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 CPUE data, and Pacific cod from assessment in Chapter 4. Euphausiid abundance (from 1967 to 1998) estimated from a trophodynamic model by Robinson and Ware (1999). solar radiation). An 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. All indices were scaled to have a mean =1. 119 5.3 Results 5.3.1 Fitting Ecosim to Time Series Data Predicted catch time series from Ecosim were consistent with 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 in Ecosim. In all cases, the relative contribution of catch time series to the total likelihood was deliberately reduced to 10% of the biomass time series. After some initial fitting runs, it became apparent that the predicted lingcod biomass time series contradicted model predictions. When the lingcod time series was included in the total likelihood, the predicted changes in herring and Pacific cod were inconsistent with observed data given the estimated sequence of primary production anomalies. As 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 all of the relative abundance data sets shown in Fig. 5.4; abundance was high in the 1970's followed by a collapse in the early 1980's (a period of low productivity apparently associated with the 1982-83 El 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 Much of the total mortality estimated for Pacific cod was a result of fishing, and Ecosim predictions correspond well with 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 mortality 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. Although Ecosim 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 with 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. During 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). Ecosim pre dicts that eulachon predation mortality declined in 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 mortality 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 mortality rates. Simulated euphausiid production was closely connected to variation in primary productivity. During mid 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 all time high, however Ecosim predicted euphausiid abundance did not increase, due to large increases in predation mortality by herring. During 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. Annual primary production anomalies were calculated in 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. Estimated mean primary production during the 1950-1977 period was 25% higher than mean primary production from 1978-2000 period (Figure 5.7). Estimated primary production reaches lows during the 1982-83 and 1997-98 El Nino years. ~n 1 1 1 1 r-1950 1960 1970 1980 1990 2000 Year Figure 5.7: Estimated 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 with large anomalies. Using the time series shown in Fig. 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). With the exception of 5 data sets, including simulated variation in primary productivity resulted in better model fits, especially for Pacific herring, Pacific cod and shrimp. For lingcod, adding variation in primary production over time increased the residual variation between predicted and observed biomass (-189% represents an increase in residual variation) catches and recruitment. Deviations in mean body weight for adult Pacific cod increased (4%) when including variation in primary production, and residual variation in shrimp mortality increased 28%. Prior to 1973, there is no information about trends in abundance for pink shrimp on 124 Table 5.3: Squared residuals for No Model (or deviations from the mean of the data), a model with no trophic interactions (assume that predation mortality and consump tion are constant), a model with 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 with variable primary productivity 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. As 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 with increases in primary productivity. From 1975 to 2000, simulated trends in shrimp biomass are consistent with combined abundance estimates from Chapter 2. In comparison with Ecosim predictions, the single species model indicates more variability in abundance, and the recent decline started later (1995) in comparison to the decline predicted by Ecosim (1989, see Fig. 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. No 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 with estimated diatom produc tion from Robinson and Ware (1999) during the 1990's. During the mid 1970's Robinson and Ware estimate an increasing trend in diatom production, whereas Ecosim predicts a declining trend in primary production (see Fig. 5.8). This discrepancy occurs because abundance estimates for shrimp, herring, and Pacific cod were all 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 Diatom 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 in Ecosim also resembles the nega tive pacific decadal oscillation (Figure 5.9). The PDO and estimated primary production are negatively correlated, thus during the negative phase of the PDO apparent primary production off the west coast of Vancouver Island is higher than average. 126 C o o ir 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 Oscillation (PDO). 5.3.2 Changes in Shrimp Mortality During the 1990's, total predation mortality on pink shrimp has declined (Fig. 5.10) cor responding with apparent declines in predators. The principal predators in the ecosystem model, juvenile and adult Pacific cod, have been declining since the mid 1970's. During the period 1950-70, mean predation mortality rates induced by Pacific cod predation were roughly 0.51~yr, or 78% of the total shrimp mortality rate. During 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 mortality 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 mortality (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: All mortality components from Ecosim compared to Z from single species model (SSM). The area between each line is the total mortality due to each predator or fishing, and the area between total mortality and fishing mortality is the unexplained mortality (1 - EE). by setting the EEi term equal to 0.95). There are additional predators that were not included in the Ecopath model, and such additional predation is assumed to be constant over time in Ecosim. In contrast, the single species assessment model in Chapter 2 treats natural mortality (M) as an arbitrary temporal variable to be estimated from the data. Overlaid on Figure 5.10 are the relative estimates of Z (F + M) from the single species model. There are two points to note about this comparison: 1) total mortality for 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 with 1982-83 and 1997-98 El Nino events (see Fig 5.11). Ecosim predicted an increase in total mortality rates during the late 1980's and 1990's. This increase in total mortality rate was a result of increased fishing mortality. During 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 Discussion The goal of this chapter has been to examine shrimp population dynamics off the west coast of Vancouver Island (WCVI) in an ecosystem context. I used an ecosystem model to look at how changes in predation mortality and variation in primary production affect the population dynamics of pink shrimp. The ecosystem model suggests that reductions in predator abundance by commercial fishing has reduced predation mortality for pink shrimp. This reduction in predation mortality has apparently not resulted in an 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 mortality (see Fig. 5.10). Both dogfish and lingeod populations have been declining during the last 2 decades as a result of by-catch (in trawl fisheries) and directed fisheries, respectively. Despite recent declines in 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 mortality that is explained by trophic interactions and fishing mortality. The unexplained proportion (1 — EE) implicitly represents other predation losses, disease, or senescence. If predation mortality is additive (i.e. modelled predators are not consuming shrimp that are going to die anyway), then unexplained loss rate (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 partition from the coast wide hake stock. The migratory hake stock is large relative to prey groups specified in the model. As discussed in Chapter 4 the spatial distribution of hake is highly variable, even within 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 mortality rates, hake should be included 132 in the model. Primary and secondary prey species for hake are euphausiids and herring (Ware and McFarlane 1995), which are also primary forage species for all large fish in 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. As noted in Chapter 4, very small increases of shrimp in the daily ration of Pacific hake could translate into very large changes in predation mortality rate. As 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 variability in Z for shrimp. The results in Table 5.3 are comparisons of 3 different models in which (1) trophic in teractions are assumed constant, (2) trophic interactions were included, and (3) both trophic interactions and variability in primary productivity are included. Assuming trophic interactions are constant over time is equivalent to using a single species model with constant natural mortality rates. By comparing a model with fixed natural mor tality rate to a model with variable natural mortality rates, we should be able to make inferences about how much of the dynamics are explained by trophic interactions. Less of the variation was explained for top predators (lingcod and Pacific cod) by including trophic interactions. In other words, much of the dynamics for lingcod and Pacific cod are determined by fishing mortality. In contrast, including trophic interactions explained a significant amount of variation in Pacific herring abundance and mortality, and variations in natural mortality are likely, as was determined by Schweigert (2000). Including trophic interactions for pink shrimp did 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 mortality model. In summary, the west coast Vancouver Island ecosystem model explains trends in 133 catch and biomass time series for pink shrimp. These trends were consistent with pre dictions generated in the single species model. The recent decline in Pacific cod explains some of the changes in total mortality rates for shrimp; however, Ecosim was unable to capture the two large natural mortality rate anomalies that occurred during the 1982-83 and 1997-98 El 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 mortality anomalies. The positive correlations in shrimp and Pacific cod abundance are best explained by variation in primary production. Estimated variation in primary production corresponds well with other independent estimates of primary production, especially in 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 optimum 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 Martell 2002). Even more complex age-structured models with hundreds of nuisance parameters require these 2 key parameters for determining optimal 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 variability (Walters and Parma 1996). However, fixed harvest rules may also lead to pathological over-fishing because of delayed responses (i.e. Walters and Kitchell 2001) or non-stationary effects, and here is where the utility 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 Variability in natural mortality 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 mortality rates using the assessment model in chapter 2. By 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 mortality rates (Z) can be inferred by tracking changes in the length com position data. Also, the addition 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. Natural mortality rates for shrimp were estimated to increase during El 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 distribution 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 with findings about pink shrimp-hake interactions in Oregon and Washington (Hannah 1995). Predicting changes in natural mortality rates for pink shrimp requires identification of the mechanisms responsible (e.g. hake predation) and establishing a functional relationship for this interaction. In chapter 5, I examine the shelf populations of pink shrimp from an ecosystem 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 mortality could be explained by including trophic in teractions and variability in primary production in the assessment process. In summary, including trophic interactions while assuming constant production degraded statistical fits to the data; however, including variable primary production explained an additional 29% and 15% of the variability 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 in the ecosystem assessment were consistent with 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 with regimes shifts documented by Beamish et al. (1999). Bibliography Allen, K. R. (1971). Relation between production and biomass. J. Fish. Res. Bd. Canada 28, 1573-1581. Anderson, P. J. and J. F. Piatt (1999). Community reorganization in the Gulf of Alaska following ocean climate regime shift. Mar. Ecol. Prog. Ser. 189, 117-123. Baranov, F. (1918). On 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. Klyashtorin, 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 multi-species models. J. Northw. Atl. Fish. Sci. 27, 69-75. Boutillier, J., I. Perry, B. Waddell, 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. Bull. Fish. Aquat. Sci. Caputi, N., J. W. Penn, L. M. Joll, and C. F. Chubb (1998). Stock-recruitment-environment relationships for invertebrate species of Western Australia. 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, NJ. 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. Wilson, M. A. Guttormsen, K. Cooke, R. Kieser, and M. E. Wilkins (1999). Status of the coastal Pacific hake/whiting stock in 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 Model Simulations of Tides and Bouyancy Currents Along 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). Multifan a likelihood-based method for estimating growth parameters and age composition from multiple length frequency data sets illustrated using data for Southern bluefin tuna (Thunnus maccoyii). Can. J. Fish. Aquat. Sci. 47, 301-317. Fournier, D., J. Sibert, and M. Terceiro (1991). Analysis of Length Frequency Samples with Relative Abundance Data for the Gulf of Maine Northern Shrimp (Pandalus borealis) by MULTIFAN 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, May 22,2002. Fu, C. and T. J. Quinn II (2000). Estimability of natural mortality 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. Carlin, H. Stern, and D. Rubin (1995). Bayesian Data Analysis. Chap man & Hall. Gillis, D. M., R. M. Peterman, and A. V. Tyler (1993). Movement dynamics in a fishery: Application of the ideal free distribution to spatial allocation of effort. Can. J. Fish. Aquat. Sci. 50, 323-333. Gotshall, 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 with impacts of salmon production. Bull. Amer. Meteorol. Soc. 78, 1069-1079. Martell, S., J. Boutillier, 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 Variability and Ecosystem Response in the Northeast Pacific. Science 10-Jul 281, 210-217. Myers, R. A., K. G. Bowen, and N. J. Barrowman (1999). Maximum reproductive rate of fish at low populations sizes. Can. J. Fish. Aquat. Sci. 56, 2404-2419. Mysak, L. (1986). El Nino, interannual variability 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 University Press. Quinn, T., C. Turnbull, and C. Fu (1998). A Length-Based Population Model 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. Sullivan, and C. I. Zhang (Eds.), Fishery Stock Assessment Models, pp. 531-556. Alaska Sea Grant College Program Report No. AK-SG-98-01, University of Alaska Fairbanks, 1998. Rexstad, E. A. and E. Pikitch (1986). Stomach contents and food consumption esti mates of pacific hake, Merluccius productus. Fish. Bull. 84, 947-956. Rinaldi, S. and M. Gatto (1978). On the determination of a commercial fishery pro duction model. Ecological Modelling 8, 165-172. Robinson, C. L. K. and D. M. Ware (1994). Modelling pelagic fish and plankton trophodynamics off Southwestern Vancouver Island, British 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. Miller (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 minimization. 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. Hill, and M. Foreman (2001). Horizontal transport of marine ogranisms resulting from interactions between diel vertical migration and tidal currents off the west caost of Vancouver Island. Can. J. Fish. Aquat. Sci 58, 736-748. Stefansson, G. (1998). Comparing different information sources in a multispecies 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 Model 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 mortality 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 with simulated data. Relative abundance data were assumed to have log-normal errors, and length frequency data were taken by randomly sampling from the true multinomial distribution in each year. An example of a length frequency sample over time is shown in 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 ° Jl—"PH—i—r5*!—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 ii—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 II i 10 15 20 25 30 35 40 Figure A.l: Example of an observed length frequency distribution and predicted frequen cies from stock assessment model. Observed frequencies were drawn at random from the true multinomial 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 mortality 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 with no observation errors, and estimated values with out penalty functions on recruitment anomalies, and natural mortality 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 will 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 distribution of all estimated parameter to look for systematic bias (Figure A.3). When length frequency samples are small and uninformative about changes in mortal ity rates, auxiliary information about relative trends in fishing mortality rates generally improves estimates of true population parameters, especially changes in 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 will over-estimate annual recruitment and under-estimate natu ral mortality, and vice versa when larger shrimp are over represented in the sample. As 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 (VJt 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 in 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: Estimated growth parameters from length and weight-at-age data from Cass et al,(1990). Description Parameter Value Asymptotic 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). This 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 lingcod biomass. Fishing mortality 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 mortality rates. The early CPUE data indicate a strong rebuilding of the stock from 1950 to 1968, and here I made the very dangerous assumption that the qualified CPUE 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 mortality 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. During the mid 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 Cl 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^. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items