MPA performance is often measured by examining gradients of various
outcomes, such as catch-per-unit-effort (CPUE) inside vs. various
distances outside the protected area. Here is an example of generating
these kinds of measurements using marlin
library(marlin)
library(tidyverse)
library(patchwork)
theme_set(marlin::theme_marlin())
resolution <- c(5,20) # resolution is in squared patches, so 20 implies a 20X20 system, i.e. 400 patches
patch_area <- 10
seasons <- 2
years <- 20
tune_type <- "depletion"
steps <- years * seasons
yft_diffusion <- 5
yft_depletion <- .2
rec_factor <- 10
mako_depletion <- 0.1
mako_diffusion <- 100
# for now make up some habitat
yft_habitat <- expand_grid(x = 1:resolution[1], y = 1:resolution[2]) %>%
mutate(habitat = 1,
habitat = habitat / max(habitat) * yft_diffusion) %>%
pivot_wider(names_from = y, values_from = habitat) %>%
select(-x) %>%
as.matrix()
mako_habitat <- expand_grid(x = 1:resolution[1], y = 1:resolution[2]) %>%
mutate(habitat = 1,
habitat = habitat / max(habitat) * mako_diffusion) %>%
pivot_wider(names_from = y, values_from = habitat) %>%
select(-x) %>%
as.matrix()
# create a fauna object, which is a list of lists
fauna <-
list(
"Yellowfin Tuna" = create_critter(
scientific_name = "Thunnus albacares",
habitat = yft_habitat, # pass habitat as lists
recruit_habitat = yft_habitat,
adult_diffusion = yft_diffusion, # cells per year
recruit_diffusion = rec_factor * yft_diffusion,
fished_depletion = yft_depletion, # desired equilibrium depletion with fishing (1 = unfished, 0 = extinct),
density_dependence = "pre_dispersal", # recruitment form, where 1 implies local recruitment
seasons = seasons,
ssb0 = 5000,
),
"Shortfin Mako" = create_critter(
scientific_name = "Isurus oxyrinchus",
habitat = list(mako_habitat), # pass habitat as lists
recruit_habitat = mako_habitat,
adult_diffusion = mako_diffusion,
recruit_diffusion = rec_factor * .1,
fished_depletion = mako_depletion,
density_dependence = "local_habitat", # recruitment form, where 1 implies local recruitment
burn_years = 200,
ssb0 = 100,
seasons = seasons,
fec_form = "pups",
pups = 2,
lorenzen_m = FALSE
)
)
fleets <- list("longline" = create_fleet(list(
`Yellowfin Tuna` = Metier$new(
critter = fauna$`Yellowfin Tuna`,
price = 100, # price per unit weight
sel_form = "logistic", # selectivity form, one of logistic or dome
sel_start = .3, # percentage of length at maturity that selectivity starts
sel_delta = .1, # additional percentage of sel_start where selectivity asymptotes
catchability = .01, # overwritten by tune_fleet but can be set manually here
p_explt = 1
),
`Shortfin Mako` = Metier$new(
critter = fauna$`Shortfin Mako`,
price = 10,
sel_form = "logistic",
sel_start = .1,
sel_delta = .01,
catchability = 0,
p_explt = 1
)),
mpa_response = "stay",
base_effort = prod(resolution),
spatial_allocation = "rpue",
resolution = resolution
))
a <- Sys.time()
fleets <- tune_fleets(fauna, fleets, tune_type = tune_type) # tunes the catchability by fleet to achieve target depletion
Sys.time() - a
#> Time difference of 3.599519 secs
# run simulations
a <- Sys.time()
spillover_sim <- simmar(fauna = fauna,
fleets = fleets,
years = years)
Sys.time() - a
#> Time difference of 0.121702 secs
proc_spillover <- process_marlin(spillover_sim, time_step = fauna[[1]]$time_step)
plot_marlin(proc_spillover, max_scale = FALSE)
mpa_locations <- expand_grid(x = 1:resolution[1], y = 1:resolution[2]) |>
arrange(x) |>
mutate(patch_name = paste(x,y, sep = "_"))
mpas <- mpa_locations$patch_name[1:round(nrow(mpa_locations) * .3)]
# mpas <- sample(mpa_locations$patch_name, round(nrow(mpa_locations) * 0.2), replace = FALSE, prob = mpa_locations$x)
mpa_locations$mpa <- mpa_locations$patch_name %in% mpas
a <- Sys.time()
mpa_spillover <- simmar(
fauna = fauna,
fleets = fleets,
years = years,
manager = list(mpas = list(locations = mpa_locations,
mpa_year = floor(years * .5)))
)
Sys.time() - a
#> Time difference of 0.03970814 secs
proc_mpa_spillover <- process_marlin(mpa_spillover, time_step = fauna[[1]]$time_step)
plot_marlin(proc_mpa_spillover, max_scale = FALSE) +
labs(y = "Spawning Stock Biomass")
plot_marlin(proc_mpa_spillover, max_scale = FALSE, plot_var = "c") +
labs(y = "Catch")
We can use marlin::get_distance_to_mpas to measure the
euclidian distance between the centroid of each cell and the nearest MPA
cell centroid.
mpa_distances <- get_distance_to_mpas(mpa_locations = mpa_locations, resolution = resolution, patch_area = patch_area)
mpa_distances |>
ggplot(aes(x,y,fill = distance_to_mpa_edge)) +
geom_tile() +
geom_tile(aes(x,y, color = mpa), size = 1.5) +
scale_fill_gradient2(low = "darkblue", high = "green", mid = "white", midpoint = 0)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
Distance of each cell to the nearest MPA edge.
mpa_distances |>
ggplot(aes(x,y,fill = total_mpa_distance)) +
geom_tile() +
geom_tile(aes(x,y, color = mpa)) +
scale_fill_viridis_c(name = "MPA gravity",option = "magma", direction = -1)
MPA gravity per cell.
We can then plot the gradients of various metrics in space and as a function of MPA distance.
We can see first off that the model generates “fishing the line” behavior where the fishing fleet is concentrated more around the borders of the protected area.
conservation_outcomes <- proc_mpa_spillover$fauna
fishery_outcomes <- proc_mpa_spillover$fleets
a = fishery_outcomes |>
filter(step == max(step), effort > 0) |>
left_join(mpa_distances, by = c("x","y","patch")) |>
ggplot(aes(distance_to_mpa_edge, effort, color = factor(step))) +
geom_point() +
facet_wrap(~critter, scales = "free_y")
b = fishery_outcomes |>
mutate(effort = if_else(effort == 0, NA, effort)) |>
filter(step == max(step)) |>
ggplot(aes(x, y, fill = effort)) +
geom_tile()
a / b
Distribution of fishing effort in space post-MPA.
The complex interactions of the life history and fishing fleet strategies mean that this “fishing the line” behavior results in a linear biomass gradient with distance from the MPA for the tuna, but actually a reverse trend for the shark where biomass is depressed right by the border due to fishing the line, and then increases with distance.
Note that we see clear “fishing the line” and spillover gradients even though this MPA has had almost no impact on total biomass or
conservation_outcomes |>
filter(step == max(step)) |>
group_by(step, x, y, patch, critter) |>
summarise(biomass = sum(b)) |>
left_join(mpa_distances, by = c("x", "y", "patch")) |>
ggplot(aes(distance_to_mpa_edge, biomass, color = factor(step))) +
geom_vline(xintercept = 0) +
geom_jitter() +
facet_wrap( ~ critter, scales = "free_y")
#> `summarise()` has grouped output by 'step', 'x', 'y', 'patch'. You can override
#> using the `.groups` argument.
Biomass of each species as a function of distance from MPA border. Negative distance means inside the MPA.