Introduction

Fisheries management is often complicated by pervasive uncertainty and competing objectives among stakeholders; however, these difficulties are particularly pronounced in Pacific salmon fisheries. Most salmon species are anadromous and highly migratory. Their complex life histories can obscure mechanisms responsible for driving interannual variability in recruitment [@Peterman2010], as well as regime shifts in population productivity [@Peterman2003]. Furthermore, although salmon populations exhibit substantial genetic diversity, they are generally harvested in multi-stock fisheries and managers have limited ability to maintain harvest rates while avoiding depleted stocks [@Walters2008]. Nevertheless, persistent declines in the abundance of many British Columbia salmon stocks have resulted in pressure to identify management strategies that can promote rebuilding and preserve biodiversity as outlined in Canada's Wild Salmon Policy [@DFO2005].

Closed-loop simulations can be used to identify management strategies that are robust to multiple dimensions of uncertainty including recruitment variation, changes in productivity regime, and deviations between target and realized catches (i.e. outcome uncertainty). Closed-loop simulations are also particularly well suited to evaluating how strategies balance conservation and socio-economic objectives. To this end closed-loop simulations can be combined with stakeholder engagement to formalize objectives and create a management strategy evaluation [@Punt2016].

Closed-loop simulations are composed of two central components. The operating model represents the ecological dynamics of the relevant stocks, as well as the fishery that harvests them. The management procedure represents the assessment process, as well as the harvest control rule that determines exploitation rates. Broadly, the models use a stock-recruit relationship to generate a cohort of recruits based on spawner abundance at the beginning of the simulation. Abundance is observed with error within the model and used to evaluate the status of the stock, which in turn determines target catch rates based on the harvest control rule. Some fraction of recruits are removed, representing harvest within the fishery, the remainder become spawners, and the entire process iterates through time.

Closed-loop simulation analyses typically involve comparing scenarios (unique combinations of operating models and management procedures) using various performance metrics. Performance metrics are selected to address the management objectives of greatest relevance to the stock aggregate and typically are composed of both conservation- and catch-based goals, at various ecological scales and time frames. For example, median catch in the first five years of the simulation may be used to evaluate short term economic objectives, while the proportion of stocks above a biological benchmark at the end of the simulation may be a proxy for conservation performance. Each scenario is simulated hundreds to thousands of times to fully incorporate stochastic processes (e.g. recruitment deviations, implementation uncertainty) into estimates of performance.

samSim - A Generic Simulation Tool

samSim is an R package containing functions to quantify the rebuilding potential for Pacific salmon populations. recoverySim is the function containing the closed-loop simulation model and is the primary function within the package. The model structure focuses on aggregates of salmon conservation units (CUs), the focal unit for assessing Canadian Pacific salmon status [@DFO2005]. Generally all CUs included within a single analysis should be managed simultaneously as a management unit (MU); however, the model can be adapted to include multiple MUs with distinct harvest rates. Although initial case studies focused on Fraser River sockeye salmon and Area 3 chum salmon, samSim is intended to be applicable to any Pacific salmon species as long as two requirements are met. First, CU-specific stock-recruit parameters and age-at-maturity values must be available to parameterize the operating model. Ideally these would be derived from observed, CU-specific time series of age-structured spawner and recruit abundance, however they can also be generated using alternative techniques (e.g. habitat-based estimates; expert opinion). Second, harvest of immature fish must be minimal because offshore fisheries are not accounted for in the model.

recoverySim uses observed time series of spawner abundance to prime the simulation so that each CU's initial status reflects the most recent assessment. Those abundances, along with externally estimated stock-recruit parameters, are used as inputs to a Ricker model (Larkin models can also be used for cyclic CUs), which generates a cohort of recruits (i.e. the total number of adult offspring produced by a given brood year of spawners). This process is stochastic, generating representative interannual variation in recruitment deviations, and can incorporate covariance among CUs, as well as temporal autocorrelation. Age-at-maturity varies among simulated recruits based on input parameters, with stochastic variability introduced via multivariate logistic error. This process creates cohorts of returns, the total number of adult offspring returning to spawn in a given year, which consist of mature fish of various ages and brood years.

Simulated returns are next subject to harvest in American fisheries, if relevant, which typically occurs during return migrations before Canadian CUs enter nearshore areas. American catches are always generated using a fixed exploitation rate passed as an input value. Canadian mixed-CU fishery harvest occurs next and catch rates are determined by one of three harvest control rules (HCRs): fixed exploitation rate, generic abundance-based, or total allowable mortality (TAM). The generic abundance-based HCR increases exploitation rates from a minimum value when return abundance exceeds user-specific reference points. The TAM rule is similar to the generic abundance-based HCR, but has additional modifications currently implemented in Fraser River sockeye salmon management [@Pestal2011]. In both the American and Canadian mixed-CU fisheries, target catches (generated using the HCR) are converted to realized catches by incorporating stochastic outcome uncertainty [@Holt2006]. CU-specific realized catch is calculated as a function of relative abundance. Once catch has been removed the remaining returns become spawners, creating the subsequent generation of recruits. Modifications to this general framwork are described below (see Additional Components section for details).

