Risk‑calibrated harvest control rules for applied fisheries under productivity shifts

A pollock‑like simulation with steepness uncertainty

1 Abstract

We evaluate whether a transparent, SPR‑based harvest control rule (HCR) can be made robust to productivity shifts by calibrating its target fishing mortality to a risk criterion. Using a pollock‑like age‑structured simulation, we compare low and high steepness scenarios and apply a single sloping HCR derived from F40% with a biomass ramp. We then choose a common Ftarget that limits the combined probability of terminal biomass below 0.2 B0 to 0.05 across both scenarios. Results show that the risk‑calibrated HCR provides consistent biomass protection under steepness uncertainty while maintaining catch levels comparable to constant‑F strategies in productive regimes. The approach offers a transparent, defensible rule that does not require time‑varying environmental covariates yet remains robust to productivity shifts.

2 Introduction

Climate‑driven changes in recruitment productivity create uncertainty in reference points and the performance of fixed harvest strategies (Szuwalski et al. 2023; Punt et al. 2024). In Alaska groundfish management, SPR‑based proxies and tiered rules provide transparency and accountability (North Pacific Fishery Management Council 2024). The question is whether a simple, transparent HCR can be calibrated to remain robust when stock productivity shifts.

This study evaluates a pollock‑like simulation under two steepness regimes (h = 0.6 and h = 0.9) and asks whether a single, risk‑calibrated HCR can balance biomass protection and yield across both regimes. The primary contribution is a calibration method that selects a common Ftarget to control terminal biomass risk across scenarios.

3 Methods

3.1 Population and recruitment

We use a 15‑age‑class model with pollock‑like growth and maturity schedules. Recruitment follows a Beverton‑Holt relationship with lognormal process error (CV = 0.7). Steepness (h) is the primary productivity uncertainty, with two scenarios representing low and high recruitment resilience (Hollowed et al. 2020).

3.2 Harvest control rule

The HCR applies F40% above 0.4 B0, ramps linearly to zero at 0.05 B0, and sets F = 0 below 0.05 B0. We compute scenario‑specific F40% values for reference but evaluate a single risk‑calibrated Ftarget shared across scenarios.

3.3 Management objectives and metrics

objectives <- data.frame(
  Objective = c(
    "Avoid low biomass",
    "Maintain long‑term yield",
    "Limit instability"
  ),
  Metric = c(
    "P(SSB_T < 0.2 B0)",
    "Mean catch",
    "CV of annual catch"
  ),
  Threshold = c(
    "≤ 0.05 (combined across scenarios)",
    "Maximize subject to risk constraint",
    "Report only"
  ),
  check.names = FALSE
)

knitr::kable(objectives)

Table 1: Management objectives and performance metrics.

Objective	Metric	Threshold
Avoid low biomass	P(SSB_T < 0.2 B0)	≤ 0.05 (combined across scenarios)
Maintain long‑term yield	Mean catch	Maximize subject to risk constraint
Limit instability	CV of annual catch	Report only

3.4 Performance metrics and risk calibration

Performance metrics include mean catch, terminal biomass, and the probability that terminal SSB falls below 0.2 B0. We calibrate Ftarget to satisfy:

\[\frac{1}{2}\left[P(SSB_T < 0.2 B_0 \mid h=0.6) + P(SSB_T < 0.2 B_0 \mid h=0.9)\right] = 0.05\]

The calibrated Ftarget is then used in the sloping HCR for both scenarios.

4 Results

4.1 Reference points and yield curves

scenario_low <- create_scenario(h = 0.6, "Low, h=0.6")
scenario_high <- create_scenario(h = 0.9, "High, h=0.9")
scenarios <- list(scenario_low, scenario_high)

comparison_df <- compare_reference_points(scenarios)
knitr::kable(comparison_df, align = c("l", "c", "c", "c", "c", "c", "c", "c"))

Table 2: Reference point comparison between low and high steepness scenarios

Scenario	Steepness	Fmsy	Bmsy (Mt)	MSY (Mt)	SSB0 (Mt)	Bmsy/SSB0	R_msy (B)
Low, h=0.6	0.6	0.3517	2.66	1.31	8.29	0.321	9.98
High, h=0.9	0.9	1.2805	1.59	2.00	8.29	0.192	12.08

plot_yield_comparison(scenarios, colors = c("blue", "red"))

Figure 1: Equilibrium yield vs. fishing mortality for both steepness scenarios. Vertical dashed lines indicate Fmsy.

