03_marlin_paper.qmd

---
format:
  pdf:
    toc: false
    fig-pos: 'H'
    mainfont: Times New Roman
    sansfont: Times New Roman
    fontsize: 12pt
    keep-tex: true
    fig-format: png
  html: default
  docx:
    reference-doc: template.docx
execute:
  echo: false
  message: false
  warning: false
  cache: false
csl: apa.csl
bibliography: references.bib
header-includes:
  - \usepackage{setspace}\doublespacing
  - \usepackage{lineno}\linenumbers
  - \usepackage{bm}
editor: 
  markdown: 
    wrap: 80
---

```{r setup}
#| echo: false
#| message: false
#| include: false

source(here::here("00_setup.R"))

theme_set(theme_classic(base_size = 10, base_family = "Helvetica") + theme(strip.background = element_blank()))

if (!file.exists(here("data","sci_to_com.csv"))){
  sci_to_com <- taxize::sci2comm(unique(marlin_inputs$scientific_name),simplify= TRUE,db='ncbi') %>% 
  bind_rows(.id = "critter") %>% 
  pivot_longer(everything(),names_to = "critter",values_to = "common_name") %>% 
  mutate(common_name = stringr::str_to_title(common_name))

  write_csv(sci_to_com, file = here("data","sci_to_com.csv"))

} else{
 sci_to_com <-  read_csv(file = here("data","sci_to_com.csv"))

}

print(results_path)
results_files <- list.files(file.path(results_path))

rdses <- results_files[str_detect(results_files,'.rds')]

walk(rdses,~assign(str_remove(.x,".rds"), read_rds(file.path(results_path,.x)), envir = .GlobalEnv)) #load and assign all .rds objects in the appropriate results folder

small_width <- 80 # figure width in mm

big_width <- 180 # width in mm


```

## Title Page

Title: Simulating Benefits, Costs, and Trade-offs of Spatial Management in
Marine Social-Ecological Systems

Running Title: Spatial Social-Ecological Simulation

Author Names: Daniel Ovando^1,2^, Darcy Bradley^3,4,5^, Echelle
Burns^3,4,5^,Lennon Thomas^3,4,5^, James Thorson^6^

Affiliations: ^1^ School of Aquatic and Fishery Sciences, University of
Washington,1122 NE Boat St Box 355020,Seattle, WA, USA ^3^ Marine Science
Institute, University of California, Santa Barbara, CA, USA ^4^ Bren School of
Environmental Science & Management, University of California, Santa Barbara, CA,
USA ^5^ Environmental Markets Lab, University of California, Santa Barbara, CA,
USA ^6^ Habitat and Ecological Processes Research Program, Alaska Fisheries
Science Center, NOAA, NMFS, NOAA, Seattle, WA, USA

Corresponding Author: Daniel Ovando
[dovando\@iattc.org](mailto:dovando@iattc.org){.email}

Present Address: \^2 Inter-American Tropical Tuna Commission,8901 La Jolla
Shores Drive, La Jolla, CA 92037, USA

Author Contributions: Methods and code written by DO and JT. Results produced by
DO. Data collected by EB and LT. All authors contributed to the
conceptualization and writing of the manuscript.

Conflict of Interest Statement: The authors declare no conflicts of interest.

\newpage

## Abstract

Designing effective spatial management strategies is challenging because marine
ecosystems are highly dynamic and opaque, and extractive entities such as
fishing fleets respond endogenously to ecosystem changes in ways that depend
upon ecological and policy context. We present a modeling framework, `marlin`,
that can be used to efficiently simulate the bio-economic dynamics of marine
systems in support of both management and research. We demonstrate `marlin’s`
capabilities by focusing on two case studies on the conservation and food
production impacts of marine protected areas (MPAs): a coastal coral reef and a
pelagic tuna fishery. In the coastal coral reef example, we explore how
heterogeneity in species distributions and fleet preferences can affect
distributional outcomes of MPAs. In the pelagic case study, we show how our
model can be used to assess the climate resilience of different MPA design
strategies, as well as the climate sensitivity of different fishing fleets. This
paper demonstrates how intermediate complexity simulation of coupled
bio-economic dynamics can help communities predict and potentially manage
trade-offs between conservation, fisheries yields, and distributional outcomes
of management policies affected by spatial bio-economic dynamics.

## Keywords

Bio-economic modeling; Conservation planning; Fisheries; Marine protected areas; Marine spatial planning; Spatio-temporal modeling

## Table of Contents

1.  Introduction
2.  Methods
    1.  Model Summary
    2.  Fleet Dynamics
    3.  Population Model
    4.  Case Studies
3.  Results
    1.  Coral Reefs
    2.  Pelagic Systems
4.  Discussion
    1.  Insights from Case Studies
    2.  Putting `marlin` in Context
    3.  General Recommendations for Use
5.  Acknowledgements
6.  Data Availability Statement
7.  References
8. Figure Legends

\newpage

## Introduction

Communities around the world are increasingly looking to spatial strategies for
managing marine ecosystems. For example, the "30x30" movement calls for
implementation of marine protected areas (MPAs), a form of spatial management,
across 30% of the world's oceans by the year 2030 [@grorud-colvert2021 and
references therein]. Recent agreements such as the United Nations (UN) agreement
on the conservation and sustainable use of marine biodiversity of areas beyond
national jurisdiction (BBNJ) calls for expanded area-based management in the
high seas. However, designing spatial management strategies to achieve desired
objectives -- which may include recovery and resilience of overfished species,
increased food production and economic well being, and equitable distribution of
benefits - is not straightforward. Different species have different resilience
to fishing, distributions in space and time, and value to fishing fleets;
fishing fleets themselves can have varying ranges of incentives and reactions to
policy structures. These dynamics pose challenges even for single-species
assessment and management strategies, which are only amplified when we consider
management policies designed around multiple species, fleets, and
spatial features in the oceans [@field2006a].

To illustrate, policies such as MPAs designed for both conservation and food
production must consider factors such as the optimal size and placement of a
protected area network given a wide array of life history, species
distributions, exploitation levels, fleet dynamics, and policy-dependent
behavior [@reimer2017], all of which may fluctuate over time, particularly given
the impacts of climate change. Efforts to effectively design spatial management
strategies such as MPAs are further constrained by a lack of empirical evidence
describing the size [@ovando2021a] and time required [@nickols2019] for MPAs to
produce substantial effects across a range of social-ecological settings
[@mcclanahan2021].

Many models have been developed to support the theory and design of spatial
management strategies (see @fulton2015 for a review).
However, these models have tended to be either highly complex and tactical
models designed for use in a specific location, or to be extremely stylized
representations intended to provide theoretical insights with less relevance for
specific applications. This lack of accessible models capable of representing reasonable amounts of complexity presents a challenge to stakeholders charged with marine spatial
planning exercises. It also presents a barrier to the scientific community,
wherein comparisons across spatial management simulation studies are clouded by
discrepancies in underlying model structures beyond differences in the core
phenomenon in question.

