The benefits and risks of incorporating climate-driven growth variation into stock assessment models, with application to Splitnose Rockfish (Sebastes diploproa)

The benefits and risks of incorporating climate-driven growth variation into stock assessment... Abstract Indices of annual growth variation are not routinely incorporated into fisheries stock assessment models, due to a lack of a general framework for deciding when to include these indices, and of a mechanistic understanding about growth drivers. Such incorporation may also not necessarily lead to improved estimation or management performance. We demonstrate a way to incorporate such an index into an assessment model (Stock Synthesis), and use risk analysis to evaluate its management-related advantages and shortcomings. We applied this method to splitnose rockfish (Sebastes diploproa), where a previously developed growth index is highly correlated with decadal-scale climate indices. We find that including a similar index in the simulated assessment increases precision and reduces bias of parameter estimates. However, not including an index or including a completely erroneous index led to highly imprecise estimates when growth was strongly climate-driven. Including this growth index when individual growth was actually constant did not lead to poorer estimation performance. The risk analysis approach can be applied to other stocks to evaluate the consequences of including an index of growth variation. Introduction Climate and environmental factors can have a large impact on population and community dynamics (Harley et al., 2006; Pinsky et al., 2013; Rowe and Terry, 2014). One of the ways in which the environment can play a major role in determining the dynamics of fish populations is through the metabolism, and ultimately growth, of individual fish (Ait Youcef et al., 2015; Sünksen et al., 2010). However, little guidance is available for whether or when to incorporate environmentally driven growth variables into the stock assessments used as the basis for advice to managers to support fisheries decision making. Such inclusion could more accurately describe the dynamics of fish populations, and potentially disentangle whether fluctuations in biomass are due to time-varying environmental conditions or fishing. Much of our current understanding of how to best incorporate environmental information into stock assessments derives from recruitment-environment relationships (e.g. Mueter et al., 2011; Schirripa et al., 2009; Thorson et al., 2013), as opposed to growth variation. Growth is an integral component of fisheries stock assessment, and the relationship between age and size is often estimated in assessments (Taylor and Methot, 2013), based on the assumption that the parameters of this relationship are constant over space and time (Lorenzen, 2016). However, it has been shown (e.g. Stawitz et al., 2015; Weatherley, 1990; Webber and Thorson, 2016) that growth in fish is inherently plastic, and that this “constant growth” paradigm needs to be challenged (Lorenzen, 2016; Thorson and Minte-Vera, 2016). Defining climate-growth relationships would allow fishery managers to better discriminate the effects of climate from those of fishing, and make more accurate predictions of future stock size, should there be a means of projecting the relationship into the future (Black et al., 2008). Black (2009) examined ring widths of rockfish (Sebastes spp) otoliths and reduced these widths to species-specific dimensionless indices that are correlated with similar dendrochronology indices from geoducks (Panopea abrupta) and trees in the same general area. These indices were also found to coincide with climate indices in the California Current System. However, these growth indices have never been incorporated in assessments due, inter alia, to the lack of a general framework for deciding whether and under what circumstances to do so. Many studies have used simulation to evaluate the benefits of incorporating indices of climate on recruitment in assessments, notably Schirripa et al. (2009) for sablefish, Hulson et al. (2013) and A’mar et al. (2009) for walleye pollock, and Haltuch and Punt (2011) for groundfish generally. Simulations are used to evaluate model performance and test the effects of different hypotheses on management outcomes (e.g. Johnson et al., 2015). However, this has yet to be done for time-varying growth. Maunder and Watters (2003) proposed a general framework for integrating environmental time series in assessments, where simulations are used to test and compare a model that includes environmental effects with the traditional approach where environment variables are ignored. Maunder and Watters (2003) focused on the ability to correctly identify the relationship between a selected variable in the model and the environmental variable, and to accurately estimate specific parameters, and not on the impact of these relationships on management outcomes. The goal of this study is to propose and illustrate a risk-analysis approach to estimating the benefits and risks of incorporating a climate index on somatic fish growth into stock assessments. Our analyses are conditioned closely upon a climate index developed by Black (2009) and the stock assessment for splitnose rockfish (Sebastes diploproa) in the Northeast Pacific Ocean (Gertseva et al., 2009), and evaluate management outcomes associated with various climate-growth scenarios. We expand on Maunder and Watters (2003) and develop a table using simulation modelling to provide guidance on whether incorporation of climate information is justified in a given case. Material and methods Overview The risk-analysis was developed using a widely used stock assessment platform—Stock Synthesis (SS; Methot and Wetzel, 2013). SS can incorporate all available data into a single analysis, and estimate parameters by maximizing the product of the likelihood function for each data type (Fournier and Archibald, 1982; Maunder and Punt, 2013). The assessment was simplified to ensure that results are widely applicable to other contexts. The numbers of fishing fleets and surveys were reduced to one of each data type. Operating models (OMs) were used to describe “true” population dynamics and quantity of data for five states of nature—one with constant growth, and four with time-varying environmentally driven growth indices that vary in terms of auto-correlation and relationship to the OM growth parameters. Each simulation was projected for 50 years using the OM (with the relevant assumptions made for each OM, i.e. constant future growth for one and time-varying future growth for the other four), with fishing occurring under a constant-catch limit obtained from the OM. Implementation error was not simulated, so catch was assumed to be taken exactly. The OM outputs based in the historical and forecast periods were considered “true” for the simulations. Observation models simulated the data-gathering process. A total of 1100 data sets were generated, with 100 data sets for each combination of OM and observation model. Each data set was analysed using two estimation methods (EMs). We examined the performance of two EMs—one that estimated time-invariant growth, and one that assumed growth was correlated with an available environmental index. Management reference points were estimated using these EMs, and used to set a catch limit that was projected forward as a constant in the simulated fishery system (as represented in the OM). The resultant forecasted outcomes were then compared to the “true” (i.e. OM) quantities to obtain performance metrics, such as lost yield and final stock biomass. A summary of the overall process can be found in Figure 1. All code for replicating this analysis is publicly available (https://github.com/lee-qi/ss3-growth-index). Figure 1. View largeDownload slide Flowchart depicting the general steps of the simulation process. Multiple operating and observation models were considered, but this is not shown here. Figure 1. View largeDownload slide Flowchart depicting the general steps of the simulation process. Multiple operating and observation models were considered, but this is not shown here. Incorporating climate indices of individual growth into assessment models As the growth model most commonly used in stock assessments (Essington et al., 2001; Francis, 2016), this study was based, for simplicity, on the von Bertalanffy growth function (Von Bertalanffy, 1957):   Lt=L∞(1-e-Kt-t0e-Kt-t0) (1) where Lt is the length of the fish at time t, L∞ is the asymptotic length of the fish and a ratio of anabolic to catabolic factors, K is the intrinsic growth rate and a catabolic factor, and t0 is a constant that adjusts the model along a time axis. Temperature is known to affect catabolism, and thus both K and L∞ (Brunel and Dickey-Collas, 2010; Lorenzen, 2016), with Kimura (2008) showing a negative correlation between the two parameters. Temporal variation in growth due to changes in environmental drivers can be incorporated in an assessment by modifying the values of the parameters of the growth function. In SS, an environmental index can be linked to the values of the growth parameters in two ways—additively or multiplicatively (Methot and Wetzel, 2013). In this study, we included an annual index ɛt that is associated with variation in growth, using a multiplicative relationship (Equation 2) for both maximum size for sex s:   L∞,s,t=L∞,s,base·eβL,sɛt (2a) and growth rate:   Ks,t=Ks,base·eβK,sɛt (2b) where β is the parameter linking the growth parameter to the environmental index, L∞,s,base and Ks,base are median values for the asymptotic length and growth rate parameters for sex s, and length-at-age in each year is calculated from length-at-age the previous year (using Equation 1). Should the index value be sufficiently large to result in a negative growth increment (e.g. if L∞,t <  Lt from Equation 1), the individuals remain at the same size for the following year and are assumed to not experience any growth, which thus avoids the occurrence of “shrinking” fish. Biological characteristics Populations were simulated based on parameters obtained from a simplified version of the assessment (Supplementary Table S1; Gertseva et al., 2009). The data sets with the longest time series for a single fishing fleet and a single fishery-independent survey were selected (Supplementary Figure S1c). In the assessment approved for management (Gertseva et al., 2009), the parameters governing natural mortality (M) and steepness (h) were pre-specified, while log initial recruitment (lnR0) and individual-growth parameters including length at age 0 (L0), L∞, K, and variability in size-at-age for fish young fish (CV1) and for old fish up to L∞ ( CV2) were estimated. Growth of splitnose rockfish is sex-specific, with males estimated to grow at a faster rate (i.e. higher Kmale), while females are estimated to attain larger sizes (i.e. higher L∞,female). This is also one of the two rockfish species for which a climate-growth index was developed by Black (2009). Simulation method Index data generation The environmental index ɛt was simulated using an autoregressive (AR) model with bias correction (Equation 3), adapted from Thorson et al. (2014), with parameters calculated for splitnose rockfish. The environmental index was unique to each replicate within each state of nature, and scaled to a mean of 0.   ɛt=ρsplitnoseɛt-1+ 1-ρsplitnose2 ωtωtρsplitnoseɛt-1+ ωtωtfor t>1for t=1  ωt ∼ Normal-σsplitnose22×1-ρsplitnose1-ρsplitnose2, σsplitnose2 (3) where ωt is independent and identically normally distributed variation for year t,  ρsplitnose is the species-specific autoregressive parameter, and σsplitnose2 is the variance of the index. As the index is treated as data, the species-specific parameters in Equation 3 ( ρsplitnose and σsplitnose2) were defined external to the OM, and fixed. Comparative runs (across OMs) used the same sequence of random numbers. Five trends in growth were investigated, based on the extent of auto-correlation (i.e. ρsplitnose) in the environmental indices and the strength of the relationship between the environmental index and the growth parameters (i.e. βL,s and βK,s; Equation 2; Table 1). Sensitivity analyses were conducted for a series of β and ρsplitnose values in the OM, using OM 2 (weakly climate-driven growth) as the base case. As the values for β increase, so too does the impact of the environmental index (and vice-versa). The series of β values tested ranged from 0 to 0.5, in 0.1 increments. The higher the value of ρsplitnose, the more that long periods of time have a series of fast or slow growth, should the index be highly correlated with growth. We tested scenarios with indices of varying levels of auto-correlation, ranging the ρsplitnose values from 0.0 (i.e. white noise) to 0.8 by 0.2 increments, and a final maximum high value of 0.95. Table 1. Summary of OMs, and their associated values. Number  OM  Strength of relationship  βL,m  βL,m  βK,f  βK,m  σp2  Degree of auto-correlation ( ρp)  1  Not Climate-Driven  No relationship  0  0  0  0  0.88  Weak (0.33)  2  Weakly Climate-Driven  Weak  −0.12  −0.083  0.13  0.22  0.88  Weak (0.33)  3  Climate-Driven (High AR)  Weak  −0.12  −0.083  0.13  0.22  0.88  Strong (0.95)  4  Climate-Driven (High Beta)  Strong  −0.3o  −0.3o  0.3o  0.3o  0.88  Weak (0.33)  5  Climate-Driven (High AR and High Beta)  Strong  −0.30  −0.3o  0.3o  0.3o  0.88  Strong (0.95)  Number  OM  Strength of relationship  βL,m  βL,m  βK,f  βK,m  σp2  Degree of auto-correlation ( ρp)  1  Not Climate-Driven  No relationship  0  0  0  0  0.88  Weak (0.33)  2  Weakly Climate-Driven  Weak  −0.12  −0.083  0.13  0.22  0.88  Weak (0.33)  3  Climate-Driven (High AR)  Weak  −0.12  −0.083  0.13  0.22  0.88  Strong (0.95)  4  Climate-Driven (High Beta)  Strong  −0.3o  −0.3o  0.3o  0.3o  0.88  Weak (0.33)  5  Climate-Driven (High AR and High Beta)  Strong  −0.30  −0.3o  0.3o  0.3o  0.88  Strong (0.95)  The values for β for OMs 2 and 3 were estimated by including the otolith index data as an index on the growth parameters, and applying the assessment to these data. The values for β for OMs 4 and 5 were fixed at a high value of 0.3, as determined by the sensitivity analysis, with positive values for K and negative values for L∞, consistent with the notion of having opposite effects. ρsplitnose and σsplitnose2 in OMs 2 and 4 were determined by fitting an AR-1 model to the otolith index obtained from Black (pers. comm.). As such, the environmental index in OM 2 had the most realistic characteristics, while OMs 3–5 were exploratory. Observation errors Observation error was added to the expected values for each data type in the OMs to generate data sets for use in the EMs (Table 2). Additional variability for each simulation was included in the form of independent and lognormally distributed yearly deviations about a Beverton-Holt stock-recruitment relationship (Methot and Wetzel, 2013). A bias-correction factor (Methot and Taylor, 2011) was applied in the EM to avoid bias in estimation of expected recruitment when there is less data. Table 2. Data types and error structures used in data generation (Gertseva et al., 2009). Data type  Years  Error distribution  Measure of uncertainty  Catch—Fishery  1900–2008  Lognormal  SE of log(value) = 0.01  Abundance indices—Survey  2003–2008  Lognormal  SE of log(value) = 0.25–0.34  Length compositions—Fishery  1978–2008  Multinomial  Effective N = 505.4  Length compositions—Fishery (Discard)  1987; 2006–2007  Multinomial  Effective N = 505.4  Length compositions—Survey  2003–2008  Multinomial  Effective N = 1339  Conditional age-at-length compositions—Survey  2003–2008  Multinomial  Effective N = 124.3  Mean body weight—Fishery  2002–2008  Lognormal  CV = 0.3  Discards—Fishery  1987; 2002–2007  Lognormal  CV = 0.5 (1987); 0.2 (2002–2007)  Environmental index—Black (2009)  1936–2006  Known exactly    Data type  Years  Error distribution  Measure of uncertainty  Catch—Fishery  1900–2008  Lognormal  SE of log(value) = 0.01  Abundance indices—Survey  2003–2008  Lognormal  SE of log(value) = 0.25–0.34  Length compositions—Fishery  1978–2008  Multinomial  Effective N = 505.4  Length compositions—Fishery (Discard)  1987; 2006–2007  Multinomial  Effective N = 505.4  Length compositions—Survey  2003–2008  Multinomial  Effective N = 1339  Conditional age-at-length compositions—Survey  2003–2008  Multinomial  Effective N = 124.3  Mean body weight—Fishery  2002–2008  Lognormal  CV = 0.3  Discards—Fishery  1987; 2002–2007  Lognormal  CV = 0.5 (1987); 0.2 (2002–2007)  Environmental index—Black (2009)  1936–2006  Known exactly    The effects of including an incorrectly observed environmental index were examined using four observation models: 158 years of environmental indices (108 historical, 50 projected) observed without error; i.e. the index used to generate the 20 years of historical environmental indices observed without error, with 50 years of projected indices, i.e. a shortened time series of data when there is insufficient historical information 158 years of auto-correlated indices generated using the same model and values for ρp and σp2 as the “true” index, but with a different starting seed, i.e. a time series of auto-correlated indices that seems plausible but is completely inaccurate 158 years of an independent time series of random errors (white noise), generated under N(0,1) These scenarios were similar to those described by Stawitz et al. (2015), where growth anomaly patterns were estimated to be in one of three predominant forms—trendless, sustained trend, and near-zero. Sensitivity of the EM to the accuracy of the environmental index was also examined, using OM 5 as a base. This was done by creating a second auto-correlated index ( ɛt,wrong) using the same AR parameters and process (Equation 3), and adding it to the true index ( ɛt,true) in varying proportions a, which ranged from 1 to 0 at 0.2 increments, to create a new index ( ɛt,new):   ɛt,new=a·ɛt,true+1-a·ɛt,wrong (4) This procedure therefore leads to a new environmental index for which a·100% of the variance is correct, and (1-a)·100% of the variance is auto-correlated noise (i.e. the correlation between ɛt,true and ɛt,new is a). The new index ɛt,wrong is calculated to have the same variance and first-order autocorrelation as the true index. Estimation methods The estimation scenarios mimicked an actual assessment process, by estimating model parameters and management reference points. Supplementary Table S1 lists the key parameters of the population dynamics model and how they are estimated in the EM. Priors for the growth parameters were set to the same values as the priors in the OMs. Two EMs were used—one where the index was included ( β is estimated), and one where the index was not ( β=0). Each combination of OM, observation model, and EM is referred to as a scenario, and 16 scenarios were evaluated (Table 3). Table 3. Summary of scenarios tested. OM  Observation of Env index  EM  Scenario label  Scenario name  1  Not Climate-Driven  NA  Constant  1-0-A  Constant growth with correct EM  Wrong  Time-varying  1-0-B  Constant growth with over-parameterized EM  2  Weakly Climate-Driven  NA  Constant  2-0-A  Weakly climate-driven growth with incorrect EM  Exact  Time-varying  2-1-B  Weakly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  2-2-B  Weakly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  3  Climate-Driven (High AR)  NA  Constant  3-0-A  Highly autocorrelated growth with incorrect EM  Exact  Time-varying  3-1-B  Highly autocorrelated growth with correct EM and index  4  Climate-Driven (High Beta)  NA  Constant  4-0-A  Highly varied growth with incorrect EM  Exact  Time-varying  4-1-B  Highly varied growth with correct EM and index  5  Climate-Driven (High AR and High Beta)  NA  Constant  5-0-A  Strongly climate-driven growth with incorrect EM  Exact  Time-varying  5-1-B  Strongly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  5-2-B  Strongly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  OM  Observation of Env index  EM  Scenario label  Scenario name  1  Not Climate-Driven  NA  Constant  1-0-A  Constant growth with correct EM  Wrong  Time-varying  1-0-B  Constant growth with over-parameterized EM  2  Weakly Climate-Driven  NA  Constant  2-0-A  Weakly climate-driven growth with incorrect EM  Exact  Time-varying  2-1-B  Weakly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  2-2-B  Weakly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  3  Climate-Driven (High AR)  NA  Constant  3-0-A  Highly autocorrelated growth with incorrect EM  Exact  Time-varying  3-1-B  Highly autocorrelated growth with correct EM and index  4  Climate-Driven (High Beta)  NA  Constant  4-0-A  Highly varied growth with incorrect EM  Exact  Time-varying  4-1-B  Highly varied growth with correct EM and index  5  Climate-Driven (High AR and High Beta)  NA  Constant  5-0-A  Strongly climate-driven growth with incorrect EM  Exact  Time-varying  5-1-B  Strongly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  5-2-B  Strongly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  Performance evaluation Performance was evaluated for each scenario by (a) comparing spawning stock biomass (SSB), recruitment, and parameters from the OM with corresponding estimates from the EM, and (b) conducting forecasts of the OM in which future catches are based on management reference points estimated by the EM and summarizing the results of the projections in terms of “lost yield” and differences between ideal and anticipated stock status when management is based on the EM. Each of these components is discussed below. The estimation performance of the EMs was summarized using relative errors by year for SSB and recruitment, using RE =  (θ^-θ)/θ, and median absolute relative errors (MARE = median( |(θ^-θ)/θ|)) for other parameter estimates (e.g. growth, R0) and derived quantities (e.g. maximum sustainable yield MSY, virgin stock biomass B0, final year biomass B2008). These median values were used as a measure of bias, as occasional outliers greatly skewed the mean values. MAREs have been used previously as a measure of precision for point parameter estimates (e.g. Johnson et al., 2015). The management implications of growth variability and how it is handled in the assessment were evaluated by assuming that management aims to follow a constant catch strategy, based on MSY, and implemented without error. This harvest strategy was selected for illustrative purposes. MSY is calculated from the base growth, the stock-recruit, and the selectivity parameters (Methot and Wetzel, 2013), as found in the first year of the simulation, and unaffected by any time-varying biology. The simulated fishery system is projected forward using this constant-catch strategy, and the consequences evaluated using selected metrics. Projecting with the MSY obtained from the OM represents the “true” state of the resource for that replicate, with all relevant hypotheses associated with that OM. Recruitment in the OM during the forecast period is simulated given recruitment variation, while management performance during the forecast period is calculated by simulating future population sizes and catches given either that future growth is constant (for OM 1) or time-varying (for OMs 2–5). The EM was also used to estimate MSY, which was then used as the basis for projections for each year of a 50-year forecast period. These “estimated” catches were provided to the OM and a second set of projections undertaken in which the catches taken from the population in the OM were set to the “estimated” catches. Thus, the OM was projected forward twice, once using the MSY from the OM (“true yield”) and once using the estimates of MSY from EM (“projected yield”). The risk for each EM given the true state of nature was calculated in terms of relative lost yield, the forecasted final year stock status, and the proportion of replicates overfished (see below). With the exception of the final metric (proportion overfished), all the metrics were calculated for each simulation replicate, and the median value and interquartile range across replicates (representing central tendency and precision, respectively) reported. Relative Lost Yield: Difference between estimated retained catch (yield from the EM) and true retained catch (yield from the OM) for the first projected year as a proportion of true catch [(projected yield—true yield)/true yield] Final Stock Status: Ratio of projected (using MSY from EM) final spawning stock biomass to true (using MSY from OM) final stock biomass (projected B2058/true B2058) Proportion Overfished: Proportion of replicates where the final biomass is less than 0.25 of B0, which is the metric by which the stock would be declared overfished {Gertseva et al., 2009; number of [(projected B2058/true B0) < 0.25]/number of runs}. The final year is used in this metric as it would be near-impossible for a fishery to recover from an overfished state under a constant-catch harvest strategy without any intervention. Results The simplified base model (Supplementary Figure S1c) had catch start in 1902, peaking in 1998 (Supplementary Figure S1a). Overall SSB increased when the environmental index was added, and the peaks and troughs in the time series appeared more pronounced (Supplementary Figure S1b). Several replicates (10/1100) failed to converge (non-positive definite Hessian matrix), and were discarded. If a replicate failed within a single scenario, the results for that replicate for all scenarios were discarded. In a true assessment, model configurations could be manually adjusted to ensure convergence, but such methods would not be feasible in a simulation study so we instead discarded these few replicates. Estimation performance Time-invariant growth in OM The two EMs—one correctly specified (1-0-A), one over-parameterized (1-0-B)—performed well when growth was constant in the OM (Figure 2). There was a slight overall positive bias (measured in median RE) for SSB and recruitment for both EMs, although the median RE for each time series was less than 7%. Over-parameterization did not result in a substantial difference in the bias of SSB and recruitment (Figure 2(ii)) because the assessment model generally estimated link parameters (β) close to zero. Base growth parameters were also estimated without bias and precisely, and were similar for the two EMs (Table 4). Table 4. MAREs for the model parameters and derived quantities.     The colours go from white to grey by order of magnitude, with 0 being white, and larger absolute values being grey. The growth parameters in these table are the base parameters of the model, i.e. before the index effects took place. Asterisks (*) denote median absolute (not relative) error as the true values for those parameters are zero. Figure 2. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the OM with constant growth. Each column represents an EM, where (i) shows the EM with constant growth, while (ii) includes a time-varying growth index. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 2. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the OM with constant growth. Each column represents an EM, where (i) shows the EM with constant growth, while (ii) includes a time-varying growth index. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Correctly specifying EM and the environmental index Most of the correctly specified scenarios (2-1-B, 3-1-B, 4-1-B, 5-1-B; Figure 3) were able to estimate SSB and recruitment accurately and precisely when growth was time-varying in the OM, with a similar bias to those from the scenarios with constant growth and either correct or over-parameterized models (Figure 2 and Table 4). The only exception was when β was set to a higher value and growth was highly varying (OM 4; Figure 3(iii)), which greatly decreased precision, even if the median SSB and recruitment estimates across simulations were fairly unbiased (median relative error <13% for all years). Figure 3. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the correctly specified EM (EM 2) given the correct index driving time-varying growth. Each column represents an OM with time-varying growth, where (i–iv) correspond to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 3. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the correctly specified EM (EM 2) given the correct index driving time-varying growth. Each column represents an OM with time-varying growth, where (i–iv) correspond to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Model misspecification (not allowing for time-varying growth) Not allowing for time-varying growth resulted in imprecise and biased estimates of SSB (Figure 4a), although recruitment was estimated precisely and without bias (Figure 4b). Growth parameters were biased, with K being negatively biased, and L∞ positively biased. The magnitude of bias in SSB and recruitment also increased with increasing β values (compare Figure 4(i) with 4(iii)). In this scenario, there was a large bias in the estimates of SSB and recruitment when the EM was misspecified (Figure 4(iii)), although the level of precision appears to be similar irrespective of whether the model was correctly specified or not (Figure 3(iii)). Misspecified EMs under OMs with high β values were the only scenarios where estimates of recruitment were biased (Figure 4b(iii)). Precision of SSB and recruitment also decreased in scenarios when autocorrelation within the environmental index increased and growth was assumed constant in the EM (Figure 4c) and Supplementary Figure S2). Figure 4. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the misspecified EM (EM 1) with time-invariant growth. Each column represents an OM with time-varying growth, where (i–iv) corresponding to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 4. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the misspecified EM (EM 1) with time-invariant growth. Each column represents an OM with time-varying growth, where (i–iv) corresponding to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Misspecifying the time-varying index We examined two ways in which the environmental index used in the EM differed from the true index driving the population dynamics. When the EM only had 20 out of the 108 years of the environmental data (observation model 2; scenarios 2-2-B and 5-2-B), the EM assumed that growth was constant up until the time series began in scenarios. Base growth parameters (K and L∞) tended to be under-estimated (Supplementary Table S2: rows 5 and 14). These errors increased when growth parameters varied more over time (comparing OM 2 to OM 5; Figure 5.1 and 5.2). SSB was generally biased over time in these cases, with both positive and negative bias during portions of the time series (Figure 5a). When growth was weakly climate-driven and only a shortened index was observed, median RE was positive at the start of the time series and became negative in the later years, and only tended to zero during the last few years of the time series (Figure 5.1a(iii)). Precision was similar to that of the correctly specified model for both SSB and recruitment (Figure 5.1(ii)). The precision of the SSB and recruitment estimates was much poorer for the scenario with strongly climate-driven growth and the shortened index data (Figure 5.2(iii)) than when growth was weakly climate-driven (Figure 5.1(iii)). SSB was over-estimated (up to 50% for some years) for this OM, although the median REs became closer to zero when the environmental data started. Figure 5. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the effects of varying the observation model. Each column represents a scenario, and the labels correspond to the scenario. Panel 1 shows results from the weakly climate-driven OM (OM 2), while panel 2 shows results from the strongly climate-driven OM (OM 5). Labels (i–v) correspond to observations of environmental index 0–4, with their respective EM. E.g. 1.i corresponds to scenario 2-0-A, 1.ii to 2-1-B, 1.iii to 2-2-B, 2.i to 5-0-B, etc. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 5. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the effects of varying the observation model. Each column represents a scenario, and the labels correspond to the scenario. Panel 1 shows results from the weakly climate-driven OM (OM 2), while panel 2 shows results from the strongly climate-driven OM (OM 5). Labels (i–v) correspond to observations of environmental index 0–4, with their respective EM. E.g. 1.i corresponds to scenario 2-0-A, 1.ii to 2-1-B, 1.iii to 2-2-B, 2.i to 5-0-B, etc. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Completely misspecifying the environmental index, either as an equally auto-correlated index or as random white noise (observation models 3 and 4), did not affect the estimation of the base growth parameters, SSB, and recruitment estimations when growth was weakly climate-driven (OM 2; Figure 5.1(iv) and 5.1 (v)). However, precision and accuracy of these derived quantity estimates greatly decreased when growth was strongly climate-driven (OM 5; Figure 5.2(iv) and 5.2(v)). The point estimates of the base growth parameters in those scenarios were biased by >10% under when the index was wrong and growth was misspecified (Table 4). In the sensitivity analysis using partially correct indices (Equation 4), bias in SSB and recruitment estimates remained relatively small when the index was correctly specified compared to being correctly specified up to 80% (a = 0.8), although precision decreased as the value of a decreased (Figure 6(ii-iii)). Thereafter, as a decreased, estimates of SSB became increasingly biased (again, both positively and negatively) and less precise than the model with a correct index (Figure 6(iv-vi)), although bias in estimates of recruitment remained relatively small until the index was completely misspecified (a = 0; Figure 6(vii)), where the estimates of SSB became highly negatively biased in the later years. Figure 6. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the results of the sensitivity analysis to the accuracy of the index for the OM with strongly climate-driven growth. Each column represents a scenario, and the labels correspond to the scenario. (i) corresponds to scenario 5-0-A, and (ii–vii) to scenarios where the a value (Equation 4) increases from 1 to 0, in 0.2 increments ((ii) also corresponds to scenario 5-1-B). The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 6. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the results of the sensitivity analysis to the accuracy of the index for the OM with strongly climate-driven growth. Each column represents a scenario, and the labels correspond to the scenario. (i) corresponds to scenario 5-0-A, and (ii–vii) to scenarios where the a value (Equation 4) increases from 1 to 0, in 0.2 increments ((ii) also corresponds to scenario 5-1-B). The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Sensitivity to high β values The magnitude of bias (in terms of median relative error) in the estimates of SSB and recruitment increased with the value for β when growth was assumed constant in the EM (Supplementary Figure S3). The pattern of bias (positive in the early years, and decreasing to become negative in the later years) remained the same across scenarios. Supplementary Figure S2 also shows that precision (measured using the inter-quantile range) of SSB and recruitment decreased as β was increased, regardless of whether allowance was made for time-varying growth in the EM. The problems in estimation found for high β values were probably caused by issues in the process of generating data. The high variance in growth caused fish to attain large sizes at young ages, and remain at that size during years where the index was at lower values, resulting in higher average sizes-at-age, as can be seen in Supplementary Table S3. Management implications Table 5 summarizes the implications of managing under a constant catch strategy, based on an estimate of MSY from monitoring data collected from the fishery and perhaps based on a misspecified assessment model. Median lost yield varied across scenarios, while median final stock status was similar, although inter-simulation variability differed. The proportion overfished differed across OMs, but was similar between EMs. Table 5. Table summary of all scenarios tested. Scenario label  Scenario name  Performance metric   Relative yield lost (interquartile range)  Est B2058/True B2058 (interquartile range)  Proportion of replicates overfished  1-0-A  Constant growth with correct EM  0.0249 (0.0682)  0.955 (0.11)  0.45  1-0-B  Constant growth with over-parameterized EM  0.0254 (0.0687)  0.955 (0.109)  0.45  2-0-A  Weakly climate-driven growth with incorrect EM  0.0666 (0.0859)  0.933 (0.106)  0.31  2-1-B  Weakly climate-driven growth with correct EM and index  0.0436 (0.0634)  0.959 (0.0837)  0.28  2-2-B  Weakly climate-driven growth with correct EM and shortened index  0.0918 (0.0818)  0.912 (0.126)  0.34  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  0.0633 (0.0919)  0.942 (0.121)  0.32  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  0.0628 (0.0922)  0.93 (0.094)  0.33  3-0-A  Highly autocorrelated growth with incorrect EM  0.0251 (0.192)  0.959 (0.293)  0.43  3-1-B  Highly autocorrelated growth with correct EM and index  0.0329 (0.0574)  0.952 (0.098)  0.47  4-0-A  Highly varied growth with incorrect EM  0.134 (0.289)  0.95 (0.171)  0.32  4-1-B  Highly varied growth with correct EM and index  0.131 (0.278)  0.941 (0.216)  0.33  5-0-A  Strongly climate-driven growth with incorrect EM  −0.0308 (0.758)  1.00 (0.555)  0.42  5-1-B  Strongly climate-driven growth with correct EM and index  0.0217 (0.1)  0.985 (0.105)  0.45  5-2-B  Strongly climate-driven growth with correct EM and shortened index  0.237 (0.577)  0.909 (0.681)  0.52  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  0.206 (1.57)  1.00 (0.844)  0.50  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  −0.097 (1.05)  1.01 (0.631)  0.47  Scenario label  Scenario name  Performance metric   Relative yield lost (interquartile range)  Est B2058/True B2058 (interquartile range)  Proportion of replicates overfished  1-0-A  Constant growth with correct EM  0.0249 (0.0682)  0.955 (0.11)  0.45  1-0-B  Constant growth with over-parameterized EM  0.0254 (0.0687)  0.955 (0.109)  0.45  2-0-A  Weakly climate-driven growth with incorrect EM  0.0666 (0.0859)  0.933 (0.106)  0.31  2-1-B  Weakly climate-driven growth with correct EM and index  0.0436 (0.0634)  0.959 (0.0837)  0.28  2-2-B  Weakly climate-driven growth with correct EM and shortened index  0.0918 (0.0818)  0.912 (0.126)  0.34  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  0.0633 (0.0919)  0.942 (0.121)  0.32  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  0.0628 (0.0922)  0.93 (0.094)  0.33  3-0-A  Highly autocorrelated growth with incorrect EM  0.0251 (0.192)  0.959 (0.293)  0.43  3-1-B  Highly autocorrelated growth with correct EM and index  0.0329 (0.0574)  0.952 (0.098)  0.47  4-0-A  Highly varied growth with incorrect EM  0.134 (0.289)  0.95 (0.171)  0.32  4-1-B  Highly varied growth with correct EM and index  0.131 (0.278)  0.941 (0.216)  0.33  5-0-A  Strongly climate-driven growth with incorrect EM  −0.0308 (0.758)  1.00 (0.555)  0.42  5-1-B  Strongly climate-driven growth with correct EM and index  0.0217 (0.1)  0.985 (0.105)  0.45  5-2-B  Strongly climate-driven growth with correct EM and shortened index  0.237 (0.577)  0.909 (0.681)  0.52  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  0.206 (1.57)  1.00 (0.844)  0.50  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  −0.097 (1.05)  1.01 (0.631)  0.47  The proportions of replicates that resulted in overfished populations ranged from 0.28 to 0.51 across the scenarios (Table 5). Some of the randomly generated recruitment deviations resulted in populations with extreme spawning stock biomass levels, both high and low, with low levels causing populations to crash. Projecting a constant catch strategy for populations with low biomass may have also contributed to the relatively high proportion of replicates being considered overfished, as shown by the 0.45 value from scenario 1-0-A, when the assessment was correctly specified and growth was time-invariant. Most EMs tended to provide MSY estimates that were median-unbiased, given the correct index data, regardless of whether time-varying growth was accounted for. However, inter-simulation variability (measured in terms of inter-quantile range), increased greatly when the EM was misspecified, i.e. growth was assumed to be constant (2-0-A, 3-0-A, 5-0-A). The exception to both these results was the OM 4, where β was set to a high level. In that case, there were negligible differences in MSY estimates between the misspecified EM with time-invariant growth and the correctly specified EM (Table 5, scenarios 4-O-A and 4-O-B). The estimates of MSY were positively biased when the EM was based on an incorrect environmental index (scenarios 2-3-B, 2-4-B, 5-3-B, 5-4-B). This resulted in up to a median of 24% of yield being lost in the OM where growth is strongly climate-driven (OM 5; Table 5). Inter-simulation variability was very high for the scenarios based on OM 5, while this was not the case for the scenarios where growth was weakly climate-driven (OM 2). Supplementary Table S4 shows an increase in inter-simulation variability of management quantities as the accuracy of the index decreased, although bias remained the same across scenarios. Model misspecification (i.e. assuming time-invariant growth when growth was actually varying over time) resulted in more biased and less precise estimates of MSY (Table 4), resulting in a larger median and wider range of lost yield across replicates, respectively (Table 5). Table 5 also shows that the inter-quantile range of lost yield (i.e. the imprecision of MSY estimates) decreased when the EM was correctly specified, particularly when growth varied over time (e.g. comparing 3-0-A vs. 3-0-B). Discussion Incorporating accurate growth indices resulted in higher precision and less bias in parameter estimates and management outcomes. These effects are more pronounced if the environmental index strongly influenced growth, was highly auto-correlated, and/or had high variance, and was accurately described. The degree of bias and inaccuracies are largely defined by the variance of the time-varying index, and the strength of its auto-correlation. When growth was either weakly time-varying or constant, including an inaccurate index of low variance (random noise) carried similar risks of obtaining inaccurate or imprecise abundance estimates to estimating time-invariant growth and not including an index at all. On the other hand, when growth varied over time with high variance, including an inaccurate index (a < 0.6) or not including one at all resulted in highly biased and inaccurate estimates of stock size. Estimating time-varying growth when the growth index is correct results in the most precise and least biased estimates. Allowing for time-varying growth in the form of an index does not risk significant degradation in estimation performance even when growth is not actually time-varying, although this scenario is unlikely (e.g. Stawitz et al., 2015; Thorson and Minte-Vera, 2016). Therefore, if the growth index under consideration is weakly auto-correlated and weakly linked to the growth parameters, inclusion of the index would likely improve the performance of the model, regardless of how accurate the index is. Based on the results for OM 2, the most “realistic” OM, and the properties of the growth index developed in Black (2009), we would recommend consideration of including this index in the assessment of splitnose rockfish. If a climate-driven growth index were to be considered for use in an assessment, we would recommend fitting an AR-1 model to the index, using these parameters to generate new simulated growth indices, and simulation-testing their inclusion in the stock assessment model. The results were obtained from only one set of life history parameters, but the simulation-based risk analysis approach could be conditioned on the data and biology of other stocks to illustrate whether to include a growth index for a different specific stock assessment. A growth index for which 60% of the variance is correct (0.77 correlation with true variation in growth; √(a=0.60)) could still lead to estimates of parameters and derived quantities that have similar bias and precision as not including the index at all, when growth is strongly time-varying (OM 5). This case study suggests that, when growth is suspected to be strongly climate-driven, environmental indices should only be used when there is high confidence in their accuracy. For example, Thorson and Minte-Vera (2016) found that the best models with a single time-varying effect could explain on average 50% of variability in weight-at-age for exploited stocks, with year-effects (i.e. annual indices like those explored here) explaining approximately 45% of variability. Stawitz et al. (2015) described a state-space model that was fit to fishery-dependent and -independent data for 29 stocks of Pacific groundfish, and found that approximately 40% of the modelled stocks exhibited time-varying growth. Whereas it is highly unlikely that it is possible to observe the true underlying mechanism by which growth varies over time and generate an exact environmental index, assuming time-invariant growth (status quo) could carry similar risks of biased and imprecise assessment outputs. A feasible method of increasing the accuracy of estimated growth variation would be by examining similarities in growth variation among species (Stawitz et al., 2015) and/or cross-validating particular years of strong or weak growth with other species (Black, 2009). Additionally, Maunder and Watters (2003) recommend an integrated approach with additional process error when including an environmental time series into stock assessment models, which should be a subject for future research. This approach could presumably be implemented using annual variation in growth, in addition to including the index multiplicatively, and the magnitude of additional variation in growth could be estimated using the Laplace approximation (Thorson et al., 2015a). Improving the accuracy of growth parameter estimates had a relatively small contribution to total predicted MSY. Furthermore, erroneous parameter and management quantity estimates had negligible impact on future stock size in this case study. This could have been an artefact of having both natural mortality M and steepness h parameters fixed in the assessment. These parameters exhibit a strong influence on reference point estimation, as shown in Mangel et al. (2013), and as such would limit the extent by which growth parameters could affect estimates of MSY. Furthermore, the studied species is long-lived and slow-growing, with a relatively high age of maturity, and a short time series of age data. MSY for long-lived species is expected to have similar sensitivity to changing growth parameters as short-lived species (Thorson et al., 2015b). Thorson et al. (2015b) found that the MSY from the three life histories tested showed similar sensitivity to time-varying biological parameters (for recruitment, size, and maturity), while sensitivities differed when natural mortality varied over time, so we do not expect MSY to be more sensitive to growth variation for shorter-lived species. Other management reference points (such as BMSY/B0) may show different sensitivities to time-varying parameters for each life history. Future work should thus consider assessments where there is more data on age and growth, and/or these parameters are not pre-specified, or based on species where growth may vary more strongly over time (e.g. California Current petrale sole Eopsetta jordani or Gulf of Alaska walleye pollock; Stawitz et al., 2015) to more fully determine the effects of incorporating a time-varying index on growth. Larger scale effects could also be tested by using alternate forms of time-varying growth, such as the sine-wave model or AR models with more lags. In addition, growth also varies at a spatial (Rahikainen and Stephenson, 2004; Thorson, 2015) and at an individual (Webber and Thorson, 2016) scale, which could increase model complexity to explain more variance in the data. Implementing constant catch exactly (i.e. without error) is unlikely to be the best strategy to use for examining the management impacts of including a time-varying growth index in assessments. The risk analysis approach can be extended to be based on alternative management strategies to perform a full management strategy evaluation (Punt et al., 2016). Conducting an MSE would allow managers to better assess the consequences of including a time-varying growth index in a management strategy that applies harvest control rules (Punt et al., 2016). Punt et al. (2014) summarize studies where this has been done for recruitment. There has been a recent push to move from single-species towards ecosystem-based fisheries management (EBFM; Pikitch et al., 2004), where environmental factors are taken into consideration to manage the ecosystem as a whole. Therefore, there is a desire to explore various methods for incorporating time-varying growth into assessments. The simulation approach described here can be used as a step towards addressing and incorporating climate effects into somatic growth for a single species. Supplementary data Supplementary material is available at the ICESJMS online version of the manuscript. Acknowledgements This work was funded by a NOAA Fisheries and the Environment (FATE) grant, project 13-01. The authors thank Timothy Essington, Melissa Haltuch, James Hastie, Michelle McClure, Felipe Hurtado Ferro, Ernesto Jardim, Kelli Johnson, Ian Taylor, Chantel Wetzel, and an anonymous reviewer for their helpful suggestions and comments throughout the working process and on the manuscript. Thanks also go to Bryan Black for sharing his otolith index data. References Ait Youcef W., Lambert Y., Audet C. 2015. Variations in length and growth of Greenland Halibut juveniles in relation to environmental conditions. Fisheries Research , 167: 38– 47. Google Scholar CrossRef Search ADS   A’mar Z. T., Punt A. E., Dorn M. W. 2009. The evaluation of two management strategies for the Gulf of Alaska walleye pollock fishery under climate change. ICES Journal of Marine Science , 66: 1614– 1632. Google Scholar CrossRef Search ADS   Black B. A. 2009. Climate-driven synchrony across tree, bivalve, and rockfish growth-increment chronologies of the northeast Pacific. The Marine Ecology Progress Series , 378: 37– 46. Google Scholar CrossRef Search ADS   Black B. A., Boehlert G. W., Yoklavich M. M. 2008. Establishing climate–growth relationships for yelloweye rockfish (Sebastes ruberrimus) in the northeast Pacific using a dendrochronological approach. Fisheries Oceanography , 17: 368– 379. Google Scholar CrossRef Search ADS   Brunel T., Dickey-Collas M. 2010. Effects of temperature and population density on von Bertalanffy growth parameters in Atlantic herring: a macro-ecological analysis. The Marine Ecology Progress Series , 405: 15– 28. Google Scholar CrossRef Search ADS   Essington T. E., Kitchell J. F., Walters C. J. 2001. The von Bertalanffy growth function, bioenergetics, and the consumption rates of fish. Canadian Journal of Fisheries and Aquatic Sciences , 58: 2129– 2138. Google Scholar CrossRef Search ADS   Fournier D., Archibald C. P. 1982. A General Theory for Analyzing Catch at Age Data. Canadian Journal of Fisheries and Aquatic Sciences , 39: 1195– 1207. Google Scholar CrossRef Search ADS   Francis R. I. C. C. 2016. Growth in age-structured stock assessment models. Fisheries. Research , 180: 77– 86. Google Scholar CrossRef Search ADS   Gertseva V. V., Cope J. M., Pearson D. E. 2009. Status of the U.S. splitnose rockfish (Sebastes diploproa) resource in 2009, Status of the Pacific Coast Groundfish Fishery through 2009 and recommended acceptable biological catches for 2011/2012: stock assessment and fishery evaluation, 2009. Pacific Fishery Management Council, Portland, OR. Haltuch M. A., Punt A. E. 2011. The promises and pitfalls of including decadal-scale climate forcing of recruitment in groundfish stock assessment. Canadian Journal of Fisheries and Aquatic Sciences , 68: 912– 926. Google Scholar CrossRef Search ADS   Harley C. D. G., Randall Hughes A., Hultgren K. M., Miner B. G., Sorte C. J. B., Thornber C. S., Rodriguez L. F., Tomanek L., Williams S. L. 2006. The impacts of climate change in coastal marine systems. Ecology Letters , 9: 228– 241. Google Scholar CrossRef Search ADS PubMed  Hulson P.-J. F., Quinn T. J., Hanselman D. H., Ianelli J. N. 2013. Spatial modeling of Bering Sea walleye pollock with integrated age-structured assessment models in a changing environment. Canadian Journal of Fisheries and Aquatic Sciences , 70: 1402– 1416. Google Scholar CrossRef Search ADS   Johnson K. F., Monnahan C. C., McGilliard C. R., Vert-pre K. A., Anderson S. C., Cunningham C. J., Hurtado-Ferro F., Licandeo R. R., Muradian M. L., Ono K., et al.   2015. Time-varying natural mortality in fisheries stock assessment models: identifying a default approach. ICES Journal of Marine Science , 72: 137–150. Kimura D. K. 2008. Extending the von Bertalanffy growth model using explanatory variables. Canadian Journal of Fisheries and Aquatic Sciences , 65: 1879– 1891. Google Scholar CrossRef Search ADS   Lorenzen K. 2016. Toward a new paradigm for growth modeling in fisheries stock assessments: embracing plasticity and its consequences. Fisheries Research , 180: 4– 22. Google Scholar CrossRef Search ADS   Mangel M., MacCall A. D., Brodziak J., Dick E. J., Forrest R. E., Pourzand R., Ralston S. 2013. A perspective on steepness, reference points, and stock assessment. Canadian Journal of Fisheries and Aquatic Sciences , 70: 930– 940. Google Scholar CrossRef Search ADS   Maunder M. N., Punt A. E. 2013. A review of integrated analysis in fisheries stock assessment. Fisheries Research , 142: 61– 74. Google Scholar CrossRef Search ADS   Maunder M. N., Watters G. M. 2003. A general framework for integrating environmental time series into stock assessment models: model description, simulation testing, and example. Fisheries Bulletin , 101: 89– 99. Methot R. D., Taylor I. G. 2011. Adjusting for bias due to variability of estimated recruitments in fishery assessment models. Canadian Journal of Fisheries and Aquatic Sciences , 68: 1744– 1760. Google Scholar CrossRef Search ADS   Methot R. D., Wetzel C. R. 2013. Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research , 142: 86– 99. Google Scholar CrossRef Search ADS   Mueter F. J., Bond N. A., Ianelli J. N., Hollowed A. B. 2011. Expected declines in recruitment of walleye pollock (Theragra chalcogramma) in the eastern Bering Sea under future climate change. Journal of Marine Science , 68: 1284– 1296. Pikitch E. K., Santora C., Babcock E. A., Bakun A., Bonfil R., Conover D. O., Dayton P., Doukakis P., Fluharty D., Heneman B., et al.   2004. Ecosystem-based fishery management. Science , 305: 346– 347. Google Scholar CrossRef Search ADS PubMed  Pinsky M. L., Worm B., Fogarty M. J., Sarmiento J. L., Levin S. A. 2013. Marine taxa track local climate velocities. Science , 341: 1239– 1242. Google Scholar CrossRef Search ADS PubMed  Punt A. E., A’mar T., Bond N. A., Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haltuch M. A., Hollowed A. B., Szuwalski C. 2014. Fisheries management under climate and environmental uncertainty: control rules and performance simulation. ICES Journal of Marine Science , 71: 2208– 2220. Google Scholar CrossRef Search ADS   Punt A. E., Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haddon M. 2016. Management strategy evaluation: best practices. Fish and Fisheries , 17: 303– 334. Google Scholar CrossRef Search ADS   Rahikainen M., Stephenson R. L. 2004. Consequences of growth variation in northern Baltic herring for assessment and management. ICES Journal of Marine Science , 61: 338– 350. Google Scholar CrossRef Search ADS   Rowe R. J., Terry R. C. 2014. Small mammal responses to environmental change: integrating past and present dynamics. Journal of Mammalogy , 95: 1157– 1174. Google Scholar CrossRef Search ADS   Schirripa M. J., Goodyear C. P., Methot R. M. 2009. Testing different methods of incorporating climate data into the assessment of US West Coast sablefish. Journal of Marine Science , 66: 1605– 1613. Stawitz C. C., Essington T. E., Branch T. A., Haltuch M. A., Hollowed A. B., Spencer P. D. 2015. A state-space approach for detecting growth variation and application to North Pacific groundfish. Canadian Journal of Fisheries and Aquatic Sciences , 72: 1316– 1328. Google Scholar CrossRef Search ADS   Sünksen K., Stenberg C., Grønkjær P. 2010. Temperature effects on growth of juvenile Greenland halibut (Reinhardtius hippoglossoides Walbaum) in West Greenland waters. Journal of Sea Research , 64: 125– 132. Google Scholar CrossRef Search ADS   Taylor I. G., Methot R. D. 2013. Hiding or dead? A computationally efficient model of selective fisheries mortality. Fisheries Research , 142: 75– 85. Google Scholar CrossRef Search ADS   Thorson J. T. 2015. Spatio-temporal variation in fish condition is not consistently explained by density, temperature, or season for California Current groundfishes. Marine Ecology Progress Series , 526: 101– 112. Google Scholar CrossRef Search ADS   Thorson J. T., Hicks A. C., Methot R. D. 2015a. Random effect estimation of time-varying factors in Stock Synthesis. Journal of Marine Science , 72: 178– 185. Thorson J. T., Jensen O. P., Zipkin E. F. 2014. How variable is recruitment for exploited marine fishes? A hierarchical model for testing life history theory. Canadian Journal of Fisheries and Aquatic Sciences , 71: 973– 983. Google Scholar CrossRef Search ADS   Thorson J. T., Minte-Vera C. V. 2016. Relative magnitude of cohort, age, and year effects on size at age of exploited marine fishes. Fisheries Research , 180: 45– 53. Google Scholar CrossRef Search ADS   Thorson J. T., Monnahan C. C., Cope J. M. 2015b. The potential impact of time-variation in vital rates on fisheries management targets for marine fishes. Fisheries Research , 169: 8– 17. Google Scholar CrossRef Search ADS   Thorson J. T., Stewart I. J., Taylor I. G., Punt A. E. 2013. Using a recruitment-linked multispecies stock assessment model to estimate common trends in recruitment for US West Coast groundfishes. The Marine Ecology Progress Series , 483: 245– 256. Google Scholar CrossRef Search ADS   Von Bertalanffy L. 1957. Quantitative laws in metabolism and growth. The Quarterly Review of Biology, 32: 217– 231. Weatherley A. H. 1990. Approaches to understanding fish growth. Transactions of the American Fisheries Society , 119: 662– 672. Google Scholar CrossRef Search ADS   Webber D. N., Thorson J. T. 2016. Variation in growth among individuals and over time: a case study and simulation experiment involving tagged Antarctic toothfish. Fisheries Research , 180: 67– 76. Google Scholar CrossRef Search ADS   © International Council for the Exploration of the Sea 2017. All rights reserved. For Permissions, please email: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png ICES Journal of Marine Science Oxford University Press

The benefits and risks of incorporating climate-driven growth variation into stock assessment models, with application to Splitnose Rockfish (Sebastes diploproa)

Loading next page...
 
/lp/ou_press/the-benefits-and-risks-of-incorporating-climate-driven-growth-RyehNi58sG
Publisher
ICES/CIEM
Copyright
© International Council for the Exploration of the Sea 2017. All rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
1054-3139
eISSN
1095-9289
D.O.I.
10.1093/icesjms/fsx147
Publisher site
See Article on Publisher Site

Abstract

Abstract Indices of annual growth variation are not routinely incorporated into fisheries stock assessment models, due to a lack of a general framework for deciding when to include these indices, and of a mechanistic understanding about growth drivers. Such incorporation may also not necessarily lead to improved estimation or management performance. We demonstrate a way to incorporate such an index into an assessment model (Stock Synthesis), and use risk analysis to evaluate its management-related advantages and shortcomings. We applied this method to splitnose rockfish (Sebastes diploproa), where a previously developed growth index is highly correlated with decadal-scale climate indices. We find that including a similar index in the simulated assessment increases precision and reduces bias of parameter estimates. However, not including an index or including a completely erroneous index led to highly imprecise estimates when growth was strongly climate-driven. Including this growth index when individual growth was actually constant did not lead to poorer estimation performance. The risk analysis approach can be applied to other stocks to evaluate the consequences of including an index of growth variation. Introduction Climate and environmental factors can have a large impact on population and community dynamics (Harley et al., 2006; Pinsky et al., 2013; Rowe and Terry, 2014). One of the ways in which the environment can play a major role in determining the dynamics of fish populations is through the metabolism, and ultimately growth, of individual fish (Ait Youcef et al., 2015; Sünksen et al., 2010). However, little guidance is available for whether or when to incorporate environmentally driven growth variables into the stock assessments used as the basis for advice to managers to support fisheries decision making. Such inclusion could more accurately describe the dynamics of fish populations, and potentially disentangle whether fluctuations in biomass are due to time-varying environmental conditions or fishing. Much of our current understanding of how to best incorporate environmental information into stock assessments derives from recruitment-environment relationships (e.g. Mueter et al., 2011; Schirripa et al., 2009; Thorson et al., 2013), as opposed to growth variation. Growth is an integral component of fisheries stock assessment, and the relationship between age and size is often estimated in assessments (Taylor and Methot, 2013), based on the assumption that the parameters of this relationship are constant over space and time (Lorenzen, 2016). However, it has been shown (e.g. Stawitz et al., 2015; Weatherley, 1990; Webber and Thorson, 2016) that growth in fish is inherently plastic, and that this “constant growth” paradigm needs to be challenged (Lorenzen, 2016; Thorson and Minte-Vera, 2016). Defining climate-growth relationships would allow fishery managers to better discriminate the effects of climate from those of fishing, and make more accurate predictions of future stock size, should there be a means of projecting the relationship into the future (Black et al., 2008). Black (2009) examined ring widths of rockfish (Sebastes spp) otoliths and reduced these widths to species-specific dimensionless indices that are correlated with similar dendrochronology indices from geoducks (Panopea abrupta) and trees in the same general area. These indices were also found to coincide with climate indices in the California Current System. However, these growth indices have never been incorporated in assessments due, inter alia, to the lack of a general framework for deciding whether and under what circumstances to do so. Many studies have used simulation to evaluate the benefits of incorporating indices of climate on recruitment in assessments, notably Schirripa et al. (2009) for sablefish, Hulson et al. (2013) and A’mar et al. (2009) for walleye pollock, and Haltuch and Punt (2011) for groundfish generally. Simulations are used to evaluate model performance and test the effects of different hypotheses on management outcomes (e.g. Johnson et al., 2015). However, this has yet to be done for time-varying growth. Maunder and Watters (2003) proposed a general framework for integrating environmental time series in assessments, where simulations are used to test and compare a model that includes environmental effects with the traditional approach where environment variables are ignored. Maunder and Watters (2003) focused on the ability to correctly identify the relationship between a selected variable in the model and the environmental variable, and to accurately estimate specific parameters, and not on the impact of these relationships on management outcomes. The goal of this study is to propose and illustrate a risk-analysis approach to estimating the benefits and risks of incorporating a climate index on somatic fish growth into stock assessments. Our analyses are conditioned closely upon a climate index developed by Black (2009) and the stock assessment for splitnose rockfish (Sebastes diploproa) in the Northeast Pacific Ocean (Gertseva et al., 2009), and evaluate management outcomes associated with various climate-growth scenarios. We expand on Maunder and Watters (2003) and develop a table using simulation modelling to provide guidance on whether incorporation of climate information is justified in a given case. Material and methods Overview The risk-analysis was developed using a widely used stock assessment platform—Stock Synthesis (SS; Methot and Wetzel, 2013). SS can incorporate all available data into a single analysis, and estimate parameters by maximizing the product of the likelihood function for each data type (Fournier and Archibald, 1982; Maunder and Punt, 2013). The assessment was simplified to ensure that results are widely applicable to other contexts. The numbers of fishing fleets and surveys were reduced to one of each data type. Operating models (OMs) were used to describe “true” population dynamics and quantity of data for five states of nature—one with constant growth, and four with time-varying environmentally driven growth indices that vary in terms of auto-correlation and relationship to the OM growth parameters. Each simulation was projected for 50 years using the OM (with the relevant assumptions made for each OM, i.e. constant future growth for one and time-varying future growth for the other four), with fishing occurring under a constant-catch limit obtained from the OM. Implementation error was not simulated, so catch was assumed to be taken exactly. The OM outputs based in the historical and forecast periods were considered “true” for the simulations. Observation models simulated the data-gathering process. A total of 1100 data sets were generated, with 100 data sets for each combination of OM and observation model. Each data set was analysed using two estimation methods (EMs). We examined the performance of two EMs—one that estimated time-invariant growth, and one that assumed growth was correlated with an available environmental index. Management reference points were estimated using these EMs, and used to set a catch limit that was projected forward as a constant in the simulated fishery system (as represented in the OM). The resultant forecasted outcomes were then compared to the “true” (i.e. OM) quantities to obtain performance metrics, such as lost yield and final stock biomass. A summary of the overall process can be found in Figure 1. All code for replicating this analysis is publicly available (https://github.com/lee-qi/ss3-growth-index). Figure 1. View largeDownload slide Flowchart depicting the general steps of the simulation process. Multiple operating and observation models were considered, but this is not shown here. Figure 1. View largeDownload slide Flowchart depicting the general steps of the simulation process. Multiple operating and observation models were considered, but this is not shown here. Incorporating climate indices of individual growth into assessment models As the growth model most commonly used in stock assessments (Essington et al., 2001; Francis, 2016), this study was based, for simplicity, on the von Bertalanffy growth function (Von Bertalanffy, 1957):   Lt=L∞(1-e-Kt-t0e-Kt-t0) (1) where Lt is the length of the fish at time t, L∞ is the asymptotic length of the fish and a ratio of anabolic to catabolic factors, K is the intrinsic growth rate and a catabolic factor, and t0 is a constant that adjusts the model along a time axis. Temperature is known to affect catabolism, and thus both K and L∞ (Brunel and Dickey-Collas, 2010; Lorenzen, 2016), with Kimura (2008) showing a negative correlation between the two parameters. Temporal variation in growth due to changes in environmental drivers can be incorporated in an assessment by modifying the values of the parameters of the growth function. In SS, an environmental index can be linked to the values of the growth parameters in two ways—additively or multiplicatively (Methot and Wetzel, 2013). In this study, we included an annual index ɛt that is associated with variation in growth, using a multiplicative relationship (Equation 2) for both maximum size for sex s:   L∞,s,t=L∞,s,base·eβL,sɛt (2a) and growth rate:   Ks,t=Ks,base·eβK,sɛt (2b) where β is the parameter linking the growth parameter to the environmental index, L∞,s,base and Ks,base are median values for the asymptotic length and growth rate parameters for sex s, and length-at-age in each year is calculated from length-at-age the previous year (using Equation 1). Should the index value be sufficiently large to result in a negative growth increment (e.g. if L∞,t <  Lt from Equation 1), the individuals remain at the same size for the following year and are assumed to not experience any growth, which thus avoids the occurrence of “shrinking” fish. Biological characteristics Populations were simulated based on parameters obtained from a simplified version of the assessment (Supplementary Table S1; Gertseva et al., 2009). The data sets with the longest time series for a single fishing fleet and a single fishery-independent survey were selected (Supplementary Figure S1c). In the assessment approved for management (Gertseva et al., 2009), the parameters governing natural mortality (M) and steepness (h) were pre-specified, while log initial recruitment (lnR0) and individual-growth parameters including length at age 0 (L0), L∞, K, and variability in size-at-age for fish young fish (CV1) and for old fish up to L∞ ( CV2) were estimated. Growth of splitnose rockfish is sex-specific, with males estimated to grow at a faster rate (i.e. higher Kmale), while females are estimated to attain larger sizes (i.e. higher L∞,female). This is also one of the two rockfish species for which a climate-growth index was developed by Black (2009). Simulation method Index data generation The environmental index ɛt was simulated using an autoregressive (AR) model with bias correction (Equation 3), adapted from Thorson et al. (2014), with parameters calculated for splitnose rockfish. The environmental index was unique to each replicate within each state of nature, and scaled to a mean of 0.   ɛt=ρsplitnoseɛt-1+ 1-ρsplitnose2 ωtωtρsplitnoseɛt-1+ ωtωtfor t>1for t=1  ωt ∼ Normal-σsplitnose22×1-ρsplitnose1-ρsplitnose2, σsplitnose2 (3) where ωt is independent and identically normally distributed variation for year t,  ρsplitnose is the species-specific autoregressive parameter, and σsplitnose2 is the variance of the index. As the index is treated as data, the species-specific parameters in Equation 3 ( ρsplitnose and σsplitnose2) were defined external to the OM, and fixed. Comparative runs (across OMs) used the same sequence of random numbers. Five trends in growth were investigated, based on the extent of auto-correlation (i.e. ρsplitnose) in the environmental indices and the strength of the relationship between the environmental index and the growth parameters (i.e. βL,s and βK,s; Equation 2; Table 1). Sensitivity analyses were conducted for a series of β and ρsplitnose values in the OM, using OM 2 (weakly climate-driven growth) as the base case. As the values for β increase, so too does the impact of the environmental index (and vice-versa). The series of β values tested ranged from 0 to 0.5, in 0.1 increments. The higher the value of ρsplitnose, the more that long periods of time have a series of fast or slow growth, should the index be highly correlated with growth. We tested scenarios with indices of varying levels of auto-correlation, ranging the ρsplitnose values from 0.0 (i.e. white noise) to 0.8 by 0.2 increments, and a final maximum high value of 0.95. Table 1. Summary of OMs, and their associated values. Number  OM  Strength of relationship  βL,m  βL,m  βK,f  βK,m  σp2  Degree of auto-correlation ( ρp)  1  Not Climate-Driven  No relationship  0  0  0  0  0.88  Weak (0.33)  2  Weakly Climate-Driven  Weak  −0.12  −0.083  0.13  0.22  0.88  Weak (0.33)  3  Climate-Driven (High AR)  Weak  −0.12  −0.083  0.13  0.22  0.88  Strong (0.95)  4  Climate-Driven (High Beta)  Strong  −0.3o  −0.3o  0.3o  0.3o  0.88  Weak (0.33)  5  Climate-Driven (High AR and High Beta)  Strong  −0.30  −0.3o  0.3o  0.3o  0.88  Strong (0.95)  Number  OM  Strength of relationship  βL,m  βL,m  βK,f  βK,m  σp2  Degree of auto-correlation ( ρp)  1  Not Climate-Driven  No relationship  0  0  0  0  0.88  Weak (0.33)  2  Weakly Climate-Driven  Weak  −0.12  −0.083  0.13  0.22  0.88  Weak (0.33)  3  Climate-Driven (High AR)  Weak  −0.12  −0.083  0.13  0.22  0.88  Strong (0.95)  4  Climate-Driven (High Beta)  Strong  −0.3o  −0.3o  0.3o  0.3o  0.88  Weak (0.33)  5  Climate-Driven (High AR and High Beta)  Strong  −0.30  −0.3o  0.3o  0.3o  0.88  Strong (0.95)  The values for β for OMs 2 and 3 were estimated by including the otolith index data as an index on the growth parameters, and applying the assessment to these data. The values for β for OMs 4 and 5 were fixed at a high value of 0.3, as determined by the sensitivity analysis, with positive values for K and negative values for L∞, consistent with the notion of having opposite effects. ρsplitnose and σsplitnose2 in OMs 2 and 4 were determined by fitting an AR-1 model to the otolith index obtained from Black (pers. comm.). As such, the environmental index in OM 2 had the most realistic characteristics, while OMs 3–5 were exploratory. Observation errors Observation error was added to the expected values for each data type in the OMs to generate data sets for use in the EMs (Table 2). Additional variability for each simulation was included in the form of independent and lognormally distributed yearly deviations about a Beverton-Holt stock-recruitment relationship (Methot and Wetzel, 2013). A bias-correction factor (Methot and Taylor, 2011) was applied in the EM to avoid bias in estimation of expected recruitment when there is less data. Table 2. Data types and error structures used in data generation (Gertseva et al., 2009). Data type  Years  Error distribution  Measure of uncertainty  Catch—Fishery  1900–2008  Lognormal  SE of log(value) = 0.01  Abundance indices—Survey  2003–2008  Lognormal  SE of log(value) = 0.25–0.34  Length compositions—Fishery  1978–2008  Multinomial  Effective N = 505.4  Length compositions—Fishery (Discard)  1987; 2006–2007  Multinomial  Effective N = 505.4  Length compositions—Survey  2003–2008  Multinomial  Effective N = 1339  Conditional age-at-length compositions—Survey  2003–2008  Multinomial  Effective N = 124.3  Mean body weight—Fishery  2002–2008  Lognormal  CV = 0.3  Discards—Fishery  1987; 2002–2007  Lognormal  CV = 0.5 (1987); 0.2 (2002–2007)  Environmental index—Black (2009)  1936–2006  Known exactly    Data type  Years  Error distribution  Measure of uncertainty  Catch—Fishery  1900–2008  Lognormal  SE of log(value) = 0.01  Abundance indices—Survey  2003–2008  Lognormal  SE of log(value) = 0.25–0.34  Length compositions—Fishery  1978–2008  Multinomial  Effective N = 505.4  Length compositions—Fishery (Discard)  1987; 2006–2007  Multinomial  Effective N = 505.4  Length compositions—Survey  2003–2008  Multinomial  Effective N = 1339  Conditional age-at-length compositions—Survey  2003–2008  Multinomial  Effective N = 124.3  Mean body weight—Fishery  2002–2008  Lognormal  CV = 0.3  Discards—Fishery  1987; 2002–2007  Lognormal  CV = 0.5 (1987); 0.2 (2002–2007)  Environmental index—Black (2009)  1936–2006  Known exactly    The effects of including an incorrectly observed environmental index were examined using four observation models: 158 years of environmental indices (108 historical, 50 projected) observed without error; i.e. the index used to generate the 20 years of historical environmental indices observed without error, with 50 years of projected indices, i.e. a shortened time series of data when there is insufficient historical information 158 years of auto-correlated indices generated using the same model and values for ρp and σp2 as the “true” index, but with a different starting seed, i.e. a time series of auto-correlated indices that seems plausible but is completely inaccurate 158 years of an independent time series of random errors (white noise), generated under N(0,1) These scenarios were similar to those described by Stawitz et al. (2015), where growth anomaly patterns were estimated to be in one of three predominant forms—trendless, sustained trend, and near-zero. Sensitivity of the EM to the accuracy of the environmental index was also examined, using OM 5 as a base. This was done by creating a second auto-correlated index ( ɛt,wrong) using the same AR parameters and process (Equation 3), and adding it to the true index ( ɛt,true) in varying proportions a, which ranged from 1 to 0 at 0.2 increments, to create a new index ( ɛt,new):   ɛt,new=a·ɛt,true+1-a·ɛt,wrong (4) This procedure therefore leads to a new environmental index for which a·100% of the variance is correct, and (1-a)·100% of the variance is auto-correlated noise (i.e. the correlation between ɛt,true and ɛt,new is a). The new index ɛt,wrong is calculated to have the same variance and first-order autocorrelation as the true index. Estimation methods The estimation scenarios mimicked an actual assessment process, by estimating model parameters and management reference points. Supplementary Table S1 lists the key parameters of the population dynamics model and how they are estimated in the EM. Priors for the growth parameters were set to the same values as the priors in the OMs. Two EMs were used—one where the index was included ( β is estimated), and one where the index was not ( β=0). Each combination of OM, observation model, and EM is referred to as a scenario, and 16 scenarios were evaluated (Table 3). Table 3. Summary of scenarios tested. OM  Observation of Env index  EM  Scenario label  Scenario name  1  Not Climate-Driven  NA  Constant  1-0-A  Constant growth with correct EM  Wrong  Time-varying  1-0-B  Constant growth with over-parameterized EM  2  Weakly Climate-Driven  NA  Constant  2-0-A  Weakly climate-driven growth with incorrect EM  Exact  Time-varying  2-1-B  Weakly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  2-2-B  Weakly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  3  Climate-Driven (High AR)  NA  Constant  3-0-A  Highly autocorrelated growth with incorrect EM  Exact  Time-varying  3-1-B  Highly autocorrelated growth with correct EM and index  4  Climate-Driven (High Beta)  NA  Constant  4-0-A  Highly varied growth with incorrect EM  Exact  Time-varying  4-1-B  Highly varied growth with correct EM and index  5  Climate-Driven (High AR and High Beta)  NA  Constant  5-0-A  Strongly climate-driven growth with incorrect EM  Exact  Time-varying  5-1-B  Strongly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  5-2-B  Strongly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  OM  Observation of Env index  EM  Scenario label  Scenario name  1  Not Climate-Driven  NA  Constant  1-0-A  Constant growth with correct EM  Wrong  Time-varying  1-0-B  Constant growth with over-parameterized EM  2  Weakly Climate-Driven  NA  Constant  2-0-A  Weakly climate-driven growth with incorrect EM  Exact  Time-varying  2-1-B  Weakly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  2-2-B  Weakly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  3  Climate-Driven (High AR)  NA  Constant  3-0-A  Highly autocorrelated growth with incorrect EM  Exact  Time-varying  3-1-B  Highly autocorrelated growth with correct EM and index  4  Climate-Driven (High Beta)  NA  Constant  4-0-A  Highly varied growth with incorrect EM  Exact  Time-varying  4-1-B  Highly varied growth with correct EM and index  5  Climate-Driven (High AR and High Beta)  NA  Constant  5-0-A  Strongly climate-driven growth with incorrect EM  Exact  Time-varying  5-1-B  Strongly climate-driven growth with correct EM and index  Exact (final 20 years only)  Time-varying  5-2-B  Strongly climate-driven growth with correct EM and shortened index  Wrong (AR-1 model)  Time-varying  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  Wrong (white noise)  Time-varying  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  Performance evaluation Performance was evaluated for each scenario by (a) comparing spawning stock biomass (SSB), recruitment, and parameters from the OM with corresponding estimates from the EM, and (b) conducting forecasts of the OM in which future catches are based on management reference points estimated by the EM and summarizing the results of the projections in terms of “lost yield” and differences between ideal and anticipated stock status when management is based on the EM. Each of these components is discussed below. The estimation performance of the EMs was summarized using relative errors by year for SSB and recruitment, using RE =  (θ^-θ)/θ, and median absolute relative errors (MARE = median( |(θ^-θ)/θ|)) for other parameter estimates (e.g. growth, R0) and derived quantities (e.g. maximum sustainable yield MSY, virgin stock biomass B0, final year biomass B2008). These median values were used as a measure of bias, as occasional outliers greatly skewed the mean values. MAREs have been used previously as a measure of precision for point parameter estimates (e.g. Johnson et al., 2015). The management implications of growth variability and how it is handled in the assessment were evaluated by assuming that management aims to follow a constant catch strategy, based on MSY, and implemented without error. This harvest strategy was selected for illustrative purposes. MSY is calculated from the base growth, the stock-recruit, and the selectivity parameters (Methot and Wetzel, 2013), as found in the first year of the simulation, and unaffected by any time-varying biology. The simulated fishery system is projected forward using this constant-catch strategy, and the consequences evaluated using selected metrics. Projecting with the MSY obtained from the OM represents the “true” state of the resource for that replicate, with all relevant hypotheses associated with that OM. Recruitment in the OM during the forecast period is simulated given recruitment variation, while management performance during the forecast period is calculated by simulating future population sizes and catches given either that future growth is constant (for OM 1) or time-varying (for OMs 2–5). The EM was also used to estimate MSY, which was then used as the basis for projections for each year of a 50-year forecast period. These “estimated” catches were provided to the OM and a second set of projections undertaken in which the catches taken from the population in the OM were set to the “estimated” catches. Thus, the OM was projected forward twice, once using the MSY from the OM (“true yield”) and once using the estimates of MSY from EM (“projected yield”). The risk for each EM given the true state of nature was calculated in terms of relative lost yield, the forecasted final year stock status, and the proportion of replicates overfished (see below). With the exception of the final metric (proportion overfished), all the metrics were calculated for each simulation replicate, and the median value and interquartile range across replicates (representing central tendency and precision, respectively) reported. Relative Lost Yield: Difference between estimated retained catch (yield from the EM) and true retained catch (yield from the OM) for the first projected year as a proportion of true catch [(projected yield—true yield)/true yield] Final Stock Status: Ratio of projected (using MSY from EM) final spawning stock biomass to true (using MSY from OM) final stock biomass (projected B2058/true B2058) Proportion Overfished: Proportion of replicates where the final biomass is less than 0.25 of B0, which is the metric by which the stock would be declared overfished {Gertseva et al., 2009; number of [(projected B2058/true B0) < 0.25]/number of runs}. The final year is used in this metric as it would be near-impossible for a fishery to recover from an overfished state under a constant-catch harvest strategy without any intervention. Results The simplified base model (Supplementary Figure S1c) had catch start in 1902, peaking in 1998 (Supplementary Figure S1a). Overall SSB increased when the environmental index was added, and the peaks and troughs in the time series appeared more pronounced (Supplementary Figure S1b). Several replicates (10/1100) failed to converge (non-positive definite Hessian matrix), and were discarded. If a replicate failed within a single scenario, the results for that replicate for all scenarios were discarded. In a true assessment, model configurations could be manually adjusted to ensure convergence, but such methods would not be feasible in a simulation study so we instead discarded these few replicates. Estimation performance Time-invariant growth in OM The two EMs—one correctly specified (1-0-A), one over-parameterized (1-0-B)—performed well when growth was constant in the OM (Figure 2). There was a slight overall positive bias (measured in median RE) for SSB and recruitment for both EMs, although the median RE for each time series was less than 7%. Over-parameterization did not result in a substantial difference in the bias of SSB and recruitment (Figure 2(ii)) because the assessment model generally estimated link parameters (β) close to zero. Base growth parameters were also estimated without bias and precisely, and were similar for the two EMs (Table 4). Table 4. MAREs for the model parameters and derived quantities.     The colours go from white to grey by order of magnitude, with 0 being white, and larger absolute values being grey. The growth parameters in these table are the base parameters of the model, i.e. before the index effects took place. Asterisks (*) denote median absolute (not relative) error as the true values for those parameters are zero. Figure 2. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the OM with constant growth. Each column represents an EM, where (i) shows the EM with constant growth, while (ii) includes a time-varying growth index. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 2. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the OM with constant growth. Each column represents an EM, where (i) shows the EM with constant growth, while (ii) includes a time-varying growth index. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Correctly specifying EM and the environmental index Most of the correctly specified scenarios (2-1-B, 3-1-B, 4-1-B, 5-1-B; Figure 3) were able to estimate SSB and recruitment accurately and precisely when growth was time-varying in the OM, with a similar bias to those from the scenarios with constant growth and either correct or over-parameterized models (Figure 2 and Table 4). The only exception was when β was set to a higher value and growth was highly varying (OM 4; Figure 3(iii)), which greatly decreased precision, even if the median SSB and recruitment estimates across simulations were fairly unbiased (median relative error <13% for all years). Figure 3. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the correctly specified EM (EM 2) given the correct index driving time-varying growth. Each column represents an OM with time-varying growth, where (i–iv) correspond to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 3. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the correctly specified EM (EM 2) given the correct index driving time-varying growth. Each column represents an OM with time-varying growth, where (i–iv) correspond to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Model misspecification (not allowing for time-varying growth) Not allowing for time-varying growth resulted in imprecise and biased estimates of SSB (Figure 4a), although recruitment was estimated precisely and without bias (Figure 4b). Growth parameters were biased, with K being negatively biased, and L∞ positively biased. The magnitude of bias in SSB and recruitment also increased with increasing β values (compare Figure 4(i) with 4(iii)). In this scenario, there was a large bias in the estimates of SSB and recruitment when the EM was misspecified (Figure 4(iii)), although the level of precision appears to be similar irrespective of whether the model was correctly specified or not (Figure 3(iii)). Misspecified EMs under OMs with high β values were the only scenarios where estimates of recruitment were biased (Figure 4b(iii)). Precision of SSB and recruitment also decreased in scenarios when autocorrelation within the environmental index increased and growth was assumed constant in the EM (Figure 4c) and Supplementary Figure S2). Figure 4. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the misspecified EM (EM 1) with time-invariant growth. Each column represents an OM with time-varying growth, where (i–iv) corresponding to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 4. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment for the misspecified EM (EM 1) with time-invariant growth. Each column represents an OM with time-varying growth, where (i–iv) corresponding to OMs 2–5. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Misspecifying the time-varying index We examined two ways in which the environmental index used in the EM differed from the true index driving the population dynamics. When the EM only had 20 out of the 108 years of the environmental data (observation model 2; scenarios 2-2-B and 5-2-B), the EM assumed that growth was constant up until the time series began in scenarios. Base growth parameters (K and L∞) tended to be under-estimated (Supplementary Table S2: rows 5 and 14). These errors increased when growth parameters varied more over time (comparing OM 2 to OM 5; Figure 5.1 and 5.2). SSB was generally biased over time in these cases, with both positive and negative bias during portions of the time series (Figure 5a). When growth was weakly climate-driven and only a shortened index was observed, median RE was positive at the start of the time series and became negative in the later years, and only tended to zero during the last few years of the time series (Figure 5.1a(iii)). Precision was similar to that of the correctly specified model for both SSB and recruitment (Figure 5.1(ii)). The precision of the SSB and recruitment estimates was much poorer for the scenario with strongly climate-driven growth and the shortened index data (Figure 5.2(iii)) than when growth was weakly climate-driven (Figure 5.1(iii)). SSB was over-estimated (up to 50% for some years) for this OM, although the median REs became closer to zero when the environmental data started. Figure 5. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the effects of varying the observation model. Each column represents a scenario, and the labels correspond to the scenario. Panel 1 shows results from the weakly climate-driven OM (OM 2), while panel 2 shows results from the strongly climate-driven OM (OM 5). Labels (i–v) correspond to observations of environmental index 0–4, with their respective EM. E.g. 1.i corresponds to scenario 2-0-A, 1.ii to 2-1-B, 1.iii to 2-2-B, 2.i to 5-0-B, etc. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 5. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the effects of varying the observation model. Each column represents a scenario, and the labels correspond to the scenario. Panel 1 shows results from the weakly climate-driven OM (OM 2), while panel 2 shows results from the strongly climate-driven OM (OM 5). Labels (i–v) correspond to observations of environmental index 0–4, with their respective EM. E.g. 1.i corresponds to scenario 2-0-A, 1.ii to 2-1-B, 1.iii to 2-2-B, 2.i to 5-0-B, etc. The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Completely misspecifying the environmental index, either as an equally auto-correlated index or as random white noise (observation models 3 and 4), did not affect the estimation of the base growth parameters, SSB, and recruitment estimations when growth was weakly climate-driven (OM 2; Figure 5.1(iv) and 5.1 (v)). However, precision and accuracy of these derived quantity estimates greatly decreased when growth was strongly climate-driven (OM 5; Figure 5.2(iv) and 5.2(v)). The point estimates of the base growth parameters in those scenarios were biased by >10% under when the index was wrong and growth was misspecified (Table 4). In the sensitivity analysis using partially correct indices (Equation 4), bias in SSB and recruitment estimates remained relatively small when the index was correctly specified compared to being correctly specified up to 80% (a = 0.8), although precision decreased as the value of a decreased (Figure 6(ii-iii)). Thereafter, as a decreased, estimates of SSB became increasingly biased (again, both positively and negatively) and less precise than the model with a correct index (Figure 6(iv-vi)), although bias in estimates of recruitment remained relatively small until the index was completely misspecified (a = 0; Figure 6(vii)), where the estimates of SSB became highly negatively biased in the later years. Figure 6. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the results of the sensitivity analysis to the accuracy of the index for the OM with strongly climate-driven growth. Each column represents a scenario, and the labels correspond to the scenario. (i) corresponds to scenario 5-0-A, and (ii–vii) to scenarios where the a value (Equation 4) increases from 1 to 0, in 0.2 increments ((ii) also corresponds to scenario 5-1-B). The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Figure 6. View largeDownload slide Time series of relative errors in estimates of (a) spawning stock biomass and (b) recruitment, showing the results of the sensitivity analysis to the accuracy of the index for the OM with strongly climate-driven growth. Each column represents a scenario, and the labels correspond to the scenario. (i) corresponds to scenario 5-0-A, and (ii–vii) to scenarios where the a value (Equation 4) increases from 1 to 0, in 0.2 increments ((ii) also corresponds to scenario 5-1-B). The black line is the median relative error, the light grey area the 50% simulation intervals, the dark grey area 95% simulation intervals. Sensitivity to high β values The magnitude of bias (in terms of median relative error) in the estimates of SSB and recruitment increased with the value for β when growth was assumed constant in the EM (Supplementary Figure S3). The pattern of bias (positive in the early years, and decreasing to become negative in the later years) remained the same across scenarios. Supplementary Figure S2 also shows that precision (measured using the inter-quantile range) of SSB and recruitment decreased as β was increased, regardless of whether allowance was made for time-varying growth in the EM. The problems in estimation found for high β values were probably caused by issues in the process of generating data. The high variance in growth caused fish to attain large sizes at young ages, and remain at that size during years where the index was at lower values, resulting in higher average sizes-at-age, as can be seen in Supplementary Table S3. Management implications Table 5 summarizes the implications of managing under a constant catch strategy, based on an estimate of MSY from monitoring data collected from the fishery and perhaps based on a misspecified assessment model. Median lost yield varied across scenarios, while median final stock status was similar, although inter-simulation variability differed. The proportion overfished differed across OMs, but was similar between EMs. Table 5. Table summary of all scenarios tested. Scenario label  Scenario name  Performance metric   Relative yield lost (interquartile range)  Est B2058/True B2058 (interquartile range)  Proportion of replicates overfished  1-0-A  Constant growth with correct EM  0.0249 (0.0682)  0.955 (0.11)  0.45  1-0-B  Constant growth with over-parameterized EM  0.0254 (0.0687)  0.955 (0.109)  0.45  2-0-A  Weakly climate-driven growth with incorrect EM  0.0666 (0.0859)  0.933 (0.106)  0.31  2-1-B  Weakly climate-driven growth with correct EM and index  0.0436 (0.0634)  0.959 (0.0837)  0.28  2-2-B  Weakly climate-driven growth with correct EM and shortened index  0.0918 (0.0818)  0.912 (0.126)  0.34  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  0.0633 (0.0919)  0.942 (0.121)  0.32  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  0.0628 (0.0922)  0.93 (0.094)  0.33  3-0-A  Highly autocorrelated growth with incorrect EM  0.0251 (0.192)  0.959 (0.293)  0.43  3-1-B  Highly autocorrelated growth with correct EM and index  0.0329 (0.0574)  0.952 (0.098)  0.47  4-0-A  Highly varied growth with incorrect EM  0.134 (0.289)  0.95 (0.171)  0.32  4-1-B  Highly varied growth with correct EM and index  0.131 (0.278)  0.941 (0.216)  0.33  5-0-A  Strongly climate-driven growth with incorrect EM  −0.0308 (0.758)  1.00 (0.555)  0.42  5-1-B  Strongly climate-driven growth with correct EM and index  0.0217 (0.1)  0.985 (0.105)  0.45  5-2-B  Strongly climate-driven growth with correct EM and shortened index  0.237 (0.577)  0.909 (0.681)  0.52  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  0.206 (1.57)  1.00 (0.844)  0.50  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  −0.097 (1.05)  1.01 (0.631)  0.47  Scenario label  Scenario name  Performance metric   Relative yield lost (interquartile range)  Est B2058/True B2058 (interquartile range)  Proportion of replicates overfished  1-0-A  Constant growth with correct EM  0.0249 (0.0682)  0.955 (0.11)  0.45  1-0-B  Constant growth with over-parameterized EM  0.0254 (0.0687)  0.955 (0.109)  0.45  2-0-A  Weakly climate-driven growth with incorrect EM  0.0666 (0.0859)  0.933 (0.106)  0.31  2-1-B  Weakly climate-driven growth with correct EM and index  0.0436 (0.0634)  0.959 (0.0837)  0.28  2-2-B  Weakly climate-driven growth with correct EM and shortened index  0.0918 (0.0818)  0.912 (0.126)  0.34  2-3-B  Weakly climate-driven growth with correct EM and wrong AR index  0.0633 (0.0919)  0.942 (0.121)  0.32  2-4-B  Weakly climate-driven growth with correct EM and wrong random index  0.0628 (0.0922)  0.93 (0.094)  0.33  3-0-A  Highly autocorrelated growth with incorrect EM  0.0251 (0.192)  0.959 (0.293)  0.43  3-1-B  Highly autocorrelated growth with correct EM and index  0.0329 (0.0574)  0.952 (0.098)  0.47  4-0-A  Highly varied growth with incorrect EM  0.134 (0.289)  0.95 (0.171)  0.32  4-1-B  Highly varied growth with correct EM and index  0.131 (0.278)  0.941 (0.216)  0.33  5-0-A  Strongly climate-driven growth with incorrect EM  −0.0308 (0.758)  1.00 (0.555)  0.42  5-1-B  Strongly climate-driven growth with correct EM and index  0.0217 (0.1)  0.985 (0.105)  0.45  5-2-B  Strongly climate-driven growth with correct EM and shortened index  0.237 (0.577)  0.909 (0.681)  0.52  5-3-B  Strongly climate-driven growth with correct EM and wrong AR index  0.206 (1.57)  1.00 (0.844)  0.50  5-4-B  Strongly climate-driven growth with correct EM and wrong random index  −0.097 (1.05)  1.01 (0.631)  0.47  The proportions of replicates that resulted in overfished populations ranged from 0.28 to 0.51 across the scenarios (Table 5). Some of the randomly generated recruitment deviations resulted in populations with extreme spawning stock biomass levels, both high and low, with low levels causing populations to crash. Projecting a constant catch strategy for populations with low biomass may have also contributed to the relatively high proportion of replicates being considered overfished, as shown by the 0.45 value from scenario 1-0-A, when the assessment was correctly specified and growth was time-invariant. Most EMs tended to provide MSY estimates that were median-unbiased, given the correct index data, regardless of whether time-varying growth was accounted for. However, inter-simulation variability (measured in terms of inter-quantile range), increased greatly when the EM was misspecified, i.e. growth was assumed to be constant (2-0-A, 3-0-A, 5-0-A). The exception to both these results was the OM 4, where β was set to a high level. In that case, there were negligible differences in MSY estimates between the misspecified EM with time-invariant growth and the correctly specified EM (Table 5, scenarios 4-O-A and 4-O-B). The estimates of MSY were positively biased when the EM was based on an incorrect environmental index (scenarios 2-3-B, 2-4-B, 5-3-B, 5-4-B). This resulted in up to a median of 24% of yield being lost in the OM where growth is strongly climate-driven (OM 5; Table 5). Inter-simulation variability was very high for the scenarios based on OM 5, while this was not the case for the scenarios where growth was weakly climate-driven (OM 2). Supplementary Table S4 shows an increase in inter-simulation variability of management quantities as the accuracy of the index decreased, although bias remained the same across scenarios. Model misspecification (i.e. assuming time-invariant growth when growth was actually varying over time) resulted in more biased and less precise estimates of MSY (Table 4), resulting in a larger median and wider range of lost yield across replicates, respectively (Table 5). Table 5 also shows that the inter-quantile range of lost yield (i.e. the imprecision of MSY estimates) decreased when the EM was correctly specified, particularly when growth varied over time (e.g. comparing 3-0-A vs. 3-0-B). Discussion Incorporating accurate growth indices resulted in higher precision and less bias in parameter estimates and management outcomes. These effects are more pronounced if the environmental index strongly influenced growth, was highly auto-correlated, and/or had high variance, and was accurately described. The degree of bias and inaccuracies are largely defined by the variance of the time-varying index, and the strength of its auto-correlation. When growth was either weakly time-varying or constant, including an inaccurate index of low variance (random noise) carried similar risks of obtaining inaccurate or imprecise abundance estimates to estimating time-invariant growth and not including an index at all. On the other hand, when growth varied over time with high variance, including an inaccurate index (a < 0.6) or not including one at all resulted in highly biased and inaccurate estimates of stock size. Estimating time-varying growth when the growth index is correct results in the most precise and least biased estimates. Allowing for time-varying growth in the form of an index does not risk significant degradation in estimation performance even when growth is not actually time-varying, although this scenario is unlikely (e.g. Stawitz et al., 2015; Thorson and Minte-Vera, 2016). Therefore, if the growth index under consideration is weakly auto-correlated and weakly linked to the growth parameters, inclusion of the index would likely improve the performance of the model, regardless of how accurate the index is. Based on the results for OM 2, the most “realistic” OM, and the properties of the growth index developed in Black (2009), we would recommend consideration of including this index in the assessment of splitnose rockfish. If a climate-driven growth index were to be considered for use in an assessment, we would recommend fitting an AR-1 model to the index, using these parameters to generate new simulated growth indices, and simulation-testing their inclusion in the stock assessment model. The results were obtained from only one set of life history parameters, but the simulation-based risk analysis approach could be conditioned on the data and biology of other stocks to illustrate whether to include a growth index for a different specific stock assessment. A growth index for which 60% of the variance is correct (0.77 correlation with true variation in growth; √(a=0.60)) could still lead to estimates of parameters and derived quantities that have similar bias and precision as not including the index at all, when growth is strongly time-varying (OM 5). This case study suggests that, when growth is suspected to be strongly climate-driven, environmental indices should only be used when there is high confidence in their accuracy. For example, Thorson and Minte-Vera (2016) found that the best models with a single time-varying effect could explain on average 50% of variability in weight-at-age for exploited stocks, with year-effects (i.e. annual indices like those explored here) explaining approximately 45% of variability. Stawitz et al. (2015) described a state-space model that was fit to fishery-dependent and -independent data for 29 stocks of Pacific groundfish, and found that approximately 40% of the modelled stocks exhibited time-varying growth. Whereas it is highly unlikely that it is possible to observe the true underlying mechanism by which growth varies over time and generate an exact environmental index, assuming time-invariant growth (status quo) could carry similar risks of biased and imprecise assessment outputs. A feasible method of increasing the accuracy of estimated growth variation would be by examining similarities in growth variation among species (Stawitz et al., 2015) and/or cross-validating particular years of strong or weak growth with other species (Black, 2009). Additionally, Maunder and Watters (2003) recommend an integrated approach with additional process error when including an environmental time series into stock assessment models, which should be a subject for future research. This approach could presumably be implemented using annual variation in growth, in addition to including the index multiplicatively, and the magnitude of additional variation in growth could be estimated using the Laplace approximation (Thorson et al., 2015a). Improving the accuracy of growth parameter estimates had a relatively small contribution to total predicted MSY. Furthermore, erroneous parameter and management quantity estimates had negligible impact on future stock size in this case study. This could have been an artefact of having both natural mortality M and steepness h parameters fixed in the assessment. These parameters exhibit a strong influence on reference point estimation, as shown in Mangel et al. (2013), and as such would limit the extent by which growth parameters could affect estimates of MSY. Furthermore, the studied species is long-lived and slow-growing, with a relatively high age of maturity, and a short time series of age data. MSY for long-lived species is expected to have similar sensitivity to changing growth parameters as short-lived species (Thorson et al., 2015b). Thorson et al. (2015b) found that the MSY from the three life histories tested showed similar sensitivity to time-varying biological parameters (for recruitment, size, and maturity), while sensitivities differed when natural mortality varied over time, so we do not expect MSY to be more sensitive to growth variation for shorter-lived species. Other management reference points (such as BMSY/B0) may show different sensitivities to time-varying parameters for each life history. Future work should thus consider assessments where there is more data on age and growth, and/or these parameters are not pre-specified, or based on species where growth may vary more strongly over time (e.g. California Current petrale sole Eopsetta jordani or Gulf of Alaska walleye pollock; Stawitz et al., 2015) to more fully determine the effects of incorporating a time-varying index on growth. Larger scale effects could also be tested by using alternate forms of time-varying growth, such as the sine-wave model or AR models with more lags. In addition, growth also varies at a spatial (Rahikainen and Stephenson, 2004; Thorson, 2015) and at an individual (Webber and Thorson, 2016) scale, which could increase model complexity to explain more variance in the data. Implementing constant catch exactly (i.e. without error) is unlikely to be the best strategy to use for examining the management impacts of including a time-varying growth index in assessments. The risk analysis approach can be extended to be based on alternative management strategies to perform a full management strategy evaluation (Punt et al., 2016). Conducting an MSE would allow managers to better assess the consequences of including a time-varying growth index in a management strategy that applies harvest control rules (Punt et al., 2016). Punt et al. (2014) summarize studies where this has been done for recruitment. There has been a recent push to move from single-species towards ecosystem-based fisheries management (EBFM; Pikitch et al., 2004), where environmental factors are taken into consideration to manage the ecosystem as a whole. Therefore, there is a desire to explore various methods for incorporating time-varying growth into assessments. The simulation approach described here can be used as a step towards addressing and incorporating climate effects into somatic growth for a single species. Supplementary data Supplementary material is available at the ICESJMS online version of the manuscript. Acknowledgements This work was funded by a NOAA Fisheries and the Environment (FATE) grant, project 13-01. The authors thank Timothy Essington, Melissa Haltuch, James Hastie, Michelle McClure, Felipe Hurtado Ferro, Ernesto Jardim, Kelli Johnson, Ian Taylor, Chantel Wetzel, and an anonymous reviewer for their helpful suggestions and comments throughout the working process and on the manuscript. Thanks also go to Bryan Black for sharing his otolith index data. References Ait Youcef W., Lambert Y., Audet C. 2015. Variations in length and growth of Greenland Halibut juveniles in relation to environmental conditions. Fisheries Research , 167: 38– 47. Google Scholar CrossRef Search ADS   A’mar Z. T., Punt A. E., Dorn M. W. 2009. The evaluation of two management strategies for the Gulf of Alaska walleye pollock fishery under climate change. ICES Journal of Marine Science , 66: 1614– 1632. Google Scholar CrossRef Search ADS   Black B. A. 2009. Climate-driven synchrony across tree, bivalve, and rockfish growth-increment chronologies of the northeast Pacific. The Marine Ecology Progress Series , 378: 37– 46. Google Scholar CrossRef Search ADS   Black B. A., Boehlert G. W., Yoklavich M. M. 2008. Establishing climate–growth relationships for yelloweye rockfish (Sebastes ruberrimus) in the northeast Pacific using a dendrochronological approach. Fisheries Oceanography , 17: 368– 379. Google Scholar CrossRef Search ADS   Brunel T., Dickey-Collas M. 2010. Effects of temperature and population density on von Bertalanffy growth parameters in Atlantic herring: a macro-ecological analysis. The Marine Ecology Progress Series , 405: 15– 28. Google Scholar CrossRef Search ADS   Essington T. E., Kitchell J. F., Walters C. J. 2001. The von Bertalanffy growth function, bioenergetics, and the consumption rates of fish. Canadian Journal of Fisheries and Aquatic Sciences , 58: 2129– 2138. Google Scholar CrossRef Search ADS   Fournier D., Archibald C. P. 1982. A General Theory for Analyzing Catch at Age Data. Canadian Journal of Fisheries and Aquatic Sciences , 39: 1195– 1207. Google Scholar CrossRef Search ADS   Francis R. I. C. C. 2016. Growth in age-structured stock assessment models. Fisheries. Research , 180: 77– 86. Google Scholar CrossRef Search ADS   Gertseva V. V., Cope J. M., Pearson D. E. 2009. Status of the U.S. splitnose rockfish (Sebastes diploproa) resource in 2009, Status of the Pacific Coast Groundfish Fishery through 2009 and recommended acceptable biological catches for 2011/2012: stock assessment and fishery evaluation, 2009. Pacific Fishery Management Council, Portland, OR. Haltuch M. A., Punt A. E. 2011. The promises and pitfalls of including decadal-scale climate forcing of recruitment in groundfish stock assessment. Canadian Journal of Fisheries and Aquatic Sciences , 68: 912– 926. Google Scholar CrossRef Search ADS   Harley C. D. G., Randall Hughes A., Hultgren K. M., Miner B. G., Sorte C. J. B., Thornber C. S., Rodriguez L. F., Tomanek L., Williams S. L. 2006. The impacts of climate change in coastal marine systems. Ecology Letters , 9: 228– 241. Google Scholar CrossRef Search ADS PubMed  Hulson P.-J. F., Quinn T. J., Hanselman D. H., Ianelli J. N. 2013. Spatial modeling of Bering Sea walleye pollock with integrated age-structured assessment models in a changing environment. Canadian Journal of Fisheries and Aquatic Sciences , 70: 1402– 1416. Google Scholar CrossRef Search ADS   Johnson K. F., Monnahan C. C., McGilliard C. R., Vert-pre K. A., Anderson S. C., Cunningham C. J., Hurtado-Ferro F., Licandeo R. R., Muradian M. L., Ono K., et al.   2015. Time-varying natural mortality in fisheries stock assessment models: identifying a default approach. ICES Journal of Marine Science , 72: 137–150. Kimura D. K. 2008. Extending the von Bertalanffy growth model using explanatory variables. Canadian Journal of Fisheries and Aquatic Sciences , 65: 1879– 1891. Google Scholar CrossRef Search ADS   Lorenzen K. 2016. Toward a new paradigm for growth modeling in fisheries stock assessments: embracing plasticity and its consequences. Fisheries Research , 180: 4– 22. Google Scholar CrossRef Search ADS   Mangel M., MacCall A. D., Brodziak J., Dick E. J., Forrest R. E., Pourzand R., Ralston S. 2013. A perspective on steepness, reference points, and stock assessment. Canadian Journal of Fisheries and Aquatic Sciences , 70: 930– 940. Google Scholar CrossRef Search ADS   Maunder M. N., Punt A. E. 2013. A review of integrated analysis in fisheries stock assessment. Fisheries Research , 142: 61– 74. Google Scholar CrossRef Search ADS   Maunder M. N., Watters G. M. 2003. A general framework for integrating environmental time series into stock assessment models: model description, simulation testing, and example. Fisheries Bulletin , 101: 89– 99. Methot R. D., Taylor I. G. 2011. Adjusting for bias due to variability of estimated recruitments in fishery assessment models. Canadian Journal of Fisheries and Aquatic Sciences , 68: 1744– 1760. Google Scholar CrossRef Search ADS   Methot R. D., Wetzel C. R. 2013. Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research , 142: 86– 99. Google Scholar CrossRef Search ADS   Mueter F. J., Bond N. A., Ianelli J. N., Hollowed A. B. 2011. Expected declines in recruitment of walleye pollock (Theragra chalcogramma) in the eastern Bering Sea under future climate change. Journal of Marine Science , 68: 1284– 1296. Pikitch E. K., Santora C., Babcock E. A., Bakun A., Bonfil R., Conover D. O., Dayton P., Doukakis P., Fluharty D., Heneman B., et al.   2004. Ecosystem-based fishery management. Science , 305: 346– 347. Google Scholar CrossRef Search ADS PubMed  Pinsky M. L., Worm B., Fogarty M. J., Sarmiento J. L., Levin S. A. 2013. Marine taxa track local climate velocities. Science , 341: 1239– 1242. Google Scholar CrossRef Search ADS PubMed  Punt A. E., A’mar T., Bond N. A., Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haltuch M. A., Hollowed A. B., Szuwalski C. 2014. Fisheries management under climate and environmental uncertainty: control rules and performance simulation. ICES Journal of Marine Science , 71: 2208– 2220. Google Scholar CrossRef Search ADS   Punt A. E., Butterworth D. S., de Moor C. L., De Oliveira J. A. A., Haddon M. 2016. Management strategy evaluation: best practices. Fish and Fisheries , 17: 303– 334. Google Scholar CrossRef Search ADS   Rahikainen M., Stephenson R. L. 2004. Consequences of growth variation in northern Baltic herring for assessment and management. ICES Journal of Marine Science , 61: 338– 350. Google Scholar CrossRef Search ADS   Rowe R. J., Terry R. C. 2014. Small mammal responses to environmental change: integrating past and present dynamics. Journal of Mammalogy , 95: 1157– 1174. Google Scholar CrossRef Search ADS   Schirripa M. J., Goodyear C. P., Methot R. M. 2009. Testing different methods of incorporating climate data into the assessment of US West Coast sablefish. Journal of Marine Science , 66: 1605– 1613. Stawitz C. C., Essington T. E., Branch T. A., Haltuch M. A., Hollowed A. B., Spencer P. D. 2015. A state-space approach for detecting growth variation and application to North Pacific groundfish. Canadian Journal of Fisheries and Aquatic Sciences , 72: 1316– 1328. Google Scholar CrossRef Search ADS   Sünksen K., Stenberg C., Grønkjær P. 2010. Temperature effects on growth of juvenile Greenland halibut (Reinhardtius hippoglossoides Walbaum) in West Greenland waters. Journal of Sea Research , 64: 125– 132. Google Scholar CrossRef Search ADS   Taylor I. G., Methot R. D. 2013. Hiding or dead? A computationally efficient model of selective fisheries mortality. Fisheries Research , 142: 75– 85. Google Scholar CrossRef Search ADS   Thorson J. T. 2015. Spatio-temporal variation in fish condition is not consistently explained by density, temperature, or season for California Current groundfishes. Marine Ecology Progress Series , 526: 101– 112. Google Scholar CrossRef Search ADS   Thorson J. T., Hicks A. C., Methot R. D. 2015a. Random effect estimation of time-varying factors in Stock Synthesis. Journal of Marine Science , 72: 178– 185. Thorson J. T., Jensen O. P., Zipkin E. F. 2014. How variable is recruitment for exploited marine fishes? A hierarchical model for testing life history theory. Canadian Journal of Fisheries and Aquatic Sciences , 71: 973– 983. Google Scholar CrossRef Search ADS   Thorson J. T., Minte-Vera C. V. 2016. Relative magnitude of cohort, age, and year effects on size at age of exploited marine fishes. Fisheries Research , 180: 45– 53. Google Scholar CrossRef Search ADS   Thorson J. T., Monnahan C. C., Cope J. M. 2015b. The potential impact of time-variation in vital rates on fisheries management targets for marine fishes. Fisheries Research , 169: 8– 17. Google Scholar CrossRef Search ADS   Thorson J. T., Stewart I. J., Taylor I. G., Punt A. E. 2013. Using a recruitment-linked multispecies stock assessment model to estimate common trends in recruitment for US West Coast groundfishes. The Marine Ecology Progress Series , 483: 245– 256. Google Scholar CrossRef Search ADS   Von Bertalanffy L. 1957. Quantitative laws in metabolism and growth. The Quarterly Review of Biology, 32: 217– 231. Weatherley A. H. 1990. Approaches to understanding fish growth. Transactions of the American Fisheries Society , 119: 662– 672. Google Scholar CrossRef Search ADS   Webber D. N., Thorson J. T. 2016. Variation in growth among individuals and over time: a case study and simulation experiment involving tagged Antarctic toothfish. Fisheries Research , 180: 67– 76. Google Scholar CrossRef Search ADS   © International Council for the Exploration of the Sea 2017. All rights reserved. For Permissions, please email: journals.permissions@oup.com

Journal

ICES Journal of Marine ScienceOxford University Press

Published: Jan 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off