This work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101000318 (SEAwise).
This work has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101000318 (SEAwise).

1. Introduction

This tutorial falls under SEAwise’s work on modelling the impact of fisheries on the ecosystem under Work Package 4: Ecological effects of fisheries.

Now that we have simulated several harvest strategies with the demersal mixed fisheries of the North Sea with FLBEIA (bioeconomic impact assessment using FLR) (Garcia et al. 2017), we want to learn how to derive certain ecological indicators from the model output. Additionally, we show on how to supplement the model output with data/assumptions from other sources to inform indicators, where an explicit incorporation into the model was not possible. The focus lies on the indicators (or descriptors) proposed by the Marine strategy framework directive (MSFD) defined to capture the concept of ‘good environmental status’ (GES). Therefore they take into account several anthropogenic pressures on the environment including fisheries.

The goals of this tutorial

  • Work with the FLBEIA model output
  • Derive important MSFD indicators directly from the output or add accompanying data
  • Visualise the results

2. R-setup

2.1 Required packages to run this tutorial

In order to execute the code in this tutorial you should have the following packages installed:

install.packages( c("FLCore", "FLFleet", "FLBEIA",
                    "FLash", "FLAssess"), 
                  repos="http://flr-project.org/R")
install.packages( c("plyr","collapse","reshape2","stringr","car","ggplot2"))

2.2 Download data

The data for this tutorial can be downloaded directly from figshare.

# automatically download input files for the tutorial from figshare

# url for the data for the tutorial on "ways to summarise env data"
url = "https://figshare.com/ndownloader/files/57395818?private_link=2a52cb4012105d77e9cb"

# specify data to download
fn <- "WP4_ecological_impacts.zip" # 670mb
# directory to store data (create data directory in the folder of each practical)
# ensure that it is correctly downloaded in the right folder
current.wd  = getwd()
current.lst.folder = gsub("(.*/\\s*(.*$))", "\\2", current.wd)
if(current.lst.folder == "WP4_ecological_impacts"){
fp <- "./"
} else{
  dir.create("./WP4_ecological_impacts",showWarnings = F)
  setwd("./WP4_ecological_impacts")
  fp <- getwd()
}

# (increase timeout for larger files)
# as the file is about ~670mb
options(timeout = max(300, getOption("timeout")))
# download only if the files do not exist yet
dir.nms = list.dirs("./",recursive = F)
if(!all(c("./input","./images","./functions","./themes") %in% dir.nms)){
# download
download.file(url = url,file.path(fp,fn),mode = "wb")
# unzip
unzip(file.path(fp,fn),exdir = fp)
# remove original zip-download
unlink(file.path(fp,fn))
# set back the wd
setwd(fp)
}

Load all the necessary packages:

library(plyr)
library(collapse)
library(reshape2)
library(stringr)
library(car)
library(FLBEIA)
library(ggplot2)

3. Prepare the FLBEIA output

At first we take the results of the FLBEIA simulation runs for 4 different management strategies (FMSY-Min,PGY,Status-quo and a specific Case-study scenario), paired with 3 climate scenarios (noCC, RCP4.5 and RCP8.5).

Management scenarios:

  • Status-quo:
    • This is a baseline scenario, simulating the effort of the recent past forward in time. Fleets are allowed to fully utilise their quota.
  • FMSY-Min:
    • This is a strict implementation of the landing obligation with Fmsy as the target. Fleets are only allowed to fish until their first quota is reached (min-scenario), leading to underutilisation of quotas and so-called choke-effects.
  • PGY-Min:
    • This scenario relaxes the constraints imposed by the landing obligation, allowing to fish each stock up to Fmsy-Upper (higher fishing mortality leading approximately to 95% of MSY) if the stock is in good condition (above the Btrigger reference point). Additionally, a TAC-constraint is in place allowing year-to-year TAC changes of only 20%.
  • Case-study:
    • This is a scenario tailored to reflect the situation in the North Sea demersal mixed fisheries with an enforced landing obligation (Min-scenario), but a more realistic choking situation. Now some flatfish stocks are not considered choke species for particular fleets. Additionally, there is an implementation-error on the Cod-advice due to TAC-overshoots for Cod in the recent past.

Climate scenarios:

  • noCC:
    • A scenario, reflecting the current level of climate change with no additional worsening in the future. The scenario is derived by removing the median-trend from the RCP4.5 scenario runs.
  • RCP4.5:
    • The Representative concentration pathway corresponding to an expected change in radiative forcing of 4.5 W/m^2 until 2100. This scenario presents an intermediary scenario, with some climate mitigation measures in place.
  • RCP8.5:
    • A “worst-case” scenario corresponding to an expected change in radiative forcing of 8.5 W/m2 until 2100.

3.1 Load/Create helper functions & lookup tables

reshape.df = function(df,id.vars,measure.vars,value.vars){
  tmp = list()
  for(i in seq(value.vars)){
    tmp[[i]] = reshape2::dcast(df,as.formula(paste(paste(id.vars,collapse = "+"),paste(measure.vars,collapse = "+"),sep = "~")),value.var = value.vars[i])
    tmp[[i]]$var = value.vars[i]
  }
  tmp = do.call(rbind,tmp)
  # order
  tmp = tmp[,c(id.vars,"var",names(tmp)[!(names(tmp) %in% c(id.vars,"var"))])]
  
  return(tmp)
}

