On the spillover effect and optimal size of marine reserves

Marine reserves are an essential component of model fisheries management. As implementing marine reserves induces an inherent tradeoff between the harvesting and conservation, to solidify the insight into fisheries management with marine reserves is fundamental for management success. Finding an optimal reserve size that improves the fishing yield is not only theoretical interest but also practically important to assess the underlying tradeoffs and to facilitate decision making. Also, since the species migration determines the degree of the spillover effect from a marine reserve, it is a key consideration to explore the performance of marine reserve. Here, we investigate an optimal reserve fraction and its management outcome under various spillover strength via a simple two-patch mathematical model, in which one patch is open to fishing, and the other is protected from fishing activities. The two-patch model is approximated by a single population dynamics when the migration rate is sufficiently larger than the growth rate of a target species. In this limit, it is shown that an optimal reserve size exists when the pre-reserve fishing is operated at the fishing mortality larger than the fMSY, the fishing mortality at the maximum sustainable yield (MSY). Also, the fishing yield with the optimal reserve size becomes as large as MSY in the limit. Numerical simulations across various migration rates between two patches suggest that the maximum harvest under the management with a marine reserve is achieved in this limit, and this contrasts with the conservation benefit in which is maximized at the intermediate migration rate.


INTRODUCTION
Marine reserves, or no-take marine protected areas (MPAs), are a central tool in modern fishery management to reduce fishing pressure and to preserve biodiversity (White & Costello, 2014;Hastings, Gaines & Costello, 2017;Sala & Giakoumi, 2018). Implementation of marine reserves leads to fishing closures, and it can cause a management tradeoff between harvest and conservation (Klein et al., 2013;Bozec et al., 2016). The spillover effect from marine reserves is essential to determine fishery and conservation benefits of reserves (Gruss et al., 2011); hence, understanding the migration effect is fundamental to management success.
In modern fishery management, in which managers seek optimal yields (Clark, 1990;Kar & Ghosh, 2013), optimal reserve sizes are widely discussed (Neubert, 2003;Hart, 2006;Gaines et al., 2010). Optimal yield is also a theoretically important concept, as it represents a baseline to assess the intensity of a management tradeoff. In models, marine reserves can reduce fishing yields when a fishery manages a stock sustainably (Holland & Brazee, 1996;Hart, 2006;Ralston & O'Farrell, 2008;White et al., 2010); hence, seeking an optimal reserve size is a key consideration. Throughout this article, I refer to the optimal reserve size in the sense of sustainable fishing yields. However, it should be noted that marine reserves can provide most fishery and conservation benefits at anywhere from short to long distances (Manel et al., 2019), including promotion of biodiversity (Almany et al., 2009;Russ & Alcala, 2011) and genetic diversity (Baskett & Barnett, 2015), improving ecological resilience (Takashina & Mougi, 2014;Barnett & Baskett, 2015) and optimal fishery profit (Sanchirico et al., 2006), mitigation of impacts from climatic change (Green et al., 2014;Roberts et al., 2017) and catastrophic events (White, Baskett & Hastings, 2019), to name just a few.
Previous studies have revealed several important properties of the optimal reserve size. Using a model of population dynamics described by a stock-recruit relationship, Hart (2006) determined that an optimal reserve size exists when the slope of the stock-recruit curve at a given fishing mortality rate is larger than the inverse of the harvest-free spawning stock biomass per recruit. Applying this principle to canary rockfish and Georges Bank sea scallop, Hart (2006) showed that a fishing mortality rate greater than f MSY is necessary to meet this condition. Takashina (2016) adopted an age-structured metapopulation model and showed that an intermediate recruitment success of eggs and a large fishing mortality rate are required for marine reserves to improve fishing yields. Bensenane, Moussaoui & Auger (2013) analyzed a bioeconomic model that describes dynamics of population size and fishing mortality rate. They derived an algebraic formulation of an optimal reserve size, composed of the market price of the target species, catchability coefficient, the carrying capacity, and the cost of fishing effort, showing that an optimal reserve size is economically feasible. Furthermore, several studies showed the existence of an optimal reserve size in models where larval dispersal attributes the spillover from marine reserves (Gaylord et al., 2005;White & Kendall, 2007;Ralston & O'Farrell, 2008;De Leo & Micheli, 2015). However, previous studies concerning optimal reserve size have avoided the effect of adult movement (Hart, 2006;De Leo & Micheli, 2015;Takashina, 2016) or have addressed only well-mixed populations (White & Kendall, 2007;Bensenane, Moussaoui & Auger, 2013).
Here, I investigate the relationship between optimal reserve size and adult spillover, and realized conservation benefit, measured by the total population size. The model is a simple two-patch model, where one patch represents a fishing ground and the other, a marine reserve. Migration connects the two patches at a species-specific rate, and affects the degree of spillover from the marine reserve. This creates source-sink dynamics between the fishing ground and marine reserve, and this is a common structure of existing spatially-explicit models (e.g., Crowder et al., 2000;Neubert, 2003;Takashina & Mougi, 2014;Ghosh et al., 2017). The two-patch model also allows us to examine different migration modes (Amarasekare, 2004), such as positive/negative density-dependent, as well as density-independent migration. This simple approach yields analytical results by introducing an aggregated model when the migration rate is sufficiently larger than the growth rate of a target species (Iwasa & Andreasen, 1987;Auger, de La Parra & Poggiale, 2008;Bensenane, Moussaoui & Auger, 2013;Takashina & Mougi, 2015;Takashina, Lee & Possingham, 2017). The aggregated model allows us to derive conditions for optimal reserve size, fishing yield, and total population size under management with the optimal reserve fraction. I show that an optimal reserve size exists when the fishing mortality is larger than f MSY , and that the fishing yield and the total population size under the management correspond to MSY and X MSY , the population size at MSY.
In addition, I numerically examine the model with various migration rates and modes when model aggregation is not valid. In addition to fishing yields, I also calculate the total population size as a measure of the conservation benefit of marine reserves and discuss the tradeoffs of using marine reserves as a management strategy. I find that an optimal reserve size exists under various migration rates and modes. Optimal reserve size tends to increase with migration rate and approaches the optimal size of the aggregated model, where the optimal size is the largest.

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: 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) 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 1 ,x 2 ) where x 1 and x 2 are the population of the fishing ground and marine reserve, respectively, Eq. (1) becomes 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 1 ,x 2 ) as 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:

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 1 +x 2 operated at the time scale of fast parameter τ = mt . Thus, Eqs. (3a) and (3b) can be approximated by a single aggregated model ( 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: and fishing yield is where, the superscript AG indicates the equilibrium of the aggregated model. By solving dȲ AG / df = 0 about f and α, respectively, I obtain the optimal fishing mortality and reserve size: 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: 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 bȳ I often use notationȲ AG * when the substitution is obvious. Similarly, I regardX AG * 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Ȳ AG * /MSY =X AG * /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 /MSY ≤Ȳ AG * /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: densitydependent (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 .

DISCUSSION
The effect of a marine reserve on fish catch is a crucial consideration because creation of a marine reserve in an existing fishing ground initially reduces fishing opportunities. By taking advantage of a simple two-patch model, I derived a necessary condition for an optimal size of the marine reserve to exist (f > f MSY ; Eq. (11)) confirming a previously reported value (Pezzey, Roberts & Urdal, 2000;Hart, 2006;Ralston & O'Farrell, 2008). I also derived a simple formula for optimal fishing mortality rate (Eq. (9)) and optimal reserve size (Eq. (10)). Numerical simulations examined various migration rates and modes and revealed several new findings. That is, the theoretically predicted fishing yield from the aggregated model gives an upper bound: a well-mixed population with a large migration rate maximizes fishing yield. In addition, even with the condition f > f MSY , a moderate migration rate is necessary, on top of previous findings (Pezzey, Roberts & Urdal, 2000;Hart, 2006;Ralston & O'Farrell, 2008), for an optimal reserve size to exist when fishing mortality is moderately high (Fig. 3B). However, this additional condition is not necessary when fishing mortality is high (Fig. 3C). I also found that different migration modes lead to qualitatively similar results, a high migration rate decreases the effect of different migration modes, and fishing yields and population sizes approach values predicted by the aggregated model. On the other hand, the model shows that a marine reserve provides a larger total population, a measure of conservation benefit when the migration ability of a target species is low or moderate, contrasting to benefits for fishing yield. This emphasizes the importance of spillover effects when considering the tradeoffs of fisheries management using marine reserves. The necessary condition derived from the aggregated model for marine reserves to improve fishing yields confirms a condition previously reported under different modeling frameworks (Pezzey, Roberts & Urdal, 2000;Hart, 2006;Ralston & O'Farrell, 2008), suggesting potentially generic characteristics of marine reserve management. If a marine reserve replaces a certain fraction of a fishing ground, the fishing mortality rate at the MSY level should increase by a factor inversely proportional to the fraction of the fishing ground, 1 − α. In practice, this corresponds to a situation where all fishermen remain in the contracted fishing ground, with proportion 1 − α, after reserve creation. Under these conditions, the fishing yield becomes as large as MSY (Hastings & Botsford, 1999;Bensenane, Moussaoui & Auger, 2013). These results suggest that the optimal reserve size spans α * ∈ (0,1), and the optimal reserve size approaches 1 as fishing mortality becomes sufficiently large. In practice, however, excluding most fishing from a region of concern may be infeasible with multiple management tradeoffs. In fact, the often suggested range of optimal reserve size to increase fishing yields and/or to meet management and conservation objectives is smaller; e.g., about 30% of the region of concern in the synthetic analysis (Gaines et al., 2010) or 20%−50% in Halpern & Warner (2003). My results suggest that the maximum harvest under management with a marine reserve is achieved when the species migration rate becomes sufficiently large (m 1). This indicates that the reduced fishing ground can be compensated by increased fishing mortality when there is sufficient migration from a marine reserve to a fishing ground. However, intensified fishing mortality to achieve the maximum fishing yield given a migration rate m often leads to a smaller total population than the MSY: X res /X MSY < 1, except for the case m 1 where X res /X MSY = 1.
On the other hand, the model shows a larger population size at an intermediate migration rate (Figs. 4B and 4C), representing the underlying tradeoff between fishing yield and population size (i.e., conservation benefit). Previous studies also showed that marine reserves may provide larger conservation benefits, such as a larger population recovery and reproductive capacity, at relatively low to moderate migration rates (Blyth-Skyrme et al., 2006;Gruss et al., 2011;Takashina & Mougi, 2014).
The finding that high migration rates (m 1) produce the highest fishing yields implies a close relevance to the results of Neubert (Neubert, 2003), where the number of marine reserves to maximize the fishing yield can become infinitely many in a finite length of one-dimensional space, leading to an infinitely large exchange rate between marine reserves and fishing grounds. It suggests that in practice, migration between marine reserves and non-reserve sites can be controlled by the configuration of marine reserves in space, and that the migration rate is not a purely biological parameter. Therefore, there may exist a marine reserve design to realize a high migration rate in the setting of my model, even for a sedentary species by, for example, creating smaller marine reserves than the home range sizes of target species. Although the two-patch model simplifies the marine environment, and the migration rate m is constant regardless of the reserve size α, investigation of the influence of migration ability (e.g., diffusivity) and mode of a target on the migration rate, m, will further improve marine reserve management. The two-patch model, however, is still a useful framework to discuss the performance of marine reserves with a simple function of migration terms, enabling us to avoid the description of complex migratory dynamics within each patch. Rosenberg et al. (1997) experimentally demonstrated that a brittle star generally shows density-dependent dispersal and estimated dispersal speed. Therefore, it may be possible to perform an experiment to estimate the parameter controlling density dependence. Also, the fact that the configuration of marine reserves influences the migration rate highlights the importance of considering the spatial scaling of management. Previously, Takashina & Baskett (2016) demonstrated that the management unit scale is decomposed from the spatial scale in which biological processes operate. They showed that fisheries management with a finer management unit scale results in larger fishing yields than management with a larger management scale, since the former can realize a fine-tuned allocation of fishing efforts across management units. Fine-tuned fishing activities across fishing grounds allows higher flexibility in management decision making than the two-patch model, where I can assign only a single fishing mortality rate. Hence, consideration of the migration rate as a function of the management unit scale may offer a more flexible way to mitigate underlying management tradeoffs.
Prediction of the equivalence in MSY and the yield from management involving marine reserves has been revised by preceding studies showing that marine reserves result in a fishing yield larger than MSY. For example, Gaylord et al. (2005) and Ralston & O'Farrell (2008) attributed excess yields to spatial structures of the model and post-dispersal density-dependent recruitment of larvae, and De Leo & Micheli (2015) demonstrated that a highly complex spatially-explicit, stage-structured stochastic model improves fishing yields. White & Kendall (2007) claimed that model complexity is not responsible for this result, and they demonstrated that only the effect of post-dispersal, density-dependent recruitment in population dynamics induces excess fishing yields. On the other hand, my simple mathematical model does not show a fishing yield larger than MSY. The lack of detailed spatial structures, stage/age-structure, and/or post-dispersal density-dependent recruitment, as in the above examples, may explain this contrasting result. Particularly, species migration and larval dispersal cause a rather different population mixing pattern. For example, post-dispersal, density-dependent recruitment causes a spontaneous reduction of the number of larvae recruited, and it is often assumed to occur at a discrete-time (i.e., seasonal breeding cycle) (Gaylord et al., 2005;White & Kendall, 2007;Ralston & O'Farrell, 2008;De Leo & Micheli, 2015), contrasting with my continuous-time model where population exchanges and density-dependent population control occur throughout the year. In marine environments connected by larval dispersal, less mobile species are more likely to increase fishing yields than highly mobile species. That is, sessile species in reserves tend to stay in the reserve longer; hence, they receive greater benefits from reserves and provide more recruits to fishing grounds (Holland & Brazee, 1996;White et al., 2010). On the other hand, high adult migration rate causes frequent population exchanges between reserves and fishing grounds, and such species are more vulnerable to fishing activities than those with lower migration rates, leading to higher fishing yields than for less mobile taxa. Provided that the optimal yield in the aggregated model is the upper bound in the two-patch model, the model serves as a basic model to assess the effect of stage/age-structure, density-dependent larval recruitment, and more complex fishing regulations under different levels of spillover. For example, size/age-based regulation is prevalent in fisheries management, and its explicit consideration with various species migration scenarios further improves our understanding of optimal reserve size.
In this article, I used the simplest possible model to discuss optimal reserve size while accommodating various adult migration rates and modes. Optimal yields provide a useful benchmark to assess the strength of management tradeoffs. There are, of course, a number of different mathematical approaches possible, particularly for models relating to marine reserves (Fulton et al., 2015). With a good understanding of adult movement of a species of concern, one natural extension of the two-patch model is the integration of metapopulation structure that realizes more complex pathways of species and enables more realistic conservation planning as discussed in Kininmonth et al. (2010) and Kininmonth et al. (2011). Multiple-species dynamics are also important components of marine ecosystem management (Pikitch et al., 2004;Takashina, Mougi & Iwasa, 2012) and discussion of the impact of species interactions on optimal reserve sizes will improve our insight. Understanding these effects on optimal reserve management, together with adult movement, will be necessary.