4.2 HCR and risk‑calibrated Ftarget

ref_low <- scenario_low$ref_points
ref_high <- scenario_high$ref_points

F40_low <- calc_F_spr(scenario_low$params, spr_target = 0.4)
F40_high <- calc_F_spr(scenario_high$params, spr_target = 0.4)

f40_table <- data.frame(
  Scenario = c("Low steepness (h=0.6)", "High steepness (h=0.9)"),
  `F40%` = round(c(F40_low, F40_high), 3),
  check.names = FALSE
)

knitr::kable(f40_table)

Table 3: F40% (SPR) values by steepness scenario.

Scenario	F40%
Low steepness (h=0.6)	0.412
High steepness (h=0.9)	0.412

target_risk <- 0.05
b0_frac <- 0.2
n_sims_risk <- 200
seed_low <- 2200
seed_high <- 3200

calc_terminal_b0_risk <- function(scenario, F_target, B0, n_sims, seed_base, b0_frac) {
  params <- scenario$params
  n_years <- params$n_years
  below <- logical(n_sims)

  for (i in 1:n_sims) {
    sim <- run_simulation_hcr(
      params,
      F_target = F_target,
      B0 = B0,
      catch_cap = NULL,
      seed = seed_base + i
    )
    below[i] <- sim$SSB[n_years] < (b0_frac * B0)
  }

  mean(below)
}

combined_b0_risk <- function(F_target) {
  risk_low <- calc_terminal_b0_risk(
    scenario_low,
    F_target,
    ref_low$SSB0,
    n_sims_risk,
    seed_low,
    b0_frac
  )
  risk_high <- calc_terminal_b0_risk(
    scenario_high,
    F_target,
    ref_high$SSB0,
    n_sims_risk,
    seed_high,
    b0_frac
  )

  list(
    risk_low = risk_low,
    risk_high = risk_high,
    risk_combined = mean(c(risk_low, risk_high))
  )
}

F_grid <- seq(0.05, 1.0, by = 0.05)
risk_grid <- lapply(F_grid, combined_b0_risk)
risk_df <- data.frame(
  F_target = F_grid,
  Risk_combined = sapply(risk_grid, function(x) x$risk_combined),
  Risk_low = sapply(risk_grid, function(x) x$risk_low),
  Risk_high = sapply(risk_grid, function(x) x$risk_high)
)

risk_min <- min(risk_df$Risk_combined)
risk_max <- max(risk_df$Risk_combined)

if (target_risk <= risk_min) {
  Ftarget <- risk_df$F_target[which.min(risk_df$Risk_combined)]
} else if (target_risk >= risk_max) {
  Ftarget <- risk_df$F_target[which.max(risk_df$Risk_combined)]
} else {
  ordered <- risk_df[order(risk_df$Risk_combined), ]
  Ftarget <- approx(ordered$Risk_combined, ordered$F_target, xout = target_risk)$y
}

risk_at_target <- combined_b0_risk(Ftarget)

ftarget_table <- data.frame(
  Metric = c(
    "Ftarget",
    "Combined P(SSB < 0.2 B0)",
    "Low steepness P(SSB < 0.2 B0)",
    "High steepness P(SSB < 0.2 B0)"
  ),
  Value = c(
    round(Ftarget, 3),
    round(risk_at_target$risk_combined, 3),
    round(risk_at_target$risk_low, 3),
    round(risk_at_target$risk_high, 3)
  )
)

knitr::kable(ftarget_table)

Table 4: Risk‑calibrated Ftarget targeting P(SSB < 0.2 B0) = 0.05 across both scenarios.

Metric	Value
Ftarget	0.596
Combined P(SSB < 0.2 B0)	0.050
Low steepness P(SSB < 0.2 B0)	0.080
High steepness P(SSB < 0.2 B0)	0.020

4.3 Monte Carlo performance

n_sims <- 100
n_years <- 50

ssb_low <- matrix(NA, n_years, n_sims)
ssb_high <- matrix(NA, n_years, n_sims)
catch_low <- matrix(NA, n_years, n_sims)
catch_high <- matrix(NA, n_years, n_sims)

for (i in 1:n_sims) {
  sim_l <- run_simulation(scenario_low$params, F_series = ref_low$Fmsy, seed = i)
  sim_h <- run_simulation(scenario_high$params, F_series = ref_high$Fmsy, seed = i)
  ssb_low[, i] <- sim_l$SSB
  ssb_high[, i] <- sim_h$SSB
  catch_low[, i] <- sim_l$Catch
  catch_high[, i] <- sim_h$Catch
}