# reformat names of scenarios and fleets
reformat.names = function(x){
  require(stringr)
  require(car)
  
  if("fleet" %in% names(x)){
    x$fleet = nameChange(x$fleet,"old")
  }
  if("metier" %in% names(x)){
    x$metier = nameChange(x$metier,"old")
  }
  if("stock" %in% names(x)){
    x$stock = nameChange(x$stock,"old")
  }
  # reformat scenario names
  climscen = str_match(x$scenario,"noCC|rcp45|rcp85")
  hcr = str_match(x$scenario,"CaseStudy|min|FmsyH20|statusQuo")

  #rename
  hcr = car::recode(hcr,"'min' = 'FMSY-Min';
                    'FmsyH20' = 'PGY-Min';
                    'statusQuo' = 'Status-quo';
                    'CaseStudy' = 'Case-study'")
  
  # build new scenario name
  scen = paste(hcr,climscen,sep = "_")
  # if string ends with underscore replace
  scen = trimws(scen, which = "right", whitespace = "\\_")
  # replace multiple underscores with only one
  scen = gsub("__","_",scen)
  x$scenario = scen
  return(x)  
}

# aggregation function to derive a data.frame with unique combinations, where those that
# are not unique are grouped together in a row seperated by a comma
agg.column = function(df,columns.unique){
  df.unique = unique(df[,columns.unique,drop = F])
  cols.not.unq = names(df)[!names(df) %in% columns.unique]
  for(i in 1:nrow(df.unique)){
    vals.not.unq = df[apply(df[,columns.unique,drop = F],1,function(x) all(x %in% df.unique[i,])),cols.not.unq]
    if(is.null(dim(vals.not.unq))){
      vals.not.unq = paste(unique(vals.not.unq),collapse = ",")
    } else{
      vals.not.unq = apply(vals.not.unq,2,function(x) paste(unique(x),collapse = ","))
    }
    df.unique[i,cols.not.unq]  = vals.not.unq
  }
  return(df.unique)
}

# load fct for name change
source("./functions/nameChange.R")

# ----------------------------------------------------------------------------- #

# read in summary tables
load("./input/Sim.summary.obj.RData")

Sc.stk = scenarios.summary$bioSum
Sc.mtsEff = scenarios.summary$metierSum

# -------------------------------------------- #
# things to calculate from the original output
# -------------------------------------------- #

# load original runs
sim.runs = list.files("./input/flbeia_stochastic_sims","*RData",
                      full.names = T)
length(sim.runs)
# get stock & average age
st = Sys.time()
Avg.age_list = list()
for(i in seq(sim.runs)){
  #print(i)
  sim.name = load(sim.runs[[i]])
  tmp = get(sim.name)
  rm(list = sim.name)
  
  #print(sim.name)
  stks = tmp$biols
  df_list = list()
  n = 1
  for(j in seq(stks)){
    #print(stk@name)
    stk = stks[[j]]
    if(grepl("NEP",stk@name)){
      next()
    }
    ages = as.numeric(rownames(stk@n))
    years = as.numeric(colnames(stk@n))
    tmp = ages*do.call(rbind,lapply(1:length(ages),function(x) (stk@mat@.Data[[1]]*stk@wt*stk@n)[x,]/apply(stk@mat@.Data[[1]]*stk@wt*stk@n,2,sum)))
    mean.age = apply(tmp,2,sum)
    df_list[[n]] = data.frame(year = years,stock = names(stks)[j],mean.age = mean.age)
    n = n+1
  }  
  # join
  mean.age_df = data.frame(sim.name = sim.name,
                           iter.ens = as.numeric(str_extract(sim.name,"(?<=iter)[0-9]*")),
                           collapse::rowbind(df_list))
  # store
  Avg.age_list[[i]] = mean.age_df  
  # remove sim to free data space (otherwise it will crash...)
  rm(tmp)
}
end = Sys.time()
time.taken = end - st
print(time.taken) # Time difference of 14.65208 mins for 1200 runs
Avg.age = collapse::rowbind(Avg.age_list)
# add scenario name
Avg.age$scenario = Avg.age$sim.name
Avg.age = reformat.names(Avg.age)

4. MSFD-Indicators

Now we calculate the MSFD-indicators: The focus is on those that can be directly derived from the FLBEIA output with a few assumptions or those where a link with other accessible data sources is possible.

  • Descriptor 1: Marine biodiversity
    • D1C1: indicators and methods for setting thresholds for the bycatch mortality rate
  • Descriptor 3: Commercial species
    • D3C3: health indicator reflecting the age and size distribution of individuals in the population of commercially-exploited species
  • Descriptor 4: Food webs
    • D4C3: Large-fish indicator (LFI)
  • Descriptor: 6: Seabed integrity
    • D6: Effort by main-demersal gear or risk for benthic communities

4.1 D1C1 - Bycatch risk/estimated bycatch for PET species

Since it was rather difficult to exactly relate the bycatch rates of Harbour porpoise (given in individuals per days at sea) to the effort represented in FLBEIA (given in kWdays), we only did a linear scaling of bycaught individuals in the year 2020 relative to the effort of bycatch-relevant gears to estimate bycaught individuals for the future scenarios. Although additional factors affecting harbour porpoise bycatch are not considered here (e.g. spatial overlap, seasonal overlap, fisher behaviour to actively avoid bycatch), this estimate can give valuable insight regarding the relative impact of each management measure on harbour porpoise bycatch.

# --------------------- #
# D1C1: Bycatch of PETs
# --------------------- #

# Here we estimate the bycatch risk/bycaught individuals of harbour porpoise per scenario

# Harbour porpoise
porpoise_bycatch = data.frame(region = "Greater North Sea",
                              metier = c("GNS","GTR","SDN","OTB"),
                              bycatch.porpoise = TRUE)

porpoise_bycatch$general.metier = ifelse(porpoise_bycatch$metier %in% c("GNS","GTR"),
                                         "Nets",
                                         ifelse(porpoise_bycatch$metier == "SDN",
                                                "Seine",
                                                ifelse(porpoise_bycatch$metier == "OTB",
                                                       "Otter", NA)))
# aggregate (so that non-dublicate values are in the same row seperated by a comma)
porpoise_bycatch = agg.column(porpoise_bycatch,c("region","bycatch.porpoise","general.metier"))

# get fleet agg-table and 
flt_agg = read.csv("./input/fleet_metier_agg.csv",sep = ";")
flt_agg = flt_agg[,c("fl","met","country_name","fishing_tech")]

flt_agg$general.metier = ifelse(grepl("DSeine|Seine",paste(flt_agg$fl,flt_agg$met,sep = "_")),
                                "Seine",
                                ifelse(grepl("Nets|GN1|GT1",paste(flt_agg$fl,flt_agg$met,sep = "_")) |
                                         grepl("DFN",flt_agg$fishing_tech),
                                       "Nets",
                                       ifelse(grepl("Otter|otter|OTB|TR1|TR2",paste(flt_agg$fl,flt_agg$met,sep = "_")),
                                              "Otter",NA)))

# merge with bycatch flag for porpoise
flt_agg = merge(flt_agg,porpoise_bycatch,
                by = c("general.metier"),all.x = T)
flt_agg[,c("general.metier","fl","met","fishing_tech","bycatch.porpoise")]

# set bycatch rates to zero if there is NA (pots and others)
flt_agg$bycatch.porpoise[is.na(flt_agg$bycatch.porpoise)] <- FALSE

flt_agg_bycatch = flt_agg[,c("fl","met","bycatch.porpoise")]
names(flt_agg_bycatch)[1:2] <- c("fleet","metier")

# merge
Sc.mets_bycatch = dplyr::inner_join(Sc.mtsEff,
                                    flt_agg_bycatch,by = c("fleet","metier"))


# do a simple bycatch scaling only taking the main gears into account
# 1627 individuals of harbour porpoise caught in 2020
# get the effort of the respective gears in 2020 and estimate in relative terms

group_mts_bycatch = fgroup_by(Sc.mets_bycatch,year,scenario,iter)
Sc_bycatch.sum = fsummarise(group_mts_bycatch,
                            effort_bycatch = sum(effort[bycatch.porpoise == TRUE]))
# derive 2020 bycatch ratio
Sc_bycatch.sum$porpoise_bycatch.ratio = 1627/unique(Sc_bycatch.sum$effort_bycatch[Sc_bycatch.sum$year == 2020])
# calculate HP bycatch
Sc_bycatch.sum$porpoise_bycatch_ind = Sc_bycatch.sum$effort_bycatch*Sc_bycatch.sum$porpoise_bycatch.ratio

# calc quantiles
group_bycatchSum = fgroup_by(Sc_bycatch.sum,scenario,year)
Sc_bycatch_quant = plyr::ddply(Sc_bycatch.sum,c("scenario","year"),summarise,
                         variable = "porpoise_bycatch_ind",
                         q025 = collapse::fquantile(porpoise_bycatch_ind,0.025,na.rm = T),
                         q50 = collapse::fquantile(porpoise_bycatch_ind,0.5,na.rm = T),
                         q975 = collapse::fquantile(porpoise_bycatch_ind,0.975,na.rm = T))
#Sc_bycatch_quant$variable = porpoise_bycatch_ind

Sc_bycatch.avg_long = pivot(Sc_bycatch_quant,ids = c("scenario","year","variable"),
                            names = list("quantile","value"),
                            how = "longer")
# calc avg. over the respective periods
D1C1_porpoise = plyr::ddply(Sc_bycatch.avg_long,c("scenario","variable","quantile"),summarise,
                      period2025_2030 = mean(value[year >= 2025 & year <= 2035]),
                      period2035_2040 = mean(value[year >= 2035 & year <= 2040]),
                      period2045_2050 = mean(value[year >= 2045 & year <= 2050]),
                      period2055_2060 = mean(value[year >= 2055 & year <= 2060]))

D1C1_porpoise = reshape.df(D1C1_porpoise,id.vars = c("scenario","quantile"),
                           measure.vars = "variable",
                           value.vars = names(D1C1_porpoise)[grep("period",names(D1C1_porpoise))])
D1C1_porpoise$var = sub("period","",D1C1_porpoise$var)
names(D1C1_porpoise)[grep("porpoise",names(D1C1_porpoise))] <- paste("D1C1",names(D1C1_porpoise)[grep("porpoise",names(D1C1_porpoise))],sep = "_")

D1C1_porpoise

MSFD_indicators = D1C1_porpoise

4.2 D3C3 - Stock health

Here we focus on average stock/spawner age (ASA) being an indicator of the age-composition within the stock. Higher values indicate a larger proportion of older individuals within the stock, which have a higher reproductive output per capita (due to higher maturity & weight), thus being beneficial for the stock reneval.

# --------------------------- #
# D3C3: Avg. spawn. age (ASA)
# --------------------------- #

Avg.age_scen = Avg.age[!grepl("NEP",Avg.age$stock),]
group_avgAge_all = fgroup_by(Avg.age_scen,scenario,year)
Avg.age_quant_all = fsummarise(group_avgAge_all,
                               q025 = collapse::fquantile(mean.age,0.025,na.rm = T),
                               q50 = collapse::fquantile(mean.age,0.5,na.rm = T),
                               q975 = collapse::fquantile(mean.age,0.975,na.rm = T))
Avg.age_quant_all$variable = "average_age"

# melt and cast again
Avg.age_quant_melt2 = pivot(Avg.age_quant_all,ids = c("scenario","year","variable"),
                            names = list("quantile","value"),how = "longer")

# cast to wide format again
Avg.age_quant_cast2 = dcast(Avg.age_quant_melt2,year + scenario + quantile ~ variable,value.var = "value")

# calc avg. over the respective periods
D3C3_ASA = ddply(Avg.age_quant_cast2,c("scenario","quantile"),summarise,
                 period2025_2030 = mean(average_age[year >= 2025 & year <= 2035]),
                 period2035_2040 = mean(average_age[year >= 2035 & year <= 2040]),
                 period2045_2050 = mean(average_age[year >= 2045 & year <= 2050]),
                 period2055_2060 = mean(average_age[year >= 2055 & year <= 2060]))
D3C3_ASA$variable = "D3C3_ASA"
D3C3_ASA = reshape.df(D3C3_ASA,id.vars = c("scenario","quantile"),
                      measure.vars = "variable",
                      value.vars = names(D3C3_ASA)[grep("period",names(D3C3_ASA))])
D3C3_ASA$var = sub("period","",D3C3_ASA$var)

head(D3C3_ASA)

# add to table
MSFD_indicators = merge(MSFD_indicators,D3C3_ASA,by = c("scenario","var","quantile"),sort = F,all.x = T)

4.3 D4C3 - Large fish indicator

The large fish indicator (LFI) defined as the proportion of fish biomass being greater than 40cm of length (in relation to the whole fish biomass of selected species) should reflect the proportion of selected species at the top of the foodweb and is considered to be responsive to the general fishing pressure (Greenstreet et al. 2010). In the North Sea it was declining due to increased fishing pressure especially from trawls, however showed a slight recovery from 2000-2012 corresponding to effort reductions occurring in the same time (Engelhard et al. 2015).

# ------------------------------- #
# D4C3: Large fish indicator
# Biomass_fish>40cm/Biomass_all
# ------------------------------- #

# length - weight relationship from fishbase.org for the respective species
lut.lw <- rbind(
  data.frame(stock = "COD-NS", a = 0.0069, b = 3.08),
  data.frame(stock = "PLE-NS", a = 0.0089, b = 3.04),
  data.frame(stock = "PLE-EC", a = 0.0089, b = 3.04),
  data.frame(stock = "HAD", a = 0.0059, b = 3.13),
  data.frame(stock = "WHG-NS", a = 0.0065, b = 3.05),
  data.frame(stock = "POK", a = 0.0074, b = 3.04),
  data.frame(stock = "SOL-NS", a = 0.0068, b = 3.10),
  data.frame(stock = "SOL-EC", a = 0.0068, b = 3.10),
  data.frame(stock = "TUR", a = 0.01122, b = 3.06),
  data.frame(stock = "WIT", a = 0.00490, b = 3.07)
)
# W = a* L^b
convert.W.to.L = function(W,a,b){
  return((W/a)^(1/b))
}

# get stock & biomass of fish larger 40cm
bio40_list = list()
for(i in seq(sim.runs)){
  print(i)
  sim.name = load(sim.runs[[i]])
  tmp = get(sim.name)
  rm(list = sim.name)
  
  print(sim.name)
  stks = tmp$biols 
  df_list = list()
  for(j in seq(stks)){
    stk = stks[[j]]
    stk.name = name(stk)
    if(grepl("NEP",stk.name)){
      bio40 = 0
      years = as.numeric(colnames(stk@n))
    } else{
      years = as.numeric(colnames(stk@n))
      # convert weight to length
      stk_wgts = stk@wt
      # convert with the coefficients per stock
      coefs = lut.lw[lut.lw$stock == nameChange(stk.name,"old"),]
      #cat(coefs$stock,stk.name,"\n")
      stk.length = convert.W.to.L(stk_wgts*1000,a = coefs[["a"]],b = coefs[["b"]])
      # select age classes that are above 40cm
      zeros = matrix(0,nrow(stk.length),ncol(stk.length))
      zeros[stk.length>40] <- (stk@wt*stk@n)[stk.length>40]
      tmp.bio40 = zeros 
      bio40 = apply(tmp.bio40,2,sum)
    }
    df_list[[j]] = data.frame(year = years,
                              stock = names(stks)[j],
                              Biomass_larger40 = as.numeric(bio40))
  }  
  # join
  bio40_df = data.frame(sim.name = sim.name,
                        iter = as.numeric(str_extract(sim.name,"(?<=iter)[0-9]*")),
                        collapse::rowbind(df_list))
  # store
  bio40_list[[i]] = bio40_df 
  # remove sim to free data space (otherwise it will crash...)
  rm(tmp)
}

biomass40 = collapse::rowbind(bio40_list)
biomass40$scenario = biomass40$sim.name
biomass40 = reformat.names(biomass40)

# sum biomass > 40 for each iteration across stokcs
group_bio40 = fgroup_by(biomass40,scenario,year,iter)
biomass40_sum = fsummarise(group_bio40,
                           Biomass_larger40_sum = sum(Biomass_larger40,na.rm = T))
group_stk = fgroup_by(Sc.stk,scenario,year,iter)
total_biomass = fsummarise(group_stk,
                           total_biomass = sum(biomass,na.rm = T))

# merge with total biomass over all fish in the model
biomass40_sum = collapse::join(biomass40_sum,total_biomass,on = c("scenario","year","iter"),how = "inner")
# calc the Large fish indicator
biomass40_sum$LFI = biomass40_sum$Biomass_larger40_sum/biomass40_sum$total_biomass

# calc quantiles over all iterations, per scenario
group_LFI = fgroup_by(biomass40_sum,scenario,year)
LFI_quant = fsummarise(group_LFI,
                       q025 = collapse::fquantile(LFI,0.025,na.rm = T),
                       q50 = collapse::fquantile(LFI,0.5,na.rm = T),
                       q975 = collapse::fquantile(LFI,0.975,na.rm = T))
LFI_quant$variable = "LFI"
LFI_quant_long = melt(LFI_quant,id.vars = c("scenario","year","variable"),variable.name = "quantile")

# calc avg. over the respective periods
D4C3_LFI = ddply(LFI_quant_long,c("scenario","variable","quantile"),summarise,
                 period2025_2030 = mean(value[year >= 2025 & year <= 2035]),
                 period2035_2040 = mean(value[year >= 2035 & year <= 2040]),
                 period2045_2050 = mean(value[year >= 2045 & year <= 2050]),
                 period2055_2060 = mean(value[year >= 2055 & year <= 2060]))
# reshape
D4C3_LFI = reshape.df(D4C3_LFI,id.vars = c("scenario","quantile"),
                      measure.vars = "variable",
                      value.vars = names(D4C3_LFI)[grep("period",names(D4C3_LFI))])
D4C3_LFI$var = sub("period","",D4C3_LFI$var)

#rename
names(D4C3_LFI)[grep("LFI",names(D4C3_LFI))] <- paste("D4C3",names(D4C3_LFI)[grep("LFI",names(D4C3_LFI))],sep = "_")

head(D4C3_LFI)

# add to table
MSFD_indicators = merge(MSFD_indicators,D4C3_LFI,by = c("scenario","var","quantile"),sort = F,all.x = T)

4.4 D6 effort by demersal gears/ RBS-indicator

To estimate the impact on the seabed one could first extract the effort used by all relevant gears that are dragged over the seabed. In the North Sea FLBEIA model this would involve beam trawls, otter trawls, and Seines or danish seines. For a first impression of bottom impact one could compare those and see which scenario has the lowest effort involved. A more refined picture can be gathered by calculating the RBS-Indicator.

Here we only describe briefly the RBS methodology, load in the results of the RBS-indicator calculations and visualize those later. RBS was previously estimated for the current fishing activities in SEAwise Task 4.3 (Van Hoey et al. 2024; Rijnsdorp et al. 2020).For calculating the RBS-indicator for the FLBEIA-output, we matched the SAR estimates of the main Benthis metiers (Eigaard et al. 2015; Rijnsdorp et al. 2020) with the fishing fleets and metiers as defined in the FLBEIA North Sea model in the year 2021. Benthis metiers were re-categorized to five groups: Nephrops otter trawls (OT_CRU, also including OT_MIX_CRU_DMF), whitefish otter trawls (OT_DMF, also including OT_MIX), flatfish-directed beam trawls (TBB_DMF), Danish seines (SDN_DMF) and fly-shooters (SSC_DMF). Three Benthis metiers (TBB_MOL, OT_MIX_DMF_PEL and OT_MIX_DMF_BEN) did not have any data in 2021, and three other Benthis metiers (DRB_MOL, OT_SPF, TBB_CRU) were not dynamically assessed in the FLBEIA model, and therefore grouped to “other” metiers and kept constant on status quo level.

The effort trajectories from FLBEIA (three climate scenarios x three management scenarios (1) Fmsy-min, (2) PGY-min and (3) a specific case study scenario) were then used to inform the RBS-indicator relative to the derived 2021-level for the future. The spatial distribution of Benthis metier-specific effort was considered constant at status-quo over the time periods and scenarios. In detail, the cumulative percentage of MSFD broad habitat types with RBS-estimates were mapped for the nine simulated scenarios, and compared those with the 2021 baseline assessment. Then the percentage of the total area with RBS values > 0.8 (indiscriminate of MSFD broad habitat types) was estimated. This latter indicator is reported in the main body of this study as a generic way to compare the effect of the simulated scenarios on benthic community biomass, and its abilities to achieve Good Environmental Status for seafloor integrity (MSFD D6).

# -------------------------------------- #
# D6: Effort by main demersal gear type
# -------------------------------------- #

Sc.flts = scenarios.summary$fltSum

# main demersal gear types:
# Beam, Otter, DSeine/Seine

Sc.flts$main.demersal.gear = str_extract(Sc.flts$fleet,"Beam|Otter|DSeine|Seine")
Sc.flts$main.demersal.gear = ifelse(Sc.flts$main.demersal.gear %in% c("DSeine","Seine"),
                                    "DSeine/Seine",Sc.flts$main.demersal.gear)
# sum avg. effort per gear type
flt.gear.type.eff = fsummarise(fgroup_by(Sc.flts,scenario,year,fleet,main.demersal.gear),
                               q025 = collapse::fquantile(effort,0.025,na.rm = T),
                               q50 = collapse::fquantile(effort,0.5,na.rm = T),
                               q975 = collapse::fquantile(effort,0.975,na.rm = T))
#flt.gear.type.eff$variable = "effort"
flt.gear.type.eff_melt = melt(flt.gear.type.eff,id.vars = c("scenario","fleet","year","main.demersal.gear"),
                              value.name = "effort",variable.name = "quantile")
flt.gear.type.eff_melt$effort = as.numeric(flt.gear.type.eff_melt$effort)
gear.type.eff_sum = fsummarise(fgroup_by(flt.gear.type.eff_melt,scenario,year,main.demersal.gear,quantile),
                               sum.eff = sum(effort,na.rm = T))


# calc avg. over the respective periods
D6.eff.per.gear = ddply(gear.type.eff_sum,c("scenario","main.demersal.gear","quantile"),summarise,
                        period2025_2030 = mean(sum.eff[year >= 2025 & year <= 2035]),
                        period2035_2040 = mean(sum.eff[year >= 2035 & year <= 2040]),
                        period2045_2050 = mean(sum.eff[year >= 2045 & year <= 2050]),
                        period2055_2060 = mean(sum.eff[year >= 2055 & year <= 2060]))
# remove NA-gear types
D6.eff.per.gear = D6.eff.per.gear[complete.cases(D6.eff.per.gear),]
# change names
D6.eff.per.gear$main.demersal.gear = paste0("D6_Eff.of.",D6.eff.per.gear$main.demersal.gear)

D6.eff.per.gear = reshape.df(D6.eff.per.gear,id.vars = c("scenario","quantile"),
                             measure.vars = "main.demersal.gear",
                             value.vars = names(D6.eff.per.gear)[grep("period",names(D6.eff.per.gear))])
D6.eff.per.gear$var = sub("period","",D6.eff.per.gear$var)
#names(D6.eff.per.gear)[which(names(D6.eff.per.gear$var) == "var")] <- "Time.period"

head(D6.eff.per.gear)

# add to table
MSFD_indicators = merge(MSFD_indicators,D6.eff.per.gear,by = c("scenario","var","quantile"),sort = F,all.x = T)
# ---------------------------------------- #
# RBS indicator (% of area with RBS < 0.8)
# ---------------------------------------- #

rbs_df = read.csv("./input/RBS.csv",header = T,sep = ";",dec = ",")

rbs_df_sub = rbs_df[rbs_df$year %in%  c("2025-2030","2055-2060"),]
unique(rbs_df_sub$scenario)

output_rbs = rbs_df_sub
names(output_rbs)[grep("Percent_RBS_higher_0.8",names(output_rbs))] <- "q50"
names(output_rbs)[grep("year",names(output_rbs))] <- "Time.period"

output_rbs$variable = "D6_RBS"

# -------------------------------------------------------- #
# more detailed look at the RBS-indicator per habitat type
# -------------------------------------------------------- #

load("./input/RBS_assessments.RData", verbose = T)
x <- prop_BHT

d <- dim(x)
n <- names(x)

y <- x |> pivot(ids = 1:2,how = "longer",names = list("Scenario","areaGES"), values = 3:d[2],factor = FALSE)
y$MSFD_BBHT[y$MSFD_BBHT == "Na"] <- NA

y$RBSthresh <- as.numeric(substr(y$Scenario,10,13))
y$Scenario <- substr(y$Scenario, start = 15, stop = 100)
y$year <- as.numeric(unlist(lapply(strsplit(y$Scenario, split = "_", fixed = T), FUN = function(x){tail(x, 1)})))
y = y[!is.na(y$year),]

y$HCR <- unlist(lapply(strsplit(y$Scenario, split = "_", fixed = T), FUN = function(x){rev(tail(x, 10))[2]}))
y$Climate <- unlist(lapply(strsplit(y$Scenario, split = "_", fixed = T), FUN = function(x){rev(tail(x, 10))[3]}))

y$yearSpan <- NA
y$yearSpan[which(y$year>=2025 & y$year<=2030)] <- "2025-2030"
y$yearSpan[which(y$year>=2035 & y$year<=2040)] <- "2035-2040"
y$yearSpan[which(y$year>=2045 & y$year<=2050)] <- "2045-2050"
y$yearSpan[which(y$year>=2055 & y$year<=2060)] <- "2055-2060"
y$yearSpan[which(y$year==2021)] <- "Current"

y <- y |> fmutate(Scen = paste(HCR, Climate, sep = "_"))

# KOBE plot
agg <- y |> fgroup_by(RBSthresh, Scen, MSFD_BBHT, yearSpan) |> 
  fsummarise(areaGES = mean(areaGES)) |> 
  fungroup()
agg

RBSthresh_lev <- 0.8
yearSpan_lev <- c("Current", "2025-2030", "2055-2060")
Scen_lev <- c("FMSY-Min_rcp45", "FMSY-Min_rcp85", "PGY-Min_rcp45", "PGY-Min_rcp85")


dfDat <- agg[agg$Scen %in% Scen_lev & 
             agg$RBSthresh %in% RBSthresh_lev & 
             agg$yearSpan %in% yearSpan_lev & 
             !is.na(agg$MSFD_BBHT),]


dfRef <- subset(agg, 
  Scen %in% "Status-quo_noCC" & 
  RBSthresh %in% RBSthresh_lev &
  yearSpan %in% yearSpan_lev & 
  !is.na(MSFD_BBHT)#
)

tmp <- dfDat |> fselect(RBSthresh, Scen, MSFD_BBHT) |> 
  unique()

tmp$yearSpan <- "Current"
tmp <- collapse::join(x = tmp, y = fselect(dfRef, RBSthresh, MSFD_BBHT, areaGES),
            how = "left")

dfDat <- rbind(dfDat, tmp)

dfDat$yearSpan <- factor(dfDat$yearSpan, levels = c("Current", "2025-2030", "2055-2060"))
unique(dfDat$yearSpan)

# rename scenarios
nested_str_replace = function(string,pattern_string,replace_string){
  for(i in 1:length(pattern_string)){
    string = str_replace(string,pattern_string[i],replace_string[i])
  }
  return(string)
}

dfDat$Scen = nested_str_replace(gsub("_"," | ",dfDat$Scen),
                                c("rcp45","rcp85"),
                                c("RCP4.5","RCP8.5"))

5. Bring everything together and visualize

# -------------------------------- #
# bring everything together & plot
# -------------------------------- #

# change var to Time.period
names(MSFD_indicators)[which(names(MSFD_indicators) == "var")] <- "Time.period"
# reformat time-period underscore to dash
MSFD_indicators$Time.period = gsub("_","-",MSFD_indicators$Time.period)

# resort columns
MSFD_indicators$Case.study = "NS"
MSFD_indicators$Model.name = "FLBEIA"

MSFD_indicators = MSFD_indicators[,c("Case.study","Model.name",
                                     names(MSFD_indicators)[!names(MSFD_indicators) %in% c("Case.study","Model.name")])]

# cast/melt
id.vars = c("Case.study","Model.name","scenario","Time.period","quantile")
output.msfd_melt = pivot(MSFD_indicators,ids = id.vars,how = "longer")
# cast
output.msfd_cast = reshape2::dcast(output.msfd_melt,Case.study + Model.name + scenario + Time.period + variable ~ quantile,value.var = "value")

# prepare data for plotting
plot_df = data.table::rbindlist(list(output.msfd_cast,output_rbs),fill = T)

plot_df = plot_df[plot_df$variable %in% c("D4C3_LFI","D1C1_porpoise_bycatch_ind","D3C3_ASA","D6_Eff.of.Beam","D6_RBS"),]
plot_df$HCR = gsub("_.*","",plot_df$scenario)
plot_df$scenMod = sapply(strsplit(plot_df$scenario,"_"),function(x) paste(x[2:length(x)],collapse = " & \n"))
plot_df = plot_df[plot_df$Time.period %in% c("2025-2030","2055-2060"),]

# order scenarios
plot_df$HCR = factor(plot_df$HCR,levels = c("FMSY-Min","PGY-Min","Case-study","Case-study-mod","Status-quo"))

num_lines <- length(unique(plot_df$HCR))
num_lines <- num_lines - 1

# number of panel rows
num_panel_rows <- ceiling(length(unique(plot_df$variable))/ 2)

# indicator lut
# select indicators that should be plotted
plot_df$variable = collapse::recode_char(as.character(plot_df$variable),
                                 D3C3_ASA = 'D3C3: Average spawner age',
                                 D1C1_porpoise_bycatch_ind = 'D1C1: Bycatch of porpoises \r\n [ind.]', 
                                 D4C3_LFI = 'D4C3: Large fish indicator \n[biomass > 40cm/total biomass]',
                                 D6_Eff.of.Beam = 'D6: Effort of beam trawls \r\n [1000kwdays]',
                                 D6_RBS = 'Relative benthic status \r\n [% of area with RBS larger 0.8]')

# remove those that have no data
plot_df = plot_df[!is.na(plot_df$q50),,]

# plot
source("./themes/ggplot_theme_Publication-2.R") 
# dark-mode ggplot themes from https://github.com/koundy/ggplot_theme_Publication

colours = c("aquamarine3","orange","orangered2")
plt.indicators = ggplot() +
  geom_point(data = plot_df, aes(x = HCR, y = q50, col = scenMod,
                                 shape = Time.period,group = paste(scenMod,Time.period)),
             position = position_dodge(width = 1)) +
  geom_errorbar(data = plot_df, aes(x = HCR, ymin = `q025`, ymax = `q975`,col = scenMod,group = paste(scenMod,Time.period)),
                width = .2, stat = "identity", position = position_dodge(width = 1),
                show.legend = F) +
  geom_vline(xintercept = 1:num_lines + .5, col = "grey30" ) +
  scale_color_manual(values = colours) +
  labs(x = "\nManagement scenario",y = "") +
  facet_wrap( .~ variable, scales = "free", labeller = label_wrap_gen(width = 29),
              nrow = num_panel_rows ) + 
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 40, hjust = 1,vjust = 1),
        legend.position = "bottom") + 
  guides(fill = guide_legend(title = "Scenario"),
         shape = guide_legend(title = "Time period")) +
  guides(color = guide_legend(title = "Scenario", override.aes = list(shape = 15, size = 5))) + 
  theme_dark_blue() + 
  theme(legend.position = "right",legend.direction = "vertical") +
  ggtitle("MSFD-indicators for the North Sea case study")
  