To help address this challenge, we present a bio-economic modeling tool
called [`marlin`](https://danovando.github.io/marlin/) that allows for efficient
simulation of spatio-temporal biological and economic dynamics. `marlin` allows
users to simulate the impacts of marine management policies across a range of
species targeted by various fleets across heterogeneous and dynamic seascapes.
This model can help users seek Pareto-optimal solutions that produce win-wins or
minimize trade-offs across multiple management objectives, such as food
production and conservation [@lester2013].

Here we present the core methods and functionality of `marlin`, and demonstrate
its use in two case-study applications; a coral reef system and a pelagic
system. Our results show how interacting ecological, economic, and design
attributes can affect the degree of benefits, costs, and trade-offs between
conservation and food production outcomes of MPAs. More broadly, this paper
demonstrates the functionality of our model and the critical need for
considering sufficiently realistic coupled ecological and economic dynamics in
policy evaluations.

## Methods

The model presented here focuses on representing the dynamics of heterogeneous
habitats and movement dynamics, along with the simultaneous impacts of fishing
fleets across multiple species (frequently called "technical interactions"). In
its general form, `marlin` simulates the behavior of populations of spatially
explicit biologically independent animal populations with age and subsequent
size structure affected by one or more fishing fleets in time steps and a
spatial resolution specified by the user.

The front end of the software package accompanying this paper,
[`marlin`](https://danovando.github.io/marlin/), is written in R
[@rcoreteam2021] to facilitate use, while the underlying population model is
written in C++, integrated through the `Rcpp` package [@eddelbuettel2018]. Users on standard computers should be able to simulate one age-structured population distributed across two-dimensional spatial surface represented by a
ten by ten grid of cells over 20 years in fractions of a second. It is very important to note that the parameters of
this model cannot be fit to data directly within the `marlin` package; users
must set model parameters themselves based on externally available data and
their best judgement.


While `marlin` does not simulate species interactions, fishing fleets in
`marlin` are capable of targeting and affecting multiple species simultaneously.
This allows the model to simulate processes such as fisheries bycatch or effort
displacement in a way that accounts for how fleet behavior might affect multiple
species in the system. The ability to track the impacts of fishing fleets across
multiple species simultaneously is particularly important as very few fishing
fleets are truly single species.

We define a few commonly used terms here. *Yields* refers to the volume of catch
from fishing activity. *Spawning stock biomass* refers to the total biomass of
reproductively mature fish in the water, a function of numbers, weight,
fecundity, and sexual maturity at age. We measure the size of the population by
the ratio of the *spawning stock biomass* (SSB) in a given time step relative to
*unfished spawning stock biomass* (SSB0), SSB/SSB0. An SSB/SSB0 value of 1 means
that the population is unfished, a SSB/SSB0 value of 0 means the population is
extinct.

As a demonstration of the use of this modeling tool, we explore two general
situations

1.  Trade-offs and distributional outcomes for food production and multi-species
    biomass in a coastal coral reef fishery
2.  Implications of climate-driven range shifts for MPA design in a pelagic
    fishery

We chose these two examples to illustrate the use of `marlin` in contrasting
systems in which spatial management strategies such as MPAs are increasingly
considered.

Below we provide a summary of the `marlin` model, as well as details of the case
studies.

### Model Summary

`marlin` simulates the dynamics of one or more species, currently best
representing fish-like species, using age-structured population dynamics. Ages
are then converted to lengths using the von Bertalanffy growth equation with
log-normally distributed variation in the length at age. Each time step, species
move throughout the simulated seascape using both diffusion and "taxis" (active
movement towards preferred habitat), experience natural and potentially fishing
mortality, and potentially spawn using one of the possible forms of
density dependence implemented in the model, with the ability to include
auto-correlated stochastic deviations in the amount of offspring produced,
generally called "recruitment deviates".

These species can be caught by fishing fleets. A fishing fleet is defined in
`marlin` by a set of fishers that have the same fishing skill, prices, and
contact selectivity [@sampson2014] for individual species (each of which we
denote as a *métier*). For example, both a longline and purse-seine fleet may
capture bigeye (*Thunnus obesus*, Scombridae) and skipjack (*Katsuwonus pelamis*, Scombridae) tunas, but the longline fleet can
be made much more likely to catch larger bigeye than skipjack, and *vice versa*.
It is important to note that contact selectivity reflects the ability of the
fishing method in question to capture fish of different sizes that come into
"contact" with the gear. The contact selectivity of each of the specified
*métiers* will then interact with the distribution of fish sizes and fishing
effort in space simulated by `marlin` to produce a net "population" selectivity,
which may differ from the individual contact selectivities of each of the
*métiers* [@waterhouse2014; @sampson2014]. Each fishing fleet distributes its
fishing effort in space according to a specified spatial allocation rule (see
*Spatial Allocation of Effort* section), for example in proportion to
profit-per-unit-effort, conditional on management policies in places such as
quotas and/or the presence of any spatial restrictions such as MPAs.

Unconstrained by management, the total amount of effort exerted by each fleet
can follow one of two dynamics: open access or constant effort (see *Calculating Total Effort*). Under constant effort, the total amount of effort of each fleet
is fixed over time, with the possible exception of attrition due to MPA
placement. Under open access, the total amount of effort in a given time step is
a function of the profitability of the fleet in the previous time step, evolving
until a bionomic equilibrium of zero total profits is reached. Profitability is
a function of the volume and price of each species caught, as well as the cost
of the total amount of fishing effort per fleet and the travel costs as a
function of distance from a port (see *Fleet Dynamics*).

The dynamics of each fleet can be modified by management in a variety of ways.
Managers can impose size limits for individual species within each fleet. They
can also impose total catch quotas for each species in each time step. When
quotas are activated, if the total catch across all fleets for a given species
under the model's effort dynamics would exceed the allowable quota for that
species, the total effort for each fleet is decreased proportional to its
contribution to the total catch until the quota is satisfied. As an alternative
to catch quotas, managers can set an effort cap per fleet, which prevents effort
from exceeding a given amount under open access dynamics, though effort may be
reduced below this cap if required by the profit equation (i.e. the fleet can
choose to fish less than the quota or effort cap). `marlin` also allows users to specify closed fishing seasons for one or more
species in the system. Lastly, managers can specify locations of no-take MPAs, which
can change in size and location if needed. When MPAs are implemented, fishing
effort that used to operate inside the MPAs can either leave the fishery, or be
redistributed outside the MPA (the default behavior).

On the biological side, at a minimum users must for each species being simulated supply the common or scientific name of the species in question, a
measure of the level of fishing intensity the species is experiencing at the
start of the simulation, and the diffusion rates for adult and larvae. Given a
common or scientific name, the model will then input default life history
parameters based on FishLife [internet connection required, @thorson2020],
though users are encouraged to check these values and supply their own life
history parameter values when possible.

On the fleet side, users must for each fleet they wish to simulate at a minimum
specify the contact selectivity curves for each species caught by the fleet, as
well relative price per unit weights greater than 0 for any target species.
Species caught by but not targeted by the fleet can be represented by prices of
0 (not targeted) or below 0 (actively avoided).

### Fleet Dynamics

Each fishing fleet (*f*) generates catches or yield (*Y*), revenues (*R*), and
costs (*C*) from fishing individual species (*s*) that it targets in a given
time step (*t*) and patch (*p*). The fleet then makes decisions around fishing
locations and intensity, subject to regulatory constraints, based on the total
profitability across all species in space and time.

Revenues for each fleet in a time step are a function of the total amount of
each species caught and its price ($\Pi$). The amount caught is a function of
the contact selectivity at age of each species for each fleet ($\alpha$), the
fishing efficiency of the fleet for that that species (also called
"catchability", *q*), the amount of fishing effort of the fleet in question in
that patch in a given time step (*E*), and the total instantaneous fishing
mortality (*u*) at age (*a*) (including all other fleets) for that particular
species in that patch and time step.

$$
 R_{t,p,f} = \sum_s^{N_s} \sum_a^{N_a} \Pi_{s,f}\frac{\alpha_{a,s,f}q_{s,f}E_{t,p,f}}{u_{t,p,s,a}} \times Y_{t,p,s,a}
$$ {#eq-revenue}

Total catches or yield *Y* are calculated through the Baranov equation
[@baranov1918], which accounts for the total instantaneous mortality (both
fishing and natural, *z*) and divides the total amount of the population biomass
(*b*) killed between natural (*m*) and fishing (*u*) sources, with the amount of
biomass killed through fishing called "yield".

$$
Y_{t,p,s,a} = \frac{u_{t,p,s,a}}{z_{t,p,s,a}} \times b_{t,p,s,a} \times (1 - e^{-z_{t,p,s,a}})
$$ {#eq-yield}

$$
u_{t,p,s,a} = \sum_{f=1}^{N_f}\alpha_{a,s,f}q_{s,f}E_{t,p,f}
$$ {#eq-u}

Contact selectivity at age $\alpha$ is modeled through either a logistic form or
a dome-shaped form. The logistic form is based on the lengths *l* at which 50%
of individuals are selected by the fishing gear, $l^{sel}$, and $\delta$ which
is the difference between the length at 50% selectivity and the length at 95%
selectivity.

$$
\alpha_{a,s,f}=\frac{1}{(1 + e^{-log(19)\times\frac{l_{a,s} - l^{sel}_{s,f}}{\delta^{sel}_{s,f}}})}
$$ {#eq-selectivity}

We approximate the dome-shaped form as a normal distribution with mean $l^{sel}$
and standard deviation $\sigma$. The normal density function is re-scaled such
that the selectivity is 1 at $l^{sel}$.

$$
\alpha_{a,s,f}= \frac{1}{\sigma_{s,f} \sqrt{2\pi}}e^{-0.5 \left( \frac{l_{a,s} - l^{sel}}{\sigma_{s,f}}  \right)^2}
$$ {#eq-dome}

Lastly *z* is total mortality, the sum of fishing mortality (*u*) and natural
mortality (*m*).

$$
 z_{t,p,s,a} = u_{t,p,s,a} + m_{s,a}
$$ {#eq-z}

This means that each species experiences a total mortality at age in a time step
in a patch, individual fractions of which are portioned off as catches and
subsequently revenues for each fleet. Given revenues *R*, we then calculate
profits $\phi$ based on revenues and costs. Costs *C* are calculated as a
function of base costs per unit effort ($\gamma$) as well as potential
additional costs per unit effort of fishing in particular patches ($\eta$) for
fleet *f*. $\beta$ allows for the cost per unit effort effort to scale
non-linearly. Travel costs ($\eta$) are calculated based on the Euclidean
distance of each patch to the nearest port of a given fleet, and users can
specify any number from zero to the number of patches of port locations. When no
ports are specified, travel costs are zero.

$$
\phi_{t,p,f} = R_{t,p,f} - C_{t,p,f}
$$ {#eq-profits}

$$
C_{t,p,f} = \gamma_f  \left( E_{t,p,f}^{\beta_f} + \eta_{f,p} E_{t,p,f} \right)
$$ {#eq-travel-costs}

#### Calculating Total Effort

`marlin` allows for two general modes of effort. The simplest is "constant
effort", in which the total effort of each fishing fleet remains constant over
time.

The more complex option is "open access". Under open access, the total amount of
effort in a given time step for fleet *f* is a function of the profitability of
that fleet in the proceeding time steps in which fishing was open, where
$\theta$ controls the responsiveness of fleet *f* to the log of the ratio of
total revenues *R* to total costs *C*, and approximates the proportional change
in effort in response to a one unit changes in the log revenue to cost ratio.

$$
E_{t+1,f} = E_{t,f} \times e^{\theta_f log(R_{t,f} / C_{t,f})} 
$$ {#eq-oa}

This is essentially a Gompertz model for fishing effort, that has been used for
other theoretical studies of fishing effort dynamics [@thorson2013].

#### Tuning Fishing Fleet Parameters

The degree of fishing pressure exerted by a given set of fishing fleets on each
species is a function of a range of parameters including the total amount of
effort (*E*), the contact selectivity ogives ($\alpha$), fishing cost (*C*), the spatial distribution of the
species affected by the fleet, the relative prices across species ($\Pi$), and
the catchability coefficient of each fleet for each species ($q$). `marlin`
provides a tuning option to help users achieve desired biological outcomes from
their fleets.

Users can tune their fleets in one of two ways, conditional on the underlying
population dynamics of the species in question. First, they can specify a target
exploitation rate *u* for each species in their simulation. Taking all the other
parameters of the model as given, `marlin` then adjusts the catchability
coefficients $q_{f,s}$ for each fleet *f* and species *s* such that the
equilibrium exploitation rate for each species matches the desired level.
Second, they can specify a target total spawning stock biomass under fishing
divided by total unfished spawning stock biomass, and the model will adjust the
catchability coefficients $q_{f,s}$ for each fleet *f* and species *s* such that
the equilibrium ratio of fished to unfished spawning biomass for each species
matches the desired level.

Users can also tune the fleet dynamics by specifying the ratio of costs to
revenues. Price data for use in the model can be obtained relatively simply
through literature reviews, market surveys, or local experts. However, cost
parameters are more complicated, as translating say cost per day of fishing into
the same representation of fishing effort used in the model is not straight
forward to accomplish. As an alternative, users can specify an equilibrium
cost-to-revenue ratio for the fleet, essentially the profit margins of the fleet
in question, and `marlin` will tune the cost parameter to achieve this desired
cost to revenue ratio given the other parameters in the model.

#### Spatial Allocation of Effort

Each fishing fleet decides how to allocate its effort in space based on one of
four possible spatial allocation strategies. The ideal free distribution (IFD)
is the standard method for distributing fishing fleets in spaces (see
@gillis2003 and references therein). However, analytical solutions to the IFD
present a number of complications for our model. The IFD for a single fleet
would commonly be modeled as a Nash-equilibrium based on the actions of each of
the individual fleet conditional on the actions of all other fleets. While
possible to implement, this would slow down our model runs to the point of
practically preventing large-scale evaluation of the spatial management
policies.

As such we explored a series of "next best" fleet distribution algorithms. While
not the IFD in any individual time step, over time they start to approximate the
IFD, as each assumes that fleets base their decisions for the current time step
on the outcomes in the prior timestep, meaning that the impacts of the actions
of other fishing fleets are eventually accounted for. The timeline required for the
fleets to reach an equilibrium distribution will vary and depend on factors such
as the population dynamics of the species in question and the number and degree of
competition across fleets. Users should explore different simulation times to
ensure that any results they wish to use are not simply a reflection of the
fleet dynamics fluctuating on their way to an equilibrium condition.

The four possible fleet distribution algorithms are

1.  Revenue per unit effort (RPUE): The fleet distributes itself in space based
    on the realized revenue per unit effort in each patch in the preceding time
    step
2.  Revenue: The fleet distributes itself in space based on the realized total
    revenue in each patch in the preceding time step
3.  Profit per unit effort (PPUE): The fleet distributes itself in space based
    on the realized profit per unit effort in each patch in the preceding time
    step
4.  Profit: The fleet distributes itself in space based on profit in each patch
    in the preceding time step

Revenue based spatial distribution is not likely to be very realistic; in
general we would expect fishing fleets to respond to profits on some core level.
However, due to the complexity of parameterizing cost functions, fleet dynamics
are often evaluated based on yield or revenue alone, and so we include those
scenarios here to allow users to evaluate the potential implications of this
choice. The decision on whether to allocate the fleet based on absolute or
relative (per unit effort) metrics is more complex. When effort represents the
actions of separate and individual fishing actors (e.g. independent fishing
vessels), a per-unit-effort strategy, in which fishers distribute themselves
based on the expected catch of their individual efforts, may be more realistic
[@hilborn1987]. Conversely, a system defined by a sole owner seeking to maximize
total profits might be better represented by a fleet model based on total
profits.

### Population Model

The underlying population model used is an age structured single-species model
in the manner of @ovando2021a. The population model requires many parameters.
However, if the user supplies either a scientific or common name for the species
in question (scientific preferred), the model will supply default values for
that species based on the values reported in FishLife [@thorson2020]. FishLife provides estimates of core life history
parameters for fish species based on a model integrating published values and
evolutionary connections. FishLife provides reasonable default values for given
species, but these default values should be not be taken as definitive and users
should check the default values and input best available parameters specific to the stock in question if they
wish to best represent the dynamics of their specific system. 

#### Movement

`marlin`'s movement dynamics are based on a continuous-time Markov chain (CTMC),
as described in @thorson2021a. Within this framework, the model allows for
movement to be broken down into three components of *advection* (drifting with
currents), *taxis* (active movement towards preferred habitat), and *diffusion*
(essentially remaining variation in movement not explained by advection or
taxis). For now, `marlin` focuses just on the diffusion and taxis components of
this model, assuming that advection is zero, though future extensions could
incorporate advection vectors from oceanographic models. In this way, `marlin`
allows users to run anything from a simple Gaussian dispersal kernel up to a
system governed by species that passively diffuse out from a core habitat
defined by a dynamic thermal range. The model currently focuses on diffusion and
taxis, which allows for general representation of animals following physical or
oceanographic features, but research on the incorporation and importance of
advective forces would be of value going forward.

We provide a brief overview of the the general CTMC method here (see
@thorson2021a for a detailed description). Under this framework, movement of
individuals from each patch to each other patch in the system in a given
timestep *t* for life stage *a* of species *s* is defined by a movement matrix
$\pmb{M}_{t,s,a}$. $\pmb{M}_{t,s,a}$ is calculated as a function of diffusion
$\pmb{D}$ and taxis $\pmb{\tau}$ matrices scaled by the width of the time step
(e.g. one year) $\Delta_{t}$ and the length of the edge of each patch (e.g. one
kilometer) $\Delta_d$ specified by the user. This parameterization allows users
to set the effective area of the spatial domain through two avenues; the number
of patches, which effectively scales the resolution of the model, and the area
of each patch, which scales the spatial extent of the simulation.

The individual components (*M*) of the movement matrix ($\pmb{M}$) are filled
based on an adjacency matrix, which defines whether two patches are both
adjacent and water (as opposed to land or another physical barrier), a diffusion
rate $D$ defined in units of area of a patch per unit of time, and a habitat
preference function *H* in units of length of a side of a patch per unit time.
For example, if we are defining the time units as years and the distance units
as kilometers, for a tuna $D$ might be 1,000 $\frac{KM^2}{Year}$. We then use
parameters $\Delta_{t}$ and $\Delta_{d}$ parameters to translate the diffusion
rate $D$ to match the time step and patch size used in a simulation. For
example, if we were to run a model on a monthly timestep given time units of
years, then $\Delta_{t} = 1/12 years$. If one square patch in the simulation has
an area of 100km^2^, then $\Delta_d = 10KM$. This "scale free" parameterization
means that appropriate value of $D$ can be identified for a species and then
set, regardless of the time step or patch size used in the simulation model
itself.

The taxis component of the movement process is a function of the difference in
habitat quality *H*. The habitat preference function itself can take any form
the user wishes. Exponentiating the difference in the habitat preference
function between patches turns the taxis matrix into a multiplier of the
diffusion rate *D*. As such, when creating habitat layers for simulation, users
can tune the scale of the habitat gradient function to result in realistic
multipliers of the diffusion rate. This parameterization ensures that the
off-diagonal elements of the movement matrix $\pmb{M}_{t,s,a}$ are all
non-negative, a requirement of the CTMC method.

$$ 
M_{p1,p2,t,s,a} = \begin{cases}
      = \frac{\Delta_{t}}{\Delta_{d}^2}De^{\frac{\Delta_t(H(p2,t,s,a) - H(p1,t,s,a))}{\Delta_d}} & \text{if p2 and p1 are adjacent}\\
     = -\sum_{p' \neq p1} M_{p1,p2,t,s,a} & \text{if p1 = p2}\\
     = 0 & \text{otherwise.}
\end{cases}
$$ {#eq-diffusion}

For both the diffusion and taxis matrices, we allow for the inclusion of
physical barriers to movement (i.e. land). Pairs of patches that are adjacent
but in which one or both patches are a barrier to movement are set as
non-adjacent. The CTMC model then produces movement dynamics that move around
barriers rather than over them. 

The movement of individuals across patches is then calculated by matrix
multiplication of the pre-movement vector of the number of individuals
($\pmb{n}$) of species *s* at age *a* in time step *t* across all patches *p*
times the matrix exponential of the movement matrix $\pmb{M}$

$$
\pmb{n}_{t+1,s,a} = \pmb{n}_{t,s,a}e^{\pmb{M}_{t,s,a}}
$$ {#eq-movement}

While this CTMC approach to movement simulation is to date not commonly seen in
the marine modeling literature, it has numerous advantages that warrant its
seeming complexity. First, the parameters of the model have interpretable
biological meaning (e.g. the diffusion rate *D*). Second, when only diffusion is
present, the model will generalize to the familiar dynamics of a Gaussian
dispersal kernel at whatever spatial and temporal resolution the simulation is
set to. Third, the taxis model allow for clearly parameterized active habitat
choices by species, allowing us to simulate preferences of species in space and
time efficiently. Lastly, the CTMC form has the advantage that its parameters
are directly estimable from real data. So, if provided with for example spatial
abundance data and a tagging study, users can estimate the diffusion and taxis
movement parameters in the same manner as @thorson2021a, and then pass their
estimated parameters to `marlin` for simulation (so long as the estimating
method uses the same functional form as the movement model in `marlin`).

#### Population Growth

For the population model, numbers *N* at time *t* for age *a* are a function of
growth, death, and recruitment

$$
N_{t,p,s,a} = \begin{cases}
      = BH(SSB_{t-1,p,s,a}) & \text{if $a = 1$}\\
     = N_{t-1,p,s,a-1}e^{-(z_{t-1,p,s,a-1})}, & \text{if $1< a < max(age)$}\\
     =  N_{t-1,p,s,a}e^{-(z_{t-1,p,s,a})} + N_{t-1,a-1}e^{-(z_{t-1,p,s,a-1})}, & \text{if $a = max(a)$}
  \end{cases}
$$ {#eq-pop}

where *BH* is the Beverton-Holt recruitment function [@beverton1957] and *SSB*
is spawning-stock-biomass. Per convention, the model allows for a "plus group",
wherein rather than tracking numbers of every possible age, individuals greater
than or equal to a given maximum age are grouped together.

Spawning stock biomass *SSB* is calculated by converting age to mean length at
age, calculating weight at age, maturity at age, and then calculating spawning
stock biomass as the sum of spawning potential at age in a given time step,
taking into account the potential for hyperallometry in the manner of
@marshall2021. Age is converted to length through the von Bertalanffy growth
equation given parameters asymptotic length ($l_\infty$), growth (*k*) and
theoretical age at length zero ($a0$) assuming log-normally distributed
variation *u* in the length at age with CV $\sigma_s$.

$$
  l_{a,s} = l_{\infty,s}\left(1 - e^{-k_s(a - a0_{s})}\right)e^{u_s}
$$ {#eq-len}

$$
u_s \sim N(0,\sigma_s)
$$ {#eq-sigma}

Users can manually supply a vector of of natural mortality at age (*m*). Or,
they can supply one value of natural mortality which is then converted into
mortality at age through one of two means. Under the default behavior, natural
mortality at age given a target mean mortality across all ages $m_s$ is
calculated using a length-inverse mortality function [@lorenzen2022].

$$
minv_{s_a} = (\frac{l_{s,a}}{l_{\infty,s}})^{-1}
$$ {#eq-minv}

$$
m_{s,a} = \frac{minv_{s,a}}{mean( minv_{s,a})}m_s
$$ {#eq-mata}

Alternatively, users can set mortality at age to be constant

$$
m_{s,a} = m_s
$$ {#eq-mata2}

Biomass *B* at age is then given by the weight at length equation governed by a
scaling coefficient $\Omega_s$ and an exponent $\Phi_s$ that controls the rate
at which volume scales with length

$$
  B_{a,s} = _s \times l_{a,s}^{wb_s}
$$ {#eq-weight}

The proportion of sexually mature individuals (*mat*) at a given age is then
calculated as a logistic function where $l_{mat}$ is the length at which on
average 50% of individuals are sexually mature, and $\delta_{mat}$ is the unit
of length beyond $l_{mat}$ at which on average 95% of fish are sexually mature.

$$
  mat_{a,s} = \frac{1}{\left(1 + e^{-log(19)\times\frac{l_{a,s} - lmat_s}{{\delta}mat_s}}\right)}
$$ {#eq-mat)}

Spawning stock biomass at time *t* is then calculated as a function of the
numbers at age, the maturity at age, and the weight at age raised by a parameter
$\gamma$. When $\gamma$ is greater than 1, the species is said to experience
hyperallometric fecundity, i.e. fecundity increases faster than weight.

$$
  SSB_{t,p,s} = \sum_{a=1}^{N_a}w_{a,s,t}^{\gamma_s} mat_{t,a} N_{t,p,s,a}
$$ {#eq-ssb}

#### Recruitment

Recruitment (i.e. the number of age 1 individuals entering the population)
follows Beverton-Holt dynamics parameterized around steepness (*h*) with
log-normally distributed recruitment deviates $\epsilon$. When steepness is one
recruitment is independent of spawning biomass. As steepness approaches 0.2
recruitment becomes a linear function of spawning biomass. `marlin` allows users
to specify a target unfinished spawning stock biomass ($SSB0$), which will be
achieved by tuning the total unfished recruitment ($r0$), given the remaining
life history parameters and independent of any characteristics of the fishing
fleets.

We allow for five variants in the timing of density dependent recruitment,
building off of @babcock2011 :

1.  Global density dependence: Density dependent recruitment is a function of
    the sum of spawning biomass across all patches, and recruits are then
    distributed according to habitat quality

$$
\begin{aligned}
N_{t,p,s,a = 1} = \left(\frac{0.8{\times}\sum_{p=1}^P{r0_{p,s}}\times{h_s}\times\sum_{p=1}^P{SSB_{t-1,p,s}}}{0.2\times{\sum_{p=1}^PSSB0_{p,s}} \times(1 - h_s)+(h_s - 0.2) \times{\sum_{p=1}^PSSB_{t-1,p,s}}}\right) \times \\ {r0_{p,s} / \sum_p^P{r0_{p,s}}\times \epsilon_{t,s}
}
\end{aligned}
$$ {#eq-dd2}

where $r0$ is is a vector of recruits under unfished conditions in a given
patch.

2.  Local density dependence: Density dependent recruitment occurs independently
    in each patch and recruits are retained in their home patch.

$$
n_{t,p,s,a = 1} = \left(\frac{0.8{\times}r0_{p,s}\times{h_s}\times{SSB_{t-1,p,s}}}{0.2\times{SSB0_{p,s}}\times(1 - h_s)+(h_s - 0.2)\times{SSB_{t-1,p,s}}}\right) \times \epsilon_{t,s}
$$ {#eq-dd1}

3.  Local density dependence then disperse: Density dependent recruitment occurs
    independently in each patch and recruits are then dispersed.

$$
  n_{t,p,s,a = 1} = \left(\frac{0.8{\times}r0_{p,s}\times{h_s}\times{SSB_{t-1,p,s}}}{0.2\times{SSB0_{p,s}}\times(1 - h_s)+(h_s - 0.2)\times{SSB_{t-1,p,s}}}\right)\times \pmb{d^l_s} \times \epsilon_{t,s}
$$ {#eq-dd3}

where **d^l^** is the recruitment movement matrix

4.  Post-dispersal density dependence: Larvae are distributed throughout the
    system, and then density dependent recruitment occurs based on the density
    of spawning biomass at the destination patch.

$$
  larv_{t,p,s} = SSB_{t-1,p,s} \times\pmb{d^l_s}
$$ {#eq-larvmove}

$$
n_{t,a = 1,p,s} = \left(\frac{0.8{\times}r0_{p,s}\times{h_s}\times{larv_{t,p,s}}}{0.2\times{SSB0_{p,s}}\times(1 - h_s)+(h_s - 0.2)\times{larv_{t,p,s}}}\right) \times \epsilon_{t,s}
$$ {#eq-dd4}

5.  Global density dependence allocated by spawning biomass: Density dependence
    is a function of the sum of spawning biomass across all patches, and
    recruits are then distributed according to the distribution of spawning
    biomass

$$
\begin{aligned}
n_{t,p,s,a = 1} = \left(\frac{0.8{\times}\sum_{p=1}^P{r0_{p,s}}\times{h_s}\times\sum_{p=1}^P{SSB_{t-1,p,s}}}{0.2\times{\sum_{p=1}^PSSB0_{p,s}}\times(1 - h_s)+(h_s - 0.2)\times{\sum_{p=1}^PSSB_{t-1,p,s}}}\right) \times \\ \frac{SSB_{t-1,p,s}}{\sum_{p=1}^PSSB_{t-1,p,s}} \times \epsilon_{t,s}
\end{aligned}
$$ {#eq-dd5}

Log-normal recruitment deviates are calculated with the potential for
autocorrelation defined with strength $\rho$

$$ 
\upsilon_{t,s} \sim \begin{cases}
       N(0,\sigma_{r,s}), & \text{if $t = 1$}\\
       \rho_s \upsilon_{t-1,s} + \sqrt{1 - \rho_s^2} N(0,\sigma_{r,s}), & \text{if $t > 1$}
  \end{cases}
$$ {#eq-recdev}

And log recruitment deviates are converted to raw units using the bias
correction factor

$$
\epsilon_{t,s} = e^{\upsilon_{t,s} - \sigma_{r,s}^2/2}
$$ {#eq-bias}

#### Reference Points

Fisheries management is often concerned with measuring stock status relative to
maximum sustainable yield (MSY) based reference points, though the exact level
of stock status relative to MSY reference points desired by societal objectives may vary widely. MSY based reference points present a
problem for a multi-fleet and spatial-temporal model such as `marlin`. MSY and
the fishing mortality rate that would produce MSY, $F_{MSY}$, are a function of
fishery selectivity. Fishery selectivity in this model can vary by fleet, and
species can be distributed unevenly in both space and time. This means that the
net effective fishing selectivity on a species can vary depending on the
dynamics at a given moment, making the definition of an equilibrium concept such
as MSY challenging [@berger2017 and references therein].

As such, we do not report MSY based reference points in the model by default.
There are many different strategies for estimating reference points in spatially
explicit systems [@kapur2021a]. We leave it to users to define and find relevant
reference points as required by their specific needs.

### Case Studies

We include two examples demonstrating how `marlin` can be used to support marine
spatial planning. In the first, we show how `marlin` can be used to compare the
total and distributional impacts of MPAs designed in a heavily fished coastal
coral reef ecosystem. In the second, we demonstrate how `marlin` can be used to
assess components of climate resilience of alternative MPA design strategies in
a pelagic system.

Each of the case studies contains too many parameters and options to be
succinctly presented in the text here. Readers should consult the accompanying code to view
the precise details of each simulation. Targeted applications must carefully
consider and document all decisions made around model parameters.

#### MPA Design Strategies

We make use of three potential "rule of thumb" MPA design strategies in our case
studies

1.  *Rate*: MPAs are placed based on the pre-MPA SSB/SSB0 weighted catch
    relative to the total catch in a patch. So, patches with high rates of catch
    of depleted species relative to total catch are prioritized.

2.  *Target Fishing*: MPAs locations are prioritized proportional to fishery
    catches. Patches with high total catches are prioritized over patches with
    low catches.

3.  *Spawning Ground*: MPAs are centered on the grounds of a known spawning
    aggregation. This strategy is only used in the coral reef case study.

In theory, the design of MPA networks can be optimized through the use of a
modeling framework, and depending on the validity of the model, this process is
likely to produce better outcomes than manually-designed strategies
[@rassweiler2012; @rassweiler2014]. However, designing optimal MPA networks
becomes increasingly difficult as the range of objectives and the complexity of
the model increase. Therefore, we focus here on the design and performance of these more rule of thumb design
strategies that may be more accessible to a wider range of users. We allow all
MPAs to be designed in a mosaic fashion in these examples, but users can easily
extend the analysis to compare outcomes between contiguous (MPA is made up of
one continuously connected block) and mosaic (MPAs can be separated in space)
MPA designs [@pons2022].

#### Coastal Coral Reef Fishery

In our coastal coral reef example, we model the dynamics of four tropical
Pacific species: a grouper (*Epinephelus fuscoguttatus*, Serranidae), a shallow-reef snapper
(*Lutjanus malabaricus*, Lutjanidae), a deep-reef snapper (*Pristipomoides filamentosus*, Lutjanidae),
and a reef shark (*Carcharhinus amblyrhynchos*, Carcharhinidae). The simulated groupers undergo
a mass migration to a spawning aggregation once per year, followed by the
sharks. Shallow-reef snappers stay in reefs closer to shore above a steep
drop-off year-round, while deep-water snappers stay in the deeper reefs past the
drop-off (@fig-coral).

These species are targeted by two different fleets. Fleet One primarily targets
the grouper and near-shore snapper populations, but will land any incidentally
captured sharks. Fleet One has a logistic selectivity pattern for all species,
as they retain any fish caught for consumption or sale. Fleet One is totally
dependent on fishing for their livelihood, meaning the local community takes
advantage of every possible opportunity to fish, and as such we model it as a
"constant effort" fishery. Due to having less efficient boats, Fleet One has a
higher cost per distance coefficient than Fleet Two. Fleet One's home port is
located near the site of the grouper spawning grounds.

Fleet Two is a more commercial fleet that primarily targets the snapper
populations. This fleet primarily sells their catch to local restaurants and
distributors where plate-sized fish are prized, and so for both snapper and
grouper Fleet Two has a dome-shaped selectivity pattern [@kindsvater2017]. While
plate-sized deep snapper are the primary target of Fleet Two, we model Fleet
Two's selectivity for deep snapper as logistic due to high levels of discard
mortality for deep-water snapper resulting from barotrauma. Fleet Two catches
groupers, though less than Fleet One, and receives no price for sharks due to
the requirements of a certification program through which they sell their
deep-water snapper. Accidental captures (bycatch) of sharks do occur, which
results in mortality. Fleet Two operates under open-access dynamics, as fishing
is not the only means of subsistence for this community; short-term effort
expands and contracts in response to profitability of the primarily
grouper-driven fishery. Fleet two's home port is located in the northwest corner
of the simulation space.

We used `marlin` to simulate the outcomes for both food production and
conservation for each of the species and both of the fleets as a function of
both MPA size and MPA design strategy. For this exercise, MPAs are placed with
perfect information and have no design constraints for continuity. We ran the
simulation for in quarter year time steps for 20 years ($\Delta_t = 1/4$) and
set the area of each patch to be 5km^2^ ($\Delta_p = \sqrt{5}$), using 144
patches for a total area of 2,000 KM^2^.

#### Pelagic Fishery

```{r}
blue_fauna <- blue_water_fauna_and_fleets$fauna

blue_fleets <- blue_water_fauna_and_fleets$fleets

blue_resolution <- sqrt(blue_fauna[[1]]$patches)

```

We model our pelagic case study loosely on the characteristics of the Western
and Central Pacific Ocean (WCPO) tuna fisheries. Note that this is an illustrative example only and simulated stock status, species distributions,  and projections presented here should not be interpreted as a indicative of the current or future state of the WCPO. We simulate trajectories of
`r n_distinct(names(blue_fauna))` species commonly caught in the region,
including both the highly abundant skipjack tuna and the heavily depleted
oceanic whitetip shark (*Carcharhinus longimanus*, Carcharhinidae) (@fig-blue-ssb). We use
publicly available data on catch-per-unit-effort of each of these species from
the WCPO as a very rough proxy for baseline habitat distributions, noting that
where possible, fishery-independent abundance indices would be preferable
(@fig-blue).

These pelagic species are caught by a longline fleet that primarily targets
large adult tunas such as bigeye and yellowfin (*Thunnus albacares*, Scombridae) for high-grade consumption, and a
purse-seine fishery that primarily targets skipjack tunas for bulk canning.
Contact selectivities were modeled as logistic for the longline fleet, and
dome-shaped for the purse-seine fleet. Both fleets operate under open-access dynamics with an effort
cap. The effort cap was set at the level of effort that resulted in the desired
levels of SSB/SSB0 for each species under open-access dynamics (@fig-blue),
intended to simulate a scenario where managers step in to prevent further
expansion of fishing effort in a fully developed fishery. For
forward-simulation, open-access dynamics can result in effort decreasing in
response to profitability, but cannot result in effort beyond the effort cap set
for each fleet.

For this exercise, we focused on using `marlin` to assess resilience of the
selected *Target Fishing* and *Rate* MPA design strategies to a climate-driven
range shift. Specifically, we simulate an extreme example where the centroid of
each population shifts northward at a rate of ~62km per year over a 20
year time horizon (@fig-blue-ssb). We designed MPA networks given the conditions
in the starting year, and then held that network constant over the years of the
experiment, running one simulation with and another without the climate-drive
range shift. We then compared the effects of this range shift on food production
and conservation outcomes from MPA networks designed based on the pre-range
shift world. We ran the pelagic simulation at a quarterly level ($\Delta_t$ =
1/4 year), and set the area of each cell to be roughly 97,000 KM^2^ across 144
patches each with a side length of roughly 311 KM, for a total area of
14e6 KM^2^, broadly commensurate with the area of the WCPO.



## Results

### Coastal Coral Reef Fishery

```{r}
#| label: coral-results

coral_mpa_experiments$mpas <- map(coral_mpa_experiments$results, "mpa")

coral_mpa_experiments$obj <- map(coral_mpa_experiments$results, "obj")


coral_results <- coral_mpa_experiments %>%
  unnest(cols = obj)

coral_fleet_results <- coral_results %>%
  pivot_longer(starts_with("fleet_"),
               names_to = "fleet",
               values_to = "fleet_yield") %>%
  group_by(prop_mpa, fleet, placement_strategy) %>%
  summarise(yield = sum(fleet_yield), biodiv = sum(unique(biodiv))) %>%
  group_by(fleet, placement_strategy) %>%
  mutate(delta_yield = yield / yield[prop_mpa == 0] - 1) |> 
  ungroup()

peak_coral_fleet_results <- coral_fleet_results |> 
  group_by(fleet) |> 
  filter(delta_yield == max(delta_yield)) |> 
  mutate(peak_yields = scales::percent(delta_yield))

tmp <- coral_results %>% 
  group_by(prop_mpa, placement_strategy) %>% 
  summarise(biodiv = sum(biodiv), yield = sum(yield)) %>% 
  group_by(placement_strategy) %>% 
  mutate(delta_yield = yield / yield[prop_mpa == 0] - 1,
         delta_biodiv = biodiv / biodiv[prop_mpa == 0] - 1) %>% 
  ungroup() %>% 
    mutate(placement_strategy= fct_relabel(placement_strategy, titler))



```

MPAs were capable of producing a range of positive and negative outcomes for
food security and conservation in the coral reef case study depending on the
design strategy used and the size of MPA implemented. Both of the MPA design
strategies were capable of increasing fisheries yield for
Fleet Two, up to a value of
`r peak_coral_fleet_results$peak_yields[peak_coral_fleet_results$fleet == "fleet_two"]`,
even when MPAs covered more than 50% of the simulated area. However Fleet One only benefited from MPAs under the *Spawning Ground*
design strategy, with a maximum increase of
`r peak_coral_fleet_results$peak_yields[peak_coral_fleet_results$fleet == "fleet_one"]`;
the *Target Fishing* design strategy produced a roughly linear decrease in
fishing yields as a function of increasing MPA size (@fig-coral-results A).

MPAs were uniformly beneficial to the spawning
biomass of all species under the *Target Fishing* design strategy, with the most
rapid increases in spawning biomass for the deep snapper population. The
*Spawning Ground* strategy primarily benefited the shallow snapper population,
producing little change in the grouper population and decreasing spawning
biomass of both the deep snapper and reef shark populations for MPA sizes
covering less than 50% of the simulated area (@fig-coral-results B).

In total, the *Spawning Ground* strategy was capable or providing net increases
in fishing yield (summed across both fleets) of up to roughly 7%, with positive
net yield impacts up to network size of nearly 60% of the area. However, even
MPAs covering 60% of the area only produced a maximum increase of roughly 5% in
total SSB/SSB0 when designed around the spawning ground. Conversely, the *Target
Fishing* design strategy was capable of producing a nearly 25% increase in total
SSB/SSB0 for the same MPA size, but at a much greater cost to the total food
production from the system's fisheries (@fig-coral-results C).

<!-- ~~These macro-level changes in total yield and conservation are the result of -->

<!-- interesting internal dynamics that marlin allows users to explore through output -->

<!-- tables at the level of individual species and fleets. The *Spawning Ground* -->

<!-- strategy increased yield to Fleet One when MPAs were relatively small, but was -->

<!-- able to maintain yields near *status quo* levels up to a very large MPA size to -->

<!-- the dependence of Fleet One on only a relatively small portion of the total -->

<!-- simulation area. However, *Spawning Ground* MPAs were unable to provide -->

<!-- increases in yield to Fleet Two. The *Spawning Ground* strategy produced the -->

<!-- most rapid conservation gains for the snapper population. Interestingly, due to -->

<!-- effort displacement from the spawning grounds onto the fishing areas inhabited -->

<!-- by the groupers the rest of the year, closure of the spawning ground had limited -->

<!-- impact on the grouper population until the MPA size exceeded 50% of the -->

<!-- simulation area. Effort displacement from the spawning grounds also resulted in -->

<!-- a small decrease in the deep snapper population for MPA sizes less than 50%. The -->

<!-- *Target Fishing* design strategy has similar yield outcomes to the *Spawning -->

<!-- Ground* strategy, but less dramatic negative impacts on Fleet Two, at the -->

<!-- expense of slightly fewer benefits for Fleet One. However, the *Target Fishing* -->

<!-- strategy provided more uniform conservation benefits to all species -->

<!-- (@fig-coral-results).~~ -->


### Pelagic Fishery

Yields of the purse-seine fleet were more resilient to the range shift under the
*Rate* design strategy, while the longline fleet had the opposite result. Under
the *Rate* strategy, longline yields were relatively stable across a large range
of MPA sizes under the status quo conditions, but declined rapidly under the
range shift conditions. For the purse-seine fleet, yields were more stable under
the range shift conditions under the *Target Fishing* design scenario, but
declined quickly as a function of MPA size under the status quo conditions
(@fig-blue-results A). The *Target Fishing* strategy produced better yield
outcomes for both fishing fleets under the range shift conditions, but the
*Rate* strategy performed best under the status quo.

The primary tuna species (bigeye, skipjack, and yellowfin) were most sensitive
to design strategy and climate scenario, with the *Target Fishing* strategies
producing for example rapid conservation gains for bigeye and skipjack under the
status quo, while the *Rate* strategy resulted in small net conservation losses
for both species until the MPA became extremely large. However, under the range
shift scenario the MPAs had little impact on many of the tuna populations until
the MPA size became extremely large, due to movement of the primary fishing
grounds outside of the current hotpots where the MPAs are placed based on status
quo conditions. The 20 year time horizon simulated here was not enough to
produce substantial gains for any of the shark species even with 100% closures
except for the more rapidly growing blue shark (@fig-blue-results B).

```{r}
#| include: false

blue_water_experiments$mpas <- map(blue_water_experiments$results, "mpa")

blue_water_experiments$obj <- map(blue_water_experiments$results, "obj")

blue_water_climate_experiments$mpas <- map(blue_water_climate_experiments$results, "mpa")

blue_water_climate_experiments$obj <- map(blue_water_climate_experiments$results, "obj")


static_results <- blue_water_experiments %>%
  unnest(cols = obj) %>%
  mutate(climate = FALSE)

climate_results <- blue_water_climate_experiments %>%
  unnest(cols = obj) %>%
  mutate(climate = TRUE)

blue_water_results <-static_results %>%
  bind_rows(climate_results) %>% 
  left_join(sci_to_com, by = "critter") %>% 
  select(-critter) %>% 
  rename(critter = common_name)


fleet_climate_results <- blue_water_results %>%
  select(placement_strategy,
         fleet_model,
         prop_mpa,
         climate,
         longline,
         critter,
         purseseine) %>%
  rename(`purse-seine` = purseseine) |> 
  pivot_longer(c(longline,`purse-seine`), names_to = "fleet", values_to = "yield") %>%
  pivot_wider(names_from = "climate", values_from =  "yield") %>%
  mutate(delta = `TRUE` - `FALSE`) %>%
  group_by(placement_strategy, prop_mpa, fleet) %>%
  summarise(delta = sum(delta))

fleet_climate_impacts <- fleet_climate_results %>%
  ggplot(aes(prop_mpa, delta, color = fleet)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_line() +
  facet_wrap(~placement_strategy) +
  scale_y_continuous(name = "Difference in Yield") +
  theme(axis.text.x = element_blank(), axis.title.x = element_blank())


critter_climate_results <- blue_water_results %>%
  select(placement_strategy,
         fleet_model,
         prop_mpa,
         climate,
         critter,
         biodiv) %>%
  pivot_wider(names_from = "climate", values_from =  "biodiv") %>%
  mutate(delta = `TRUE` - `FALSE`)


critter_climate_impacts <- critter_climate_results %>%
  ggplot(aes(prop_mpa, delta, color = critter)) +
  geom_line() +
  facet_wrap(~placement_strategy) +
  scale_x_continuous(name = "MPA",labels = scales::label_percent()) +
  scale_y_continuous(name = "Difference in SSB/SSB0")

fleet_climate_impacts / critter_climate_impacts

```


## Discussion

Marine ecosystems are driven by complex social-ecological dynamics. Communities
must often make decisions on how to manage these systems based on limited
empirical evidence. Modeling tools such as the one presented here can help users
answer scientific questions and design marine management policies informed by a
better understanding of sensitivities to key uncertainties.

### Insights from Case Studies

Our coastal coral reef example illustrates both the potential for MPAs to benefit conservation and
food production in these systems, and the potential for the same MPA to benefit
some fleets and species while harming others when fishing fleets affected by an
MPA do not share the same objectives and species are not uniformly distributed.
Our pelagic case study illustrates how `marlin` can be used to assess the
climate resiliency of spatial management strategies, and identify strategies
that best meet the needs of both current and future conditions. Results from
this type of work could be used to help prioritize communities and species at
particular risk to climate change impacts.

Holding constant other social-ecological variables in the coral reef example, the *Spawning Ground* MPA network
was able to provide more equitable yield outcomes across the two fishing
communities, whereas the *Target Fishing* MPA network only benefited Fleet Two (@fig-coral-results A). This is because Fleet One has two
primary fishing grounds; the spawning grounds, and the offshore area where the
deep snapper live. Fleet Two primarily fishes in the northern portions of the
simulation grid. Under the *Spawning Ground* strategy, while Fleet One quickly
losses fishing grounds on the spawning grounds, it can compensate for this by
moving offshore and fishing the deep snapper population harder, while also being
able to fish some of the spillover of snapper from the spawning ground closure.
This displacement of Fleet One's fishing effort is why biomass of deep snapper
actually declines under the *Spawning Ground* strategy.

Conversely, most of the catch in the coral reef fishery comes from the nearshore and deep
water snapper populations, which overlaps with both of Fleet One's fishing
grounds. The *Target Fishing* strategy then begins closing both the nearshore
and offshore snapper fishing grounds, resulting in too much of a loss in fishing
grounds for Fleet One to be offset by spillover from the MPAs. As a result of
primarily placing MPAs on Fleet One's fishing grounds, Fleet Two gains spillover
benefits at little cost to their fishing grounds until the MPAs reach their
fishing grounds once protection nears 100%.

By protecting both the spawning ground and the offshore areas, the *Target
Fishing* network provides conservation benefits to all of the species, in
contrast to the *Spawning Ground* strategy that only provides meaningful
conservation gains to the snapper population up until very large MPA sizes. Our
result that protection of a dedicated spawning ground did not provide
substantially greater conservation outcomes for the species using that spawning
ground (groupers and sharks) than an alternative design strategy is supported by
other modeling studies that show that the impacts of spawning ground protection
on conservation and yields may be highly variable, and that displacement of high
levels fishing effort from the spawning grounds can offset potential
conservation gains of the protection [@grüss2014]. That being said there is also
evidence for the benefits of spawning aggregation protection [@erisman2015] ,
and our results are not nearly sufficiently resolved to provide any general
statements as to the relative value of spawning aggregation protection relative
to other design strategies. Further research could for example alter both the distribution of species in space and the susceptibility of the species to fishing gear when aggregated. 

The two design strategies produced very different outcomes in terms of total
changes in yield and conservation outcomes for the coral reef case study
(@fig-coral-results C). Both MPA network design strategies were capable of
producing net "win-win" outcomes in which both food security and conservation. However, in general for the same size MPA network the *Spawning Ground* design strategy produced better food security outcomes but worse
conservation outcomes, and *vice* *versa* for the *Target Fishing* strategy. The
coral reef case study shows how the modeling framework presented here can help
stakeholders explore how different management strategies affect outcomes both in
total (@fig-coral-results C) and across species (@fig-coral-results A) and
fishing communities (@fig-coral-results B).

Turning to our Pelagic Ecosystem case study, the yield outcomes of the
purse-seine fleet were much more sensitive to the presence of a range shift than
the longline fleet, particularly under the *Target Fishing* design strategy.
This is because the purse-seine fleet primarily targets skipjack tuna, which in
this simulation are concentrated in a relatively narrow latitudinal band
(@fig-blue), and the purse-seine fleet makes up a large portion of the total
catch in the simulated fishery. So, the "Target Fishing" strategy starts by
closing off the main purse-seine fishing grounds, which while only small part of
the spatial domain of the model represents a large portion of the purse-seine
fleet's fishing grounds, resulting both in more rapid conservation gains and
fishery losses under the status quo species distributions. The *Rate* strategy
places more MPAs in areas that are of lesser importance fo the purse-seine fleet
but overlap more with species such as oceanic whitetip shark. Conversely, the
purse-seine fleet appears to do have better MPA yield outcomes under the range
shift scenario, not because of rebuilding of the skipjack population, but due to
future fishing grounds being essentially unprotected by smaller MPAs targeting
the current skipjack distribution ( @fig-blue-results ).

The conservation outcomes of the skipjack, bigeye, and yellowfin tunas were
among the most sensitive to range shifts, whereas even complete closure of the
region was not sufficient to see significant increases in severely depleted
species like oceanic whitetip sharks within the 20-year timeline of the
simulation. Future research could be conducted then to see what sorts of
timelines might be needed to see recovery of these species, and which design
strategies result in the fastest recovery at the lowest cost to other objectives
such as food security.


These case studies are intended to illustrate the capabilities of the model and
the importance of considering the bio-economic dynamics represented in `marlin`;
it is beyond the scope of this paper to provide broader conclusions about the
performance of MPAs under different contexts or a comprehensive comparison of
simulation and empirical results of MPAs around the world. That being said,
the kinds of dynamics resulting in our case studies are well supported by both
modeling and empirical studies. Our results support conceptual [@hilborn2004a;
@gaines2010] and empirical [@ban2019] evidence that under the right conditions
no-take MPAs can benefit fisheries and conservation. @rassweiler2014
demonstrated that the kinds of design choices presented in our case studies can
greatly drive MPA outcomes.

MPAs based around coral-reef style ecosystems have been extensively studied
around the world, with much of the empirical evidence of their performance
centered on demonstrating higher metrics such as biomass densities inside
protected areas relative to fished reference areas ["response ratios",
@lester2009]. Our model predicts a similar rapid increase of the simple ratio of
mean biomass inside MPAs relative to outside (Fig. S5, acknowledging that
designing a proper response ratio would require controlling for habitat
characteristics and MPA design criteria). However, despite producing clear
response ratios, the net conservation and fishery outcomes of our simulated
coral-reef ecosystem MPAs were not nearly as large or clear as the response
ratio results (@fig-coral-results); this result is supported by works such as
@ferraro2018 an @ovando2021a showing that response ratios alone may be a poor
indicator of the net causal impacts of MPAs at the population scale.

Studies such as @gilman2020, @abbott2012, @davies2018 and @pons2022 provide
empirical support for our case study results showing that MPAs can produce
trade-offs across the conservation of different species and fleets when the
spatial-temporal distributions and life history traits of affected species are
heterogeneous, and MPAs result in some degree of effort displacement or
concentration in the remaining fishing grounds. @hampton2023 supports our
pelagic case study result showing that both conservation and fishery impacts of
MPAs on the highly mobile species of the open oceans can be limited unless
protected areas are very large. @brown2018 supports our result showing how
heterogeneity in the behavior and objectives of fishing fleets sharing an
ecosystem can affect the magnitude and equity of fishery reform outcomes.
@davies2017 discusses the importance of considering the potential of climate
driven range shifts when assessing spatial management policies; our work builds
on this by allowing users to not just simulate species distributions but also
distribution of biomass and age composition of fish in space and time under climate
change.

### Putting `marlin` In Context

The role of this paper is not to conduct a review of the many tools for
spatial-temporal modeling available in the literature, each of which provide
useful functionality for specific applications. However, we highlight here the
specific gaps that we feel `marlin` fills in the modeling literature and in the
policy support toolbox using some selected publications. More end-to-end models
such as ATLANTIS [@audzijonyte2019] , OSMOSE [@shin2001], Ecopath with Ecosim
[@christensen2004], POSEIDON [@bailey2019], DISPLACE [@bastardie2013], or
SEAPODYM [@lehodey2008] are capable of representing tremendous amounts of
complexity, but can be time consuming to design and run. `marlin`
allows some of the realism of these more end-to-end models while being simpler
and faster to parameterize and run.

Tools such as virtualspecies [@leroy2016], STEPS [@visintin2020], RangeShifter
[@bocedi2014], and
[SMS](https://www.aquaveo.com/software/sms-surface-water-modeling-system-introduction)
can be efficiently constructed to model the dispersal and distribution of
species in space and time as a function of environmental covariates. However,
these typically require specifying a covariate-response curve without explicitly
acknowledging how this arises from habitat-specific movement and demography
(virtualspecies), or model movement using a dispersal kernel (RangeShifter,
STEPS) or a least-cost path algorithm (SMS, STEPS). By contrast, `marlin` uses a
continuous-time Markov chain movement model, which integrates multiple paths
(including their path-dependent probability based on intervening habitat types
and species preferences) while using scale-free parameters that can be measured
experimentally in laboratory or tagging studies. In addition, `marlin` allows
for simulating not only the distribution of the species but also biomass, age,
and length structure in space and time.

Simplified bio-economic models such as those used in @hastings2017, @sala2021,
and @cabral2019 can be applied at scale and provide analytically tractable
results, but as a result must abstract over many bio-economic dynamics that can
be important for more tactical applications. Our results demonstrate how nuances
in fleet dynamics and species distributions can dramatically impact MPA
outcomes. `marlin` allows for more realistic representations of spatial
social-ecological systems while maintaining processing speed.

`marlin` does not represent trophic interactions. Anthropgenic changes in
species abundance can result in trophic cascades.
However, studies such as @gilman2020, @ovando2021a, @bruno2019 , and
@malakhoff2021 found no clear signs of MPA driven trophic cascades within the
first decades of protection. Signals of trophic cascades may be masked by
variations in the direction and strength of species interactions driven by
environmental context [@liu2022], or may simply take longer to develop
detectable effects than the coverage of many time series of MPAs. While `marlin`
does not incorporate trophic interactions, what empirical evidence we have does
not suggest that management-mediated trophic cascades are so common and clear
that they must be incorporated into any credible multi-species simulation model. However,
research on the trophic impacts of spatial-temporal management actions is
clearly of value and models like Atlantis [@audzijonyte2019], Ecopath with
Ecosim [@christensen2004], and EASI-Fish [@griffiths2019] can help users explore
those kinds of trophic processes, in the manner of @baskett2007.

### General Recommendations for Use

The model presented here is designed to help users explore the impacts of different variables on policy
outcomes. But, that freedom means that users have a large
number of options at their disposal that they must decide on. Our recommended
strategy is for users to try and narrow down a list of parameters that they feel
are sufficiently "known" and another list of parameters that are highly
uncertain and / or contentious which the users feel may impact results. Where
possible, parameters from locally estimated stock assessment models can be used
to provide a foundation around which sensitivity analyses around specific
uncertainties of interest [@berger2017], though care should be taken interpreting population selectivity curves estimated from stock assessment as contact selectivity curves required as inputs to `marlin` [@sampson2014]. 

For example, a community seeking to a model a well-studied coastal coral reef
ecosystem might leave as fixed the habitat distribution (as represented by reef
locations) and general life history of the species in question (growth rates,
age at maturity, etc). From there, users may wish to test
the sensitivity of proposed MPA networks to key unknowns such as adult and
larval dispersal rates or the the economic incentives and contact selectivity of
the fishing fleets. As an example of this, we ran an alternative version of our
coral reef case study in which Fleet Two was assigned logistic selectivity for
the snapper and grouper species, rather than the dome-shaped selectivity
presented in our main results. We found that the impacts of MPAs on conservation
and food security were relatively insensitive to the form of contact selectivity
specified for Fleet Two (Fig.S1). This does not mean that knowing the "correct"
form of contact selectivity may not be extremely important in providing accurate
assessment results [@waterhouse2014], but rather that in this case misspecifying
the contact selectivity curve is projected to have little impact on the
simulated outcomes of MPAs, conditional on holding other parameters constant.

Given a set of model runs, the choice of whether to use those results
"strategically" or "tactically" will depend on the needs of the user. Since `marlin` is not fit to data directly it is more easily applied to strategic questions. The extent to which users are comfortable using the
outputs of `marlin` tactically will depend on the confidence they have in their
parameterization of the model relative to the precision policy-makers
require in order to make a decision. Other models may be better suited to address specific forms of
complexity. Particularly for more tactical applications, we would encourage
users to explore multiple modeling frameworks to help design policies that are
likely to be robust to many different kinds of complexities.

At this time, `marlin` assumes that all evaluated policies are perfectly
implemented; e.g. that MPAs are 100% no-take and perfectly enforced, that quotas
and closed seasons are respected, that there is no discard mortality, etc. In
reality no policy is perfectly implemented, and users should consider the extent
to which the policies they simulate are actually feasible to implement. There is
value though in being able to simulate and compare the outcomes of perfectly
implemented policies to isolate the concept of the policy itself from its
implementation.

Fish populations often exhibit variation in demographic traits across dimensions such as space, time, and sex. Models
such as Stock Synthesis [@methotjr.2013] capture these processes through the use
of "morphs". `marlin` does not currently allow for these kinds of processes
explicitly. Users should be
cautious interpreting simulation results from `marlin` for species in which
these sorts of dynamics are likely to be particularly prevalent. For strongly
sexually dimorphic species, we would recommend picking the sex most likely to
drive the outcomes of management policies, which is often females given the
general prevalence of eggs as the limiting reproductive material in marine
ecosystems. Users should also proceed with caution using `marlin` for species
with more complex reproductive biology, such as sex-changing fish
[@kindsvater2017].

## Conclusions

`marlin` complements the existing spatial marine modeling literature by allowing
scientists, decision makers, and other stakeholders to efficiently examine the
impacts of realistic bio-economic dynamics on academic and applied problems. We envision
`marlin` being applicable to research on dynamic ocean management, range shifts,
management strategy evaluation, policy interactions, and spatial stock
assessment. `marlin` can help researchers generate data for further testing of
the performance and design of spatially explicit integrated population models,
in the manner of @bosley2022. In addition, the process-based movement model used
in `marlin` can directly use empirical estimates of movement dynamics derived in
the manner of @thorson2021a, providing a link between empirical and simulation
based approaches to marine resource management that has been challenging to
implement in spatial simulations [@berger2017].

Fisheries models, assessments, management have often abstracted away many of the
spatial-temporal complexities of marine social-ecological systems [@ovando2021;
@berger2017]. The modeling framework described here can help facilitate the
science and application of spatial fisheries management by supporting the
simulation of different spatial-temporal dynamics to aid in testing of various
aspects of the marine resource management process.

Even the most complex marine model is a stylized cartoon of the true dynamics of
ocean ecosystems. However, for all their limitations, models
can help users understand factors that drive the performance of marine
management strategies. The goal of this tool is to empower people to design
policies based on evaluation of key uncertainties and trade-offs, and in doing
so support more effective and equitable marine resource management.

## Acknowledgements

This project was funded by the Waitt Foundation. We thank Dr. Alexa Fredston and
two anonymous reviewers for helpful revisions and comments of this manuscript.

## Data Availability Statement

The code, novel code, data, and materials needed to reproduce the results and
all aspects of this manuscript are publicly available at
https://github.com/DanOvando/marlin-paper, with supporting novel code available
at https://github.com/DanOvando/marlin/. Upon acceptance, all data, materials,
code, and novel code needed to reproduce the results and all aspects of this
manuscript will be publicly available via GitHub at
https://github.com/DanOvando/marlin-paper and
https://github.com/DanOvando/marlin/ and through figshare at
https://figshare.com/articles/preprint/marlin-paper/21843582

## References

::: {#refs}
:::

## Figure Legends


```{r}
#| label: fig-coral
#| fig-cap: "Distribution of fish spawning biomass (cell color, yellow = higher, blue = lower) each season (columns) under unfished conditions.  The x-axis represents longitude, y-axis latitude. Numbers show the approximate location of each fleet's (one or two) port. Columns indicate quarterly seasons that repeat each year. "

coral_fauna <- coral_fauna_and_fleets$fauna

coral_fleets <- coral_fauna_and_fleets$fleets

resolution <- sqrt((coral_fauna[[1]]$patches))

grid <- expand_grid(x = 1:resolution, y= 1:resolution) %>%
  mutate(patch = 1:nrow(.))


mpa_locations <- expand_grid(x = 1:resolution, y = 1:resolution) %>%
  mutate(mpa = TRUE)

coral_sim <- simmar(
  fauna = coral_fauna,
  fleets = coral_fleets,
  manager = list(mpas = list(
    locations = mpa_locations,
    mpa_year = 1
  )),
  years = 2
)




patch_biomass <-
  map_df(coral_sim, ~ map_df(.x, ~tibble(biomass = rowSums(.x$ssb_p_a), patch = 1:nrow(.x$ssb_p_a)), .id = "critter"), .id = "step") %>%
  mutate(step = marlin::clean_steps(step)) %>%
  left_join(grid, by = "patch")

patch_biomass <- bind_cols(patch_biomass,marlin::process_step(patch_biomass$step))


 fig_1 <- patch_biomass %>%
  group_by(critter, step) %>%
  mutate(sbiomass = biomass / max(biomass)) %>%
  ungroup() %>%
  filter(between(year,1,2)) %>%
  mutate(Season = season) %>% 
  ggplot()+
  geom_tile(aes(x,y, fill = sbiomass),show.legend = FALSE) +
  geom_text_repel(data = coral_ports, aes(x, y, label = fleet), color = "red") +
  facet_grid(critter~Season, labeller = labeller(critter = titler, Season = label_both)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  scale_fill_viridis_c() + 
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.line = element_blank())
 
 ggsave(
   filename = file.path(results_path, "fig_1.pdf"),
   fig_1,
   width = small_width,
   units = "mm"
 )
 
  ggsave(
   filename = file.path(results_path, "fig_1.png"),
   fig_1,
   width = small_width,
   units = "mm"
 )
 
 
 fig_1
 
```

```{r}
#| label: fig-blue-ssb
#| fig-cap: "Spawning biomass divided by unfished spawning biomass in the time period prior to implementation of MPAs for the pelagic case study. Simulated species include skipjack tuna, yellowfin tuna, bigeye tuna, shortfin mako (*Isurus oxyrinchus*, Lamnidae), swordfish (*Xiphias gladius*, Xiphiidae), albacore tuna (*Thunnus alalunga*, 	Scombridae), blue shark (*Prionace glauca*, Carcharhinidae), silky shark (*Carcharhinus falciformis*, Carcharhinidae), and oceanic whitetip shark."

grid <- expand_grid(x = 1:blue_resolution, y= 1:blue_resolution) %>%
  mutate(patch = 1:nrow(.))


mpa_locations <- expand_grid(x = 1:blue_resolution, y = 1:blue_resolution) %>%
  mutate(mpa = TRUE)

blue_sim <- simmar(
  fauna = blue_fauna,
  fleets = blue_fleets,
  years = 50
)
critter_biomass <-
  map_df(blue_sim, ~ map_df(.x, ~tibble(ssb = rowSums(.x$ssb_p_a), patch = 1:nrow(.x$ssb_p_a)), .id = "critter"), .id = "step") 

critter_biomass <- bind_cols(critter_biomass,marlin::process_step(critter_biomass$step)) %>%
  mutate(step = marlin::clean_steps(step)) %>%
  left_join(grid, by = "patch") %>% 
  group_by(critter, year, season,time_step) %>% 
  summarise(ssb = sum(ssb)) %>% 
  group_by(critter) %>% 
  mutate(depletion = ssb / ssb[time_step == min(time_step)])

fig_2 <- critter_biomass %>%
  ungroup() %>% 
  left_join(sci_to_com, by = "critter") %>%
  select(-critter) %>%
  rename(critter = common_name) %>%
  filter(time_step == max(time_step)) %>% 
  ggplot(aes(reorder(critter, depletion), depletion)) +
  geom_col() +
  scale_y_continuous(limits = c(0, 1), name = "SSB/SSB0", expand = expansion(mult = c(0.0,0.05))) +
  labs(x = "") + 
  coord_flip()

 ggsave(
   filename = file.path(results_path, "fig_2.pdf"),
   fig_2,
   width = small_width,
   units = "mm"
 )
 
  ggsave(
   filename = file.path(results_path, "fig_2.png"),
   fig_2,
   width = small_width,
   units = "mm"
 )
 

fig_2

```

```{r}
#| label: fig-blue
#| fig-cap: "Distribution of unfished spawning biomass in space under status quo (left column) and range-shifted (right column) conditions for blue-water simulation for each species. The x-axis represents longitude, the y-axis latitude."
#| fig-height: 8
#| fig-width: 6

blue_climate_sim <- simmar(
  fauna = blue_fauna,
  fleets = blue_fleets,
  habitat = future_habitat,
  years = length(future_habitat[[1]]),
  manager = list(mpas = list(
    locations = mpa_locations,
    mpa_year = 1
  ))
)

patch_biomass <-
  map_df(blue_climate_sim, ~ map_df(.x, ~tibble(biomass = rowSums(.x$ssb_p_a), patch = 1:nrow(.x$ssb_p_a)), .id = "critter"), .id = "step") %>%
  mutate(step = marlin::clean_steps(step)) %>%
  left_join(grid, by = "patch")




sq_habitat <- patch_biomass %>%
  left_join(sci_to_com, by = "critter") %>%
  select(-critter) %>%
  rename(critter = common_name) %>%
  filter(step == min(step)) %>%
  group_by(critter, step) %>%
  mutate(sbiomass = biomass / max(biomass)) %>%
  ungroup() %>%
  ggplot() +
  geom_tile(aes(x, y, fill = sbiomass),
            linewidth = .1,
            show.legend = FALSE) +
  facet_wrap( ~ critter, ncol = 1) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    legend.position = "top",
    axis.line = element_blank(),
    panel.spacing.y = unit(0, "lines")
  ) +
  scale_fill_viridis() +
  labs(title = "Status Quo")

rs_habitat <- patch_biomass %>%
  left_join(sci_to_com, by = "critter") %>% 
  select(-critter) %>% 
  rename(critter = common_name) %>% 
  filter(step == max(step)) %>% 
  group_by(critter, step) %>%
  mutate(sbiomass = biomass / max(biomass)) %>%
  ungroup() %>%
  ggplot()+
  geom_tile(aes(x,y, fill = sbiomass), linewidth = .1, show.legend = FALSE) +
  facet_wrap(~critter, ncol = 1) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme(axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        legend.position = "top",
        axis.line = element_blank(),
        panel.spacing.y = unit(0, "lines")) + 
  scale_fill_viridis() + 
  labs(title = "Range Shift")

range_shift_plot <- fit_3 <- sq_habitat | rs_habitat


 ggsave(
   filename = file.path(results_path, "fig_3.pdf"),
   fit_3,
   width = big_width,
   units = "mm"
 )
 
  ggsave(
   filename = file.path(results_path, "fig_3.png"),
   fit_3,
   width = big_width,
   units = "mm"
 )
 
 
range_shift_plot



```

```{r}
#| label: fig-coral-results
#| fig-cap: "Percent change in yield per fleet (A) and SSB/SSB0 (spawning biomass divided by unfished spawning biomass) per species (B) as a function of MPA size and MPA design strategy for the coastal coral reef case study. Percent change in total yield across both fleets (y-axis) and total SSB/SSB0 across all species (x-axis) and placement strategies (line color) (C). Color of points along each line in panel C indicates the percent of the simulation area in an MPA. All results reflect outcomes after 20 years of simulated MPA protection."
#| fig-height: 4
#| fig-width: 6


coral_fleet_frontier <- coral_results %>%
  pivot_longer(starts_with("fleet_"), names_to = "fleet", values_to = "fleet_yield") %>%
  group_by(prop_mpa, fleet, placement_strategy) %>%
  summarise(yield = sum(fleet_yield),biodiv = sum(unique(biodiv))) %>%
  ggplot(aes(biodiv, yield, color = placement_strategy)) +
  geom_line() +
  scale_x_continuous(name = "Change in SSB/SSB0",limits = c(0, NA)) +
    scale_y_continuous(name = "Total Yield",limits = c(0, NA)) +
  facet_wrap(~ fleet, scales = "free_y")


frontier_labels <- tmp %>%
  filter(prop_mpa < 0.7) |>
  mutate(delta_50 = (delta_yield - -0.5)^2,
         delta_0 = delta_yield^2) %>%
  group_by(placement_strategy) %>%
  filter(prop_mpa > 0) %>%
  filter(yield == max(yield) | delta_0 == min(delta_0) | delta_50 == min(delta_50)) %>%
  ungroup() %>%
  mutate(pmpa = scales::percent(prop_mpa))

coral_frontier <- tmp %>%
  filter(prop_mpa < 0.7) |>
  ggplot(aes(delta_biodiv, delta_yield)) +
  geom_hline(aes(yintercept = 0), linetype = 2) +
  geom_path(alpha = 0.8, aes(color = placement_strategy),size = 1.1) +
  geom_point(aes(fill = prop_mpa), shape = 21, size = 3) +
  # geom_text_repel(data = frontier_labels, aes(label = pmpa), min.segment.length = 0, box.padding = 0.5) +
  scale_x_continuous(name = "Change in Total SSB/SSB0", labels = scales::label_percent(accuracy = 1), guide = guide_axis(n.dodge = 2), limits = c(0,0.3)) +
  scale_y_continuous(name = "Change in Yield", labels = scales::label_percent(accuracy = 1)) +
  scale_color_manual(values = c("tomato", "steelblue"),name = '') +
  # guides(fill = guide_legend(order = 1, override.aes = list(size = 4))) +
  scale_fill_viridis(name = "% MPA", labels = scales::label_percent(accuracy = 1),guide = guide_colorbar(ticks.colour = "black",frame.colour = "black"), option = "cividis") +
  labs(tag = "C")  +
  theme(
    legend.position = c(0.7, .73),
    legend.background = element_rect(fill = "transparent"),
    legend.text = element_text(size = 10)
  )
coral_fleet_yield <- coral_fleet_results |> 
  mutate(fleet= fct_relabel(fleet, titler)) %>% 
  ggplot(aes(prop_mpa, delta_yield, color = fleet)) +
    geom_hline(aes(yintercept = 0), linetype = 2) +
  geom_line() +
  facet_wrap( ~ placement_strategy, labeller = labeller(placement_strategy = titler)) +
  scale_x_continuous(name = "MPA Size",
                     labels = scales::label_percent(accuracy = 1)) +
  scale_y_continuous(name = "Yield",
                     labels = scales::label_percent(accuracy = 1)) + 
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 8),
        legend.margin = margin())  + 
    scale_color_discrete(name = '') + 
    labs(tag = "A") 




coral_critter_bio <- coral_results %>%
  group_by(prop_mpa, critter, placement_strategy) %>%
  summarise(biodiv = sum(unique(biodiv))) %>%
  group_by(critter, placement_strategy) %>%
  mutate(delta_bio = biodiv) %>%
  mutate(critter = fct_relabel(critter, titler)) %>% 
  # filter(prop_mpa < 0.4) %>% 
  ggplot(aes(prop_mpa, pmin(Inf,delta_bio), color = critter)) +
  geom_line() +
  facet_wrap( ~ placement_strategy, labeller = labeller(placement_strategy = titler)) +
  scale_x_continuous(name = "MPA Size",
                     labels = scales::label_percent(accuracy = 1),
                     guide = guide_axis(n.dodge = 2)) +
  scale_y_continuous(name = "SSB/SSB0") + 
  scale_color_brewer(name = '', type = "qual", palette = "Dark2") + 
    labs(tag = "B") + 
  theme(axis.text = element_text(size = 8), 
        axis.text.x = element_text(size = 8),
        legend.margin = margin())






fig_4 <- ((coral_fleet_yield / coral_critter_bio) | (coral_frontier + plot_layout(guides = "keep")) ) & theme(strip.text = element_text(size  = 6), plot.margin = margin(.05, .05, .05, .05, "cm"), panel.spacing = unit(0.25, "lines"), legend.text = element_text(size = 8))


 ggsave(
   filename = file.path(results_path, "fig_4.pdf"),
   fig_4,
   width = big_width,
   units = "mm"
 )
 
 ggsave(
   filename = file.path(results_path, "fig_4.png"),
   fig_4,
   width = big_width,
   units = "mm"
 )
 
 
fig_4

```

```{r}
#| label: fig-blue-results
#| fig-cap: "Change in yield (A) and SSB/SSB0 (spawning biomass divided by unfished spawning biomass) (B) as a function of MPA size and design strategy by fishing fleet and species. Blue lines indicate impacts of MPAs under status quo habitat, red impacts under climate-driven range shift. Results reflect the outcome of 20 years of simulated MPA protection."
#| fig-height: 6
#| fig-width: 6

linesize <- 0.5

# blue_water_results <- blue_water_results %>% 
#   rename(`purse seine` = purseseine)

blue_fleet_frontier <- blue_water_results %>%
  pivot_longer(names(blue_fleets), names_to = "fleet", values_to = "fleet_yield") %>%
  group_by(prop_mpa, fleet, placement_strategy, climate) %>%
  summarise(yield = sum(fleet_yield),biodiv = sum(unique(biodiv))) %>%
  ggplot(aes(biodiv, yield, color = placement_strategy, linetype = climate)) +
  geom_line() +
  scale_x_continuous(name = "Change in Total SSB/SSB0",limits = c(0, NA)) +
    scale_y_continuous(name = "Total Yield",limits = c(0, NA)) +
  facet_wrap(~ fleet, scales = "free_y")

blue_frontier <- blue_water_results %>% 
  filter(climate == FALSE) %>% 
  group_by(prop_mpa, placement_strategy, climate) %>% 
  summarise(biodiv = sum(biodiv), yield = sum(yield)) %>% 
  group_by(placement_strategy) %>% 
  mutate(delta_yield = yield / yield[prop_mpa == 0] - 1,
         delta_biodiv = biodiv / biodiv[prop_mpa == 0] - 1) %>% 
  ggplot(aes(delta_biodiv, delta_yield, color = placement_strategy, linetype = climate)) + 
  geom_vline(aes(xintercept = 0), linetype = 2) +
  geom_hline(aes(yintercept = 0), linetype = 2) +
  geom_line(size = linesize) + 
  scale_x_continuous(name = "Change in Conservation", labels = scales::label_percent(accuracy = 1)) +
  scale_y_continuous(name = "Change in Yield", labels = scales::label_percent(accuracy = 1)) + 
  theme(legend.position = "top")


blue_fleet_yield <- blue_water_results %>%
  filter(placement_strategy %in% c("target_fishing", "rate")) %>%
  pivot_longer(names(blue_fleets),
               names_to = "fleet",
               values_to = "fleet_yield") %>%
  mutate(fleet = ifelse(fleet == "purseseine","purse-seine",fleet)) |> 
  group_by(prop_mpa, fleet, placement_strategy, climate) %>%
  summarise(yield = sum(fleet_yield), biodiv = sum(unique(biodiv))) %>%
  group_by(fleet, placement_strategy) %>%
  mutate(delta_yield = yield / yield[prop_mpa == 0] - 1) %>%
  ggplot(aes(prop_mpa, delta_yield, color = climate)) +
  geom_hline(aes(yintercept = 0), linetype = 2) +
  geom_line(size = linesize, show.legend = FALSE) +
  facet_grid(fleet ~ placement_strategy, scales = "free_y", labeller = labeller(fleet = titler, placement_strategy = titler)) +
  scale_x_continuous(name = "MPA Size",
                     labels = scales::label_percent(accuracy = 1)) +
  scale_y_continuous(name = "Change in Yield",
                     labels = scales::label_percent(accuracy = 1)) +
  scale_color_manual(
    values = c("#63B8FF", "#FF4500"),
    labels = c("Status Quo", "Range Shift"),
    name = ''
  ) + 
  theme(strip.text = element_text(size = 8)) + 
  labs(tag = "A")



blue_critter_bio <- blue_water_results %>%
    filter(placement_strategy %in% c("target_fishing", "rate")) %>% 
  group_by(prop_mpa, critter, placement_strategy, climate) %>%
  summarise(biodiv = sum(unique(biodiv))) %>%
  group_by(critter, placement_strategy) %>%
  mutate(delta_bio = biodiv) %>%
  ungroup() %>% 
  mutate(placement_strategy = fct_relabel(as.factor(placement_strategy), titler)) %>% 
  ggplot(aes(prop_mpa, pmin(1,delta_bio), color = climate, linetype = placement_strategy)) +
  geom_line(size = linesize) + 
  facet_wrap(~critter) +
  scale_color_manual(values = c("#63B8FF","#FF4500"), labels = c("Status Quo","Range Shift"), name = '') +
  scale_x_continuous(name = "MPA Size",
                     labels = scales::label_percent(accuracy = 1), guide = guide_axis(n.dodge = 2)) +
  scale_y_continuous(name = "SSB/SSB0",limits = c(0, NA)) +
  scale_linetype(name = "MPA Strategy") + 
    labs(tag = "B")+
  theme(strip.text = element_text(size = 7),
        axis.text.y = element_text(size = 6))


# blue_frontier / blue_critter_bio + plot_layout(heights = c(1,2))
fig_5 <- (blue_fleet_yield / 
blue_critter_bio) + plot_layout(heights = c(1.5,2))

 ggsave(
   filename = file.path(results_path, "fig_5.pdf"),
   fig_5,
   width = big_width,
   units = "mm"
 )

  ggsave(
   filename = file.path(results_path, "fig_5.png"),
   fig_5,
   width = big_width,
   units = "mm"
 )
 
fig_5
```