ssb_low_df <- data.frame(
  Year = rep(1:n_years, times = n_sims),
  Sim = rep(1:n_sims, each = n_years),
  SSB = as.vector(ssb_low) / 1e9,
  Scenario = "h=0.6"
)

ssb_high_df <- data.frame(
  Year = rep(1:n_years, times = n_sims),
  Sim = rep(1:n_sims, each = n_years),
  SSB = as.vector(ssb_high) / 1e9,
  Scenario = "h=0.9"
)

ssb_df <- rbind(ssb_low_df, ssb_high_df)

median_df <- rbind(
  data.frame(Year = 1:n_years, SSB = apply(ssb_low, 1, median) / 1e9, Scenario = "h=0.6"),
  data.frame(Year = 1:n_years, SSB = apply(ssb_high, 1, median) / 1e9, Scenario = "h=0.9")
)

bmsy_df <- data.frame(
  Scenario = c("h=0.6", "h=0.9"),
  Bmsy = c(ref_low$Bmsy, ref_high$Bmsy) / 1e9
)

scenario_colors <- c("h=0.6" = "blue", "h=0.9" = "red")

p_mc <- ggplot(ssb_df, aes(x = Year, y = SSB, group = interaction(Scenario, Sim), color = Scenario)) +
  geom_line(alpha = 0.1, size = 0.4, show.legend = FALSE) +
  geom_line(
    data = median_df,
    aes(x = Year, y = SSB, color = Scenario),
    size = 1.1,
    inherit.aes = FALSE
  ) +
  geom_hline(
    data = bmsy_df,
    aes(yintercept = Bmsy, color = Scenario),
    linetype = "dashed",
    size = 0.8,
    inherit.aes = FALSE
  ) +
  facet_wrap(~ Scenario, ncol = 2, scales = "free_y") +
  scale_color_manual(values = scenario_colors) +
  labs(x = "Year", y = "Female SSB (Mt)", title = "Monte Carlo SSB trajectories") +
  ggthemes::theme_few() +
  guides(color = "none")

p_mc

Figure 2: Ensemble of 100 simulations for each scenario showing SSB variability. Thin lines are individual simulations, thick solid lines are medians, and dashed lines indicate Bmsy.

5 Discussion

A single risk‑calibrated HCR can maintain biomass protection across low and high productivity regimes without relying on environmental covariates. This reduces dependence on uncertain climate‑productivity linkages while preserving a transparent decision rule. The approach does not replace full‑feedback MSE, but it provides a practical, defensible rule that can be communicated and implemented within current tier‑based frameworks.

5.1 Limitations

This analysis omits observation error, assessment feedback, and fleet behavior, so realized management performance may differ from the idealized operating‑model results. The population dynamics are stylized and do not represent spatial structure, ecosystem interactions, or time‑varying selectivity. Future work should embed the risk‑calibrated HCR within a full‑feedback MSE framework and test robustness to assessment bias, implementation error, and alternative operating models.

6 Conclusions

Risk calibration of a transparent SPR‑based HCR provides a simple method to manage productivity uncertainty. In this pollock‑like case, a single Ftarget achieves consistent biomass protection across steepness regimes while preserving yield in more productive states.

7 Data and code availability

All code and inputs for the simulations and manuscript are available in the public repository: https://github.com/jimianelli/HCR_paper

References

Hollowed, A. B., J. N. Ianelli, P. A. Livingston, B. Planque, and W. S. Wooster. 2020. “Changed States of the North Pacific Ecosystem.” ICES Journal of Marine Science 77: 886–97.

North Pacific Fishery Management Council. 2024. Fishery Management Plan for Groundfish of the Bering Sea and Aleutian Islands Management Area. North Pacific Fishery Management Council.

Punt, A. E., C. S. Szuwalski, J. N. Ianelli, and A. B. Hollowed. 2024. “Evaluating Harvest Control Rules Under Ecosystem Change.” ICES Journal of Marine Science 81 (1): 45–62.

Szuwalski, C. S., J. N. Ianelli, A. E. Punt, and P. D. Spencer. 2023. “Climate-Informed Recruitment Dynamics: Implications for Harvest Control Rules.” Fisheries Research 265: 106745.