print(plt.indicators)

# ------------------------- #
# visualise RBS per habitat
# with 0.8 threshold
# ------------------------- #

# visualise
RBSthresh_filter <- 0.8
dfsub <- dfDat[dfDat$RBSthresh == RBSthresh_filter,]
habitat.cols = colorRampPalette(c("#984ea3","#00798c", "#fb9a99","#d1495b","#660000","#66a182","#93c572","#EFDFBB","#feca01"))(length(unique(dfsub$MSFD_BBHT)))
p1 <- ggplot(dfsub) + aes(x = yearSpan, y = areaGES, color = MSFD_BBHT, group = MSFD_BBHT) +
  facet_wrap(~Scen, ncol = 2) +
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0.75,
    fill = 'orange',alpha = 0.5)  + 
  annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 0.75, ymax = Inf,
    fill = 'cyan3',alpha = 0.5)  + 
  geom_point() +
  geom_line() +
  scale_color_manual(values = habitat.cols) + 
  scale_y_continuous(limits = c(0.5,1), expand = expansion(mult = c(0, 0), add = c(0, 0.02))) +
  labs(y = "area GES [frac.]", x = element_blank(), title = paste("RBS threshold =", RBSthresh_filter)) +
  theme_dark_blue() +
  theme(legend.position = "right",legend.direction = "vertical") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
