[]

Introduction

This tutorial is based on the work carried out within Bay of Biscay (Western Waters) Case study of SeaWise – Work Package 6: “Evaluation of fisheries management strategies in an Ecosystem context”. This tutorial demonstrates how to simulate temporal closures and annual effort reduction using FLBEIA modeling framework (Garcia et al. 2017). Originally the code was developed to simulate the temporal closure in the Bay of Biscay to protect Dolphins that was implemented in the region in 2024. In the original work a more detailed fleet-metier description was used, in the tutorial fleet data has been aggregated at gear level.

Objective

This tutorial demonstrates how to simulate a temporal closure scenario and a limitation on the effort of trawlers using FLBEIA modeling framework. It guides users through the conditioning and configuration of the model and analysis of the results.

You will learn to simulate two scenarios:

Temporal closure: The implementation of a temporary ban on fishing activity is based on a policy measure in force since 2024 in the Bay of Biscay, where most fisheries undergo a one-month closure (22 January to 20 February) to reduce common dolphin bycatch (EU 2024; ICES 2023a).
Annual effort reduction: This hypothetical scenario simulates a 50% reduction in the annual fishing effort of trawlers, to mitigate ecosystem impacts.

By the end of this tutorial, you will be able to modify FLBEIA configuration to reflect changes in fleet behavior derived from these two management actions, run your own scenario-based simulations, and interpret the outcomes.

Case study overview

The tutorial is focused on demersal mixed-fisheries in the Bay of Biscay embedded in the Western Waters case study. More details about the case study are available on the SeaWise Project wedpage and deliverables.

Model overview:

Time coverage: 2000-2021 (history), 2022-2050 (projection).
Time step: Annual.
Uncertainty: Introduced through Monte Carlo simulation, 100 model replicates that include variability in both stock dynamics and fleet behavior during the projection period.

Operating model:

Stocks dynamics: 16 stocks, 9 age-structured and 7 in biomass.

All stocks with analytical assessment (ICES category 1) - except nephrops - are incorporated using an age-structured approach (Figure 1), conditioned with the ICES assessment model carried out in 2022 using data up to 2021. For Northern hake, a density-dependent individual growth model is applied, with recruitment influenced by chlorophyll levels to reflect the impact of global change.
Stocks with no analytical assessment (ICES Category 3) and nephrops are modeled using a constant CPUE (Catch Per Unit Effort) approach (Figure 1). Population biomass is assumed constant over time and the catch production is proportional to effort.

Figure 1: Stocks included in the FLBEIA.

Fleet dynamics: 3 fleet-metier segments with a single metier (MIS).

To avoid confidentiality issues with metier level data, the original set of 33 metier segments (ICES 2023b) were aggregated into three fleet segments distinguished by the gear used: TRW (trawlers), GNS (gillnetters), and LLS (longliners).
For each stock one ghost fleet is included to account for the catches occurring outside the Bay of Biscay.

Covariates:

Economic covariates are included to calculate economic indicators, the data are obtained form the STECF Annual Economic Report (STECF 2022). Chlorophyll data to simulate hake recruitment was taken from NEMO-MEDUSA for RCP 8.5 (Yool, Popova, and Coward 2015).

Management Procedure

A short cut approach is used in the management procedure, in the observation model biological and catch data are observed perfectly. No assessment model is used, and population numbers and exploitation level at age is taken from the operating model.

The TAC is simulated using the ICES Harvest Control Rule in the case of age structured stocks and for biomass-agregated stocks the TAC in the whole projection is uniformly distributed with mean equal to the last observed TAC and a coefficient of variation of 25%.

In all the scenarios landing obligation is fully implemented and hence fleets stop fishing when the first of the quotas, of a predefined set of stocks (those subject to landing obligation), is consumed.

Required packages

You need the following R packages:

CRAN: ggplot2, dplyr, reshape2, tidyr
FLR: FLCore, FLAssess, FLash, FLBEIA, FLFleet, ggplotFL

Use 64-bit R on Windows, as 32-bit may cause issues.

# Install packages
install.packages( c("ggplot2","dplyr", "reshape2", "tidyr"))
install.packages( c("FLCore", "FLFleet", "FLBEIA", "ggplotFL", 
                    "FLash", "FLAssess"), 
                  repos="http://flr-project.org/R")
remotes::install_github("flr/FLBEIAshiny")

Important: Install FLCore, FLFleets and FLBEIA in this order to avoid compatibility issues.

# Load packages
library(FLBEIA)
library(FLBEIAshiny)
library(ggplot2)
library(dplyr)
library(reshape2)
library(tidyr)

Load base case data and settings

# Clear environment
rm(list=ls())

# Set working directory
myWD <- "C:/USE/GitHub/seawise/demersal_upd22"
#myWD <- "C:/use/OneDrive - AZTI/0_Proiektuak/SeaWise/WP6/Course"

setwd(myWD)

# Start time
init <- Sys.time() 

