In Bai-Li-NOAA/IFA4EBFM: Interdisciplinary Forecasting Approach for Ecosystem-Based Fisheries Management

knitr::opts_chunk$set(echo = TRUE)

devtools::load_all()

required_pkg <- c(
  "RColorBrewer", "ggplot2", "tidyverse",
  "reshape2", "here"
)

pkg_to_install <- required_pkg[!(required_pkg %in%
  installed.packages()[, "Package"])]
if (length(pkg_to_install)) install.packages(pkg_to_install, repos = "http://cran.us.r-project.org")

invisible(lapply(required_pkg, library, character.only = TRUE))

options(digits = 2)

Simulation scenarios

We want to incorporate climate, environmental, or biological drivers to the Ecopath with Ecosim (EwE) model. Base on Buchheister et al. (2016) paper, we explored the following scenarios:

devtools::load_all()

scenarios <- read.csv(here::here("data", "ewe", "7ages", "ecosim_with_environmental_driver", "ecosim_scenarios.csv"))
functional_groups = c(
  "StripedBass0",
  "StripedBass2_5",
  "StripedBass6",
  "AtlanticMenhaden0",
  "AtlanticMenhaden1",
  "AtlanticMenhaden2",
  "AtlanticMenhaden3",
  "AtlanticMenhaden4",
  "AtlanticMenhaden5",
  "AtlanticMenhaden6",
  "SpinyDogfish",
  "BluefishJuvenile",
  "BluefishAdult",
  "WeakfishJuvenile",
  "WeakfishAdult",
  "AtlanticHerring0_1",
  "AtlanticHerring2",
  "Anchovies",
  "Benthos",
  "Zooplankton",
  "Phytoplankton",
  "Detritus"
)

bratio <- matrix(NA, nrow=nrow(scenarios), ncol=length(functional_groups))
colnames(bratio) <- functional_groups
rownames(bratio) <- scenarios$driver_code

for (id in 1:nrow(scenarios)){
  ewe_output <- read_ewe_output(
    file_path = here::here("data", "ewe", "7ages", "ecosim_with_environmental_driver", scenarios$driver_code[id]),
    file_names = "biomass_annual.csv",
    skip_nrows = 8,
    plot = FALSE,
    figure_titles = NULL,
    functional_groups =functional_groups,
    figure_colors = NULL
  )

  for (i in 2:ncol(ewe_output[[1]])){
    bratio[id, i-1] <- ewe_output[[1]][nrow(ewe_output[[1]]),i]/ewe_output[[1]][1,i]
  }

}

scenarios <- cbind(scenarios, bratio)
row.names(scenarios) <- NULL
write.csv(scenarios, here::here("data", "ewe", "7ages", "ecosim_with_environmental_driver", "ecosim_scenarios.csv"), row.names = F)

subsenarios_id <- 1:(nrow(scenarios)-1)

knitr::kable(scenarios[, c("scenario_id", "driver_code", "driver_name", "driver_category", "justification", "data_source", "sum_of_squared_differences_menhaden0", "sum_of_squared_differences_total")])

Relative biomass

To check if any functional groups have relatively large decrease, we calculated relative biomass ($B_{2017}/B_{1985}$) for each functional group and each scenario.

heatmap(bratio[subsenarios_id, ], Rowv = NA, Colv = NA,
        scale="none", xlab="", ylab="Simulation Scenarios", main="B2017/B1985", 
        col=colorRampPalette(brewer.pal(8, "Blues"))(20))
legend(x="bottomright", legend=c("min", "max"), 
       fill=colorRampPalette(brewer.pal(8, "Blues"))(2))

Key findings:

• Base scenario 1: EwE without any environmental drivers and the sum of squared differences was 2339.
• Scenario 5: Base EwE with incorporation of precipitation data had the lowest sum of squared differences of 2298.
• Scenario 3: Base EwE with incorporation of Atlantic Multidecadal Oscillation Index (i.e., lag of one year) had a sum of squared differences of 2343.
• Scenario 10: EwE with incorporation of both AMO index and precipitation data had a sum of squared differences of 2290.
• The minimum $B_{2017}/B_{1985}$ was 0.16 (Anchovies from pp33 scenario) and other $B_{2017}/B_{1985}$ values were above 0.2.

Temporal trends of climate and environmental drivers from key scenarios (1, 3, 5, 10)

driver_dir <- file.path(here::here(), "data", "ewe")
sim_years <- 1985:2017
amo <- read.csv(file.path(driver_dir, "amo_unsmooth_lag1.csv"))
pcp <- read.csv(file.path(driver_dir, "precipitation.csv"))

data <- data.frame(
  Time = rep(sim_years, each=12), 
  Month = rep(1:12, times=length(sim_years)),
  AMO = amo$scaled_value,
  PCP = pcp$scaled_value
)