samSim Outputs

library(samSim)
library(tidyverse)

simPar <- read.csv(here::here("data", "manProcScenarios",
                        "fraserMPInputs_techReport.csv"),
                   stringsAsFactors = F)
scenNames <- unique(simPar$scenario)
dirNames <- sapply(scenNames, function(x) paste(x, unique(simPar$species),
                                                sep = "_"))
agVarsToPlot <-  c("medRecRY", "medSpawners", "ppnCUUpper", "ppnCUExtinct",
                   "medCatch")

#set global figure options (links to .tex file)
# knitr::knit_hooks$set(plot = function(x, options)  {
#   hook_plot_tex(x, options)
# })

The simulation model produces a range of outputs automatically. A PDF of diagnostics shows simulated stock-recruit figures, as well as time series of abundance, population parameters, and various performance metrics. Automatically generated R data and .csv files contain arrays of CU-specific time series, matrices of aggregate time series, and data frames of aggregate or CU-specific performance metrics. Each file contains all the Monte Carlo trials for a specific scenario.

The output files are used by additional samSim functions to create various summary figures. Each figures incorporates percentile intervals that are calculated across Monte Carlo trials to provide an estimate of uncertainty. The plotContTradeoffs function generates double-y axis line plots that demonstrate how spawner abundance, stock status, and extinction probability change as exploitation rates increase within a single scenario (Fig. 1).

```rChanges in aggregate conservation- and catch-based performance metrics (escapement - purple; catch - blue; proportion of CUs above benchmark - green; extinction risk - red) as fixed exploitation rates increase. Dark lines represent median values and polygons the 90th percentile interval across Monte Carlo trials."} agDat <- buildDataAgg(dirNames, agVars = agVarsToPlot, keyVarName = "expRate") refAgDat <- agDat %>% filter(om == "ref", hcr == "fixedER") plotContTradeOffs(refAgDat, keyVar = "expRate", double = TRUE)

The `plotAgDot` and `plotCUDot` functions provide summaries of multiple performance metrics (aggregate or CU-specific respectively) across different management procedures or different operating models. For example, we can visualize declines in return size and the proportion of CUs above their biological benchmark as exploitation rates increase, as well as how such patterns differ among two different productivity regimes (Fig. 2).

```rPerformance across three metrics as a function of fixed exploitation rate (x-axis) and two productivity operating models (color). Points represent median values and whiskers the 90th percentile interval across Monte Carlo trials."}
prodDat <- agDat %>%
  filter(om %in% c("ref", "lowProd"),
         expRate %in% c("0.1", "0.3", "0.5", "0.7", "0.9"),
         hcr == "fixedER", 
         var %in% c("medRecRY", "ppnCUUpper", "medCatch", "med")) %>% 
  mutate(var = recode(factor(var), "medRecRY" = "Median Return",
                      # "medSpawners" = "Median Escapement",
                      "ppnCUUpper" = "Ppn. CUs Above\nBenchmark",
                      # "ppnCUExtinct" = "Ppn. CUs Extinct",
                      "medCatch" = "Median Catch"),
         om = recode(factor(om), "ref" = "Reference", "lowProd" = "Low"))

plotAgDot(prodDat, group = "om", legendLab = "Productivity",
          xLab = "Exploitation Rate", yLab = NULL, plotTitle = NULL, 
          axisSize = 10, dotSize = 3, lineSize = 1, legendSize = 11)
png(here::here("outputs", "figs", "tech_rep", "ribb_plot.png"), height = 4.5, 
    width = 5.5, units = "in", res = 300)
plotContTradeOffs(refAgDat, keyVar = "expRate", double = TRUE)
dev.off()

png(here::here("outputs", "figs", "tech_rep", "agg_dot.png"), height = 3.5, 
    width = 6, units = "in", res = 300)
plotAgDot(prodDat, group = "om", legendLab = "Productivity",
          xLab = "Exploitation Rate", yLab = NULL, plotTitle = NULL, 
          axisSize = 10, dotSize = 3, lineSize = 1, legendSize = 11)
dev.off()

The plotAgTradeoff and plotCUTradeoff functions plot a conservation metric on the x-axis and a catch metric on the y-axis to visualize how tradeoffs between objectives vary among management procedures or operating models. Such figures can be particularly useful in identifying subsets of management procedures that meet pre-specified objectives (e.g. minimum median return sizes and catches), while readily incorporating different operating models. For example, intermediate fixed exploitation rates or an abundance-based (TAM) harvest control rule may lead to optimal outcomes under reference productivity, but fail to meet threshold objectives when productivity is low (Fig. 3).

