29 Reference Points
This chapter gives the full mathematical specification of how reference points are calculated in openMSE. For a conceptual overview see Chapter 22.
Reference points are calculated per stock complex (see Chapter 20) using a per-recruit framework combined with the stock-recruitment relationship. Results are stored in Hist@Reference and accessed via Reference(Hist).
29.1 Per-Recruit Framework
The per-recruit framework tracks the expected number, biomass, and spawning production of a single recruit as it ages under specified mortality conditions.
29.1.1 Numbers Per Recruit
Let \(N_a^{PR}\) denote the cumulative survival from recruitment to age class \(a\), expressed per unit recruit. Survival accumulates over age classes, accounting for total mortality \(Z_a = M_a + F_a^{\text{dead}}\) and spawn timing \(\omega\) (see Chapter 26 for the spawn-timing survival formulation):
\[N_1^{PR} = \exp(-Z_1 \cdot \omega)\]
\[N_a^{PR} = N_{a-1}^{PR} \cdot \exp\!\bigl(-Z_{a-1}(1 - \omega) - Z_a \cdot \omega\bigr) \cdot (1 - \gamma_{a-1}), \quad a > 1\]
where \(\gamma_{a-1}\) is semelparous mortality at age \(a-1\). The unfished numbers-per-recruit \(N_a^{PR,0}\) is obtained by setting \(F_a^{\text{dead}} = 0\).
For the plus group, survival is adjusted for ongoing mortality beyond the maximum age:
\[N_{A}^{PR} = \frac{N_{A-1}^{PR} \cdot \exp(-Z_{A-1}(1-\omega) - Z_A \cdot \omega)}{1 - \exp(-Z_A)}\]
29.1.2 Spawning Production Per Recruit and SPR
The unfished spawning production per recruit \(\Phi_0\) is:
\[\Phi_0 = \sum_a N_a^{PR,0,SP} \cdot \psi_a\]
where \(N_a^{PR,0,SP}\) is the survival to spawning time under unfished conditions and \(\psi_a\) is fecundity-at-age. This is stored in Hist@Reference@SPR0 and computed by CalcSPR0().
At fishing mortality \(F\), the fished spawning production per recruit is:
\[\Phi_F = \sum_a N_a^{PR,F,SP} \cdot \psi_a\]
The spawning potential ratio is:
\[\text{SPR}(F) = \frac{\Phi_F}{\Phi_0}\]
29.1.3 Yield Per Recruit
The removals per recruit attributable to fleet \(f\) are:
\[Y_f^{PR} = \sum_a N_a^{PR,F} \cdot \bigl(1 - \exp(-Z_a)\bigr) \cdot \frac{F_{a,f}^{\text{dead}}}{Z_a} \cdot \bar{w}_{a,f}\]
where \(\bar{w}_{a,f}\) is the fleet-specific mean weight of removed fish. Landed yield per recruit replaces \(F_{a,f}^{\text{dead}}\) with \(F_{a,f}^{\text{ret}}\) in the numerator.
CalcPerRecruit() returns a perrecruit object containing numbers-, spawning-, and biomass-per-recruit, yield-per-recruit, and SPR evaluated at one or more user-specified apical \(F\) values.
29.2 Fleet Allocation
Within a complex, apical \(F\) is distributed across fleets in proportion to effort-weighted catchability:
\[\alpha_f = \frac{E_f \cdot q_f}{\sum_{f'} E_{f'} \cdot q_{f'}}\]
where \(E_f\) is fleet effort and \(q_f\) is fleet catchability (Fleet@Catchability@Efficiency), both evaluated at the specified year(s). The fleet-specific fishing mortality-at-age is then:
\[F_{a,f}^{\text{int}} = F^{\text{apical}} \cdot \alpha_f \cdot v_{a,f}\]
with retained, discarded, and dead components computed as in Chapter 27. Fleet allocations reflect the historical effort pattern at the specified year(s).
29.3 Equilibrium Recruitment
The per-recruit framework gives yield and spawning production per unit recruit. To compute absolute equilibrium quantities, recruitment must be scaled using the stock-recruitment relationship (SRR).
29.3.1 From the SRR to equilibrium recruitment
At equilibrium, the stock-recruitment relationship and the population dynamics must be simultaneously satisfied. The SRR expresses recruitment as a function of spawning production: \(R = g(SP)\). At equilibrium, spawning production is itself the product of recruitment and the fished spawning production per recruit: \(SP = R_{eq} \cdot \Phi_F\). Substituting gives the fixed-point equation
\[R_{eq} = g(R_{eq} \cdot \Phi_F)\]
For standard stock-recruitment models this can be solved analytically for \(R_{eq}\) as a function of \(\Phi_F\) (or equivalently of \(\text{SPR} = \Phi_F / \Phi_0\)). The result can always be written as
\[R_{eq}(F) = R_0 \cdot \text{RelRec}\!\bigl(\text{SPR}(F)\bigr)\]
where \(\text{RelRec}(\text{SPR}) = R_{eq}/R_0\) is a dimensionless function of SPR only. By construction \(\text{RelRec}(1) = 1\) (unfished) and \(\text{RelRec}(0) = 0\) (recruitment collapses when spawning production is eliminated). The shape of \(\text{RelRec}\) between these limits is determined by the SRR steepness and functional form.
Both models express \(\text{RelRec}\) in terms of the compensation ratio \(CR\) (recruits per spawner at the origin implied by steepness), floored at zero to rule out negative recruitment:
For the Beverton-Holt SRR with steepness \(h\):
\[\text{RelRec}(\rho) = \max\!\left(0,\; \frac{CR\cdot\rho - 1}{(CR-1)\,\rho}\right), \qquad CR = \frac{4h}{1-h}, \qquad \rho = \frac{\Phi_F}{\Phi_0}\]
For the Ricker SRR with steepness \(h\):
\[\text{RelRec}(\rho) = \max\!\left(0,\; \frac{\log(CR\cdot\rho)}{\log(CR)}\right), \qquad CR = (5h)^{1.25}\]
For the Hockey-Stick SRR with hinge point \(S_\text{hinge}\) (a fraction of \(S_0\), Shinge), equilibrium recruitment is a step function of SPR rather than a smooth curve:
\[\text{RelRec}(\rho) = \begin{cases} 1 & \rho \geq S_\text{hinge} \\ 0 & \rho < S_\text{hinge} \end{cases}\]
Equilibrium abundance, biomass, and yield follow by multiplying the per-recruit quantities by \(R_{eq}(F)\). Absolute equilibrium quantities across a grid of apical \(F\) values are returned by CalcEquilibrium().
For multi-stock complexes, when stock \(j\) provides the spawning production driving the SRR of stock \(i\) (Stock@SRR@SPFrom), the relative recruitment for stock \(i\) is evaluated using the SPR of stock \(j\).
29.4 MSY Reference Points
\(F_{\text{MSY}}\) is found by numerical optimisation:
\[F_{\text{MSY}} = \arg\max_{F > 0} \; Y_{eq}(F)\]
where \(Y_{eq}(F) = R_{eq}(F) \cdot \sum_f Y_f^{PR}(F)\) is total equilibrium yield across all fleets in the complex. The type argument to CalcMSY() controls whether yield is defined as removals (landings plus dead discards; default) or landings only. The optimisation searches over a range of plausible \(F\) values up to OM@maxF.
Once \(F_{\text{MSY}}\) is found, the following quantities are computed and stored in Hist@Reference@MSY:
| Slot | Description |
|---|---|
FMSY |
Apical fishing mortality at MSY |
BMSY |
Total biomass at MSY |
SBMSY |
Spawning biomass at MSY |
SPMSY |
Spawning production at MSY |
SPRMSY |
Spawning potential ratio at MSY |
MSYLandings |
Landed catch at MSY |
MSYDiscards |
Dead discards at MSY |
All slots have dimensions Sim × Stock × Year. Reference points can be recalculated at any time using CalcMSY().
29.4.1 Multi-Stock Complexes
For complexes with more than one stock, the optimisation finds the single apical \(F\) maximising total yield across all stocks, with fleet allocation as above. Stock-level quantities (BMSY, SBMSY, etc.) are then evaluated at the complex-level \(F_{\text{MSY}}\), which will generally differ from the single-stock optimum for any individual stock.
29.5 Per-Recruit and Replacement-Line Reference Points
CalcRefPoints() computes \(F_{0.1}\), \(F_\text{Max}\), \(F_{x\%\text{SPR}}\), \(F_\text{Med}\), \(F_\text{crash}\), and \(\text{SPR}_\text{crash}\), storing them in Hist@Reference@F01, @FMax, @FSPR, @FMed, @FCrash, and @SPRcrash. It is called automatically by Simulate() (SimControl(RefPoints = TRUE), the default).
All six quantities are evaluated on a common log-spaced grid of nF apical \(F\) values (default 100) spanning approximately 0 to OM@maxF, using a single call to CalcPerRecruit() over the grid rather than re-optimising per reference point.
29.5.1 \(F_\text{Max}\) and \(F_{0.1}\)
\(F_\text{Max}\) is the grid point maximising yield per recruit:
\[F_\text{Max} = \arg\max_{F \in \text{grid}} \; \sum_f Y_f^{PR}(F)\]
For asymptotic selectivity, yield per recruit has no interior maximum and increases monotonically with \(F\); in that case \(F_\text{Max}\) is reported as the upper end of the grid (OM@maxF) rather than left undefined.
\(F_{0.1}\) is the apical \(F\) at which the slope of the YPR curve is 10% of the slope at the first grid interval:
\[\left.\frac{d}{dF}\sum_f Y_f^{PR}(F)\right|_{F = F_{0.1}} = 0.1 \cdot \left.\frac{d}{dF}\sum_f Y_f^{PR}(F)\right|_{F \to 0}\]
found by computing the discrete slope between consecutive grid points and linearly interpolating to the target slope. Neither \(F_\text{Max}\) nor \(F_{0.1}\) requires a stock-recruitment relationship.
29.5.2 \(F_{x\%\text{SPR}}\)
For each target in SPRTarget (default c(0.3, 0.4, 0.5)), \(F_{x\%\text{SPR}}\) is the apical \(F\) at which \(\text{SPR}(F)\) equals that target, found by linear interpolation of \(\text{SPR}(F)\) over the grid. Results are stored with an additional Target dimension, one value per requested target.
29.5.3 \(F_\text{Med}\)
\(F_\text{Med}\) is the apical \(F\) at which the per-recruit replacement line equals the median historical recruits-per-spawner ratio. For each stock, the historical ratio of recruits (numbers at the first age class) to spawning biomass is computed for every historical year, and its median across years gives the target. \(F_\text{Med}\) is then the \(F\) at which \(1/\Phi_F\) equals that target, found by interpolation over the grid.
29.6 SPR-Based Reference Points
\(F_\text{crash}\) and \(\text{SPR}_\text{crash} = \text{SPR}(F_\text{crash})\) are the apical \(F\) (and corresponding SPR) at which the per-recruit replacement line \(1/\Phi_F\) crosses a compensation-ratio threshold \(\alpha\) derived from stock-recruitment steepness. They are only defined for Beverton Holt and Ricker stock-recruitment models; other models return NA.
\(\alpha\) is the compensation ratio \(CR\) divided by a fixed reference unfished spawning-per-recruit \(\Phi_0^{\text{ref}}\), evaluated at the first historical year rather than the reference year used elsewhere, matching how the SRR’s fixed \(\alpha\)/\(\beta\) parameters were originally calibrated:
\[\alpha = \frac{CR}{\Phi_0^{\text{ref}}}, \qquad CR = \begin{cases} \dfrac{4h}{1-h} & \text{BevertonHolt} \\[4pt] (5h)^{1.25} & \text{Ricker} \end{cases}\]
\(F_\text{crash}\) is found by interpolating \(1/\Phi_F\) against \(\alpha\) over the same \(F\) grid used for \(F_\text{Max}\) and \(F_{0.1}\):
- If \(1/\Phi_F \geq \alpha\) across the entire grid (the stock is already unable to replace itself even at the lowest grid \(F\)), \(F_\text{crash} = 0\) and \(\text{SPR}_\text{crash} = 1\).
- If \(1/\Phi_F \leq \alpha\) across the entire grid (replacement is never exceeded within the searched range), \(F_\text{crash}\) is reported at the upper grid bound and \(\text{SPR}_\text{crash}\) is extrapolated (floored at 0.01).
- Otherwise, both are found by linear interpolation at the point where \(1/\Phi_F\) crosses \(\alpha\).
29.7 Reference Yield
The reference yield is a simulation-based yield benchmark: the highest yield achievable under a fixed effort policy over the projection period.
For each simulation, a log-scalar \(\lambda\) is applied uniformly to the effort from the last complete historical year:
\[E_{f,t}^{\text{proj}} = E_{f,t_0} \cdot e^{\lambda}\]
The scaled effort is held constant across the projection period and the population is projected forward using the standard population dynamics equations. The mean yield over the last \(n_\text{avg}\) calendar years of the projection is recorded (default \(n_\text{avg} = 5\), OM@Control$RefYield$lastnTS; internally converted to time steps via OM@Seasons so seasonal models average whole years). The \(\lambda\) that maximises this mean yield is found by one-dimensional optimisation:
\[\lambda^* = \arg\max_\lambda \; \overline{Y}(\lambda)\]
and the resulting maximum mean yield is the reference yield. For seasonal models, effort is extracted for all \(S\) seasons of the last historical year and tiled across the projection, preserving the seasonal effort distribution (see Section 29.11).
29.8 Mean Generation Time
MGT is computed by CalcMGT() and stored in Hist@Reference@MGT. It is calculated automatically by Simulate() (SimControl(MGT = TRUE), the default):
\[\text{MGT} = \frac{\sum_a a \cdot l_a \cdot m_a}{\sum_a l_a \cdot m_a}\]
where \(l_a\) is survivorship to age \(a\) under natural mortality only, with \(l_1 = 1\) and \(l_{a+1} = l_a \cdot \exp(-M_a)\), and \(m_a\) is maturity-at-age. Ages \(a\) are in years; for seasonal models they are in decimal-year units so MGT is returned in years regardless of seasonal resolution.
29.9 \(B_\text{Low}\)
\(B_\text{Low}\) is a rebuilding-based limit reference point used in DFO Canada’s Precautionary Approach framework, computed by CalcBLow() and stored in Hist@Reference@BLow. For each simulation and stock, numbers-at-age in the last historical year are scaled by a factor \(k\) and projected forward under zero fishing effort (a fixed zero-effort MP passed to Project()) for \(\text{round}(\text{HZN} \times \text{MGT})\) years. \(k\) is found by one-dimensional optimisation so that the projected spawning biomass at the end of that horizon equals the target:
\[k^* = \arg\min_{k} \; \Bigl(\log SB_{\text{HZN}\times\text{MGT}}(k) - \log(B_\text{frac} \cdot SB_\text{MSY})\Bigr)^2\]
\[B_\text{Low} = k^* \cdot SB_{\text{now}}\]
Default values are \(B_\text{frac} = 0.5\) and \(\text{HZN} = 2\). SB_MSY is taken from Hist@Reference@MSY (computed via CalcMSY() if not already present) and MGT from CalcMGT(). The rebuilding horizon is capped at the number of available projection years, with a warning, if it would otherwise exceed them.
29.10 Depletion and Stock Status
Current stock status relative to unfished and MSY reference levels is expressed as ratios of historical biomass, spawning biomass, spawning production, or fishing mortality to their respective reference values. These ratios have dimensions Sim × Stock × Year and are accessible via the helper functions listed in Chapter 22. Unfished reference values (\(B_0\), \(SB_0\), \(SP_0\)) come from Hist@Unfished; MSY-based reference values (\(B_\text{MSY}\), \(SB_\text{MSY}\), \(SP_\text{MSY}\), \(F_\text{MSY}\)) come from Hist@Reference@MSY.
29.11 Seasonal Models
When OM@Seasons > 1, age classes correspond to season-steps rather than annual years and the per-recruit calculation must track cohorts through the full seasonal cycle. All reference point quantities are returned on an annual scale. For a conceptual overview see Section 22.8.
29.11.1 Notation
| Symbol | Meaning |
|---|---|
| \(S\) | Number of seasons per year |
| \(s_0 \in \{1,\ldots,S\}\) | Birth season of a cohort |
| \(\sigma(a,s_0)\) | Calendar season experienced by a cohort born in \(s_0\) at age class \(a\) |
| \(\pi_{s_0}\) | Fraction of annual unfished recruitment entering in season \(s_0\) |
| \(Z_{a,s}\) | Total mortality (M + F dead) for age class \(a\) in season \(s\) |
| \(\psi_{a,s}\) | Fecundity-at-age \(a\) in season \(s\) (seasonally adjusted by AdjustSeasonalFecundity) |
| \(\omega\) | Spawn-timing fraction within a season (SRR@SpawnTimeFrac) |
The calendar-season function: \[\sigma(a, s_0) = (s_0 + a - 2) \bmod S + 1\] maps birth season \(s_0\) and age class \(a\) to the calendar season that cohort occupies at age \(a\).
29.11.2 Year specification
Reference points are evaluated using parameters from one or more complete calendar years. The Years argument to CalcMSY(), CalcPerRecruit(), and CalcEquilibrium() accepts integer calendar years; all \(S\) season-steps within each requested year are used automatically. By default the last complete historical year is used.
29.11.3 Birth-season cohort tracking
For each birth season \(s_0\), the numbers per recruit are computed by tracking the cohort through seasons using season-specific mortality:
\[N_{1,s_0}^{PR} = \exp\!\bigl(-Z_{1,\,\sigma(1,s_0)}\cdot\omega\bigr)\]
\[N_{a,s_0}^{PR} = N_{a-1,s_0}^{PR} \cdot \exp\!\Bigl(-Z_{a-1,\,\sigma(a-1,s_0)}(1-\omega) - Z_{a,\,\sigma(a,s_0)}\,\omega\Bigr) \cdot \bigl(1 - \gamma_{a-1,\,\sigma(a-1,s_0)}\bigr), \quad a > 1\]
where \(\gamma_{a,s}\) is semelparous post-spawning mortality in season \(s\) (zero for iteroparous species). The unfished variants \(N_{a,s_0}^{PR,0}\) and \(N_{a,s_0}^{PR,0,SP}\) are obtained by setting \(F = 0\) (i.e., \(Z_{a,s} = M_{a,s}\)).
For the plus group, the current implementation applies a single-season correction using the mortality in the season the cohort occupies at the maximum age class:
\[N_{n_\text{age},s_0}^{PR} \;\leftarrow\; \frac{N_{n_\text{age},s_0}^{PR}}{1 - \exp(-Z_{n_\text{age},\,\sigma(n_\text{age},s_0)})}\]
Spawning-time numbers \(N_{a,s_0}^{PR,SP}\) are computed with the same recurrence but using \(Z^{\text{SP}}_{a,s} = Z_{a,s}\cdot\omega\) for the initial step and \(Z_{a-1,s}(1-\omega) + Z_{a,s}\,\omega\) for subsequent steps, as in Section 29.1.
29.11.4 Annual per-recruit aggregation
Annual per-recruit quantities aggregate over birth seasons weighted by the seasonal recruitment proportions \(\pi_{s_0}\):
\[\Phi_0 = \sum_{s_0=1}^{S} \pi_{s_0} \sum_{a=1}^{n_\text{age}} N_{a,s_0}^{PR,0,SP} \cdot \psi_{a,\,\sigma(a,s_0)}\]
\[\Phi_F = \sum_{s_0=1}^{S} \pi_{s_0} \sum_{a=1}^{n_\text{age}} N_{a,s_0}^{PR,F,SP} \cdot \psi_{a,\,\sigma(a,s_0)}\]
\[\text{SPR}(F) = \frac{\Phi_F}{\Phi_0}\]
AdjustSeasonalFecundity() scales \(\psi_{a,s}\) by season in proportion to the seasonal recruitment pattern \(\pi_{s_0}\), so that \(\Phi_0\) and \(\Phi_F\) measure annual spawning production per recruit weighted by the timing and magnitude of spawning. Seasons where \(\pi_{s_0} = 0\) have \(\psi_{a,s} = 0\); otherwise fecundity is scaled toward seasons of higher recruitment.
Unlike \(\Phi_0\), \(\Phi_F\), and yield, which are flow quantities summed over a cohort’s full lifetime, biomass and spawning biomass are stock quantities reported as a snapshot at one point in the year. Rather than averaging over all \(S\) within-year snapshots, \(B^{PR}\) and \(SB^{PR}\) are a cross-sectional snapshot at one or more reference season(s) OM@RefSeason: for each reference season \(s_\text{ref}\), the age classes that each birth-season cohort currently occupies at \(s_\text{ref}\) are summed, weighted by that cohort’s recruitment fraction \(\pi_{s_0}\):
\[B^{PR} = \sum_{s_0=1}^{S} \pi_{s_0} \sum_{a:\, \sigma(a,s_0) = s_\text{ref}} N_{a,s_0}^{PR,F} \cdot \bar{w}_{a,s_\text{ref}}\]
with \(SB^{PR}\) computed identically but weighting by \(\bar{w}_{a,s_\text{ref}} \cdot \text{mat}_{a,s_\text{ref}}\) instead. When RefSeason is NULL (the default), the reference season(s) are auto-detected per simulation as the season(s) in which fecundity (or maturity, if fecundity isn’t set) is nonzero (i.e. the spawning season(s)) falling back to an equal weighting across all seasons if none is detected. If more than one reference season applies, \(B^{PR}\) and \(SB^{PR}\) are averaged equally across them.
29.11.5 Annual yield per recruit
Annual removals per recruit (summed over fleets \(f\)):
\[Y^{PR} = \sum_{s_0=1}^{S} \pi_{s_0} \sum_{a=1}^{n_\text{age}} \sum_f N_{a,s_0}^{PR,F} \cdot \bigl(1 - \exp(-Z_{a,\,\sigma(a,s_0)})\bigr) \cdot \frac{F_{a,f,\,\sigma(a,s_0)}^{\text{dead}}}{Z_{a,\,\sigma(a,s_0)}} \cdot \bar{w}_{a,f,\,\sigma(a,s_0)}\]
Landed yield per recruit replaces \(F_{a,f,s}^{\text{dead}}\) with \(F_{a,f,s}^{\text{ret}}\) in the numerator.
29.11.6 Seasonal fleet allocation
In a seasonal model, effort and catchability vary by season. The fleet allocation array \(\alpha_{f,s}\) gives the relative fishing pressure of fleet \(f\) in season \(s\):
\[\alpha_{f,s} = E_{f,s} \cdot q_{f,s}\]
where \(E_{f,s}\) and \(q_{f,s}\) are effort and catchability efficiency for fleet \(f\) in season \(s\) from the specified year(s). The scalar apical \(F\) multiplies this pattern so that the most heavily fished season carries exactly the apical rate:
\[F_{a,f,s}^{\text{int}} = F^{\text{apical}} \cdot \frac{\alpha_{f,s}}{\displaystyle\max_s \sum_{f'} \alpha_{f',s}} \cdot v_{a,f,s}\]
The internal optimisation scalar therefore corresponds to the peak-season apical fishing mortality. The seasonal effort pattern determines how \(F\) is distributed across seasons; only the overall scale is optimised.
29.11.7 MSY optimisation and reporting
The optimisation and equilibrium recruitment steps are identical to the annual case (Section 29.4, Section 29.3), operating on the annually-aggregated per-recruit quantities above. The optimisation searches over the peak-season scalar; \(F_\text{MSY}\) is then expressed as the annual apical fishing mortality (sum of seasonal instantaneous rates at the apical age), directly comparable to \(F_\text{MSY}\) from an annual model.
The same peak-season-scalar optimisation and annual-apical-\(F\) reporting applies to \(F_{0.1}\), \(F_\text{Max}\), \(F_{x\%\text{SPR}}\), and \(F_\text{crash}\) (Section 29.5, Section 29.6), since all four are evaluated over the same annually-aggregated per-recruit grid used for \(F_\text{MSY}\).
29.12 Spatial Structure
The per-recruit framework tracks a single cohort through time without spatial structure. When a stock occupies multiple areas with movement, the \(F\) experienced by a cohort depends on which area it occupies at each time step, which a scalar per-recruit calculation cannot resolve.
In the current implementation, parameters used for reference point calculations (M, selectivity, weight, catchability) do not vary spatially. Reference points for spatially structured models should therefore be interpreted as population-level averages and used cautiously when effort or closures are strongly area-specific.
A spatially-explicit equilibrium approach, iterating the full spatial dynamics to equilibrium at each candidate \(F\), is planned as a user-controlled option via OM@Control.