par(mar = c(5,5,2,5))
with(data, plot(1:nrow(data), AMO, type="l", col="black", 
                xlab="Time",
             ylab="scaled AMO_lag1",
             ylim=c(0,3)))
par(new = T)
with(data, plot(1:nrow(data), PCP, axes=F, xlab=NA, ylab=NA, type="l", lty=2, col="blue"))
axis(side = 4, col="blue")
mtext(side = 4, line = 3, 'Scaled PCP', col="blue")
legend("topleft",
       legend=c("AMO_lag1", "PCP"),
       lty=c(1,2), col=c("black", "blue"))

Biomass over time from key scenarios

scenario_code <- c("base", "pcp", "amo_lag1", "amo_pcp")
labels <- c(
  "Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+",
  "Atlantic Menhaden 0", "Atlantic Menhaden 1",
  "Atlantic Menhaden 2", "Atlantic Menhaden 3",
  "Atlantic Menhaden 4", "Atlantic Menhaden 5",
  "Atlantic Menhaden 6+",
  "Spiny Dogfish", "Bluefish juv", "Bluefish adult",
  "Weakfish juv", "Weakfish adult",
  "Atlantic Herring 0-1", "Atlantic Herring 2+",
  "Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus"
)

for (id in 1:length(scenario_code)){
  ewe_output <- read_ewe_output(
    file_path = here::here("data", "ewe", "7ages", "ecosim_with_environmental_driver", scenario_code[id]),
    file_names = "biomass_annual.csv",
    skip_nrows = 8,
    plot = FALSE,
    figure_titles = NULL,
    functional_groups =functional_groups,
    figure_colors = NULL
  )
  names(ewe_output[[1]])[2:ncol(ewe_output[[1]])] <- labels
  ewe_output_melt <- reshape2::melt(ewe_output[[1]], id.vars = "Year")
  ewe_output_melt$Scenario <- scenario_code[id]

  if (id==1) {
    ewe_data <- ewe_output_melt
  } else{
    ewe_data <- rbind(ewe_data, ewe_output_melt)
  }

}

ggplot(ewe_data, aes(x = Year, y = value, color = Scenario)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  labs(
    title = "Comparison between scenarios",
    x = "Year",
    y = "Biomass (million mt)"
  ) +
  theme_bw()

Key findings:

• The climate and environmental drivers were linked with age 0 of menhaden-like species only, only menhaden, zooplankton, phytoplankton, and detritus showed fluctuations over time. Other functional groups showed relatively smooth trends.
• The sharp decrease in age 0 menhaden biomass was caused by the incorporation of AMO index. The AMO index decreases over time in 1985.

Bratio and biomass over time from scenario 11: link precipitation with primary production and link AMO with age-0 fish

Bratio:

knitr::kable(bratio[c(1, nrow(bratio)), ])

Biomass over time:

scenario_code <- c("base", "pcp_pp_amo_am")
labels <- c(
  "Striped Bass 0", "Striped Bass 2-5", "Striped Bass 6+",
  "Atlantic Menhaden 0", "Atlantic Menhaden 1",
  "Atlantic Menhaden 2", "Atlantic Menhaden 3",
  "Atlantic Menhaden 4", "Atlantic Menhaden 5",
  "Atlantic Menhaden 6+",
  "Spiny Dogfish", "Bluefish juv", "Bluefish adult",
  "Weakfish juv", "Weakfish adult",
  "Atlantic Herring 0-1", "Atlantic Herring 2+",
  "Anchovies", "Benthos", "Zooplkanton", "Phytoplankton", "Detritus"
)

for (id in 1:length(scenario_code)){
  ewe_output <- read_ewe_output(
    file_path = here::here("data", "ewe", "7ages", "ecosim_with_environmental_driver", scenario_code[id]),
    file_names = "biomass_annual.csv",
    skip_nrows = 8,
    plot = FALSE,
    figure_titles = NULL,
    functional_groups =functional_groups,
    figure_colors = NULL
  )
  names(ewe_output[[1]])[2:ncol(ewe_output[[1]])] <- labels
  ewe_output_melt <- reshape2::melt(ewe_output[[1]], id.vars = "Year")
  ewe_output_melt$Scenario <- scenario_code[id]

  if (id==1) {
    ewe_data <- ewe_output_melt
  } else{
    ewe_data <- rbind(ewe_data, ewe_output_melt)
  }

}

ggplot(ewe_data, aes(x = Year, y = value, color = Scenario)) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  labs(
    title = "Comparison between scenarios",
    x = "Year",
    y = "Biomass (million mt)"
  ) +
  theme_bw()

Bai-Li-NOAA/IFA4EBFM documentation built on April 21, 2024, 7:58 p.m.