On the spillover effect and optimal size of marine reserves for sustainable fishing yields

Basic model

The starting point is the commonly used Schaefer model (Gordon, 1954; Schaefer, 1954) in which population dynamics of the target species x are described with growth rate r, carrying capacity K, and fishing mortality rate f as follows: (1)${\displaystyle \frac{dx}{dt}=rx\left(1-\frac{x}{K}\right)-fx.}$ In this model, MSY, the fishing mortality rate at MSY, f_MSY, and the population size at MSY, X_MSY, are described as Clark (1990) (2)${\displaystyle \mathrm{MSY}=\frac{rK}{4},f_{MSY}=\frac{r}{2},\mathrm{and}{X}_{MSY}=\frac{K}{2}.}$ To investigate the effect of a marine reserve, I employ a common two-patch model (e.g.,  Takashina, Mougi & Iwasa, 2012; Takashina & Mougi, 2014; Ghosh et al., 2017) where one patch is open to fishing (i = 1; with fraction 1 − α), the other patch (i = 2; with fraction α) is protected from fishing (i.e.,  f = 0), and migration connects the two patches (Fig. 1). With the migration function M(x₁, x₂) where x₁ and x₂ are the population of the fishing ground and marine reserve, respectively, Eq. (1) becomes (3a)${\displaystyle \frac{d{x}_1}{dt}=r{x}_1\left(1-\frac{x_1}{\left(1-\alpha \right)K}\right)-f{x}_1+M\left(x_1,x_2\right),}$ (3b)${\displaystyle \frac{d{x}_2}{dt}=r{x}_2\left(1-\frac{x_2}{\alpha K}\right)-M\left(x_1,x_2\right),}$ I consider the situation where migrations between the patches are either random or positively/negatively density-dependent. With the migration rate between the marine reserve and the fishing ground, m, and the parameter controlling the migration mode, s, I describe the migration function M(x₁, x₂) as (4)${\displaystyle M\left(x_1,x_2\right)=m\left\{{\left(\frac{x_2}{\alpha K}\right)}^s\left(1-\alpha \right)x_2-{\left(\frac{x_1}{\left(1-\alpha \right)K}\right)}^s\alpha x_1\right\},}$ where, the function represents the negative density-dependent migration when −1 < s < 0, random migration when s = 0, and density-dependent migration when s > 0 (Amarasekare, 2004).

I use the following notations to describe the fishing yield and the total population size under management with a marine reserve: (5)${\displaystyle Y_{res}=f{x}_1,\mathrm{and}{X}_{res}=x_1+x_2.}$

Model aggregation

When the migration rate is sufficiently larger than the growth rate (m ≫ r), there are fast and slow dynamics operating at different time scales (Iwasa & Andreasen, 1987; Auger, de La Parra & Poggiale, 2008). Then, the migration term has a negligible effect on the total population X = x₁ + x₂ operated at the time scale of fast parameter τ = mt. Thus, Eqs. (3a) and (3b) can be approximated by a single aggregated model (Takashina & Mougi, 2015) (6)${\displaystyle \frac{dX}{d\tau }=rX\left(1-\frac{X}{K}\right)-f\left(1-\alpha \right)X.}$ Hereafter, I discuss analytical aspects of the aggregated model Eq. (6), and perform numerical calculations across the migration rate m, including a situation where the model aggregation is not valid.

Analysis of the aggregated model Eq. (6)

The aggregated model allows us to derive the explicit form of the equilibrium as follows: (7)${\displaystyle {\bar{X}}^{AG}=K\left(1-\frac{f\left(1-\alpha \right)}{r}\right),}$ and fishing yield is (8)${\displaystyle {\bar{Y}}^{AG}=f\left(1-\alpha \right)K\left(1-\frac{f\left(1-\alpha \right)}{r}\right),}$ where, the superscript AG indicates the equilibrium of the aggregated model. By solving $d{\bar{Y}}^{AG}∕$ df = 0 about f and α, respectively, I obtain the optimal fishing mortality and reserve size: (9)${\displaystyle f^{AG\ast }=\frac{r}{2\left(1-\alpha \right)}=\frac{f_{MSY}}{1-\alpha },}$ (10)${\displaystyle {\alpha }^{AG\ast }=1-\frac{r}{2f}.}$ Equation (9) suggests that one needs to increase the fishing mortality rate at the MSY level with a rate inversely proportional to the fraction of fishing ground 1 − α after establishment of a marine reserve. It becomes infinitely large when the fraction of the marine reserve approaches unity (Fig. 2). On the other hand, Eq. (10) shows that for an optimal reserve size to exist, the fishing mortality rate should be larger than MSY: (11)${\displaystyle f>\frac{r}{2}=f_{MSY}.}$ In other words, if fishing mortality is smaller than MSY, there is no optimal reserve size that can improve fishing yield. Also, the optimal reserve size approaches 1 (i.e., complete fishing ban) as fishing mortality becomes large (Fig. 2).

