All coding ws done using R [@R] and FLR [@kell2007flr] in Rmarkdown [@Rmdbook], and ran as vignettes. The information required to run the analysis, is available at ???, including the ICES assessment data set, and the inputs and outputs from the ecosystem models.

Code for processing the input data, running the MSE, and summarising the results is contained in the vignettes like this one, and all code can be found in the ‘Rmd’ files. Only Key steps in the analysis are shown, e.g. code used for plotting is not echoed in the output but can be found in the ‘Rmd’ file.

Introduction

The Operating Model is conditioned on the ICES stock assessment [@ICES2023]. A Reference Case Operating Model was conditioned, then to evaluate robustness a range of hypotheses were used to condition a Robustness Set. These are a limited set of scenarios which plausible are likely to to have major impacts on the performance of the management strategies.

To model environmental factors and predation natural mortality was split into predation (M2) and background (M1) mortality, based on hypotheses from Strategic Ecosystem Models. The assessment was then refitted, and and stock recruitment relationship estimated under a variety of hypotheses, including alternative hypotheses about steepness using priors from FishLife [@thorson2017predicting], and depensation.

This vignette is divided into sections for Installing the required libraries, sourcing bespoke functions, and loading the data, before the Operating Condition is conducted, based on hypotheses about Natural Mortality and the Stock Rcruitment relationship.

Details on Funding, Software Versions andReferences can be found at the end of the document.

Installation

Libraries (i.e. packages) can be installed from gihtub using the remotes package:

install.packages("remotes")

library(remotes)

install.packages("ggplot2")
install.packages("plyr")
install.packages("dplyr")
install.packages("reshape")

remotes::install_github("flr/FLCore")
remotes::install_github("flr/ggplotFL")
remotes::install_github("flr/FLBRP")
remotes::install_github("flr/FLasher")
remotes::install_github("flr/FLAssess")
remotes::install_github("flr/FLife")
remotes::install_github("flr/mydas")
remotes::install_github("lauriekell/FLCandy")

remotes::install_github("flr/FLSRTMB")

install.packages("patchwork")

remotes::install_github("james-thorson/FishLife")

remotes::install_github('fishfollower/SAM/stockassessment', ref='components')
remotes::install_github("flr/FLfse/FLfse")

The libraries used are

library(FLCore)
library(ggplotFL) 
library(FLBRP)
library(FLasher)
library(FLAssess)

library(ggplot2)
library(plyr)
library(dplyr)
library(reshape)

library(FLife)
library(FLBRP)
library(FLCandy)

library(FLSRTMB)
library(patchwork)

library(stockassessment)
library(FLfse)

library(FishLife)

theme_set(theme_bw(16))

Data

The data are the ICES inputs and historical estimates [@ICES2023] generated by “1_ICES” which re-runs the ICES assessment

Bespoke functions

Functions are used to run key parts of the analysis, this makes it easier to document the analysis, debug code and replicate it. The ICES assessment is run using the ‘runSam’ function, so that the assessment can be redone in a comparable way when new data become available, or based on alternative hypotheses

source("../R/runSam.R")  
source("../R/stars.R")
source("../R/eqFn.R")
source("../R/mFn.R")

Back to Top

Natural Mortality

In the ICES assessment natural mortality (\(M\)) is assumed to be 0.15 for a ll ages in all year. However, mass-at-age of mackerel has varied over time.

m=read.csv("../data/inputs/ecosystem/Mackerel_M.csv")

Figure1 Predation mortality (M2) estimates for Mackerel from ecosystem models of the Northern Barents Sea (NorBar), West Coast of Scotand (WCofS) and North Sea. The dashed horizontal line at M=0.15 represents the assumed M in the ICES stock assessment for North Atlantic mackerel.

Natural mortality (M) has been shown to be a function of size-at-age [@lorenzen2022size], e.g. \(M=aW^b\)

\(a\) was chosen so that M at age-at-maturity (\(A_{50}\)) is the same as in the ICES assessment (e.g. 0.15) and M1 was assumed to be 0.025, based on ecosystem models, giving M2=M-M1. M-at-age was then re-estimated based on the stock weights-at-age which have varied over time.

ddM<-function(wt,par){ 
  par["m1"]%*%(wt%^%par["m2"])}

M1Fn=function(x,y=0.025) FLQuant(y,dimnames=dimnames(m(x))) 
M2Fn=function(x) m(x) - M1Fn(x)

par=FLPar(m1=c(0.15/(0.300^-0.288)),
          m2= -0.288)

There are no values for the weights-at-age in 0-group in the stock so these were replaced with the catch values.

load("../data/om/icesAssessment.RData")