print(p1)

And this is it. Now you know how to derive relevant MSFD-indicators from an MSE-type simulation in FLBEIA and estimate wider impacts in an ecosystem context. Happy coding. 😄

6. References

Eigaard, Ole R., Francois Bastardie, Mike Breen, Grete E. Dinesen, Niels T. Hintzen, Pascal Laffargue, Lars O. Mortensen, et al. 2015. “Estimating Seabed Pressure from Demersal Trawls, Seines, and Dredges Based on Gear Design and Dimensions.” ICES Journal of Marine Science 73 (suppl_1): i27–43. https://doi.org/10.1093/icesjms/fsv099.
Engelhard, Georg H., Christopher P. Lynam, Bernardo García-Carreras, Paul J. Dolder, and Steven Mackinson. 2015. “Effort Reduction and the Large Fish Indicator: Spatial Trends Reveal Positive Impacts of Recent European Fleet Reduction Schemes.” Environmental Conservation 42 (3): 227–36. https://doi.org/10.1017/S0376892915000077.
Garcia, D., S. Sanchez, R. Prellezo, A. Urtizberea, and M. Andres. 2017. “FLBEIA: A Simulation Model to Conduct Bio-Economic Evaluation of Fisheries Management Strategies.” SoftwareX 6: 141–47.
Greenstreet, Simon P. R., Stuart I. Rogers, Jake C. Rice, Gerjan J. Piet, Emma J. Guirey, Helen M. Fraser, and Rob J. Fryer. 2010. “Development of the EcoQO for the North Sea Fish Community.” ICES Journal of Marine Science 68 (1): 1–11. https://doi.org/10.1093/icesjms/fsq156.
Rijnsdorp, A D, J G Hiddink, P D van Denderen, N T Hintzen, O R Eigaard, S Valanko, F Bastardie, et al. 2020. “Different Bottom Trawl Fisheries Have a Differential Impact on the Status of the North Sea Seafloor Habitats.” ICES Journal of Marine Science 77 (5): 1772–86. https://doi.org/10.1093/icesjms/fsaa050.
Van Hoey, Gert, Luke Batts, Esther Beukhof, Pierluigi Carbonara, Asbjørn Christensen, Jochen Depestele, Josefine EGEKVIST, et al. 2024. “SEAwise Report on Changes in the Spatiotemporal Benthic Effects of Fishing on Benthic Habitats Relative to Suggested Reference Levels, Both with Respect to Area Impacted and Impact Intensity in Response to Spatial Management.”

7. Software Versions

R version 4.4.2 (2024-10-31 ucrt)

car: 3.1.3

collapse: 2.0.19

FLash: 2.5.11

FLAssess: 2.6.3

FLBEIA: 1.16.0.0

FLCore: 2.6.20.9335

FLFleet: 2.6.1

ggplot2: 3.5.2

plyr: 1.8.9

reshape2: 1.4.4

stringr: 1.5.1

Compiled: 2025-Aug-25