29 Reference Points
This chapter gives the full mathematical specification of how reference points are calculated in openMSE. For a conceptual overview see Chapter 21.
Reference points are calculated per stock complex (see Chapter 19) 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.
For the Beverton-Holt SRR with steepness \(h\):
\[\text{RelRec}(\rho) = \frac{4h\rho - (1-h)}{(5h-1)\rho}, \qquad \rho = \frac{\Phi_F}{\Phi_0}\]
For the Ricker SRR with steepness \(h\):
\[\text{RelRec}(\rho) = \rho \cdot \exp\!\Bigl(\kappa \bigl(1 - \rho\bigr)\Bigr)\]
where \(\kappa = \log(5h)/(1 - 0.2^{\kappa})\) … [exact steepness parameterisation to be confirmed against the MSEtool SRR implementations and documented here].
The exact analytical forms of \(\text{RelRec}\) for each SRR model (Beverton-Holt, Ricker, hockey-stick) should be confirmed against the MSEtool SRR class implementations and documented here. The BH form above is the standard steepness parameterisation; the Ricker form requires verifying the steepness-to-\(\kappa\) mapping used internally.
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 Yield-Per-Recruit Reference Points
\(F_\text{Max}\) and \(F_{0.1}\) are not yet calculated automatically. Code to compute and store these reference points is under development.
\(F_\text{Max}\) is the apical fishing mortality that maximises yield per recruit:
\[F_\text{Max} = \arg\max_{F > 0} \; \sum_f Y_f^{PR}(F)\]
\(F_{0.1}\) is the apical \(F\) at which the slope of the YPR curve is 10% of its value at \(F = 0\):
\[\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 = 0}\]
Neither reference point requires a stock-recruitment relationship.
29.6 SPR-Based Reference Points
\(F_\text{crash}\) and \(\text{SPR}_\text{crash}\) are not yet calculated automatically. Code to compute and store these reference points is under development.
\(F_\text{crash}\) is the fishing mortality at which equilibrium recruitment falls to zero, and \(\text{SPR}_\text{crash} = \text{SPR}(F_\text{crash})\) is the corresponding spawning potential ratio. The existence and value of \(F_\text{crash}\) depend on the SRR.
For the Beverton-Holt SRR, equilibrium recruitment approaches zero asymptotically as \(F \to \infty\). There is no finite \(F_\text{crash}\); by convention \(\text{SPR}_\text{crash} = 0\).
For the Ricker SRR, the critical condition is that each unit of spawning production produces exactly one recruit, i.e. the per-capita productivity at low density \(\alpha\) satisfies \(\alpha \cdot \Phi_F = 1\). The analytical result is:
\[\Phi_{F_\text{crash}} = \frac{1}{\alpha} \qquad \Longrightarrow \qquad \text{SPR}_\text{crash} = \frac{1}{\alpha \cdot \Phi_0}\]
where \(\alpha\) is recoverable from the steepness \(h\) and \(\Phi_0\) via the standard Ricker steepness parameterisation.
For the hockey-stick SRR, \(F_\text{crash}\) corresponds to the SPR at which the linear segment of the SRR intersects zero recruitment; the exact value follows from the slope and breakpoint parameters.
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. It is calculated separately for landings and removals by CalcRefYield() and stored in Hist@Reference@Landings and Hist@Reference@Removals. It is not calculated by default; enable it via Simulate(DoRefLandings = TRUE).
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\)). 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.8).
29.8 Mean Generation Time
CalcMGT() exists but is not yet called automatically during Simulate(). It should be added to the default computation sequence.
MGT is computed by CalcMGT() and stored in Hist@Reference@MGT:
\[\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 not calculated by default due to its computational cost. It must be enabled as an argument to Simulate(). Code for this is under development.
\(B_\text{Low}\) is a rebuilding-based limit reference point used in DFO Canada’s Precautionary Approach framework. It is the largest current spawning biomass \(SB_{\text{now}}\) from which the population, projected forward under zero fishing, fails to reach \(B_\text{frac} \times SB_\text{MSY}\) within \(\text{HZN} \times \text{MGT}\) years:
\[B_\text{Low} = \sup\Bigl\{SB_{\text{now}} : \max_{t \leq \text{HZN} \times \text{MGT}} SB_t\bigl(SB_{\text{now}},\, F=0\bigr) < B_\text{frac} \cdot SB_\text{MSY}\Bigr\}\]
Default values are \(B_\text{frac} = 0.5\) and \(\text{HZN} = 2\). The projection uses the operating model’s stock-recruitment relationship and biological parameters from the terminal historical year. \(B_\text{Low}\) is stored in Hist@Reference@BLow.
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 21. 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.
For seasonal models, the depletion accessor functions (B_BMSY(), SB_SBMSY(), SP_SPMSY(), etc.) return ratios where the numerator is historical biomass at seasonal time-step resolution while the denominator (\(B_\text{MSY}\) etc.) is a mean-seasonal-snapshot annual quantity. This temporal mismatch means the ratios fluctuate within a year at seasonal frequency rather than varying only between years. The accessor functions should aggregate historical biomass to calendar-year means before computing the ratio. This is not yet implemented.
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 21.9.
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)})}\]
This is an approximation. The correct denominator for a seasonal model should use the annual mortality of the plus-group cohort cycling through all \(S\) seasons: \(1 - \exp(-\sum_{s=1}^{S} Z_{n_\text{age},s})\). The current single-season correction underestimates the denominator by a factor of approximately \(S\), leading to overestimation of plus-group abundance. This will be addressed in a future update.
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.
Annual biomass per recruit is:
\[B^{PR} = \frac{1}{S}\sum_{s_0=1}^{S} \pi_{s_0} \sum_a N_{a,s_0}^{PR,F} \cdot \bar{w}_{a,\,\sigma(a,s_0)}\]
The \(1/S\) factor converts the cohort-weighted sum to a mean seasonal snapshot biomass, directly comparable to start-of-year biomass in an annual model. The same correction applies to spawning biomass per recruit \(SB^{PR}\) (using \(\bar{w}_{a,s} \cdot \text{mat}_{a,s}\)). Annual spawning production \(\Phi_F\) and yield \(Y^{PR}\) are flow quantities that sum correctly over seasons without this correction.
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.
Reported quantities at \(F_\text{MSY}\):
| Quantity | Definition |
|---|---|
| \(F_\text{MSY}\) | Annual apical \(F\); sum of instantaneous rates at the apical age across all seasons |
| \(B_\text{MSY}\) | Mean seasonal snapshot total biomass at MSY \(= R_{eq} \cdot B^{PR}\) |
| \(SB_\text{MSY}\) | Mean seasonal snapshot spawning biomass at MSY \(= R_{eq} \cdot SB^{PR}\) |
| \(SP_\text{MSY}\) | Annual spawning production at MSY \(= R_{eq} \cdot \Phi_{F_\text{MSY}}\) |
| MSY | Annual yield \(= R_{eq} \cdot Y^{PR}\) |
| \(\text{SPR}_\text{MSY}\) | \(\Phi_{F_\text{MSY}} / \Phi_0\) |
29.11.8 Reference yield in seasonal models
CalcRefYield() extracts the full effort pattern from the last complete historical year (all \(S\) seasons) as the baseline. The scalar \(e^\lambda\) is applied to all fleets and seasons simultaneously, preserving the relative seasonal effort distribution. The scaled pattern is then tiled across the projection period so the seasonal effort cycle repeats each year.
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.