om=ices 

stock.wt(om)["0"]=catch.wt(om)["0"]

m(om)=ddM(stock.wt(om),par)

Figure2 Time series of mass-at-age.

Figure3 Time series of natural mortality at age.

M-at-age decreases with size.

Figure4 Plot of mass verse M.

Figure5 Barents Sea: Time series of natural mortality at age based on Norwegian and Barents Seas predation mortalities.

Figure6 North Sea: Time series of natural mortality at age based on Norwegian and Barents Seas predation mortalities.

Scenarios

There are four M scenarios, the ICES assessment (M=0.15), the **Reference Case** where M is a function of mass-at-age, and two sensitivity tests where M varied based on the Barents Sea data (Environmental Driver) and where M was increased (High M).

m=t(m(om)[drop=T])
m=rbind(m,m[rep(41,2),])

cat(c("Natural Mortality", 
      "1 5",
      "1980 2022",
      "0 12", 
      "1"), sep="\n",file="../data/inputs/ices/nm-1.dat")
write(t(m),file="../data/inputs/ices/nm-1.dat",ncol=13,append=T)

m=t(m(oms[["Barents Sea"]])[drop=T])
m=rbind(m,m[rep(41,2),])

cat(c("Natural Mortality",
      "1 5",
      "1980 2022",
      "0 12", 
      "1"), sep="\n",file="../data/inputs/ices/nm-2.dat")
write(t(m),file="../data/inputs/ices/nm-2.dat",ncol=13,append=T)

m=t(m(om)[drop=T])
m=rbind(m,m[rep(41,2),])*1.25

cat(c("Natural Mortality",
      "1 5",
      "1980 2022",
      "0 12", 
      "1"), sep="\n",file="../data/inputs/ices/nm-3.dat")
write(t(m),file="../data/inputs/ices/nm-3.dat",ncol=13,append=T)

m=t(m(oms[["North Sea"]])[drop=T])
m=rbind(m,m[rep(41,2),])

cat(c("Natural Mortality",
      "1 5",
      "1980 2022",
      "0 12", 
      "1"), sep="\n",file="../data/inputs/ices/nm-4.dat")
write(t(m),file="../data/inputs/ices/nm-4.dat",ncol=13,append=T)

Figure 7 M-at-age 4

Once M has been modified then it is necessary to rerun the assessment using ‘runSam()’ to re-estimate numbers and fishing mortality-at-age.

Historical Time Series

res=list("M=0.15"               =runSam("nm"  ,dir="../data/inputs/ices"),
         "Reference"            =runSam("nm-1",dir="../data/inputs/ices"),
         "Barents Sea"          =runSam("nm-2",dir="../data/inputs/ices"),
         "High M"               =runSam("nm-3",dir="../data/inputs/ices"),
         "North Sea"            =runSam("nm-4",dir="../data/inputs/ices")) 

## function returns both assessment outputs and an FLStock object
## FLStock object is used for running the MSE
oms =FLStocks(llply(res, function(x) x[[1]]))
sams=llply(res, function(x) x[[2]])

save(oms, file="../data/om/oms.RData")  
save(sams,file="../data/om/sams.RData")  
[1] "C:/flrpapers/erp/Rmd"

Figure8 Historical time series.

(#top)

Stock recruitment relationship

Three hypotheses were run for the stock recruitment relationship which was fitted post-hoc, i.e. after numbers-at-age were estimated within the assessment.

  • Prior for steepness
  • Steepness = 0.9, i.e. recruitment environmentally driven
  • Depensation

Prior for steepness

     [,1]           [,2]         
[1,] "K"            "M"          
[2,] "Winfinity"    "Loo"        
[3,] "tmax"         "tm"         
[4,] "Lm"           "Temperature"
[5,] "ln_margsd"    "rho"        
[6,] "logitbound_h" "ln_r"

Figure 9 Steepness prior from FishBase

Figure 10 Stock recruitment relationship fits

Figure 11 Stock recruitment relationship fits

Figure 12 Production functions

Figure 13 Recruitment residuals, with regimes indicated by box.

Based on the residuals increase the expected recruitment

Hindcast Projections

The MSE was run in the form of a hindcast, i.e. a projection was run from 2001 forward using the observed weights-at-age

Initial Conditions

In order to represent the uncertainty on the starting conditions (i.e. uncertainty on all estimates from the assessment, including parameters and states), 1000 replicates of the stock were generated by resampling all the SAM parameters and all the states (abundances and fishing mortalities at age) from a multivariate normal distribution of mean and variance-covariances taken from the SAM output.

Back to Top

Software Version

Back to Top

Funding

Back to Top

References

Back to Top