# Base case data
load('./data/FLR_Objs_tutorial_100iter.RData') 
load('./data/ctrls_tutorial_100iter.RData') 
load('./data/monthly_catch_effort_tutorial.RData')

# Other settings
load('./data/settings_tutorial_100iter.RData') 
stknms <- c(dyn.stk, cnt.stk) 
stks.ntac <- setdiff(stknms, stks.tac)
fltnm <- names(fleets)[which(substr(names(fleets),1,2) != 'OT')]

Set the names and the location of the files to save the results

FLR_Objs        <- './output/out_FLR_Objs_tutorial_100iter.RData'
ctrls           <- './output/out_ctrls_tutorial_100iter.RData'
out_file        <- './output/out_results_tutorial_100iter.RData' 
out_sum_file    <- './output/out_results_sum_tutorial_100iter.RData' 
out_sum_file_qs <- './output/out_results_sumQ_tutorial_100iter.RData'

The conditioning of the scenarios

The following sections detail how to condition FLBEIA to simulate the temporal closure and the effort reduction in some of the fleets.

Temporal closure

In this scenario, we assume that the activity of trawlers and gillnetters is subject to a one-month closure measure, while longliners continue operating uninterrupted, unaffected by the closure measure. This implementation of FLBEIA is annual, hence the closure is simulated implicitly through a change in catchability and reduction in the maximum effort (capacity) the fleet can perform along the year. Seasonal catch per unit of effort (CPUE) data are analysed outside the model to re-estimate catchability after the one month closure.

Plot: monthly CPUE for impacted fleets

The CPUE trends, based on catch and effort data from the last three historical years, are in Figure 2 and Figure 3. These plots provide insight into monthly fleet performance, showing the productivity of the fleet along the year. This information is useful for understanding if the fishery will be able to compensate the catches lost in February (the month when the closure happens) fishing more in the rest of the year.

For example in trawlers (Figure 2), the CPUE for horse mackerel, mackerel and blue whiting is higher during the first quarter of the year, when the stocks migrate to Bay of Biscay. As a result, the February closure is likely to negatively affect their catches, because the trawlers will need more effort in later months than at the beginning of the year to obtain the same catches, or even with more effort it will not be possible to reach the February level of catches.

cls_fleets  <- c("TRW","GNS")

# Calculate monthly catch-to-effort ratio for the impacted fleets (TRW, GNS)
cat_eff <- monthly_catch %>%
  filter(fleet %in% cls_fleets) %>%
  left_join(monthly_effort, by = c("year", "month", "fleet", "metier")) 

cat_eff_rc <- cat_eff %>% group_by(month, fleet, metier, stock) %>%
  summarise(rC = sum(catch) / sum(effort), .groups = "drop")

# Plot catch/effort by stock, month and fleet
plotCPUE <- function(data, fleet_name) {
  ggplot(filter(data, fleet == fleet_name), aes(stock, rC)) +
    geom_bar(stat = "identity") +
    facet_wrap(~month) +
    labs(title = fleet_name, y = "CPUE") +
    theme_bw() +
    theme(
      legend.position = "none",
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
    )
  }

plotCPUE(cat_eff_rc, "TRW")

plotCPUE(cat_eff_rc, "GNS")

Plot: monthly catch share

The monthly catch share -the proportion of total annual catch per month- calculated over the last three years are in Figure 4 and Figure 5. This figure reveals when catches are most concentrated for each fleet and stock, and helps evaluate how the February closure measure may impact in fleets depending on the stock.

We can see that catches of horse mackerel, mackerel and blue whiting are concentrated in the first quarter of the year. Species such as black and white anglerfish, hake, megrim and pollack show approximately 25% higher catches during the first four month of the year compared to the remainder months of the year. In contrast, seabass and sole have lower catch volumes in the early months, with catches increasing later in the year. Meanwhile, species like neprhops, the rays (RJC, RJN and RJU) and smooth-hounds and whiting have a uniform catch distribution throughout the year.

So for example, as commented earlier, in the case of widely distributed pelagic stocks (horse mackerel, mackerel and blue whiting), as the catches are predominantly concentrated in February, a closure during this critical period necessitates significantly higher catch levels in other months to maintain annual catches.

mean_monthly_share <- cat_eff %>%
  group_by(month, fleet, stock) %>%
  summarise(monthly_catch = sum(catch), .groups = "drop") %>%
  group_by(fleet, stock) %>%
  mutate(annual_catch = sum(monthly_catch)) %>%
  ungroup() %>%
  mutate(catch_share = monthly_catch / annual_catch)

plotCatchShare <- function(data, fleet_name) {
  ggplot(filter(data, fleet == fleet_name), aes(x = factor(month), y = catch_share)) +
    geom_bar(stat = "identity", position = "dodge") +
    facet_wrap(~ stock) +
    labs(x = "Month", y = "Annual catch proportion") +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
}

plotCatchShare(mean_monthly_share, "TRW")

plotCatchShare(mean_monthly_share, "GNS")

Calculate annual catchability with and without temporal closure

