Dynamics equations
The surplus production model uses the Fletcher (1978) formulation. The biomass \(B_t\) in year \(t\) is \[B_t = B_{t-1} + P_{t-1} - C_{t-1},\] where \(C_t\) is the observed catch and \(P_t\) is the surplus production given by: \[P_t = \gamma \times MSY \times \left(\dfrac{B_t}{K}-\left[\dfrac{B_t}{K}\right]^n\right), \] where \(K\) is the carrying capacity, \(MSY\) is the maximum sustainable yield, and \(n\) is the parameter that controls shape of the production curve, and \(\gamma\) is \[\gamma = \dfrac{1}{n-1}n^{n/(n-1)}.\]
The fishing mortality \(F_t\) is conditioned on the observed catch \(C_t\). If fishing occurs as a pulse fishery, then \[\hat{F}_t = \dfrac{C_t}{\hat{B}_t},\] where \(F_t\) is effectively an exploitation rate.
Continuous surplus production
Smaller time steps are used in the model to approximate continuous production and fishing. With the biomass at the start of the year, the model iteratively solves for the year-specific fishing mortality rate that produces the observed annual catch over sub-annual time steps.
For each seasonal time step \(t'\) within year,
\[ \begin{align} P_{t'} &= \Delta \times \gamma \times MSY \times \left(\dfrac{B_{t'}}{K}-\left[\dfrac{B_{t'}}{K}\right]^n\right)\\ C_{t'} &= \Delta \times F_t \times B_{t'}\\ B_{t'+1} &= B_{t'} + P_{t'} - C_{t'}\\ C_t &= \sum_{t'} C_{t'} \end{align} \]
By default, four sub-annual time steps \(\Delta = 0.25\) with \(t'=1,2,3,4\) are used to solve for \(F_t\).
Estimated and derived parameters
The model is configured with \(F_{MSY}\) and \(MSY\) as leading parameters to be estimated.
The biomass \(B_{MSY}\) at \(MSY\) is \[B_{MSY} = \dfrac{MSY}{F_{MSY}},\] the carrying capacity \(K\) is \[K = n^{1/(n-1)} \times B_{MSY} ,\] and the intrinsic rate of population increase \(r\) is \[ r = n \times F_{MSY}.\]
Schaefer model
The production parameter \(n\) is typically fixed, and the Schaefer model with a symmetric productive curve (\(B_{MSY}/K = 0.5\)) is produced when \(n = 2\) (default).
Fox model (SP_Fox)
The Fox model is the limiting case of the Fletcher parameterization as \(n \rightarrow 1\), where
\[K = e \times B_{MSY}\] \[r = F_{MSY}\] \[P_t = -e \times MSY \times \dfrac{B_t}{K} \times \log\left(\dfrac{B_t}{K}\right) \] where \(e\) is Euler’s number.
State-space version (SP_SS)
In the state-state version, annual biomass deviates are estimated as random effects. Similar to Meyer and Millar (1999), the biomass \(B_t\) in year \(t\) is \[B_t = (B_{t-1} + P_{t-1} - C_{t-1})\exp(\delta_t - 0.5 \tau^2),\] where \(\delta_t \sim N(0, \tau^2)\) are biomass deviations in lognormal space and \(\tau\) is the standard deviation of the biomass deviations.
Likelihood
By conditioning the model on observed catch, the model is fitted to the indices of abundance \(I_{s,t}\). The likelihood of the observed indices, using assuming a lognormal distribution, is \[L^I = \sum_s \lambda_s \sum_t \left(-\log(\sigma_{s,t}) - 0.5\left[\dfrac{\log(I_{s,t}) - \log(\hat{q}_s \hat{B}_t)}{\sigma_{s,t}}\right]^2\right)\] where \(lamba_s\) is the weighting coefficient of index \(s\), \(q_s\) is the cathabilitiy coefficient, and the circumflex denotes an estimate.
Log-biomass deviations \(\hat{\delta}_t\) are typically estimated as penalized parameters towards the likelihood, with the penalty: \[L^{\delta} = \sum_t \left(-\log(\tau) - 0.5\left[\dfrac{\delta_t}{\tau}\right]^2\right).\]
Prior for MSY
A lognormal prior (a penalty to the likelihood) can be provided to the model:
\[L^{\textrm{MSY}} = - 0.5\left(\dfrac{\log(\widehat{\textrm{MSY}}) - \log(\mu_{\textrm{MSY}})}{\sigma_{\textrm{MSY}}}\right)^2\]
Prior for r
A lognormal prior for the \(r\) can be provided to the model:
\[L^r = - 0.5\left(\dfrac{\log(\hat{r}) - \log(\mu_r)}{\sigma_r}\right)^2,\] Users can explicitly specify \(\log(\mu_r)\) and \(\sigma_r\).
Alternatively, the Euler-Lotka method can be used to generate the prior for the intrinsic rate of increase, by sampling natural mortality \(M\) and steepness \(h\) from a distribution. Natural mortality is modelled to be age-invariant and is sampled from a lognormal distribution. With either a Beverton-Holt or Ricker stock-recruit relationship, steepness is sampled from a transformed beta or transformed lognormal distribution, respectively. From these samples, \(\log(\mu_r)\) and \(\sigma_r\) are calculated.
For each pair of sampled M and h values, the corresponding value of \(r\) is obtained by solving a modified Euler-Lotka equation: \[\alpha \sum_a l_a f_a \exp(-r \times a) = 1.\]
Equation 16 is modified to include the \(\alpha\) term from the stock-recruit relationship (Stanley et al. 2009). In this way, the recruits-per-spawner at low stock sizes, i.e., as spawning biomass approaches zero, is considered for calculating \(r\).
The numbers-per-recruit at age \(a\) is \[ l_a = \begin{cases} 1 & a = 1\\ l_{a-1} \exp(-M_{a-1}) & a = 2, \ldots, A-1\\ \dfrac{l_{a-1} \exp(-M_{a-1})}{1 - \exp(-M_a)} & a = A \\ \end{cases}. \]
Fecundity at age \(f_a\) is \[m_a = w_a\left[1 + \exp\left(-\log(19) \dfrac{a - \tilde{a}_{50}}{\tilde{a}_{95} - \tilde{a}_{50}}\right)\right]^{-1}, \] where \(\tilde{a}_{50}\) and \(\tilde{a}_{95}\) are the ages of 50% and 95% maturity, respectively.
Weight-at-age \(w_a\) is \[w_a = W_{\infty}(1 - \exp[K\{a-a_0\}])^b.\]
The recruits per spawner at the origin of the stock-recruit relationship \(\alpha\) is \[\alpha = \dfrac{4h}{(1-h)\phi_0}, \] or \[\alpha = \dfrac{(5h)^{1.25}}{\phi_0},\] for a Beverton-Holt and Ricker stock-recruit relationship, respectively. Unfished recruits-per-spawner \(\phi_0\) is \[\phi_0 = \sum_a l_a f_a.\]
References
Fletcher, R.I. 1978. On the restructuring of the Pella-Tomlinson system. Fishery Bulletin 76:515-521.
Meyer, R., and Millar, R.B. 1999. BUGS in Bayesian stock assessments. Canadian Journal of Fisheries and Aquatic Science 56:1078-1086.
Stanley, R.D., M. McAllister, P. Starr and N. Olsen. 2009. Stock assessment for bocaccio (Sebastes paucispinis) in British Columbia waters. DFO Can. Sci. Advis. Sec. Res. Doc. 2009/055. xiv + 200 p.