```rTrade-offs between aggregate catch and escapement for fixed exploitation rate and total allowable mortality harvest control rules. Points represent median values and whiskers the 90th percentile interval across Monte Carlo trials."} prodDat2 <- agDat %>% filter(om %in% c("ref", "lowProd"), expRate %in% c("0.1", "0.3", "0.5", "0.7", "0.9", "TAM"), var %in% c("medCatch", "medSpawners")) %>% mutate(om = recode(factor(om), "ref" = "Reference Productivty", "lowProd" = "Low Productivity"))

plotAgTradeoff(prodDat2, consVar = "medSpawners", catchVar = "medCatch", facet = "om", shape = "hcr", showUncertainty = TRUE, mainLab = "", legendLab = "Exploitation\nRate", xLab = "Median Catch", yLab = "Median Spawners", axisSize = 12, legendSize = 10, scaleAxis = "fixed") ```

Additional Components

The basic recoverySim function is intended to provide an indication of how changes in exploitation rate will influence the long term status of an aggregate of CUs. However additional components of the simulation model can be adjusted to explore how these patterns change under different ecological assumptions or, alternatively, how more nuanced harvest control rules perform. The principal optional components are described below.

  1. Changes in productivity

A large number of Pacific salmon populations show evidence of declines in productivity (i.e. mean number of recruits produced per spawner; [@Peterman2012; @Dorner2017]), which may limit the potential for rebuilding efforts even if precautionary management strategies are in place. These impacts can be incorporated into the simulation process in a number of ways. First, the input file contains an argument prodRegime that can be used to set productivity at the long term average (default), apply a scalar multiplier, a linear decline over a specified number of years, generate divergent trends among CUs, or cause a single CU to increase/decrease. Second, an alternative dataset of parameters generated externally can be passed to the model. For example, the default dataset may represent a stationary stock-recruit model and the alternative a time-varying.

  1. Alternative forms of recruitment deviation

The simplest interpretation of a stock-recruit model assumes that deviations from expected recruitment are log-normally distributed, are temporally independent, and are independent among CUs. The impact of these assumptions can be explicitly tested in forward simulations. First, the distribution that generates recruitment deviations can be changed to a student-t distribution, which increases the probability of extreme recruitment events, or a skewed distribution, which shifts the median of the distribution away from zero so that deviations are biased high/low. Second, temporal autocorrelation in recruitment deviations can be added via an autoregressive lag-1 process, which causes deviations in a given year to be correlated with those in the previous year. Third, covariance among CUs can be added by specifying a non-zero correlation coefficient. As this value increases, the simulated dynamics of CUs change from independent to synchronized. The impact of changes in the magnitude of recruitment deviations can also be explicitly considered by applying a scalar to the variance parameter in the stock-recruit relationship.

  1. En-route mortality

Pacific salmon populations, particularly those with long-distance freshwater migrations, can exhibit substantial mortality between marine fisheries and the spawning grounds [@Cooke2004]. These losses are commonly referred to as en-route mortality and can severely impact dynamics by driving spawner abundance below management targets, particularly when exploitation rates are relatively high. recoverySim can incorporate mean en-route mortality rates, as well as normally distributed interannual variability via a SD term, as CU-specific input parameters. In the case of Fraser River sockeye salmon en route mortality was parameterized using observed difference between in-river estimates of spawner abundance and estimates of abundance from spawning grounds.

  1. Single-CU fisheries

Given that salmon fisheries often harvest a mix of abundant and depleted CUs, there is growing interest in evaluating how changes in fishery selectivity could increase the probability of rebuilding populations. While the precise structure of single-CU fisheries will vary among regions and species, recoverySim includes an option to incorporate a generic single-CU HCR intended to reduce harvest on CUs at low abundance. Briefly, a TAC for all Canadian fisheries is calculated based on the specified mixed-CU fishery HCRs described above. A set proportion of the Canadian TAC is allocated to mixed- and single-CU fisheries based on user-defined inputs. The mixed-CU fishery occurs as previously described, followed by the single-CU fishery that has its own unique outcome uncertainty. To increase conservation performance, catch targets for single-CU fisheries are set to zero for any CUs with estimates of abundance below a biological benchmark, generally resulting in a decline in realized catches for a given harvest rate as single-CU fishery allocations increase. Although selectivity could theoretically be increased via spatial or temporal closures, the recoverySim framework basically represents a mixed-CU marine fishery and a terminal single-CU fishery. Since such a fishery may be influenced by en-route mortality, the simulation model also allows users to specify whether en-route mortality occurs before or after single-CU harvest.

  1. Observation component

recoverySim includes model components that simulate the observation and assessment process, which can be used to evaluate how changes in observation error influence deviations between "true" and estimated status. In short, simulated spawner and catch abundances are converted into observed abundances by applying log-normally distributed error. A simple run reconstruction within the model uses these data to generate age-specific estimates of recruitment, once again after applying observation error. These values are then used to construct observed spawner-recruit time series for each CU within the model. The stock-recruit data are used to estimate each CU's status relative to stock-recruit or percentile benchmarks [@Holt2009]. Althugh not currently implemented within recoverySim, the observation model could also be directly linked to an HCR resulting in feedback between the observation process and harvest rates.

Literature Cited



CamFreshwater/samSim documentation built on Sept. 25, 2023, 10:22 a.m.