In this section, the annual CPUE for each fleet and stock with and without temporal closure are calculated. To do the calculations without February the effort and catch in that month are set to zero.

We then compute the relative change in CPUE (alpha = closure CPUE / original CPUE), which is later used to adjust catchability.

# Loops through each fleet and stock to compute average CPUE:
# - `cpue_org`: original CPUE: full-year baseline CPUE
# - `cpue_cls`: closure CPUE : CPUE after simulating February closure
# - `alpha`: relative change ratio (used to adjust catchability later)

q.cls.dat <- list()

for (fl in cls_fleets) {
  df_fl <- cat_eff %>% filter(fleet == fl)
  
  for (mt in unique(df_fl$metier)) {
    df_mt <- df_fl %>% filter(metier == mt)
    
    for (sp in unique(df_mt$stock)) {
      df_sp <- df_mt %>% filter(stock == sp)
      
      # whole year
      cpue_org <- df_sp %>%
        group_by(year) %>%
        summarise(cpue_org = sum(catch) / sum(effort), .groups = "drop") %>%
        summarise(cpue_org = mean(cpue_org)) %>%
        pull(cpue_org)
      
      # removing February (month == '2')
      cpue_cls <- df_sp %>%
        mutate(effort = if_else(month == "2", 0, effort),
               catch  = if_else(month == "2", 0, catch)) %>%
        group_by(year) %>%
        summarise(cpue_cls = sum(catch) / sum(effort), .groups = "drop") %>%
        summarise(cpue_cls = mean(cpue_cls)) %>%
        pull(cpue_cls)
      
      q.cls.dat[[length(q.cls.dat)+1]] <- data.frame(
        fleet = fl, metier = mt, stock = sp, alpha = cpue_cls / cpue_org
      )
    }
  }
}

q.cls.dat <- bind_rows(q.cls.dat) %>% drop_na()

Modify fleet capacity and catchability

Here, we adjust the fleet object to reflect change in catchability derived from the February closure. Specifically, capacity is decreased by 1/12 to represent the one-month closure, and catchability is updated using previously calculated alpha values.

# Clone original data to preserve base case
fleets.cl <- fleets

# Reduce fleet capacity by 1/12 to reflect one-month closure
# Adjust catch.q (catchability) using alpha computed earlier
for (fl in cls_fleets) {
  fleets.cl[[fl]]@capacity@.Data <- fleets.cl[[fl]]@capacity@.Data * 11/12
  new_q.fl <- subset(q.cls.dat, fleet == fl)
  for (mt in unique(new_q.fl$metier)) {
    new_q.mt <- subset(new_q.fl, metier == mt)
    for (sp in unique(new_q.mt$stock)) {
      alpha <- subset(new_q.mt, stock == sp)$alpha
      fleets.cl[[fl]]@metiers[[mt]]@catches[[sp]]@catch.q[, ac((yrd + 1):yrf)] <-
        fleets.cl[[fl]]@metiers[[mt]]@catches[[sp]]@catch.q[, ac((yrd + 1):yrf)] * alpha
    }
  }
}

Plot: alpha values

In Figure 6 alpha values are ploted by stock. An alpha value below 1 indicates a reduction in catchability, while a value above 1 suggests an increase in catchability. An alpha equal to 1 implies that catchability remained unchanged after closure. This figure highlights the stocks most affected in each fleet by the management action.

For example, as inferred from Figure 2 and Figure 4, the catchabilities of horse mackerel, mackerel, and blue whiting are reduced which means the catch produced per unit of effort will be lower.

q.cls.dat$fleet    <- factor(q.cls.dat$fleet, levels = c("TRW", "GNS"))

ggplot(q.cls.dat, aes(x = stock, y = alpha)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
  geom_segment(aes(x = stock, xend = stock, y = 1, yend = alpha), color = "gray") +
  geom_point(size = 3) +
  facet_wrap(~ fleet) +
  labs(
    x = "Stock",
    y = "Alpha (CPUE after closure / original CPUE)"
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1),
    legend.position = "bottom"
  )

Figure 6: Alpha (CPUE after closure / original CPUE)

Annual effort reduction

In this scenario, the activity of trawlers is reduced by 50%. Thus, the capacity of the fleet was reduced accordingly to limit the effort that the fleet was able to perform annually. This would potentially result in a reduction in the use of quota-share of each of the stock. As the aim of this management measure is not to reduce the use of fishing opportunities but just to reduce the trawling activity to limit its impact in the ecosystem, the likely unused fishing opportunities are distributed among the rest of the gears. The distribution is done proportionally to their current quota shares. But, how do we calculate the potential decrease in the use of quota-share? In an scenario with constant effort, we project the system until equilibrium reducing the total effort of trawlers in 50%, while keeping the effort of the rest of the fleets equal to the mean of last three historical years, and calculate the catch share of the trawlers in equilibrium. Then, as fleet dynamics are constrained by quota shares in this implementation of FLBEIA we update the quota share in the model with this new value.

