2 Description of the Operating Model (OM)
This section describes the base configuration of the population dynamics equations used in the simulation framework. Please be aware that this configuration can be tailored to specific applications. An operating model is a mathematical representation of the biology of a fish stock, the fishery that operates on the fish stock, and ways in which the data is collected. The operating model should also include a sub-model that reflects how the management regulations are implemented and are adhered to in practice (Punt et al. 2016).
2.1 Population dynamics
The operating model (OM) describes the age-structured population dynamics of a fish population or stock. Population dynamics are specified for a single sex, with dynamics operating on an annual time step, and allowing for migration (of all age classes) between multiple areas. Abundance at age, \(a\), at the start of year \(t\), and in area \(i\), \(N_{a,t,i}\) is given by:
\[ N_{0,t,i} = \rho_i R_t, \tag{eq. 1} \label{eq:first} \]
And,
\[ N_{a+1,t+1,i} = S_{a,t,i} \left(N_{a,t,i}\theta_{i \rightarrow i}+\sum_{j \neq i} N_{a,t,j}\theta_{j \rightarrow i}\right). \tag{eq. 2} \label{eq:second} \]
Equation \(\ref{eq:first}\) is the fraction, \(\rho_i\), of total age-0 recruits, \(R_t\), at the start of the year that are added to area \(i\). In equation \(\ref{eq:second}\), the term in parenthesis describes the fraction of abundance-at-age that does not emigrate from area \(i\), \(\theta_{i \rightarrow i}\), plus the summation of abundance-at-age arriving from all other areas, \(j\), \(\theta_{a,j \rightarrow i}\). In this form, Equation \(\ref{eq:second}\) specifies migration to occur at the beginning of the year, followed by survival, \(S_{a,t,i}\), from age \(a\) to age \(a+1\).
The multi-area model is implemented in matrix form, analogous to a Leslie matrix, as described by Quinn and Deriso (1999). For brevity, the matrix form is summarized using a two-area example, reflecting the description provided by Quinn and Deriso (1999), but note that the operating model is generalized as a multi-area model. To account for migration between two areas, \(i\) and \(j\), abundance-at-age in the matrix form is:
\[ \mathbf{N}_{t+1} = \mathbf{P}\mathbf{N}_t, \tag{eq. 3} \label{eq:third} \]
where \(\mathbf{N}_t\) for ages 0 to maximum age, \(A\), is written:
\[ \mathbf{N}_t = \begin{pmatrix} \mathbf{N}_{t,i} \\ \mathbf{N}_{t,j} \end{pmatrix} = \begin{pmatrix} N_{1,t,i} \\ \vdots \\ N_{A,t,i} \\ N_{1,t,j} \\ \vdots \\ N_{A,t,j} \end{pmatrix}. \tag{eq. 4} \label{eq:fourth} \]
The projection matrix \(\mathbf{P}\) is
\[ \mathbf{P} = \begin{pmatrix} P_{i,j} & P_{i,j} \\ P_{i,j} & P_{i,j} \end{pmatrix}, \tag{eq. 5} \label{eq:fifth} \]
where each element, \(P_{i,j}\), is a matrix that accounts for movement and survival. The matrices, \(P_{i,j}\), are populated as:
\[ P_{i,i} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ S_{1,t,i}\theta_{1,i\rightarrow i} & 0 & 0 & 0 \\ 0 & S_{2,t,i}\theta_{2,i\rightarrow i} & 0 & 0 & 0 \\ 0 & 0 & \ddots & 0 & 0 \\ 0 & 0 & 0 & S_{A-1,t,i}\theta_{A-1,i\rightarrow i} & S_{A,t,i}\theta_{A,i\rightarrow i} \end{pmatrix}, \tag{eq. 6} \label{eq:sixth} \]
And,
\[ P_{i,j} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ S_{1,t,i}\theta_{1,j\rightarrow i} & 0 & 0 & 0 \\ 0 & S_{2,t,i}\theta_{2,j\rightarrow i} & 0 & 0 & 0 \\ 0 & 0 & \ddots & 0 & 0 \\ 0 & 0 & 0 & S_{A-1,t,i}\theta_{A-1,j\rightarrow i} & S_{A,t,i}\theta_{A,j\rightarrow i} \end{pmatrix}. \tag{eq. 7} \label{eq:seventh} \]
The operating model specifies a plus group in position (\(A,A\)).
Total annual recruitment is calculated according to the Beverton and Holt (1957) stock-recruitment function, which is parameterized using steepness \((h)\):
\[ R_t = \left(\frac{0.8R_0hB_t}{0.2B_0\left(1-h\right) + \left(h-0.2\right)B_t}\right) \exp{\left(d_t - \frac{{\sigma_R}^2}{2}\right)}, \tag{eq. 8} \label{eq:eighth} \]
where \(B_0\) is unfished reproductive output, \(R_0\) is unfished recruitment, and \(d\) is the annual recruitment deviation, where \(\sigma_R\) is recruitment variability. Annual reproductive output, \(B_t\), is aggregate output (e.g., spawning biomass) for all areas combined, excluding any reproductive contribution of age-0 recruits created at the beginning of the year. Calculated total recruits are distributed to each area in proportion to the area-specific fraction, \(\rho_i\), with \(\sum\rho_i=1\). An overall log-scale recruitment deviation signal (e.g., a surrogate for an environmental condition affecting recruitment success) is generated from a normal distribution with mean zero and standard deviation, \(\sigma_R\), (e.g., \(\sigma_R\) = 0.6). Inter-annual autocorrelation in recruitment can be specified, producing 1-year lagged correlation in log-scale recruitment deviations.
2.2 Life history
Length is calculated at the start of the year according to the von Bertalanffy growth curve, where \(L_\infty\) is asymptotic length, \(K\) is the growth rate, and parameter \(t_0\):
\[ L_a = L_\infty \left(1- \exp{\left(-K\left(a - t_0\right)\right)}\right). \tag{eq. 9} \label{eq:nineth} \]
Weight-at-age, \(W_a\), with parameters \(a\) and \(b\) is specified as an exponential function:
\[ W_a=\alpha L_a^\beta. \tag{eq. 10} \label{eq:tenth} \]
Maturation follows a logistic function with parameters \(L50\) and \(L95\), reflecting the lengths at which 50% and 95% of the population are mature, respectively. Optionally, species can be specified as protogynous hermaphroditic species, with proportion of male in the population following an increasing logistic function with parameters \(H50\) and \(H95\), reflecting the lengths at which 50% and 95% of the population are male, respectively. For gonochoristic species, a 50:50 sex ratio is assumed at all lengths or ages. Total reproductive output, \(B_t\), is a summation of mature biomass across age classes and areas:
\[ B_t=\sum_{i}\sum_{a}{N_{a,t,i}W_a\mathit{mat}_a{propFemale}_a}, \tag{eq. 11} \label{eq:eleventh} \]
Where \(propFemale\) is proportion of the population female at age, with a value of 0.5 for all ages for gonochoristic species and values of 1-\(propMale_a\) when sexual transition from female to male is specified.
Both natural mortality and maximum age can be specified. When maximum age, \(A\), is specified, this quantity is used in constructing abundance matrices. Maximum must be equal to or greater than 2, as this modeling framework is not well suiting to species with very fast life histories. When maximum age is not specified, the age to which 1% the population survives in an unfished system is used to calculate maximum age, using the formula:
\[ A=\ ceiling\left(-\frac{\log{\left(0.01\right)}}{M}\right) \tag{eq. 12} \label{eq:twelftth} \]
Uncertainty in life history can be accounted for by specifying parameter ranges, rather than point estimates for most life history parameters. Each iteration will produce a unique set of life history parameters based on independent draws from uniform distributions that correspond to the specified minimum and maximum for each parameter.
2.3 Initial conditions
This modeling framework was developed to create historical dynamics of fish stocks that begin (i.e., year 0) in a fished state, meaning that fishing mortality (and consequently fishing effort) are greater than zero in the initial equilibrium year (year 0). Thus, the modeling framework is not suitable for circumstances for initializing the model in an unfished or pre-fishing state. Accordingly, initial depletion (spawning biomass relative to unfished spawning biomass) should always be less than 1.0. An alternative formulation is available where initial ‘depletion’ can instead be specified as initial Spawning Potential Ratio (SPR). SPR is spawning biomass per recruit relative to unfished spawning biomass per recruit). Subsequently, the population is initialized as follows. First, equilibrium age-structure is determined for the given depletion or SPR level assuming area-specific recruitment fractions, \(\rho_i\), but no movement, resulting in an equilibrium fishing mortality rate and equilibrium abundance scaled relative to the specified \(R_0\). Second, a burn-in period is used to project the population forward for \(A\)x4 years at the estimated equilibrium fishing mortality rate, allowing a stable age distribution between areas to be obtained through migration.
2.4 Fishery dynamics
Survival (\(S\)) consists of natural mortality (\(M\) year\(^{-1}\), fishing mortality (\(F\) year\(^{-1}\)):
\[ S_{a,t,i}=exp{\left(-M-{\mathit{Removal}}_aF_{t,i}\right)}, \tag{eq. 13} \label{eq:thirteenth} \]
where \(Removal\) is a component of fishery selectivity and is covered in a subsequent section.
Landings in numbers (\(C^{N}\) is:
\[ C_{a,t,i}^N=\frac{{\mathit{Keep}}_aF_{t,i}}{\left(M_s+{\mathit{Removal}}_aF_{t,i}\right)}\left(1-S_{a,t,i}\right)N_{a,t,i}, \tag{eq. 14} \label{eq:fourteenth} \]
and \(keep\) is a component of fishery selectivity. Landings in weight (\(C^B\)) is:
\[ C_{a,t,i}^B=C_{a,t,i}^NW_a \tag{eq. 15} \label{eq:fifteenth} \]
Fishery selectivity is defined as follows. Vulnerability to the fishing gear includes several options. Retention, \(Ret\), can be specified as ‘full’, resulting in full retention across all vulnerable size classes, ‘logistic’ or ‘slot limit’. Additionally, the maximum level of retention (e.g., a quantity between 0 and 1) can be specified for any of the above stated retention types. Finally, a discard mortality proportion, \(D\), (e.g. quantity between 0 and 1) can be specified to affect the fate of discards. These inputs are used in calculating the following components of fishery selectivity.
Note: vulnerability and retention are specified as functions of length. These components of fishery selectivity are converted to age using corresponding length-at-age.
\(Keep\) - the resulting probability of being landed \[ Keep_a=Vul_aRet_a \tag{eq. 16} \label{eq:sixteenth} \]
\(Dead\,discards\) - deaths resulting from vulnerable abundance that is not retained \[ Dead\,discards_a=Vul_a(1-Ret_a)D \tag{eq. 17} \label{eq:seventeenth} \]
\(Total\,removals\) – probability of removal from the population via landing or discard \[ Removal_a=Vul_a(Ret_a+(1-Ret_a)D) \tag{eq. 18} \label{eq:eighteenth} \]
Uncertainty in vulnerability and retention can be accounted for by specifying parameter ranges, rather than point estimates. Each iteration will produce a unique selectivity and/or retention characteristics based on independent draws from uniform distributions that correspond to the specified minimum and maximum for each parameter.
2.5 Observation and monitoring
Depending on the specific application of the model, an observation model may be required. This requirement is driven by whether a harvest control rule is used in decision-making. That is, whether some form of data collection will inform year-to-year adjustments in a total allowable catch or total fishery effort. If so, then simulation of the data observation processes (including imperfect observation) is required. Observation models tend to be tailored to the type of data collection program that is in place. This component should be developed as needed for each unique application.
2.6 Growth-type group
To account for individual variation in growth trajectories, the population is divided into growth-type groups (GTGs) or cohorts Hordyk et al. (2016). GTGs describe the variation in growth of a fish population through the creation of \(G\) groups and with dimensions indexed \(g=1,2…G\). Each group differs in terms of its \(L_\infty\). Given a mean value, \(\bar{L}_\infty\) and a standard deviation (calculated as a coefficient of variation, \(CV_L\) times \(\bar{L}_\infty\)), the range of \(L_\infty\) ± two standard deviations is divided into \(G\) equal increments. A default value of \(CV_L\) = 0.1 and \(G\) = 13 groups, results in \(g\) = 7 representing \(\bar{L}_\infty\). Cohort size assigned to each group is determined by distributing a fraction of annual recruits in each group, \(p\), where \(\sum{p}=1\). Vector \(p\) is determined as proportional to the \(G\) increments in \(L_\infty\) along the probability density function of a normal distribution with mean \(\bar{L}_\infty\) and variance \(\left(CV_L{\bar{L}}_\infty\right)^2\). Functionality is also included such that if \(G\) = 1, the model collapses to a simpler age-based model.