From Eq. (6), maximum fishing yields coincide with MSY of Eq. (1). In fact, substituting either Eq. (9) or Eq. (10) into Eq. (8) recovers Eq. (2). I denote this by (12)${\displaystyle {\bar{Y}}^{AG\ast }{|}_{f^{\ast }}={\bar{Y}}^{AG\ast }{|}_{{\alpha }^{\ast }}=\mathrm{MSY}.}$ I often use notation ${\bar{Y}}^{AG\ast }$ when the substitution is obvious. Similarly, I regard ${\bar{X}}^{AG\ast }$ as the population size when fishing yield is given by Eq. (12).

Numerical investigation for general situation

When the migration rate is not large, the aggregated model is not valid. Nonetheless, I will show that the analytical results above provide a maximum fishing yield, and this becomes a benchmark to discuss the performance of an introduced marine reserve.

Here, I numerically solve Eqs. (3a) and (3b) to find the fishing yield and the optimal reserve size α^∗, as well as the total population size under various reserve sizes and migration rates. To compare these quantities with those of management without a marine reserve, I introduce the following two normalized quantities: the fishing yield normalized by MSY, Y_res/MSY, and the total population size normalized by X_MSY, X_res∕X_MSY. Note from the analysis above, the aggregated model under the optimal reserve size, α^AG∗, gives the values of normalized fishing yield and population size ${\bar{Y}}^{AG\ast }∕\mathrm{MSY}={\bar{X}}^{AG\ast }∕X_{MSY}=1$.

Figures 3A–3C show the normalized fishing yield under density-independent migration (s = 0). As expected, the optimal reserve size α^∗, if any, approaches that of the aggregated model α^AG∗ as the migration rate m becomes large. Also, the normalized fishing yield of the aggregated model gives the upper bound: ($Y_{res}∕\mathrm{MSY}\le {\bar{Y}}^{AG\ast }∕\mathrm{MSY}$). These calculations also suggest that the condition for the positive optimal reserve size to exist (Eq. (11)) still holds for various migration rates. Therefore, the improvement of fishing yield by introducing a marine reserve occurs only when the initial fishing mortality exceeds MSY, f_MSY. Although this condition is necessary for a marine reserve to increase the harvest when the migration rate is not high enough (i.e., the model aggregation is not valid), it does not guarantee the existence of an optimal reserve size. For example, when fishing mortality is moderately high (f = 0.75; Fig. 3B), a spillover effect from the marine reserve is necessary for an optimal reserve size to exist. A high fishing mortality rate (f = 1) relaxes this requirement (Fig. 3C) and an optimal reserve size exists for all migration rates above m > 0.01. Figures 3D–Figures 3F show a slice of the heat map shown in the top panels at m = 0.1, 1, or 10. The fishing mortality rate f corresponds to the top panel. These demonstrate that the fishing yield tends to be higher when the migration rate m is high. Also, there are peaks when fishing mortality is larger than f_MSY = 0.5, and these correspond to the optimal reserve size α^∗.

On the other hand, the bottom three panels of Fig. 4 show the normalized total population size. This value represents the conservation benefit of the marine reserve. Management with an optimal marine reserve size tends to give a smaller population size than X_MSY (i.e., X_res∕X_MSY < 1). However, if fishing mortality is larger than f_MSY, this value approaches 1, and the prediction of the aggregated model, as the migration rate m becomes sufficiently large (Figs. 4B, 4C). Also, increasing the reserve size provides a higher normalized population size, and the population size becomes larger at low to moderate migration rates (m is about 1 in Fig. 4) at a given reserve size. The bottom three panels in Fig. 4 show more explicitly this relationship when the migration rate is m = 0.1, 1, or 10. Unlike fishing yields, the population size, a measure of conservation benefit, is larger when the migration rate is moderate (m = 1) rather than high (m = 10). This suggests a mismatch between optimal harvesting and the conservation benefits.

I found qualitatively similar trends in two alternative migration modes: density-dependent (s = 1; Fig. 5) and negatively density-dependent (s =  − 0.5; Fig. 6) migrations between patches where I show only the heat maps. Some quantitative differences appear when the migration rate m is around 1, but these approach the prediction of the aggregated model as the migration rate becomes sufficiently large, and optimal reserve size appear when fishing mortality is larger than f_MSY.