Adjust capacity, effort, and number of vessels for trawlers

To simulate reduced activity, we reduced by half the effort of the fleet. The number of vessels is recalculated based on the adjusted effort and maximum allowable fishing days, and the fleet capacity is updated based on those.

# Clone original data to preserve base case
fleets.50 <- fleets
covars.50 <- covars
fleets.ctrl.fixedfixedEff <- fleets.ctrl.fixedEff
advice.50 <- advice

# Apply 50% reduction to TRW fleet effort
fleets.50[["TRW"]]@effort@.Data[1, ac((yrd+1):yrf), 1, 1, 1, 1:Nit] <- fleets.50[["TRW"]]@effort@.Data[1, ac((yrd+1):yrf), 1, 1, 1, 1:Nit]*0.5

# Update number of vessels based on effort and max days
covars.50[["NumbVessels"]]@.Data["TRW", ac((yrd+1):yrf), 1, 1, 1, 1:Nit]  <- ceiling(fleets.50[["TRW"]]@effort@.Data[1, ac((yrd+1):yrf), 1, 1, 1, 1:Nit]/covars.50 [["MaxDays"]]@.Data["TRW", ac((yrd+1):yrf), 1, 1, 1, 1:Nit])

# Recalculate fleet capacity based on new number of vessels and max days
fleets.50[["TRW"]]@capacity@.Data[1, ac((yrd+1):yrf), 1, 1, 1, 1:Nit] <- covars.50[["NumbVessels"]]@.Data["TRW", ac((yrd+1):yrf), 1, 1, 1, 1:Nit] * covars.50 [["MaxDays"]]@.Data["TRW", ac((yrd+1):yrf), 1, 1, 1, 1:Nit]

Run simulations with fixed effort

FLBEIA was runned under two fixed effort scenarios to calculate the excedent of quota share that could be shared with the rest of the gears once the effort of the trawlers is limited by a 50% reduction. Scenarios: (1) historical fixed effort, and (2) historical fixed effort reduced by 50%.

# Set effort model to fixed effort for all fleets
for(fl in names(fleets.ctrl.fixedfixedEff)[-c(1,2)]){
  fleets.ctrl.fixedfixedEff[[fl]]$effort.model  <- "fixedEffort"
  fleets.ctrl.fixedfixedEff[[fl]][["effort.restr"]] <- "max"}

# Run simulations with and without 50% TRW reduction
fix.qs    <- FLBEIA(biols = biols, SRs = SRs, BDs = NULL, fleets = fleets, covars = covars,
                    indices = NULL, advice = advice, main.ctrl = main.ctrl, biols.ctrl = biols.ctrl, fleets.ctrl = fleets.ctrl.fixedfixedEff, covars.ctrl = covars.ctrl, obs.ctrl = obs.ctrl, assess.ctrl = assess.ctrl, advice.ctrl = advice.ctrl.msy)

fix.qs.50 <- FLBEIA(biols = biols, SRs = SRs, BDs = NULL, fleets = fleets.50, covars = covars.50,
                    indices = NULL, advice = advice, main.ctrl = main.ctrl, biols.ctrl = biols.ctrl, fleets.ctrl = fleets.ctrl.fixedfixedEff, covars.ctrl = covars.ctrl, obs.ctrl = obs.ctrl, assess.ctrl = assess.ctrl, advice.ctrl = advice.ctrl.msy)

