Hierarchical stock assessment methods improve management performance in multi-species, data-limited fisheries
HHierarchical stock assessment methods improve management performance inmulti-species, data-limited fisheries.
Samuel D. N. Johnson a, ∗ , Sean P. Cox a a School of Resource and Environmental Management, Simon Fraser University, 8888 University Drive, BC, V5K 1S6, Canada
Abstract
Managers of multi-species fisheries aim to balance harvest of multiple interacting species that vary in abundanceand productivity. Technical interactions are a defining characteristic of multi-species fisheries, and are often ignoredwhen setting catch limits or target exploitation rates, despite their effect on the ability to balance multi-speciesharvest. Balancing harvest is more difficult when stock assessments are impossible to perform due to data-limitations,involving inter alia short time series of fishery independent observations, poor catch monitoring, limited biologicaldata, and spatial mismatching of harvest decisions to system structure. In this paper, we asked whether usinghierarchical, data-limited stock assessment methods to set catch limits improved management performance indata-limited, multi-species fisheries. Management performance of five alternative stock assessment methods wasevaluated by using them to set harvest levels targeting multi-species maximum yield in a multi-species flatfish fishery,including single-species and hierarchical multi-species models, and methods that pooled data across species andspatial strata, with catch outcomes of each method under three data scenarios compared to catch under an omniscientmanager simulation. Operating models included technical interactions between species intended to produce chokeeffects often observed in output controlled multi-species fisheries. Hierarchical multi-species models outperformed allother methods under data-poor and data-moderate scenarios, and outperformed single-species models under thedata-rich scenario. Hierarchical models were least sensitive to prior precision, sometimes improving in performancewhen prior precision was reduced. Choke effects were found to both positive and negative effects, sometimes leadingto underfishing of non-choke species, but at other times preventing overfishing of non-choke species. We highlightthe importance of including technical interactions in multi-species assessment models and management objectives,how choke species can indicate mismatches between management objectives and system dynamics, and recommendhierarchical multi-species models for multi-species fishery management systems.
Keywords:
Data-limited fisheries management; Multi-species fisheries management; technical interactions;Management Strategy Evaluation; hierarchical Multi-species Stock Assessment; Choke species. ∗ Corresponding author
Email addresses: [email protected] (Samuel D. N. Johnson), [email protected] (Sean P. Cox)
Preprint submitted to Fisheries Research June 26, 2020 a r X i v : . [ q - b i o . P E ] J un . Introduction Managers of multi-species fisheries aim to balance harvest of multiple interacting target and non-target speciesthat vary in abundance and productivity. Among-species variation in productivity implies variation in single-speciesoptimal harvest rates, and, therefore, differential responses to exploitation. Single-species optimal harvest rates (e.g.,the harvest rate associated with maximum sustainable yield) typically ignore both multi-species trophic interactionsthat influence individual demographic rates (Gislason, 1999; Collie and Gislason, 2001), and technical interactionsthat make it virtually impossible to simultaneously achieve the optimal harvest rate for all species (Pikitch, 1987).Technical interactions among species that co-occur in non-selective fishing gear are a defining characteristic ofmulti-species fisheries (Pikitch, 1987; Punt et al., 2002) and, therefore, play a central role in multi-species fisheriesmanagement outcomes for individual species (Ono et al., 2017; Kempf et al., 2016). Catch limits set for individualspecies without considering technical interactions subsequently lead to sub-optimal fishery outcomes (Ono et al., 2017;Punt et al., 2011a, 2020). For example, under-utilization of catch limits could occur when technically interactingquota species are caught at different rates (i.e., catchability) by a common gear, leading to a choke constraint inwhich one species quota is filled before the others (Baudron and Fernandes, 2015). Choke constraints are generallyconsidered negative outcomes for multi-species fishery performance, because they reduce harvester profitability asincreasingly rare quota for choke species may limit access to fishing grounds, as well as driving quota costs above thelanded value of the choke species (Mortensen et al., 2018).Setting catch limits for individual species in any fishery usually requires an estimate of species abundance, whichcontinues to be a central challenge of fisheries stock assessment (Hilborn and Walters, 1992; Quinn, 2003; Maunderand Piner, 2015), especially when species data are insufficient for stock assessment models. Where such datalimitations exist, data pooling is sometimes used to extend stock assessments to complexes of similar, interactingstocks of fish (Appeldoorn, 1996). A couple of examples include pooling data for a single species across multiplespatial strata when finer scale data are unavailable or when fish are believed to move between areas at a sufficientlyhigh rate (Benson et al., 2015; Punt et al., 2018), and pooling data for multiple species of the same taxonomic groupwithin an area when data are insufficient for individual species or during development of new fisheries (DeMartini,2019). Data-pooled estimates of productivity represent means across the species complex; therefore the general costof such pooling approaches is that catch limits will tend to overfish unproductive species and underfish productiveones (Gaichas et al., 2012). Hierarchical stock assessment models may help reduce these problems by representing acompromise between uninformative single-species single-species approaches and problematic data-pooling approaches.In particular, hierarchical models allow data sharing between multiple species within a complex via hierarchicalshrinkage prior distributions on model parameters that represent biological similarities or technical interactions(Thorson et al., 2015). Shared priors shrink species-specific parameters towards an estimated complex level meanvalue, improving model convergence for data-poor species and basing estimates more on available data than strong apriori assumptions like fixed (i.e. known) parameters, identical parameters among species, or strongly informativepriors (Jiao et al., 2009, 2011; Punt et al., 2011b).Although hierarchical stock assessments produce better estimates of species biomass and productivity thansingle-species methods in data-limited contexts, it remains unclear whether such improved statistical performancetranslates into better management outcomes (Johnson and Cox, 2018). On one hand, improving assessment modelstatistical performance may sacrifice future fishery benefits of adaptive learning (Walters, 1986), as large assessmenterrors caused by poor quality data in the present may generate high contrast data that are more informative in thefuture. Alternatively, adaptive learning may fail to occur in a reasonable time in the presence of large errors, leadingto stock collapse or harvesters divesting from the fishery because of over- or under-exploitation. Improving the rateof adaptive learning may be possible via hierarchical modeling of spatially replicated groundfish stocks, becausethe shared priors take lessons learned from responses to disturbances in one area and use them to improve stockassessments in other areas (Collie and Walters, 1991).In this paper, we investigated whether hierarchical stock assessment models improved management performance ina simulated multi-species, data-limited fishery. The simulated fishery was modeled on a complex of Dover sole, Englishsole, and southern rock sole - three species of right-eyed flounders (
Pleuronectidae spp. ) - Dover sole (
Microstomuspacificus ), English sole (
Parophrys vetulus ), and southern rock sol (
Lepidopsetta bilineata ) - fished off the coast ofBritish Columbia, Canada, in three spatial management areas. Closed-loop feedback simulation was used to estimatefishery outcomes when catch limits were set based on estimates of biomass from single-species, data-pooling, andhierarchical state-space surplus production models under three data scenarios. Assessment models were either fitto species specific data as single-species or hierarchical multi-species models, or fit to data pooled spatially across2anagement units, pooled across species within a spatial management unit, or totally aggregated across both speciesand spatial management units. Management performance of each assessment approach was measured as cumulativeabsolute loss in catch, measured against optimal catch trajectories generated by an omniscient manager, who couldset annual effort to maximize total multi-species/multi-stock complex yield given perfect knowledge of all futurerecruitments (Walters, 1998; Martell et al., 2008).
2. Methods
British Columbia’s multi-species complex of right-eyed flounders is a technically interacting group of flatfishesmanaged across the whole BC coast. Although there are several right-eyed flounders in BC waters, we focus on Doversole, English sole, and Rock sole. Taken together, these species comprise a multi-stock complex (DER complex),managed as part of the BC multi-species groundfish fishery.The DER complex is managed in three spatially distinct stock areas (Figure 1) (Fisheries and Oceans, Canada,2015). From north to south, the first stock area - Hecate Strait/Haida Gwaii (HSHG) - corresponds to BC groundfishmanagement area 5CDE, extending from Dixon Entrance and north of Haida Gwaii, south through Hecate Strait.The second stock area - Queen Charlotte Sound (QCS) - corresponds to BC groundfish management area 5AB,extending from the southern tip of Haida Gwaii to the northern tip of Vancouver Island. Finally, the third area -West Coast of Vancouver Island (WCVI) - corresponds to groundfish management area 3CD, extending from thenorthern tip of Vancouver Island south to Juan de Fuca Strait.Each area had two commercial catch rate series and at least one survey biomass index time series for eachspecies (Table 1). The two commercial fishery catch-per-unit-effort series span 1976 - 2016, split between a historicaltrawl fishery (1976 - 1995) and a modern trawl fishery (1996 - 2016), with the split corresponding to pre- andpost-implementation of an at-sea-observer program, respectively. The fishery independent trawl survey biomassindices were the biennial Hecate Strait Assemblage survey in the HSHG area (1984 - 2002), and the GroundfishMulti-species Synoptic Survey, with separate biennial legs in all three stock areas (2003 - 2016).
Closed-loop simulation is often used to evaluate proposed feedback management systems. In fisheries, closed-loopsimulation evaluates fishery management system components, such as stock assessment models or harvest decisionrules, by simulating repeated applications of these components, along with propagating realistic errors in monitoringdata, stock assessment model outputs, and harvest advice (de la Mare, 1998; Cox and Kronlund, 2008; Cox et al.,2013). At the end of the simulation, pred-defined performance metrics are used to determine the relative performanceof the system components being tested.Our closed-loop simulation framework included a stochastic operating, representing the stock and fishery dynamics,as well as an observation model, and an assessment model component that estimated stock biomass from simulatedobservations and fishery catches. The operating model simulated population dynamics of a spatially stratified,multi-species flatfish complex in response to a multi-species trawl fishery in each of the three stock areas. Althoughtotal fishing effort was not restricted across the entire area, effort in each individual area was allocated such that nospecies-/area-specific catch exceeded the species-/area-specific total allowable catch (TAC). Within an area, fishingeffort was allowed to increase until at least one species-/area-specific TAC was fully caught.The simulation projected population dynamics for all nine stocks forward in time for 32 years, with annualsimulated assessments, harvest decisions, and catch removed from the population from 2017 - 2048. The followingfour steps summarise the closed-loop simulation procedure for each projection year t :1. Update stochastic population dynamics and generate new realized catch C s,p,t ( E p,t ) in each area from effort(eq 1);2. Generate new observation data I s,p,f,t (eqs 4-11)3. Apply an assessment model (defined below) to estimate the spawning biomass for the upcoming year ˆ B s,p,t +1 (eq 3);4. Apply a harvest rate to generate a total allowable catch T AC s,p,t +1 (eqs 12 - 16);5. Allocate effort E p,t +1 to fully realise at least one TAC in each stock area (eq 2).3 .2.1. Operating Model The operating model (OM) was a multi-species, multi-stock age- and sex-structured population dynamics model(Appendix A). Population life-history parameters for the operating model were estimated by fitting a hierarchicalage- and sex-structured model to data from the real DER complex (Johnson, Cox and Knowler, in Prep).Fishing mortality for individual stocks was scaled to commercial trawl effort via species-specific catchabilityparameters, i.e., F s,p,t = q Fs,p · E p,t , (1)where F s,p,t is the fishing mortality rate applied to species s in stock-area p by fleet f in year t , and q Fs,p is thecommercial catchability coefficient scaling trawl effort E p,t in area p to fishing mortality (Table 2) [Johnson andCox, in prep].The above effort model implies a multi-species maximum yield for the trawl fishery in each area, depending onindividual species productivity, commercial trawl catchability, and relative catchabilities among species in the area(Pikitch, 1987; Punt et al., 2011a), i.e., M SY
MS,p = max E Y MS (cid:0) E, q F ,p , q F ,p , q F ,p , Q p (cid:1) , where E is the total commercial trawl fishing effort, q s,p is the commerical trawl catchability coefficient scalingfishing effort to fishing mortality for species s in area p , and Q p is the set of life-history parameters for all DERcomplex species in area p . The function Y MS is the total multi-species yield in area p as a function of fishing effort E p (Figure 2, heavy black lines), defined as the sum of species-specific equilibrium yields Y MS,p (cid:0) E p , q F ,p , q F ,p , q F ,p , Q p (cid:1) = X s =1 Y SS,s,p ( E p , q s,p , Q s,p ) , where Q s,p is the subset of Q p containing life history parameters for species s in area p , and the function Y SS,s,p isthe traditional single-species yield curve expressed as a function of effort rather than fishing mortality (Figure 2,coloured lines).The level of fishing effort producing multi-species maximum yield
M SY
MS,p in area p is defined as E MSY,p = arg max E Y MS (cid:0) E, q F ,p , q F ,p , q F ,p , Q p (cid:1) , with arguments defined above (Figure 2, vertical dashed line). The fishing effort level E MSY,p has an associatedequilibrium biomass B MSY,MS,s,p and yield
M SY
MS,s,p = Y SS,s,p ( E MSY,p , q
Fs,p , Q s,p )for each species s (Figure 2, three lower horizontal dashed lines showing individual species yield), which differ fromthe single-species B MSY and
M SY when commercial catchability scalars differ between species.For simulated fishing in the projection time period, we allocated maximum fishing effort to each area, defined asthe effort required to fully utilize the TAC of at least one species, while never exceeding the TAC of any one of thethree species, i.e., E p,t +1 = max { E | C s,p,t +1 ( E ) ≤ T AC s,p,t +1 ∀ s } , (2)where C s,p,t +1 ( E ) is the catch of species s when effort E is applied in area p in year t + 1 (the method for determiningTACs is explained below). We used maximum effort instead of an explicit effort dynamics model because the formerwas adequate for understanding the relative management performance of the assessment models that we tested.Furthermore, the maximum effort model may be simplified over reality, but it did an adequate job of capturingchoke effects observed in the real BC groundfish fishery system. At each time step t , simulated annual assessments were used to estimate the expected future spawning biomassestimate ˆ B t +1 via a state-space Schaefer production model (Schaefer, 1957; Punt et al., 2002), modified to betterapproximate the biomass-yield relationship underlying the age-/sex-structured operating model (Pella and Tomlinson,1969; Winker et al., 2018). We extended the Johnson and Cox (2018) hierarchical state-space model to fit a4ulti-species, spatially stratified complex, as well as fit a single-stock model to data from individual or data-pooledstocks, via the biomass equation B s,p,t +1 = " B s,p,t + U MSY,s,p · m s,p m s,p − · B s,p,t · − m s,p · (cid:18) B s,p,t B MSY,s,p (cid:19) m s,p − ! − C s,p,t e ζ s,p,t , (3)where B MSY (optimal biomass) and U MSY (optimal harvest rate) are the leading biological model parameters, m s,p is the Pella-Tomlinson parameter controlling skew in the biomass/yield relationship (derived from operating modelyield curve), and ζ s,p,t are annual process error deviations.In total, we defined the following five potential assessment model configurations:1. Total Aggregation (1 stock);2. Species Pooling (3 stocks, independent);3. Spatial Pooling (3 species, independent);4. Single-stock (9 stocks, independent);5. Hierarchical Multi-stock (9 stocks, sharing information);where the number of management units is shown in parentheses (delete subscripts in eq. 7 for species or stock-areaas appropriate, e.g., Spatial Pooling models have no p subscript).Prior distributions on optimal biomass B MSY,s,p , optimal harvest rate U MSY,s,p , catchability q s,p,f , and processerror deviations ζ s,p,t were defined for each assessment model (Table 3). For all non-hierarchical models, prior meanvalues were derived from the operating model biomass/yield relationship, with prior mean B MSY values µ B MSY,s,p = B MSY,s,p ,µ B MSY,s = X p B MSY,s,p ,µ B MSY,p = X s B MSY,s,p ,µ B MSY = X s,p B MSY,s,p , for single-stock, spatial pooled, species pooled, and totally aggregated configurations, where B MSY,s,p is theequilibrium biomass at which single-species maximum yield is achieved. Similarly, log-normal prior mean log U MSY values were given by µ log U MSY,s,p = log(
M SY s,p /B MSY,s,p ) ,µ log U MSY,s = log( X p M SY s,p / X p B MSY,s,p ) ,µ log U MSY,p = log( X s M SY s,p / X s B MSY,s,p ) ,µ log U MSY = log( X s,p M SY s,p / X s,p B MSY,s,p ) , for single-stock, spatial-pooled, species-pooled, and totally aggregated configurations, respectively, where M SY s,p is the single-species maximum yield. Log-normal survey biomass and commercial CPUE catchability had priormeans set to OM values for the single-stock model, and prior means that were averaged OM values over the sourcesof pooled data under the data-pooled models. Prior standard deviations for catchability were set to s f = 1 . s f = 1 . s f = 0 . .
05 for allstocks, species, and model configurations.The hierarchical multi-stock model and the single-stock model used the same priors for B MSY and catchabilityfor the commercial CPUE and Hecate Strait Assemblage survey biomass indices. Log-normal hierarchical shrinkagepriors were applied to area-specific Synoptic survey catchabilities q s,p,f within each species, and in two levels to5 MSY,s,p (as a proxy for intrinsic growth rates) across the whole complex. Both hierarchical prior distributionshad fixed standard deviation parameters σ U MSY = τ q = 0 .
05 for differences between stocks within a species forcatchability q s,p,f and U MSY , and the same standard deviation for the log-normal prior on species average U MSY,s values. Hierarchical species mean catchability for the Synoptic survey q s,f had a log-normal hyperprior with mean m log q s,f and standard deviation s f = .
5, and DER complex mean U MSY had a log-normal hyperprior with mean m log U MSY and standard devation of s U MSY = . Time series of catch and biomass indices were simulated in the historical and projection periods for fittingassessment models. Simulated data matched the time periods of the real DER complex biomass index data sources,with the exception of the Synoptic trawl survey, which was also simulated in the projection, period (Table 1).Observation error precision for each data source is given in Table A.1 (Appendix A).Biomass indices for individual stocks depended on the nature of the index. Assemblage and Synoptic Surveybiomass indices were defined as trawlable biomass¯ I s,p,f,t = q s,p,f · B s,p,f,t , (4)where ¯ I s,p,f,t is the index without observation error, q s,p,f is the survey trawl efficiency for species s and stock-area p , and B s,p,f,t is the biomass of species s in stock-area p vulnerable to survey f in year t . Commercial CPUE indiceswere defined as ¯ I s,p,f,t = C s,p,f,t E p,f,t , (5)where C s,p,f,t was commercial catch by fleet f of species s from area p at time t , and E p,f,t was commercial fishingeffort in area p from fleet f at time t , with both catch and effort required to be positive.The method for generating pooled catch and biomass index data depended on the data type. Catch datawere pooled by summation, and the index data were pooled according to the stock-specific definition above. ForAssemblage and Synoptic survey biomass indices, pooled indices without observation error were defined as¯ I pooleds,f,t = X p I ( I s,p,f,t > · q s,p,f · B s,p,f,t (6)¯ I pooledp,f,t = X s I ( I s,p,f,t > · q s,p,f · B s,p,f,t (7)¯ I pooledf,t = X s,p I ( I s,p,f,t > · q s,p,f · B s,p,f,t (8)where ¯ I pooleds,f,t is the spatially pooled index for species s , ¯ I pooledp,f,t is a species pooled index for area p , ¯ I pooledf,t is thetotally aggregated index (all without error), and I ( I s,p,f,t >
0) is the indicator function that takes value 1 whensurvey leg for fleet f in area p was running in year t . Pooled commercial CPUE indices were simulated as¯ I pooleds,f,t = P p C s,p,f,t P p E p,f,t (9)¯ I pooledp,f,t = P s C s,p,f,t E p,f,t (10)¯ I pooledf,t = P s,p C s,p,f,t P p E p,f,t (11)where C s,p,f,t is catch, and E p,f,t is commercial trawl effort, with subscripts as defined above. Simulated harvest decision rules applied a constant target harvest rate to generate TACs from one-year aheadbiomass forecasts obtained from each assessment model, i.e.,
T AC s,p,t +1 = U s,p · ˆ B s,p,t +1 , (12)6here U s,p is the target harvest rate for species s in stock-area p , and ˆ B s,p,t +1 is the year t + 1 biomass forecast fromthe assessment model. Target harvest rates were derived from the multi-species maximum yield relationship via U s,p = M SY
MS,s,p B MSY,MS,s,p , (13)where M SY
MS,s,p and B MSY,MS,s,p are the yield and spawning biomass, respectively, associated with E MSY,MS,p for species s in area p .Finally, inter-annual increases in TAC were limited to 20% for all individual stocks as a practical precautionarymeasure, i.e., T AC s,p,t +1 = min { T AC s,p,t +1 , . ∗ T AC s,p,t } , (14)where T AC s,p,t +1 is the proposed TAC determined above, and T AC s,p,t is the previous year’s TAC. The limitedincrease reduces overfishing as a result of optimistic assessment errors, potentially requiring several years of increasesto realise target harvest rates as set by assessment model estimates. Conversely, no limit was applied to TACreductions, so that steep declines in biomass forecasts would lead to steep declines in catch. Note that this constrainton inter-annual TAC changes is reflected as a constraint on inter-annual changes in effort in the objective functionused in omniscient manager solutions described below.Pooled TACs were set analogously to the stock-specific case above, with pooled target harvest rates applied tobiomass projections from pooled assessments. For a spatially pooled assessment of species s , we defined the spatiallypooled target harvest rate as U s = P p M SY
MS,s,p P p B MSY,MS,s,p , (15)where the notation is as defined above. Species pooled harvest rates and total aggregation harvest rates were definedanalogously.Pooled TACs were split within an area or across spatial strata proportional to Synoptic trawl survey indices forthe individual stocks. For example, if the TAC for area p is set by a species pooled assessment, then the proposedTAC for species s is defined as T AC s,p,t +1 = ¯ I s,p, P s ¯ I s ,p, T AC p,t +1 , (16)where ¯ I s,p, is the 2-year running average of individual biomass indices from the Synoptic survey for species s inarea p . The 2-year average is used because the synoptic survey alternates the legs each year, so some individualstock indices are missing. We ran a total of 15 simulation experiments comprising five assessment models and three data quality scenarios.Simulations integrated over the stochastic processess by running at least 100 random replicates of each combination,where each simulation was initialized with the same set of random seeds to eliminate random effects amongcombinations of assessments and data scenarios. Simulations stopped when they either reached 100 replicates wheresimulated assessment models converged in at least 95% of projection time steps for all stocks, or when 140 randomreplicates were attempted. We defined convergence as a positive definite Hessian matrix and a maximum gradientcomponent less than 10 − in absolute value. A lower threshold of 95% convergence was chosen given that fittingmodels becomes more difficult as data quality is deliberately reduced, and a simulated assessment can not always betuned like a real assessment performed by a real-life analyst. Any OM/assessment model/stock combinations thatcould not reach 100 replicates meeting the 95% convergence criterion over a total of 140 attempts (approximately a70% success rate) were considered non-significant and were indicated as such in the results. The operating modelwas run for two Dover sole generations (32 years; Seber, 1997), because this species had the longest generation time.Operating model population dynamics and index data were identical among replicates for each stock during theoperating model historical period, except for the last few years near the end, where the operating model simulatesrecruitment process errors because the conditioning assessment was unable to estimate them. Simulated log-normalobservation and process errors in the projection were randomly drawn with the same standard deviations as the7rrors used in the historical period, and bias corrected so that asymptotic medians matched their expected values,i.e., I s,p,f,t = ¯ I s,p,f,t · exp( τ s,p,f · δ s,p,f,t − . τ s,p,f ) (17) R s,p,t = R s,p,t · exp( σ s,p · (cid:15) s,p,t − . σ s,p ) (18)where ¯ I s,p,f,t is the index without error defined above, τ s,p,f is the log-normal observation error standard deviation, δ s,p,f,t is the annual standard normal observation error residual, R s,p,t is the equilibrium recruitment from theBeverton-Holt stock-recruitment curve, σ s,p is the recruitment process error standard deviation, (cid:15) s,p,t is the annualstandard normal recruitment process error, and subscripts s, p, f, t are for species, stock, fleet and year, respectively.Error is added to biomass indices for pooled data independently of the error added to individual indices, i.e., I s,f,t = ¯ I s,f,t · exp( τ s,f · δ s,f,t − . τ s,f ) (19) I p,f,t = ¯ I p,f,t · exp( τ p,f · δ p,f,t − . τ p,f ) (20) I f,t = ¯ I f,t · exp( τ f · δ f,t − . τ f ) (21)where τ s,f , τ p,f , τ f were averaged over the components of the pooled index. The three data quality scenarios range from relatively data-rich to data-poor by successively removing commercialCPUE index series from the full set, i.e.,1. Data-
Rich : Historical CPUE, Modern CPUE, Assemblage survey, Synoptic survey;2. Data-
Mod erate: Modern CPUE, Assemblage survey, Synoptic survey;3. Data-
Poor : Assemblage survey, Synoptic survey.To improve convergence, the Hierarchical Multi-stock and Single-stock assessment models were initialised laterunder the
Mod and
Poor data scenarios, with the starting year of the assessments set to the first year with indexdata, which was 1984 in HSHG for both scenario, and 1997 or 2003 for other areas under the
Mod and
Poor scenarios, respectively.
Assessment model performance was measured against a simulatedomniscient fishery manager who is aware of all the future consequences of harvest decisions and is, therefore, able toadapt the management to meet specific quantitative objectives under any process error conditions (Walters, 1998).Omniscient manager solutions were used rather than equilibrium based metrics (Punt et al., 2016) because moststocks were in a healthy state in 2016 (i.e. above single-species B MSY , Table 2) and, therefore, the time-path offishery development was important (Walters et al., 1988).The omniscient manager was implemented as an optimisation of future fishing effort by area (Appendix B), withthe objective function defined as O = "X s,p − log( ¯ C s,p, · ) + P diff X p E p, · ! + P init X p E p, ! , (22)where − log ¯ C s,p, · is the negative log of total future catch for species p in area p over the projection period (equivalentto maximising catch). Penalty functions P (eq. B.1) were applied for annual changes in total effort across all threeareas being above 20% ( P diff ) to match the TAC smoother in stochastic experiments, and differences greater than10% between the last year of historical effort and the first year 2017 of simulated effort ( P init ).An omniscient manager solution was obtained for each stochastic trajectory in the stochcastic managementsimulations. Each replicate was run for 80 years to reduce end effects, such as transient dynamics at the beginningof the projection, or a lack of consequences for overfishing at the end of the projection.8 .3.2.2. Cumulative catch loss. For each stochastic trajectory, the cumulative absolute loss in catch was calculatedas (Walters, 1998): L s,p = T X t = T | C s,p,t,sim − C s,p,t,omni | , (23)where the C s,p,t, · values were commercial trawl catch for species s and stock p from stochastic simulations ( sim ) orthe omniscient manager simulation ( omni ) simulation. When repeated over all random seed values, the loss functionsgenerated a distribution of cumulative catch loss, which were then used to determine relative performance of eachassessment model under the three data scenarios. Cumulative loss was calculated for the ten year period T = 2028to T = 2037, chosen in the middle of the projection period because dynamics in the earlier time were dominated bythe smoothers on effort and catch for the omniscient manager and TACs, respectively, and after 2028 the omniscientmanager’s median effort has reached a stable state near the multi-species optimum. Biomass loss was also calculated,but inferences about assessment method performance in preliminary simulations were not qualitatively differentbetween the two metrics, so catch loss was reported only, with biomass risk for each stock indicated in reference to acritically overfished level, defined as 40% of their individual single-species B MSY value (DFO, 2006).Cumulative absolute catch loss was used to calculate the relative rank of each assessment model for eachspecies/stock/OM scenario combination, where lower loss ranked higher. Rank distributions across species andstocks were then used to calculate an average, minimum, and maximum rank of each assessment model to determinethe overall performance of each assessment model in each OM scenario.
Parameter prior distributions area a key feature of most contemporary stock assessment models, even in data-richcontexts. Moreover, prior distributions are a defining feature of hierachical multi-species stock assessment models.Therefore, we focused sensitivity analyses on fixed prior standard deviations for leading parameters B MSY and U MSY , and the hierarchical shrinkage prior SDs τ q and σ U MSY (Table 4). Analyses were run with 50 simulationreplicates for only data-rich and data-poor scenarios.
3. Results
As expected, the omniscient manager was able to achieve the theoretical multi-species optimal yield in thepresence of technical interactions (Figure 3, blue closed circle). Median biomass, catch, and fishing mortality reachthe equilibrium after a transition period of about 20 years. During the transitionary period, effort is slowly rampedup in each area from the end of the historical period, stabilising around area-specific E MSY after about 12 years(Figure 4, blue closed circles).The multi-species optimal effort for the complex results in overfishing Dover sole stocks between 53% and 92%of area-/single-species B MSY to increase fishery access to English and Rock sole. Despite this tendency towardoverfishing Dover sole, very few optimal solutions risk severe overfishing of Dover sole below 40% of B MSY , indicatingthat lost Dover sole yield from further overfishing relative to single-species optimal levels is not compensated byincreased English and Rock sole yield (Table 5). The probability of being critically overfished in the period 2028 -2037 was 0% for most stocks, with two Dover sole stocks having 3% (HSHG) and 1% (QCS).Although DER stocks begin the simulations in an overall healthy state, the omniscient manager reduced fishingeffort to near 0 in HSHG and QCS areas early in the projection period in some replicates (Figure 4, HSHG andQCS, 2016-2020). In these cases, anticipatory feedback control by the omniscient manager reduced fishing effort toavoid low spawning stock biomasses, thus and ensuring higher production in later time steps where recruitmentprocess error deviations were sustained at low levels.
Hierarchical Multi-stock assessment models ranked highest, on average, under the Poor and Mod data qualityscenarios (Table 6). Under the Poor data quality scenario, the Hierarchical model had the lowest cumulative absolute9atch loss, on average, and ranked in the top three assessment methods for all stocks. As data quantity increasedfor the Mod scenario, the average absolute rank of the Hierarchical model degraded slightly, but still came first inaverage rankings. It was only under the Rich scenario that the Hierarchical model average rank fell to fourth place,although it still ranked above the Single-stock model.Under the Rich data quality scenario, the Species Pooling and Total Aggregation methods ranked best, neverranking below fourth for all stocks (Table 6, Rich). As data was removed for the Mod data quality scenario, theaverage ranks of both Species Pooling and Total Aggregation methods degraded, placing them second and third.The Single-stock model had the worst average rank under all data-quality scenarios, with its worst performanceobserved under the Mod data quality scenario where it ranked 5th across all stocks.Cumulative catch loss distributions for the Hierarchical model and three data-pooling methods were fairly closelyclustered for the Rich and Mod data quality scenarios for most stocks (Figure 5). Only the Single-stock model wasqualitatively different, but this was not always the case (e.g., HSHG Dover, Figure 5). Relative performance of eachmethod across data quality scenarios was stock dependent, with no consistent trends across scenarios and stocks.Technical interactions had both a positive and negative effect on meeting the multi-species objectives andminimising catch loss. For example, over the course of the projection period, the Hierarchical Multi-stock model wasable to bring median biomass levels close to the multi-species optimal level B MSY,MS,s,p for six out of the nine stocksunder the Poor data quality scenario (Figure 6, blue dots). The three stocks that did not meet the multi-speciesoptimal biomass were all in the HSHG area, indicating that technical interactions may have been a factor in theinability to meet the target harvest rates, choking off TACs for some species in HSHG. For the HSHG area, catch ofDover and Rock soles are choked by English sole TAC (Figure 7, HSHG, unfilled TAC bars), which is low relative tothe target level because the Hierarchical assessment model viewed HSHG English sole as a larger and less productivestock than it truly was, resulting in a persistent negative bias in the stock assessment biomass. Furthermore, thereis evidence that a large perturbation is needed to improve the assessment for HSHG English sole and overcomethe choke effect, as the assessment model believes the TACs are appropriately scaled and biomass is approaching B MSY,MS . For QCS and WCVI areas, all stocks were able to meet the target harvest objectives and minimise lossdirectly because technical interactions protected against overfishing. In QCS, English sole had close to unbiasedassessments, and was able to limit catches for Dover and Rock soles despite large positive assessment errors for both(Figure 7, QCS). Similarly, WCVI English sole catch was choked by TACs for Dover and Rock sole, which both hadunbiased assessments (Figure 7, WCVI).
Distributions of catch and biomass relative to
M SY MS and B MSY,MS , respectively, were produced for the timeperiod 2028 ≤ t ≤ B MSY,SS towards the multi-species optimallevel B MSY,MS , asymptotically approaching it as biomass equilibrates to a harvest rate well below the single speciesoptimal rate U MSY,SS .Catch-biomass trade-offs in the HSHG area indicated that Hierarchical models performed more similar tothe omniscient manager under the Rich data quality scenario (Figure 8, HSHG). In the same area, the SpeciesPooling method (ranked highest under the Rich scenario) produced median biomass between 140% and 210% of themulti-species optimal level (Figure 8, HSHS, Purple diamond), despite the superior performance with respect tocatch loss for Dover and English soles (indicated by proximity to the horizontal median catch line).As mentioned above, several stocks were pushed into a critically overfished state by different methods. TheHierarchical Multi-stock model (QCS, Rich data quality scenario) and Spatial Pooling method (WCVI, Poor dataquality scenario) produced median biomass levels that were below 40% of B MSY,SS for two out of three species ineach area (Figure 8, QCS and WCVI). These methods avoided pushing the QCS English and WCVI Rock soles,respectively, into critically overfished states thanks largely to their relatively low commercial catchability scalars.Low catchability means that their multi-species optimal biomass B MSY,MS was well above the single-species optimalbiomass level B MSY,SS , providing plenty of room to overshoot the optimal biomass level while still avoiding thecritically overfished level. Even so, median biomass is still well outside the omniscient manager’s distribution forboth stocks, with catches far above the omniscient manager’s median catch, ranging between 150% and 180% ofmaximum yield under multi-species E MSY . Finally, Total Aggregation and Spatial Pooling assessment methodshad higher than 10% probability of being below 40% of B MSY,SS under the Poor data scenario for (i) HSHG Doversole (37 % Total Aggregation, 42% Species Pooling), (ii) QCS Dover sole (37 % Total Aggregation, 11% SpeciesPooling), (iii) HSHG English (12 % Total Aggregation, 11 % Species Pooling), and (iv) QCS Rock sole (11% TotalAggregation).
We summarised average model sensitivities by fitting linear regressions to the distributions of median cumulativeloss. To remove the effect of absolute catch scales on the regressions, median loss distributions were standardisedacross assessment models, stratified by species, stock-area, and data scenario. Regressions with positive slopes haveincreasing catch loss with increasing prior uncertainty, and negative or zero slopes indicate a decrease or no changein catch loss with increasing prior uncertainty.The Hierarchical Multi-stock model was least sensitive to increasing B MSY prior CVs in both Rich and Poor dataquality scenarios. Cumulative catch loss actually decreased on average with increasing B MSY prior CVs when catchlimits were set by the Hierarchical model under the Poor data quality scenario, decreasing by about 0.3 standarddeviations over the range of CVs tested (Figure 9, left column, Poor). While cumulative catch loss from Hierarchicalmodels increased on average with increasing B MSY prior CVs under the Rich data quality scenario, the slope of theregression line was the lowest among methods (equal with the Species Pooling and Total Aggregation methods),implying the least sensitivity to prior specifications.Cumulative catch loss under the Hierarchical model was also insensitive to the complex mean U MSY hyperpriorstandard deviation, with zero slope regression lines under both Rich and Poor data quality scenarios over the rangeof tested log-normal prior SDs (Figure 9, middle column). The lack of sensitivity of the Hierarchical model to the U MSY hyperprior was partially because target harvest rates were based on theoretically optimal levels, and notestimated from the assessment model; however, this implies that the biomass estimates were also insensitive to thehyperprior standard deviation, which may be because the effect of the hyperprior was soaked up by the two levels ofshrinkage priors between the complex mean and stock-specific estimates of U MSY .Hierarchical Multi-stock assessment models were more able to suitably scale catch limits to operating modelbiomass under relaxed catchability and U MSY hierarchical shrinkage priors. As σ U MSY and τ q increased from 0.1to 0.5, average catch loss dropped by about 0.4 standard deviations under the Rich data quality scenario, and 0.8standard deviations for the Poor data quality scenario (Figure 9, third column). Improvements stemmed from achange in assessment model errors in the HSHG area, in particular (Figures C.1, C.2, Appendix C). The change inassessment performance switched the choke species from English to Rock sole in HSHG, and while most assessmentswere biased and had strong retrospective patterns, the combination of higher catch limits and favourable technical11nteractions produced catches that were more closely aligned with the omniscient manager’s and the multi-speciesmaximum yield.As expected, cumulative catch loss increased with decreasing prior precision under all other assessment methods.The Single-stock model was most sensitive to B MSY
CVs under both Rich and Poor data quality scenarios, and to U MSY log-normal SDs under the Poor data quality scenario, indicating that the single-stock model required themost prior knowledge for scaling biomass estimates correctly. The Total Aggregation and Spatial Pooling methodswere most sensitive to the U MSY prior standard deviation, which may indicate a that the aggregate prior meanvalues we calculated for those assessment models were not well supported by the aggregated data.
4. Discussion
In this paper, we demonstrated that hierarchical stock assessment models may improve management performancein a data-limited, multi-species flatfish fishery. When available data quality was moderate or poor (indicatedhere by time-series length), biomass estimates from hierarchical stock assessment models resulted in catches thatwere closer to an omniscient manager’s optimal reference series compared to catch limits derived from single-stockand data-pooling assessment methods. Under high data quality scenarios, data-pooling methods outperformedhierarchical models, but the latter still outperformed single-stock assessment methods. Improved performancerelative to single-stock models under lower data quality conditions is consistent with our previous study, wherestatistical performance of hierarchical multi-stock assessments improved with decreasing data quantity and quality(Johnson and Cox, 2018). This suggests that hierarchical assessment methods would generally be a better approachthan conventional single-species methods under typical fisheries data quality conditions.Our results arise from models that are necessarily a simplification of the real stock-management system. Theharvest rules applied to DER complex species were relatively simple and may require more detail or complexityfor practical applications. First, the harvest rules were all constant target harvest rates, which do not includeprecautionary “ramping-down” of catch towards a limit biomass level (DFO, 2006; Cox et al., 2013). Including aramped harvest rule may have reduced the probability of some stocks being critically overfished in some cases, butprobably at some further cost of choke effects. Second, catch limits for the simulated DER complex were set based onfixed target harvest rates that were derived a priori from multi-species yield curves, and not estimated as part of theassessment models that we tested. Incorporating multi-species yield curve calculations based on assessment modelouput into the harvest decision would be simple to do, but would require either a model of increased complexity tolink fishing effort to single-species yield, or an extra assumption linking effort to surplus production model yieldcalculations, which would likely increase assessment model errors. Finally, the TAC allocation model for data-pooledmethods was only one example from a large set of potential options. Understanding the relative risks of data-poolingwould require testing alternative allocation methods, which was beyond the scope of this paper.We only considered multi-species technical interactions, which although an important part of exploited systemdynamics, are not the entire story. Although there is limited evidence for ecological interactions among DERcomplex species (Pikitch, 1987; Wakefield, 1984), what does exist may influence the multi-species yield relationshipwith fishing effort or, as with technical interactions, inhibit the ability of the management system to meet targetcatch levels. For example, individual survival or growth may change in response to varied fishing pressure due tounmodeled linkages (Collie and Gislason, 2001). Yet, including such ecological interactions would imply a highlydata rich scenario, which is counter to our focus on data-limited, multi-species fisheries. Furthermore, accountingfor potential ecological interactions would require multiple OMs to test performance against a range of plausiblehypotheses since ecological uncertainties are much broader in complexity and scope than technical interactions alone.Nevertheless, future work combining technical interactions with minimum realistic models for ecological interactionscould help determine the extent to which assessment approaches affect these more complex multi-species fisheriesoutcomes (Punt and Butterworth, 1995). For example, while diet overlap between the three species is small off thecoast of Oregon, the major Rock sole prey was recently settled pleuronectiform fishes, which may include Dover andEnglish sole young and therefore shift the complex equilibrium as fishing pressure is applied, reducing predationmortality for Dover and English sole and prey availability for Rock sole (Wakefield, 1984; Collie and Gislason, 2001).Our effort model applied to the DER complex was also a simplification of reality, where effort was limited onlyby the TACs in each area. Limiting by TACs was intended to reflect the management of the real BC groundfishfishery in which harvester decisions drive TAC utilisation among target species (via increasing catchability; Puntet al., 2011a), and non-target or choke species (via decreasing catchability; Branch and Hilborn, 2008). Changing12atchability for targeting or avoidance could be simulated as a random walk in the projections, with correlationand variance based on the historical period, or perhaps simulated via some economic sub-model that accounted forex-vessel prices and variable fishing costs. These economic factors could affect targeting and avoidance behaviouramong species (Punt et al., 2011a, 2020), as well as effort allocation among stock-areas (Hilborn and Walters,1987; Walters and Bonfil, 1999); however, it is not clear that our median results would be significantly differentgiven the potential magnitude of assessment model errors in data-limited scenarios. Impacts of a detailed effortdynamics sub-model would probably be more important in more extreme data-limited scenarios that relied solelyon fishery CPUE as an index of abundance. In fact, it would be interesting to determine whether the hierarchicalinformation-sharing approach would exacerbate assessment model errors in such a (common) context where fisheryCPUE is the main abundance index.Despite the limitations above, our results indicate that even in fisheries with long time series of catch and effortdata, hierarchical multi-species assessment models may be preferable over typical single-species methods. The poorperformance of the single-species models in all scenarios highlights the difference between data-rich (i.e., a higherquantity of data) and information-rich (i.e., data with higher statistical power) fisheries. The data-rich scenariodiffered from data-moderate and data-poor scenarios by the inclusion of a historical series of fishery dependent CPUE,which was quite noisy and subject to the effects of changing harvester behaviour like targeting (variable catchability),and therefore, additional historical CPUE data had little effect on cumulative catch loss under the single-speciesmodels. In contrast, the data-pooling procedures all ranked higher than single-species and multi-species modelsunder the data-rich scenario, as they were able to leverage additional statistical power from the historical CPUE byeffectively increasing the sample size through data aggregation. The superior performance of the hierarchical modelover the single-species model under the data-rich scenario indicate that shared priors partially compensate for lowstatistical power, but not as much as data-pooling.Although the data-pooled methods performed better under the rich scenario, they were more sensitive to priors,so those results may be optimistically biased in all scenarios. Data-pooled observation errors were simulatedas independent of the observation errors in the component indices, using the average standard deviation of thecomponents. If aggregate indices pooled errors from each component index, then the resulting observation errorvariance would be additive in the components, especially if those errors were positively correlated, which may be thecase under a common survey or fishery.A dual effect of control was observed under catch limits set by the hierarchical models in the data-poor scenario, inwhich lower contrast in assessment model outputs reduced the statistical power of assessment model data, sacrificinglong term adaptability of the management system in favour of short term stability in catch. Moreover, there waslimited evidence that adaptive learning was facilitated through the shared hierarchical prior distributions (Collieand Walters, 1991). Uninformative catch series resulting from a lack of large perturbations (sometimes caused bylarge catch errors) resulted in hierarchical multi-species assessment model estimates that viewed HSHG Englishsole as a larger and less productive stock, causing negatively biased assessment estimates of biomass to approachthe multi-species B MSY,MS value for that stock, and indicating that the assessments believed the TACs wereappropriately scaled to the target harvest rates. Relaxing hierarchical prior standard deviations in the sensitivityanalysis removed this persistent bias and all stocks approached biomass levels associated with multi-species maximumyield; however, information sharing was not wholly responsible for the improved behaviour, which appeared to alsorely on a favourable combination of assessment errors and technical interactions as in the other stock areas. Adaptivelearning catalysts like large stock perturbations may have been limited by a low recruitment process error standarddeviation of σ s,p = . .1. Conclusion Hierarchical multi-species assessment models can out-perform single-species assessment models in meetingmulti-species harvest objectives across data-rich, data-moderate, and data-poor scenarios. As expected, biomassestimation performance of hierarchical models improved relative to other methods as data quantity was reduced,and - as hoped - this translated into improved management performance across the multi-species flatfish fishery.We therefore recommend that assessment and management of multi-species fisheries include technical interactionswhen designing harvest strategies and management procedures aimed at achieving strategic objectives. Otherwise,error-prone single-species approaches may give a misleading picture of the expected performance of multi-speciesfishery management.
5. Acknowledgements
Funding for this research was provided by a Mitacs Cluster Grant to S.P. Cox in collaboration with the CanadianGroundfish Research and Conservation Society, Wild Canadian Sablefish, and the Pacific Halibut ManagementAssociation. We thank S. Anderson and M. Surry at the Fisheries and Oceans, Canada Pacific Biological Station forfulfilling data requests. Further support for S.P.C. and S.D.N.J. wass provided by an NSERC Discovery Grant toS.P. Cox. 14 . Tables
Table 1: Fishery dependent and independent indices of biomass available for stock assessments of DER complex species.
Series Extent DescriptionHistorical Commercial 1976 - 1995 Historical period commercial CPUEModern Commercial 1996 - 2016 Modern period commercial CPUEHS Assemblage 1984 - 2002 Hecate Strait Assemblage trawl survey biomass index, biennialSynoptic 2003 - 2016 Multi-species Synoptic trawl survey biomass index, biennial able 2: Unfished biomas B , single-species MSY based reference points B MSY,SS and U MSY,SS , stock status as absolute biomass in2016 B , depletion relative to single-species optimal biomass B /B MSY,SS , and commercial trawl catchability scalar q F for allDER complex stocks in 2016. Biomass quantities are given in kilotonnes, and depletion and harvest rates are unitless. Single-species Reference Points Stock Status Comm. CatchabilityStock B B MSY,SS U MSY,SS B B /B MSY,SS q F Dover sole
HSHG 17.50 6.06 0.17 6.06 1.00 0.022QCS 5.98 1.99 0.15 1.74 0.88 0.017WCVI 15.20 5.03 0.19 4.05 0.81 0.039
English sole
HSHG 9.98 3.42 0.31 5.95 1.74 0.024QCS 0.57 0.19 0.29 0.38 2.04 0.011WCVI 0.90 0.30 0.29 0.38 1.26 0.042
Rock sole
HSHG 16.35 5.75 0.23 8.40 1.46 0.009QCS 5.57 1.90 0.21 1.63 0.86 0.014WCVI 1.72 0.63 0.23 0.68 1.09 0.009 able 3: Prior distributions used for each assessment model. For single stock and hierarchical multi-stock assesssment models, s f is 1 forcommercial indices, and 0.5 for survey indices. Prior
Totally aggregated AM B MSY ∼ N ( µ B MSY , . · µ B MSY )log U MSY ∼ N ( µ log U MSY , . q f ∼ N ( µ log q f , . Species-Pooled AM B MSY,p ∼ N ( µ B MSY,p , . · µ B MSY,p )log U MSY,p ∼ N ( µ log U MSY,p , . q p,f ∼ N ( µ log q p,f , . Spatial-Pooled AM B MSY,s ∼ N ( µ B MSY,s , . · µ B MSY,s )log U MSY,s ∼ N ( µ log U MSY,s , . q s,f ∼ N ( µ log q s,f , . Single-stock AM B MSY,s,p ∼ N ( µ B MSY , . · µ B MSY )log U MSY,s,p ∼ N ( µ log U MSY , . q s,p,f ∼ N ( µ log q s,p,f , s f ) Hierarchical multi-stock AM B MSY,s,p ∼ N ( µ B MSY , . · µ B MSY )log U MSY,s,p ∼ N (log U MSY,s , σ U MSY )log U MSY,s ∼ N (log U MSY , σ U MSY )log U MSY ∼ N ( µ log U MSY , . q s,p,f ∼ N (log q s,f , τ q )log q s,f ∼ N ( µ log q s,f , s f ) All AMs ζ s,p,t ∼ N (0 , . able 4: Summary of sensitivity analyses, showing the total number of experiments, the factor being varied, the levels of that factor, andthe data scenarios and AMs included in the analysis. N Factor Levels Scenarios AMs30 B M SY prior CV 0 . , . , . U M SY prior SD 0 . , . , . τ q , σ U MSY . , . , . able 5: Probability of being overfished and experiencing overfishing with respect to single-species reference points, and catching less than the historical miniumum during thetime period 2028 - 2037 for all nine DER complex stocks when managed by the omniscient manager. Prob. of being overfished Prob. of overfishing Prob. of low catchStock P ( B t < . B MSY,SS ) P ( B t < . B MSY,SS ) P ( C t > M SY SS ) P ( F t > F MSY,SS ) P ( C t < min C ) Dover sole
HSHG 0.03 0.70 0.74 0.97 0QCS 0.01 0.46 0.66 0.86 0WCVI 0.00 0.23 0.54 0.63 0
English sole
HSHG 0.00 0.12 0.64 0.49 0QCS 0.00 0.00 0.24 0.00 0WCVI 0.00 0.03 0.53 0.23 0
Rock sole
HSHG 0.00 0.01 0.50 0.06 0QCS 0.00 0.04 0.53 0.29 0WCVI 0.00 0.00 0.01 0.00 0 able 6: Summary of AM rankings with respect to cumulative absolute catch loss between 2026 and 2035 under each scenario. The rankcolumn shows the average rank of the AMs over all species/areas, and the range of that distribution in parentheses. AMs are ordered bymean rank within a scenario. AM Average Loss Rank (range)
Rich
Species Pooling 1.67 (1, 3)Total Aggregation 2.33 (1, 4)Spatial Pooling 3.00 (2, 4)Hierarchical Multi-stock 3.22 (1, 5)Single Stock 4.78 (4, 5)
Mod
Hierarchical Multi-stock 2.00 (1, 4)Species Pooling 2.11 (1, 4)Total Aggregation 2.44 (1, 4)Spatial Pooling 3.44 (2, 4)Single Stock 5.00 (5, 5)
Poor
Hierarchical Multi-stock 1.44 (1, 3)Spatial Pooling 2.67 (1, 5)Total Aggregation 3.33 (2, 5)Single Stock 3.67 (1, 5)Species Pooling 3.89 (2, 5)20 . Figures over English Rock kg < B t r a w l < kg kg < B t r a w l < kg kg < B t r a w l < kg kg < B t r a w l < kg B t r a w l > kg Figure 1: Mininum trawlable survey biomass B trawl estimates for DER complex species on the BC coast, aggregated to a 10km squaregrid. Estimates are produced by scaling average trawl survey ( kg/m ) density values in each grid cell by the cell’s area in m . Locationsthat do not show a coloured grid cell do not have any survey blocks from which to calculate relative biomass. Survey density for eachgrid cell is calculated from data for the Hecate Strait Assemblage Survey and the BC Groundfish Trawl Synoptic Survey, stored in theGFBio data base maintained at the Pacific Biological Station of Fisheries and Oceans, Canada. Thick black lines delineate the majorstatistical areas 3CD and 5ABCDE used for groundfish management BC, while the dashed grey lines makr out latitude and longitude,indicating the rotation of the coordinates to save space. The full colour figure is available in the online version of the article. .00.51.01.52.0 H S H G Q C S W C V I DoverEnglishRockComplex
Commercial Trawl Effort (1000 hrs) E qu ili b r i u m Y i e l d ( k t ) Figure 2: Operating model equilibrium yield curves for a given fishing effort for the nine DER complex management units, as well ascomplex yield curves within each stock area. Each panel shows the three individual DER complex species yield curves within a givenstock area in different colours, and the complex yield curve found by summing the three yield curves as a thick black line. ange(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Dover range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) English range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Rock H S H G range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Q C S range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) W C V I B t B Year
Figure 3: Spawning biomass depletion and relative catch simulation envelopes for all nine DER complex management units from the omniscient manager simulations. Medianbiomass is shown by the thick black line, with the grey region showing the central 95% of the distribution of spawning biomass, and thin black lines showing three randomlyselected simulation replicates. Catch is shown as grey bars in the historical period, which represent median catch in the projection, with thin vertical line segments showing thecentral 95% of the catch distribution. The depletion level associated with the traditional single species optimal biomass B MSY,SS is shown as a dashed horizontal green line,while the depletion level associated with our derived complex level multi-species optimal yield is shown as a blue horizontal dashed line. ange(yrs) c ( , m a x ( E _qp ft[, p , , ], na . r m = T )) H S H G range(yrs) c ( , m a x ( E _qp ft[, p , , ], na . r m = T )) Q C S range(yrs) c ( , m a x ( E _qp ft[, p , , ], na . r m = T )) W C V I C o mm e r c i a l T r a w l E ff o r t ( h r s ) Figure 4: Commercial fishing effort simulation envelopes for each stock area. Historical and median simulated effort in the projectionperiod are shown by a thicker black line, while the central 95% of the distribution of simulated effort is shown as grey shaded region inthe projection period. Single simulation replicates in the projection period are shown as thinner black lines. (0.5, nScen + 0.5) c ( , l o ss R ange ) Dover c(0.5, nScen + 0.5) c ( , l o ss R ange ) c ( , l o ss R ange ) Rich Mod Poor012345 c(0.5, nScen + 0.5) c ( , l o ss R ange ) English c(0.5, nScen + 0.5) c ( , l o ss R ange ) c ( , l o ss R ange ) Rich Mod Poor0.00.10.20.30.40.50.6 c(0.5, nScen + 0.5) c ( , l o ss R ange ) Rock H S H G c(0.5, nScen + 0.5) c ( , l o ss R ange ) Q C S c(0.5, nScen + 0.5) c ( , l o ss R ange ) Rich Mod Poor0.00.20.40.60.8 W C V I C u m u l a t i v e A b s o l u t e C a t c h Lo ss ( k t ) Data Scenario
SS HMS SpecPool SpatPool TA
Figure 5: Distributions of cumulative absolute loss in catch (kt) for the projection years 2026 to 2036 under each assessment model and OM data scenario (x-axis labels). Thepoints show the median loss, while line segments show the central 90% of cumulative loss distributions. Panels show the cumulative loss distributions by species (columns) andstocks (rows). Each assessment model point type and colour is shown in the legend at the bottom. Loss distributions for assessment models that couldn’t reach 100 replicateswith satisfactory convergence are shown as a dashed line. ange(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Dover c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) English c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Rock H S H G range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Q C S range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) W C V I B t B Year
Figure 6: Spawning biomass depletion and relative catch simulation envelopes for all nine DER complex management units when assessed by the Hierarchical Multi-stockassessment model under the Poor data-quality scenario. Median biomass is shown by the thick black line, with the grey region showing the central 95% of the distribution ofspawning biomass, and thin black lines showing three randomly selected simulation replicates. Catch is shown as grey bars in the historical period, which represent mediancatch in the projection, with thin vertical line segments showing the central 95% of the catch distribution. Coloured circles on the right hand vertical axis show the biomassdepletion level associated with the multi-species (closed blue circle) and single-species (open green circle) maximum sustainable yield. ange(yrs) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) Dover c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) English c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) Rock H S H G range(yrs) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) Q C S range(yrs) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) W C V I Year S p a w n i ng B i o m ass , T AC , a nd C a t c h ( k t) Figure 7: Operating model spawning stock biomass (red line), commercial trawl vulnerable biomass (grey dotted line), retrospective assesment model estimates of spawningstock biomass (thin grey/purple lines), and catch and TACs (grey bars) from the first simulation replicate in the Poor data-quality scenario and under the HierarchicalMulti-stock assessment model. Catch bars show realised catch in grey for the whole simulation period, and unfilled bars in the projection period show the difference betweenMP set TACs and realised catch. Coloured circles on the right hand vertical axis show the biomass level associated with the multi-species (closed blue circle) and single-species(open green circle) maximum sustainable yield. ange(B_qSAsp[2, , , s, p], omniB_qsp[, s, p], omniC_qsp[, s, p]) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) Dover range(B_qSAsp[2, , , s, p], omniB_qsp[, s, p], omniC_qsp[, s, p]) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) English range(B_qSAsp[2, , , s, p], omniB_qsp[, s, p], omniC_qsp[, s, p]) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) H S H G Rock range(B_qSAsp[2, , , s, p], omniB_qsp[, s, p], omniC_qsp[, s, p]) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) Q C S range(B_qSAsp[2, , , s, p], omniB_qsp[, s, p], omniC_qsp[, s, p]) r ange ( C _q SA s p [ , , , s , p ], o m n i B _q s p [, s , p ], o m n i C _q s p [, s , p ] ) W C V I B B
MSY C M SY SS HMS SpePool SpaPool TA Rich Poor
Figure 8: Tradeoff between catch and biomass during the 2028 - 2037 period implied by switching between different assessment models under Rich and Poor data-qualityscenarios. Panels are gridded by species (columns) and stocks (rows), with biomass relative to B MSY,MS,s,p on the horizontal axis, and catch relative to
MSY
MS,s,p on thevertical axis. Distributions of biomass and catch under the omniscient manager are shown by the black crosshair, with points indicating optimal biomass and yield for singlespecies maximum yield (open green circles) and multi-species maximum yield (closed blue circles). The biomass level at which a stock is critically overfished is shown as avertical red dashed line. Coloured point symbols show median biomass and catch for over all replicates for different assessment models, with assessment models under the samedata-quality scenario joined by a solid line (Rich) or dashed line (Poor). .0 0.2 0.4 0.6 0.8 1.0 1.2−1012 0.0 0.2 0.4 0.6 0.8 1.0 1.2 CV ( B MSY ) −1012 0.0 0.2 0.4 0.6 0.8 1.0 1.2−10123 0.0 0.2 0.4 0.6 0.8 1.0 1.2 SD ( U MSY ) −1012 0.0 0.1 0.2 0.3 0.4 0.5 0.6−0.6−0.4−0.20.00.20.40.6 R i c h t q and s U MSY −0.6−0.4−0.20.00.20.40.6 P oo r S t d . d i ff f r o m m ean c u m u l a t i v e l o ss Single−stock Hierarchical Species Pooling Spatial Pooling Total Aggregation
Figure 9: Regressions showing the average sensitivity of cumulative catch loss to the prior standard deviations under each assessment model (colours, line types) for thedata-rich (left column) and data-poor (right column) scenarios. The horizontal axis on each plot shows the prior standard deviation (CV for B MSY ), while the vertical axisshows the standardised difference between median cumulative loss for an assessment model and the mean of median cumulative loss values over AMs, stratified by species andarea. See the online version of the journal for a full colour version of the plot. ppendix A. The operating model The operating model was a standard age- and sex-structured operating model, with additional structure formulti-species and multi-stock population dynamics. DER complex species and stocks were simulated assuming noecological interactions or movement between areas. The lack of movement may be unrealistic, especially for DoverSole given their extent, but this is how the the DER complex stocks are currently managed in practice. The lackof ecological interactions is more realistic for Dover and English soles, as although both species are benthophagus,there is evidence that they belong to different feeding guilds (Pikitch, 1987).DER complex abundance N a,x,s,p,t for age a , sex x , species s and stock p at the start of year t was given by N a,x,s,p,t = . R s,p,t a = 1 ,N a − ,x,s,p,t − · e − Z a − ,x,s,p,t − < a < A,N a − ,x,s,p,t − · e − Z a − ,x,s,p,t − + N a,x,s,p,t − · e − Z a,x,s,p,t − a = A s , where R t is age-1 recruitment in year t , Z a,x,s,p,t is the instantaneous total mortality rate, and A s is the plus groupage for species s .Numbers-at-age were scaled to biomass-at-age by sex/species/area- specific weight-at-age. Weight-at-age was anallometric function of length-at-age w a,x,s,p = α x,s,p · L β x,s,p a,x,s,p where α x,s,p scaled between cm and kg, β x,s,p determined the rate of allometric growth, and L a,x,s,p was the lengthin cm of a fish of age a , sex x , species s and stock p . Length-at-age was given by the following Schnute formulationof the von-Bertalanffy growth curve (Schnute, 1981; Francis et al., 2016) L a = L A − ( L A − L A ) · (cid:18) e − kA − e − ka ) e − kA − e − kA (cid:19) where A and A are well spaced reference ages, L A and L A are the mean lengths in cm of fish at ages A and A ,and k is the growth coefficient. Note that in the growth model we dropped the sex, species and stock subscripts forconcision.The maturity-at-age ogive was modelled as a logistic function m a,s,p = e − ln 19( a − amat ,s,p ) amat ,s,p − amat ,x,s,p − , where m a,s,p was the proportion of age- a female fish of species s in stock p that were mature, and a mat ,s,p and a mat ,s,p are the ages at which 50% and 95% of fish of age- a , species s and stock p were mature.Female spawning stock biomass was calculated as B s,p,t = X a N a,x ,s,p,t m a,s,p w a,x ,s,p , where x denotes female fish only. Spawning stock biomass was used to calculate expected Beverton-Holt recruitment,which then had recruitment process errors applied R s,p,t +1 = R s,p, · h s,p · B s,p,t B s,p, · (1 − h s,p ) + (5 h s,p − · B s,p,t · e (cid:15) s,p,t +1 − . σ R,s,p , where R s,p, is unfished equilibrium recruitment, B s,p,t is the spawning stock biomass at time t , B s,p, is unfishedspawning stock biomasss, h s,p is stock-recruit steepness (average proportion of R s,p, produced when B s,p,t = . B s,p, ),and (cid:15) s,p,t is the recruitment process error with standard deviation σ R,s,p .The operating model was initialised in 1956 at unfished equilibrium for all species s and areas p , with numbers-at-age in 1956 given by N a,x,s,p, = . R s,p, a = 1 ,N a − ,x,s,p, · e − M x,s,p < a < A,N a − ,x,s,p, · e − Mx,s,p − e − Mx,s,p a = A, F a,x,s,p,f,t = S a,x,s,p,f · F s,p,f,t , where F s,p,f,t is the fully selected fishing mortality rate for fleet f at time t , and S a,x,s,p,f is the selectivity-at-age a for sex x in species s and area p by fleet f . Selectivity-at -age was modeled as a logistic function of length-at-age S a = (cid:18) (cid:18) − ln 19( L a − l sel ) l sel − l sel (cid:19)(cid:19) − , where L a is length-at-age, defined above, and l sel and l sel are the length-at-50% and length-at-95% selectivity,respectively; stock, species and fleet subscripts are left off for concision. Catch-at-age was then found via the Baranovcatch equation C a,x,s,p,f,t = (1 − e − Z a,x,s,p,f,t ) · N a,x,s,p,t w a,x,s,p F a,x,s,p,f,t Z a,x,s,p,f,t , where total mortality-at-age is defined as Z a,x,s,p,f,t = M x,s,p + S a,x,s,p,f · F a,x,s,p,f,t . Appendix A.1. Observation error standard deviations
Operating model observation error standard deviations were derived from estimates from fitting a hierarchicalage-structured model to DER complex data [Johnson and Cox, in prep]. To improve convergence in the simulations,we multipled all estimates by 0 . Table A.1: Log-normal observation error standard deviations for all DER complex biomass indices
Observation Error SDStock Historical Modern HS Ass. Syn
Dover sole
HSHG 0.174 0.158 0.328 0.138QCS 0.187 0.155 0.149WCVI 0.221 0.161 0.092
English sole
HSHG 0.170 0.163 0.255 0.191QCS 0.225 0.171 0.203WCVI 0.216 0.174 0.136
Rock sole
HSHG 0.164 0.161 0.254 0.187QCS 0.176 0.184 0.205WCVI 0.198 0.233 0.21432 ppendix B. Omniscient Manager Optimisation
We defined penalty functions so that inside their respective desired regions the penalty was zero, and otherwisethe penalty grew as a cubic function of distance from the desired region. For example, a penalty designed to keep ameasurement x above a the desired region boundary (cid:15) is of the form P ( x, (cid:15) ) = (cid:26) x ≥ (cid:15), | x − (cid:15) | x < epsilon. (B.1)This form has a several advantages over simple linear penalties, or a logarithmic barrier penalty (Srinivasan et al.,2008). First, the cubic softens the boundary threshold (cid:15) , effectively allowing a crossover if doing so favours anotherportion of the objective function. Second, unlike lower degree polynomials, cubic functions remain closer to the x -axis when | x − (cid:15) | <
1. Third, zero penalty within in the desirable region stops the objective function fromfavouring regions far from the boundaries of penalty functions. In contrast, a logarithmic function would favouroverly conservative effort series to keep biomass far from a lower depletion boundary. Finally, the cubic penaltyfunction and its first two derivatives are continuous at every point x , allowing for fast derivative-based optimisationmethods.We used a cubic spline of effort in each area to reduce the number of free parameters in the optimisation. Foreach area, 9 knot points were distributed across the full 40 year projection, making them spaced by 5 years. Wepadded the omniscient manager simulations by an extra eight years over the stochastic simulations to avoid anypossible end effects of the spline entering the performance metric calculations. Effort splines were constrained to bebetween 0 and 120 times the operating model E MSY,p , by replacing any value outside that range with the closestvalue inside the range (i.e. negative values by zero , large values by 120 E MSY p ).33 ppendix C. Hierarchical model performance under relaxed shrinkage prior SDs ange(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Dover c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) English c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Rock H S H G range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) Q C S range(yrs) c ( , m a x ( SB _q s p t[, s , p , ], na . r m = T )) W C V I B t B Year
Figure C.1: Spawning biomass depletion and relative catch simulation envelopes for all nine DER complex management units when assessed by the hierarchical multi-stockasssessment model under the Poor OM data quality scenario, when hierarchical shrinkage prior SDs were σ U MSY = τ q = 0 .
5. Median biomass is shown by the thick black line,with the grey region showing the central 95% of the distribution of spawning biomass, and thin black lines showing three randomly selected simulation replicates. Catch isshown as grey bars in the historical period, which represent median catch in the projection, with thin vertical line segments showing the central 95% of the catch distribution.Coloured circles on the right hand vertical axis show the biomass depletion level associated with the multi-species (closed blue circle) and single-species (open green circle)maximum sustainable yield. ange(yrs) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) Dover c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) English c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) Rock H S H G range(yrs) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) Q C S range(yrs) c ( , m a x ( VB _ s p t[ s , p , ], SB _ s p t[ s , p , ], na . r m = T )) W C V I Year S p a w n i ng B i o m ass , T AC , a nd C a t c h ( k t) Figure C.2: Operating model spawning stock biomass (red line), commercial trawl vulnerable biomass (grey dotted line), retrospective assesment model estimates of spawningstock biomass (thin grey/purple lines), and catch and TACs (grey bars) from the first simulation replicate in the Poor data quality scenario under the hierarchical multi-stockassessment model when hierarchical shrinkage prior SDs were σ U MSY = τ q = 0 .
5. Catch bars show realised catch in grey for the whole simulation period, and unfilled bars inthe projection period show the difference between MP set TACs and realised catch. Coloured circles on the right hand vertical axis show the biomass level associated with themulti-species (closed blue circle) and single-species (open green circle) maximum sustainable yield. eferences Richard S Appeldoorn. Model and method in reef fishery assessment. In
Reef fisheries , pages 219–248. Springer,1996.Alan R Baudron and Paul G Fernandes. Adverse consequences of stock recovery: European hake, a new “choke”species under a discard ban?
Fish and Fisheries , 16(4):563–575, 2015. doi: 10.1111/faf.12079. URL https://doi.org/10.1111/faf.12079.Ashleen J Benson, Sean P Cox, and Jaclyn S Cleary. Evaluating the conservation risks of aggregate harvestmanagement in a spatially-structured herring fishery.
Fisheries Research , 167:101–113, 2015.Trevor A Branch and Ray Hilborn. Matching catches to quotas in a multispecies trawl fishery: targeting andavoidance behavior under individual transferable quotas.
Canadian Journal of Fisheries and Aquatic Sciences , 65(7):1435–1446, 2008. doi: 10.1139/F08-065. URL https://doi.org/10.1139/F08-065.Jeremy S Collie and Henrik Gislason. Biological reference points for fish stocks in a multispecies context.
CanadianJournal of Fisheries and Aquatic Sciences , 58(11):2167–2176, 2001. doi: 10.1139/f01-158. URL https://doi.org/10.1139/f01-158.Jeremy S Collie and Carl J Walters. Adaptive management of spatially replicated groundfish populations.
CanadianJournal of Fisheries and Aquatic Sciences , 48(7):1273–1284, 1991. doi: 10.1139/f91-153. URL https://doi.org/10.1139/f91-153.Sean P Cox and Allen Robert Kronlund. Practical stakeholder-driven harvest policies for groundfish fisheries inBritish Columbia, Canada.
Fisheries Research , 94(3):224–237, 2008. doi: 10.1016/j.fishres.2008.05.006. URLhttps://doi.org/10.1016/j.fishres.2008.05.006.Sean P Cox, Allen R Kronlund, and Ashleen J Benson. The roles of biological reference points and operationalcontrol points in management procedures for the Sablefish (
Anoplopoma fimbria ) fishery in British Columbia,Canada.
Environmental Conservation , 40(4):318–328, 2013. doi: 10.1017/S0376892913000271. URL https://doi.org/10.1017/S0376892913000271.William K de la Mare. Tidier fisheries management requires a new MOP (management-oriented paradigm).
Reviewsin Fish Biology and Fisheries , 8(3):349–356, 1998. doi: 10.1023/A:1008819416162. URL https://doi.org/10.1023/A:1008819416162.Edward E DeMartini. Hazards of managing disparate species as a pooled complex: A general problem illustrated bytwo contrasting examples from hawaii.
Fish and Fisheries , 20(6):1246–1259, 2019.DFO. A Harvest Strategy Compliant with the Precautionary Approach. Technical Report 2006/023, DFO Can. Sci.Advis. Rep., 2006.Fisheries and Oceans, Canada. Pacific Region Integrated Fisheries Management Plan: Groundfish, 2015.RIC Chris Francis, Alexandre M Aires-da Silva, Mark N Maunder, Kurt M Schaefer, and Daniel W Fuller. Estimatingfish growth for stock assessments using both age–length and tagging-increment data.
Fisheries research , 180:113–118, 2016.Sarah Gaichas, Robert Gamble, Michael Fogarty, Hugues Benoît, Tim Essington, Caihong Fu, Mariano Koen-Alonso, and Jason Link. Assembly rules for aggregate-species production models: simulations in support ofmanagement strategy evaluation.
Marine Ecology Progress Series , 459:275–292, 2012. doi: 10.3354/meps09650.URL https://doi.org/10.3354/meps09650.Henrik Gislason. Single and multispecies reference points for baltic fish stocks.
ICES Journal of Marine Science , 56(5):571–583, 1999. doi: 10.1006/jmsc.1999.0492. URL https://doi.org/10.1006/jmsc.1999.0492.Ray Hilborn and Carl J Walters. A general model for simulation of stock and fleet dynamics in spatially heterogeneousfisheries.
Canadian Journal of Fisheries and Aquatic Sciences , 44(7):1366–1369, 1987. doi: 10.1139/f87-163. URLhttps://doi.org/10.1139/f87-163. 37ay Hilborn and Carl J Walters.
Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty/Bookand Disk . Springer Science & Business Media, 1992.Yan Jiao, Christopher Hayes, and Enric Cortés. Hierarchical Bayesian approach for population dynamics modellingof fish complexes without species-specific data.
ICES Journal of Marine Science: Journal du Conseil , 66(2):367–377, 2009. doi: 10.1093/icesjms/fsn162. URL https://doi.org/10.1093/icesjms/fsn162.Yan Jiao, Enric Cortés, Kate Andrews, and Feng Guo. Poor-data and data-poor species stock assessment using abayesian hierarchical approach.
Ecological Applications , 21(7):2691–2708, 2011. doi: 10.1890/10-0526.1. URLhttps://doi.org/10.1890/10-0526.1.Samuel D N Johnson and Sean P Cox. Evaluating the role of data quality when sharing information in hierarchicalmulti-stock assessment models, with an application to dover sole.
Can. J. Fish. Aquat. Sci , 2018. doi: 10.1139/cjfas-2018-0048. URL https://doi.org/10.1139/cjfas-2018-0048.Alexander Kempf, John Mumford, Polina Levontin, Adrian Leach, Ayoe Hoff, Katell G Hamon, Heleen Bartelings,Morten Vinther, Moritz Staebler, Jan Jaap Poos, et al. The MSY concept in a multi-objective fisheries environment–Lessons from the North Sea.
Marine Policy , 69:146–158, 2016. doi: 10.1016/j.marpol.2016.04.012. URLhttps://doi.org/10.1016/j.marpol.2016.04.012.Steven JD Martell, Carl J Walters, and Ray Hilborn. Retrospective analysis of harvest management performance forbristol bay and fraser river sockeye salmon (oncorhynchus nerka).
Canadian Journal of Fisheries and AquaticSciences , 65(3):409–424, 2008. doi: 10.1139/f07-170. URL https://doi.org/10.1139/f07-170.Mark N Maunder and Kevin R Piner. Contemporary fisheries stock assessment: many issues still remain.
ICESJournal of Marine Science , 72(1):7–18, 2015.Lars O Mortensen, Clara Ulrich, Jan Hansen, and Rasmus Hald. Identifying choke species challenges for an individualdemersal trawler in the north sea, lessons from conversations and data analysis.
Marine Policy , 87:1–11, 2018. doi:10.1016/j.marpol.2017.09.031. URL https://doi.org/10.1016/j.marpol.2017.09.031.Kotaro Ono, Alan C Haynie, Anne B Hollowed, James N Ianelli, Carey R McGilliard, and André E Punt. Managementstrategy analysis for multispecies fisheries, including technical interactions and human behavior in modellingmanagement decisions and fishing.
Canadian Journal of Fisheries and Aquatic Sciences , (999):1–18, 2017. doi:10.1139/cjfas-2017-0135. URL https://doi.org/10.1139/cjfas-2017-0135.Jerome J Pella and Patrick K Tomlinson. A generalized stock production model.
Inter-American Tropical TunaCommission Bulletin , 13(3):416–497, 1969. URL http://aquaticcommons.org/id/eprint/3536.Ellen K Pikitch. Use of a mixed-species yield-per-recruit model to explore the consequences of various managementpolicies for the oregon flatfish fishery.
Canadian Journal of Fisheries and Aquatic Sciences , 44(S2):s349–s359,1987. doi: 10.1139/f87-336. URL https://doi.org/10.1139/f87-336.AE Punt and DS Butterworth. The effects of future consumption by the cape fur seal on catches and catch rates ofthe cape hakes. 4. modelling the biological interaction between cape fur seals arctocephalus pusillus pusillus andthe cape hakes merluccius capensis and m. paradoxus.
South African Journal of Marine Science , 16(1):255–285,1995.André E Punt, Anthony DM Smith, and Gurong Cui. Evaluation of management tools for australia’s south eastfishery. 1. modelling the south east fishery taking account of technical interactions.
Marine and FreshwaterResearch , 53(3):615–629, 2002. doi: 10.1071/MF01007. URL https://doi.org/10.1071/MF01007.André E Punt, Roy Deng, Sean Pascoe, Catherine M Dichmont, Shijie Zhou, Éva E Plagányi, Trevor Hutton,William N Venables, Rob Kenyon, and Tonya Van Der Velde. Calculating optimal effort and catch trajectories formultiple species modelled using a mix of size-structured, delay-difference and biomass dynamics models.
FisheriesResearch , 109(1):201–211, 2011a. doi: 10.1016/j.fishres.2011.02.006. URL https://doi.org/10.1016/j.fishres.2011.02.006. 38ndrè E Punt, David C Smith, and Anthony DM Smith. Among-stock comparisons for improving stock assessmentsof data-poor stocks: the “Robin Hood” approach.
ICES Journal of Marine Science: Journal du Conseil , 68(5):972–981, 2011b. doi: 10.1093/icesjms/fsr039. URL https://doi.org/10.1093/icesjms/fsr039.André E Punt, Doug S Butterworth, Carryn L Moor, José A A De Oliveira, and Malcolm Haddon. Managementstrategy evaluation: best practices.
Fish and Fisheries , 2016. doi: 10.1111/faf.12104. URL https://doi.org/10.1111/faf.12104.André E Punt, Daniel K Okamoto, Alec D MacCall, Andrew O Shelton, Derek R Armitage, Jaclyn S Cleary, Ian PDavies, Sherri C Dressel, Tessa B Francis, Phillip S Levin, et al. When are estimates of spawning stock biomassfor small pelagic fishes improved by taking spatial structure into account?
Fisheries Research , 206:65–78, 2018.André E Punt, Michael G Dalton, and Robert J Foy. Multispecies yield and profit when exploitation rates varyspatially including the impact on mortality of ocean acidification on north pacific crab stocks.
Fisheries Research ,225:105481, 2020. doi: 10.1016/j.fishres.2019.105481. URL https://doi.org/10.1016/j.fishres.2019.105481.Terrance J Quinn. Ruminations on the development and future of population dynamics models in fisheries.
NaturalResource Modeling , 16(4):341–392, 2003.Milner B Schaefer. Some considerations of population dynamics and economics in relation to the management of thecommercial marine fisheries.
Journal of the Fisheries Board of Canada , 14(5):669–681, 1957. doi: 10.1139/f57-025.URL https://doi.org/10.1139/f57-025.Jon Schnute. A versatile growth model with statistically stable parameters.
Canadian Journal of Fisheries andAquatic Sciences , 38(9):1128–1140, 1981.GA Seber.
Estimation of Animal Abundance . Oxford University Press, 2 edition, 1997.Balasubramaniam Srinivasan, Lorenz T Biegler, and Dominique Bonvin. Tracking the necessary conditions ofoptimality with changing set of active constraints using a barrier-penalty function.
Computers & ChemicalEngineering , 32(3):572–579, 2008.James T Thorson, Jason M Cope, Kristin M Kleisner, Jameal F Samhouri, Andrew O Shelton, and Eric J Ward.Giants’ shoulders 15 years later: lessons, challenges and guidelines in fisheries meta-analysis.
Fish and Fisheries ,16(2):342–361, 2015. doi: 10.1111/faf.12061. URL https://doi.org/10.1111/faf.12061.Willard Waldo Wakefield. Feeding relationships within assemblages of nearshore and mid-continental shelf benthicfishes off oregon. 1984.Carl Walters. Adaptive management of renewable resources. 1986.Carl Walters. Evaluation of quota management policies for developing fisheries.
Canadian Journal of Fisheries andAquatic Sciences , 55(12):2691–2705, 1998. doi: 10.1139/f98-172. URL https://doi.org/10.1139/f98-172.Carl J Walters and Ramón Bonfil. Multispecies spatial assessment models for the British Columbia groundfish trawlfishery.
Canadian Journal of Fisheries and Aquatic Sciences , 56(4):601–628, 1999. doi: 10.1139/f98-205. URLhttps://doi.org/10.1139/f98-205.Carl J Walters, Jeremy S Collie, and Timothy Webb. Experimental designs for estimating transient responses tomanagement disturbances.
Canadian Journal of Fisheries and Aquatic Sciences , 45(3):530–538, 1988.Henning Winker, Felipe Carvalho, and Maia Kapur. JABBA: just another Bayesian biomass assessment.