# Combine results
fltStk  <- rbind((fltStkSum(fix.qs, scenario = 'fix.qs', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                 (fltStkSum(fix.qs.50, scenario = 'fix.qs.50', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)))

Calculate new quota shares after annual effort reduction

Using catch data from the previous two simulations, we calculate the quota shares for each stock across fleets with and without reduction in the effort of trawlers.

Qshare.data <- fltStk %>%
  filter(year %in% (yrf-20):yrf, indicator == "catch") %>%
  select(-indicator)

# Create an empty data frame
df_template <- data.frame(
  fleet = rep(names(fleets), times = Nit),
  iter = rep(1:Nit, each = length(names(fleets))),
  scenario = NA,
  matrix(NA, nrow = length(names(fleets)) * Nit, ncol = length(stknms))
)
colnames(df_template)[-(1:3)] <- stknms

# Loop through scenarios
qshare.trwl <- data.frame()

for (sc in unique(Qshare.data$scenario)) {
  data.sc <- Qshare.data %>% filter(scenario == sc)
  df.sc <- df_template
  df.sc$scenario <- sc
  
  for (it in 1:Nit) {
    for (fl in names(fleets)) {
      for (st in stknms) {
        # Filter relevant data
        data.fl.st <- data.sc %>% filter(fleet == fl, stock == st, iter == it)
        data.st <- data.sc %>% filter(stock == st, iter == it)
        
        if (nrow(data.fl.st) > 0 && nrow(data.st) > 0) {
          # Calculate yearly ratios
          ratios <- sapply(unique(data.fl.st$year), function(yr) {
            val.fl <- sum(data.fl.st$value[data.fl.st$year == yr])
            val.st <- sum(data.st$value[data.st$year == yr])
            if (val.st > 0) val.fl / val.st else NA
          })
          # Assign mean ratio
          df.sc[df.sc$fleet == fl & df.sc$iter == it, st] <- mean(ratios, na.rm = TRUE)
        }
      }
    }
  }
  # Append scenario result
  qshare.trwl <- bind_rows(qshare.trwl, df.sc)
}

# Replace NA with 0
qshare.trwl[is.na(qshare.trwl)] <- 0

Adjust advice object for new quota shares

Comparing the original and new quota shares, we calculate the quota share that originally belonged to trawlers and with the reduction in effort, in order to avoid losing fishing opportunities, it can be shared among other fleets. The remnant quota from trawlers is redistributed across other gears based on their historical quota shares. We create a new advice.50 object equal to the original advice object and update it with the new quota shares.

# Identify TRW and non-TRW fleets
trw_flt   <- "TRW"
notrw_flt <- setdiff(fltnm, trw_flt)

# Calculate difference in quota share due to reduction
qshare.sc  <- filter(qshare.trwl, scenario == "fix.qs.50")
qshare.fix <- filter(qshare.trwl, scenario == "fix.qs")

# Redistributed TRW quota share values
redistributed_qs <- data.frame(stock = stknms, qs = NA, stringsAsFactors = FALSE)

for (it in 1:Nit){
  for (st in stknms) {
    trw_diff <- qshare.fix[qshare.fix$fleet == trw_flt & qshare.fix$iter == it, st] - qshare.sc[qshare.sc$fleet == trw_flt & qshare.sc$iter == it, st]
    if (it == 1) {
      redistributed_qs[redistributed_qs$stock == st, "qs"] <- trw_diff
    }
    for (fl in fltnm) {
      if (st %in% catchNames(fleets.50[[fl]])) {
        if (fl == trw_flt) {
          advice.50$quota.share[[st]][fl, ac((yrd+1):yrf),1,1,1,it] <- advice$quota.share[[st]][fl, ac(yrf),1,1,1,it] - trw_diff
        } else { # Reallocate TRW quota among other fleets proportionally
          share_total <- sum(advice$quota.share[[st]][dimnames(advice$quota.share[[st]])$fleet %in% notrw_flt, ac(yrf),1,1,1,it])
          share_fl    <- sum(advice$quota.share[[st]][fl, ac(yrf),1,1,1,it])
          redistributed <- trw_diff * (share_fl / share_total) + share_fl
          advice.50$quota.share[[st]][fl, ac((yrd+1):yrf),1,1,1,it] <- redistributed
        }
      } else {
        advice.50$quota.share[[st]][fl,ac((yrd+1):yrf),1,1,1,it] <- 0
      }
    }
  }
}

# # Check for quota balance
# lapply(advice.50$quota.share, quantSums)

Plot: excess quota in trawler fleet

In Figure 7, we show the excess of quota generated by the trawlers as a result of effort reduction. This quota surplus is redistributed to other fleets proportionally, based on their historical catch records.

ggplot(redistributed_qs, aes(x = stock, y = qs)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    x = "Stock",
    y = "Excess quota share of trawler"
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
    legend.position = "bottom"
  )

Figure 7: Trawler quota share excess after its effort reduction

Plot: fleet catch contribution

The following two plots (Figure 8 and Figure 9) provide the relative catch contribution of different fleets across stocks. The first plot focuses on the three selected fleets (gillnetters, longliners and trawlers, longliners, trawlers) for all the stocks, while the second includes all fleets but only for the stocks with proper population dynamics.

These plots are useful to understand how the excess quota from trawlers is redistributed among other fleets. For example, the excess quota of sole from trawlers is primarily reallocated to gillnetters because the catches of longliners represent a very really low proportions of the total catches, whereas the excess quota of horse mackerel is mostly transferred to longliners (Figure 8). In the case of neprhops, as only trawlers exploit it the excess quota is lost.

# Considering the three fleets (GNS, LLS, TRW) for all the stocks
flt_stk <- monthly_catch %>% 
  filter (fleet %in% fltnm, stock %in% stknms) %>%
  group_by(fleet, stock) %>%
  summarise(annual_fl_catch = sum(catch), .groups = "drop") %>%
  group_by(stock) %>%
  mutate(annual_tot_catch = sum(annual_fl_catch)) %>%
  ungroup() %>%
  mutate(catch_prc = annual_fl_catch / annual_tot_catch) %>%
  group_by(fleet, stock) %>%
  summarise(fl_tot_catch = mean(catch_prc), .groups = "drop")

ggplot(flt_stk, aes(x = factor(fleet), y = fl_tot_catch)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ stock) +
  labs(x = "fleet", y = "Proportion of catch") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Figure 8: Catch Share by Fleet and Stock

# Considering all the fleets for dynamic stocks
allflt_dynstk <- monthly_catch %>%
  filter (stock %in% dyn.stk) %>%
  group_by(fleet, stock) %>%
  summarise(annual_fl_catch = sum(catch), .groups = "drop") %>%
  group_by(stock) %>%
  mutate(annual_tot_catch = sum(annual_fl_catch)) %>%
  ungroup() %>%
  mutate(catch_prc = annual_fl_catch / annual_tot_catch) %>%
  group_by(fleet, stock) %>%
  summarise(fl_tot_catch = mean(catch_prc), .groups = "drop")

ggplot(allflt_dynstk, aes(x = factor(fleet), y = fl_tot_catch)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ stock) +
  labs(x = "fleet", y = "Proportion of total catch") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Figure 9: Fleet contribution to total catch of each stock

Run the simulations

Additionally to the temporal closure and effort reduction scenario a base case scenario was run in order to compare the results of the management measures with the current situation. In this base case scenario MSY advice for all stocks with real population dynamics was combined with landing obligation. The other two scenarios were simulated on top of this base case.

# Base case
bc    <- FLBEIA(biols = biols, SRs = SRs, BDs = NULL, fleets = fleets, covars = covars, indices = NULL, advice = advice, main.ctrl = main.ctrl, biols.ctrl = biols.ctrl, fleets.ctrl = fleets.ctrl.min, covars.ctrl = covars.ctrl, obs.ctrl = obs.ctrl, assess.ctrl = assess.ctrl, advice.ctrl = advice.ctrl.msy)

# Temporal closure
sc <- FLBEIA(biols = biols, SRs = SRs, BDs = NULL, fleets = fleets.cl, covars = covars, indices = NULL, advice = advice, main.ctrl = main.ctrl, biols.ctrl = biols.ctrl, fleets.ctrl = fleets.ctrl.min, covars.ctrl = covars.ctrl, obs.ctrl = obs.ctrl, assess.ctrl = assess.ctrl, advice.ctrl = advice.ctrl.msy)

# Annual Effort Reduction
er <- FLBEIA(biols = biols, SRs = SRs, BDs = NULL, fleets = fleets.50, covars = covars.50, indices = NULL, advice = advice.50, main.ctrl = main.ctrl, biols.ctrl = biols.ctrl, fleets.ctrl = fleets.ctrl.min, covars.ctrl = covars.ctrl, obs.ctrl = obs.ctrl, assess.ctrl = assess.ctrl, advice.ctrl = advice.ctrl.msy)

Make summaries

By iteration

brp   <- data.frame(stock = rep(names(biols)[1:9], each = 1), iter = rep(1:Nit, 9),
                    Bpa  = rep(sapply(names(biols)[1:9], function(x) advice.ctrl.msy[[x]]$ref.pts['Btrigger',1]), each = 1),
                    Blim = rep(sapply(names(biols)[1:9], function(x) advice.ctrl.msy[[x]]$ref.pts['Blim',1]), each = 1),
                    Ftarget = rep(sapply(names(biols)[1:9], function(x) advice.ctrl.msy[[x]]$ref.pts['Fmsy',1]), each = 1),
                    Fmsy = rep(sapply(names(biols)[1:9], function(x) advice.ctrl.msy[[x]]$ref.pts['Fmsy',1]), each = 1),
                    Btarget = NA,
                    Fpa = NA,
                    Flim = NA,
                    Fupp = rep(sapply(names(biols)[1:9], function(x) base::ifelse("Fupp" %in% rownames(advice.ctrl.msy[[x]]$ref.pts), advice.ctrl.msy[[x]]$ref.pts['Fupp',1], NA)), each = 1),
                    Flow = rep(sapply(names(biols)[1:9], function(x) base::ifelse("Flow" %in% rownames(advice.ctrl.msy[[x]]$ref.pts), advice.ctrl.msy[[x]]$ref.pts['Flow',1], NA)), each = 1))

bio   <-  rbind((bioSum(bc, scenario = 'bc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE, brp= brp)),
                (bioSum(sc, scenario = 'sc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE, brp= brp)),
                (bioSum(er, scenario = 'er', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE, brp= brp)))

flt    <- rbind((fltSum(bc, scenario = 'bc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (fltSum(sc, scenario = 'sc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (fltSum(er, scenario = 'er', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)))

fltStk <- rbind((fltStkSum(bc, scenario = 'bc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (fltStkSum(sc, scenario = 'sc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (fltStkSum(er, scenario = 'er', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)))

mt     <- rbind((mtSum(bc, scenario = 'bc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (mtSum(sc, scenario = 'sc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (mtSum(er, scenario = 'er', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)))

mtStk  <- rbind((mtStkSum(bc, scenario = 'bc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (mtStkSum(sc, scenario = 'sc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (mtStkSum(er, scenario = 'er', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)))

adv    <- rbind((advSum(bc, scenario = 'bc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (advSum(sc, scenario = 'sc', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)),
                (advSum(er, scenario = 'er', years = ac(yr0:main.ctrl[[1]][2]), long = TRUE)))

Calculte quantiles

BIOQ     <- bioSumQ(bio, prob = c(0.05, 0.5, 0.95))
FLTQ     <- fltSumQ(flt, prob = c(0.05, 0.5, 0.95))
FLTSTKQ  <- fltStkSumQ(fltStk, prob = c(0.05, 0.5, 0.95))
MTQ      <- mtSumQ(mt, prob = c(0.05, 0.5, 0.95))
MTSTKQ   <- mtStkSumQ(mtStk, prob = c(0.05, 0.5, 0.95))
ADVQ     <- advSumQ(adv, prob = c(0.05, 0.5, 0.95))

Save updated objects, results and summaries

save(biols, advice, fleets, SRs, covars,
     fleets.cl, q.cls.dat,
     advice.50, fleets.50, covars.50, 
     fix.qs, fix.qs.50, qshare.trwl,
     file = FLR_Objs)

save(main.ctrl, biols.ctrl, covars.ctrl, obs.ctrl, assess.ctrl,
     fleets.ctrl.fixedEff, fleets.ctrl.min, fleets.ctrl.prev, fleets.ctrl.fixedfixedEff,
     advice.ctrl.msy,
     file = ctrls)

save(bc, sc, er,
     file = out_file)

save(brp, bio, flt, fltStk, mt, mtStk, adv,
     file = out_sum_file)

save(BIOQ, FLTQ, FLTSTKQ, MTQ, MTSTKQ, ADVQ, file =  out_sum_file_qs)

Results

The SSB for the stocks showed minimal differences across the three management scenarios (Figure 10). The SSB of horse mackerel, mackerel, and blue whiting were only marginally affected by changes in the management of demersal mixed-fisheries (Figure 10), as these stocks are primarily caught outside the Bay of Biscay (Figure 9).

Across all fleets, the temporal closure scenario resulted in catch levels similar to the base case, though with some variations (Figure 12, Figure 13, and Figure 14). However, effort dynamics varied significantly between fleets under these two scenarios (Figure 11). In the trawl fleet, the temporal closure scenario led to higher overall effort than the base case. Horse mackere, mackerel, and blue whiting are mainly caught during the first quarter of the year (Figure 4), when CPUE is also highest (Figure 2). With February closed to fishing, trawlers had to increase effort in the remaining months to maintain catch levels, resulting in an overall rise in effort (Figure 11). In contrast, the impact of temporal closure on gillnetters was minimal. In this fleet catch is more evenly distributed throughout the year (Figure 5), and monthly CPUE differences are less pronounced than in trawlers (Figure 3). Consonantly, effort remained relatively stable between the temporal closure and base case scenarios in GNS (Figure 11).

On the other hand, the Annual effort reduction scenario (50% reduction) produced clear differences in both catch and effort across all fleets. As expected, in the trawl fleet, catch and effort decreased by approximately 50% compared to the base case (Figure 12 and Figure 11). Meanwhile, in the gillnet and longline fleets, the reduction in trawl effort led to an increase in catch across all stocks (Figure 12), along with a rise in overall effort (Figure 11). This reflects a redistribution of fishing pressure from TRWs to these fleets.

add_scenario_name <- function(df, code_column) {
  df$scenario_name <- ifelse(df[[code_column]] == "bc", "Base Case",
                      ifelse(df[[code_column]] == "sc", "Temporal Closure",
                      ifelse(df[[code_column]] == "er", "Annual Effort Reduction", NA)))
  return(df)
}

BIOQ    <- add_scenario_name(BIOQ, "scenario")
FLTQ    <- add_scenario_name(FLTQ, "scenario")
FLTSTKQ <- add_scenario_name(FLTSTKQ, "scenario")


# Color palette for scenarios
scnms  <- unique(BIOQ$scenario_name)
pal.sc <- setNames(c("#A6CEE3", "#B2DF8A", "#FDBF6F"), scnms)

# Reordering factor levels
BIOQ$scenario_name    <- factor(BIOQ$scenario_name, levels = c("Base Case", "Temporal Closure", "Annual Effort Reduction"))
FLTQ$scenario_name    <- factor(FLTQ$scenario_name, levels = c("Base Case", "Temporal Closure", "Annual Effort Reduction"))
FLTSTKQ$scenario_name <- factor(FLTSTKQ$scenario_name, levels = c("Base Case", "Temporal Closure", "Annual Effort Reduction"))

ref.pts  <- brp %>% select(-iter, -Btarget, -Fpa, -Flim) %>% distinct()

BIOQ_ssb <- BIOQ %>%
  filter(stock %in% dyn.stk, indicator == "ssb")

ggplot(BIOQ_ssb, aes(x = year, y = q50, color = scenario_name, fill = scenario_name)) +
  geom_line(size = 1.2) +
  facet_wrap(~stock, scales = "free_y", ncol = 3) +
  scale_color_manual(values = pal.sc) +
  scale_fill_manual(values = pal.sc) +
  geom_hline(data = ref.pts, aes(yintercept = Blim), linetype = "dashed", color = "red") +
  geom_hline(data = ref.pts, aes(yintercept = Bpa), linetype = "dashed", color = "gold1") +
  theme_bw(base_size = 12) +
  theme(strip.text = element_text(size = 12), legend.position = "bottom") +
  labs(y = "SSB (tons)", x = NULL, fill = NULL, color = NULL)

Figure 10: Spawning stock biomass time series for the modelled stocks (Figure 1) by scenario. Horizontal dashed red and yellow lines correspond to Blim and Bpa, respectively.

FLTQ_eff <- FLTQ %>%
  filter(year > 2021, fleet %in% fltnm, indicator == "effort")

FLTQ_eff$fleet    <- factor(FLTQ_eff$fleet, levels = c("TRW", "GNS", "LLS"))

ggplot(FLTQ_eff, aes(x = year, y = q50, color = scenario_name, fill = scenario_name)) +
  geom_line(size = 1.2) +
  facet_wrap(~fleet, scales = "free_y", ncol = 3) +
  scale_color_manual(values = pal.sc) +
  scale_fill_manual(values = pal.sc) +
  theme_bw(base_size = 12) +
  theme(strip.text = element_text(size = 12), legend.position = "bottom") +
  labs(y = "Effort (days at sea)", x = NULL, fill = NULL, color = NULL)

Figure 11: Effort time series for the modelled fleets by scenario.

FLTSTKQ_catch   <- FLTSTKQ %>%
  filter(year > 2021, fleet %in% fltnm, stock %in% dyn.stk, indicator == "catch")

plot_fleet_catch <- function(data, fleet_name, pal.sc) {
  df <- subset(data, fleet == fleet_name)
  
  ggplot(df, aes(x = year, y = q50, color = scenario_name, fill = scenario_name)) +
    geom_line(size = 1.2) +
    facet_wrap(~ stock, scales = "free_y", ncol = 3) +
    scale_color_manual(values = pal.sc) +
    scale_fill_manual(values = pal.sc) +
    theme_bw(base_size = 12) +
    theme(
      strip.text = element_text(size = 12),
      legend.position = "bottom",
      plot.title = element_text(face = "bold", vjust = 0.5, hjust = 1)
    ) +
    labs(
      title = fleet_name,
      y = "Catch (tons)",
      x = NULL,
      fill = NULL,
      color = NULL
    )
}

plot_fleet_catch(FLTSTKQ_catch, "TRW", pal.sc)

Figure 12: Catch time series of trawler fleet for the modelled stocks (Figure 1) by scenario.

plot_fleet_catch(FLTSTKQ_catch, "GNS", pal.sc)

Figure 13: Catch time series of gillneter fleet for the modelled stocks (Figure 1) by scenario.

plot_fleet_catch(FLTSTKQ_catch, "LLS", pal.sc)

Figure 14: Catch time series of longliner fleet for the modelled stocks (Figure 1) by scenario.

Explore results in FLBEIA Shiny

flbeiaApp(bio = BIOQ, flt = FLTQ, fltStk = FLTSTKQ, adv = ADVQ, mt= MTQ, mtStk=MTSTKQ)

Time needed to complete this tutorial

fin <- Sys.time()
duration_TC <- fin - init

# print(duration_TC)
print("Time difference of 15.52843 hours")

[1] "Time difference of 15.52843 hours"

References

EU. 2024. “Commission Delegated Regulation (EU) 2024/3089 of 30 September 2024 Amending Regulation (EU) 2019/1241 as Regards Measures to Reduce Incidental Catches of Common Dolphin (Delphinus Delphis) and Other Small Cetaceans in the Bay of Biscay.” http://data.europa.eu/eli/reg_del/2024/3089/oj.

Garcia, D., S. Sanchez, Raul Prellezo, A. Urtizberea, and M. Andres. 2017. “FLBEIA: A Simulation Model to Conduct Bio-Economic Evaluation of Fisheries Management Strategies.” SoftwareX 6.

ICES. 2023a. “EU Additional Request on Mitigation Measures to Reduce by-Catches of Common Dolphin (Delphinus Delphis) in the Bay of Biscay and Iberian Coast.” sr.2023.01. ICES Advisory Committee; https://doi.org/10.17895/ices.advice.21946634.

———. 2023b. “Working Group on Mixed Fisheries Advice Methodology (WGMIXFISH-METHODS).” 105. Vol. 5. ICES Scientific Reports; https://doi.org/10.17895/ices.pub.24496048.

STECF. 2022. “The 2022 Annual Economic Report on the EU Fishing Fleet (STECF 22-06).” Luxembourg: Publications Office of the European Union; https://doi.org/10.2760/120462.

Yool, A., E. E. Popova, and A. C. Coward. 2015. “Future Change in Ocean Productivity: Is the Arctic the New Atlantic?” Journal of Geophysical Research: Oceans 120: 7771–90. https://doi.org/10.1002/2015JC011167.