Marine reserves and optimal dynamic harvesting when fishing damages habitat

Marine fisheries are a significant source of protein for many human populations. In some locations, however, destructive fishing practices have negatively impacted the quality of fish habitat and reduced the habitat’s ability to sustain fish stocks. Improving the management of stocks that can be potentially damaged by harvesting requires improved understanding of the spatiotemporal dynamics of the stocks, their habitats, and the behavior of the harvesters. We develop a mathematical model for both a fish stock as well as its habitat quality. Both are modeled using nonlinear, parabolic partial differential equations, and density dependence in the growth rate of the fish stock depends upon habitat quality. The objective is to find the dynamic distribution of harvest effort that maximizes the discounted net present value of the coupled fishery-habitat system. The value derives both from extraction (and sale) of the stock and the provisioning of ecosystem services by the habitat. Optimal harvesting strategies are found numerically. The results suggest that no-take marine reserves can be an important part of the optimal strategy and that their spatiotemporal configuration depends both on the vulnerability of habitat to fishing damage and on the timescale of habitat recovery when fishing ceases.


Introduction
It is not uncommon for modern fishing gear to impact the habitats of the fish being harvested. Trawls and dredgesfishing gear that is dragged along the ocean's floor-can be particularly harmful to the benthic habitat upon which many commercially valuable species rely. Indeed the effects of these gears in some locations have been compared to forest clearcutting (Watling and Norse 1998). These effects are well known; reviews of fishing gear impacts can be found in papers by Dayton et al. (1995) and Chuenpagdee et al. (2003) and Grabowski et al. (2014). Recent analyses suggest that the extent of trawling, particularly in the deep sea, has been underestimated (Victorero et al. 2018).
Given the collateral damage that fishing may impose on essential habitat, it is reasonable to ask whether marine reserves-areas that are closed to fishing-might be a useful part of management. After all, reserves would protect both stocks and their habitats. Of course, closures might also negatively affect yields or rents, by closing off a portion of the stock to harvest as well as by changing the distribution of harvesters (Kaiser et al. 2002;Kellner et al. 2007). The analysis of bioeconomic models is useful for understanding this tradeoff.
Spatially explicit models are necessary to fully address the question of marine reserve utility, since reserves are a form of spatial management (Herrera and Lenhart 2010). Previous analyses of spatial models (that ignore the potential for habitat damage) have found that reserves may be necessary as part of a strategy designed to maximize yield or rent depending on the ecological and economic circumstances (e.g., Neubert 2003;Joshi et al. 2009;Ding and Lenhart 2009;Kelly et al. 2015). Building on these results, Neubert (2013, 2015) have analyzed spatially explicit models with the object of evaluating the role of marine reserves in optimal (i.e., rent maximizing) harvesting when fishing has negative habitat effects. In Moeller and Neubert (2015), they analyzed a two-patch system and showed that reserves were necessary to maximize sustainable profit when dispersal between the patches was sufficiently high and habitat was especially vulnerable to damage. They also showed that spatially heterogeneity was not necessary for the result to hold. In Moeller and Neubert (2013), fishing had one of two negative effects: it either reduced the intrinsic growth rate of the fish or increased the effects of density dependence. They found that in the latter case, when fishing intensified the negative effects of density dependence, marine reserves were optimal under similar circumstances (that is, when habitat was vulnerable to fishing damage).
Neither of these models explicitly track the amount or quality of habitat. Rather, they assume that changes in fishing effort directly and instantaneously change the vital rates of the fish. However, in reality, habitat may recover slowly when harvesting is stopped. In addition, neither model accounted for the other benefits that habitat might provide (e.g., existence values or ecosystem services). In this article, we will construct a system of two reaction-diffusion partial differential equations (PDEs) that incorporates the deleterious effects of harvesting on habitat, and that explicitly accounts for habitat dynamics. It also will incorporate an existence value for the habitat to account, for example, for non-extractive ecosystem services that intact habitat may provide. Using that model, we derive a harvesting strategy that maximizes the net present value of the stock and explore the circumstances under which marine reserves are part of that optimal strategy.
In "Model and objective functional," we formulate the coupled stock-habitat model and describe the optimization problem. In "Derivation of the optimality system," we derive the optimality system-partial differential equations and associated boundary conditions that must be satisfied by the optimal state and adjoint variables with the corresponding optimal control characterization. We solve the optimality system numerically in "Numerical simulations" for several scenarios and show how habitat dynamics and habitat vulnerability (i.e., the degree to which fishing degrades habitat) affect the optimal harvest. We finish with a brief discussion. Following Joshi et al. (2009) and Ding and Lenhart (2009) and Kelly et al. (2015), we consider a stock living in an onedimensional spatial domain over the time interval [0, T ].

Model and objective functional
At the initial time (t = 0), we assume that the stock density, u(x, t), is known. i.e., on × (t = 0). For t > 0, we assume that changes to the stock density derive from four processes: population growth due to reproduction and natural mortality (as described by the function f ), spatial movement in the form of diffusion (with strictly positive diffusion coefficient a 1 (x, t)) and/or advection (with advection speed b 1 (x, t)), and harvesting at a rate that is proportional to the effort h(x, t). Combining these processes gives the parabolic partial differential equation The "catchability coefficient" c measures the effect of a unit of effort on fishing mortality; without loss of generality, we take c = 1 for the remainder of the paper. We imagine that the stock cannot survive outside of the habitat and so adopt the Dirichlet boundary condition An important component of model (2) is the nonlinear growth term f (u, k). This term depends both on the stock density, u(x, t), and on k(x, t), which we take as a measure of habitat quality. The function k(x, t) increases with improved habitat quality. Therefore, larger values of k(x, t) reduce the effects of density-dependent mortality and increase the maximum density that the stock can attain at location x. These effects can be captured by assuming that f (u, k) has the form where r 1 is the intrinsic growth rate of the stock (Fogarty 2005;Foley et al. 2012). We will require k(x, t) to be nonnegative. In the absence of harvesting or dispersal, the stock would grow to density M when k = 0, so we think of M as the carrying capacity for the stock when its habitat is most degraded. Next, we must account for the dynamics of habitat quality k(x, t), which we would like to have certain properties. We expect that (1) habitat quality is nonnegative; (2) habitat quality will, in the absence of harvesting, increase up to a limit; (3) if habitat quality is low in some location, it can be "filled in" by diffusive spread from adjacent locations; and (4) habitat quality is reduced at a rate that is proportional to the fishing effort employed at that location. A model that has these properties is similar to the PDE (2): which holds on on × (0, T ), with boundary and initial conditions k(x, t) = 0 on ∂ × (0, T ), and (6) Here, a 2 (x, t) (strictly positive) and b 2 (x, t) are the diffusion and advection coefficients for habitat quality. We assume nonlinear growth of habitat quality with g(k) given by The negative effect of harvesting on habitat quality is given by the term σ kh(x, t); the parameter σ measures the vulnerability of habitat to harvesting. We expect the diffusion coefficient to be smaller for habitat quality k than for the stock u. For the reader's convenience, we have summarized the parameters, and their units, in Table 1. We now take the view of a social planner whose objective is to maximize the present value of the stream of future rent. Rent comes from selling the harvest (at a fixed unit price P N ) less the costs of effort, and from the existence value of quality habitat, which we take as constant P K . The cost of effort, in turn, has two components. The first can be thought of as the wage that must be paid to harvesters W 0 . The second is a congestion cost W 1 h that accounts for fishermen getting in each other's way when effort is concentrated in space. We assume that W 0 ≥ 0 and W 1 > 0.
If rent is discounted at the rate δ, then the net present value is The social planner seeks to maximize J (h) by choosing the harvest effort in space and time from the control set The optimal harvest h * (x, t) thus satisfies

Derivation of the optimality system
Our state system consists of the PDEs for u and k and their corresponding initial and boundary conditions. Given a control h, there exists a unique weak solution to the state system in L ∞ (Q) ∩ L 2 (0, T : H 1 0 ( )) with u t , k t in L 2 (0, T : H −1 ( )) using the standard techniques (Kelly et al. 2015;De Silva et al. 2017). The solutions of the state system satisfy a priori estimates in our weak solution space, which leads to the existence of an optimal control (using compactness results from Simon (1987)).
To characterize an optimal control, we need to differentiate the map h → J (h) with respect to the control. Since the state solutions depend on the control h and are explicitly in the integrand of J , we also must differentiate the map h → (u, k), and the derivatives of this control-to-states map are called the sensitivity functions. The limits in these derivative calculations can be justified by the appropriate a priori estimates (Evans 2010). We show the framework for finding the sensitivity functions and then the corresponding adjoint functions needed to characterize an optimal control (Lenhart and Workman 2007;De Silva et al. 2017).

Sensitivity functions
To derive the optimality system, we first differentiate the map h → (u, k), and those derivatives are the sensitivity functions, which satisfy a linearized version of the state system. Denote an optimal control by h * , the sensitivity functions are defined by the limits, and where h = h * + l is another control in H . Note that we are using directional derivatives in the direction l. The state systems corresponding the controls h and h * are given by We form difference quotients, dividing by , to obtain the following PDE: Extra consideration is required for the growth terms for the population density of the stock and for the habitat quality. For the growth term of the stock population density, by adding and subtracting appropriate terms, we have Similarly for the habitat quality, the following PDE results for the difference quotients: By taking the limits as → 0, we have the PDE system with corresponding boundary and initial conditions for the sensitivity functions: with and and where f u , f k , and g k are partial derivatives with respect to u and k.

Adjoint Functions
We use the operators in the sensitivity PDEs (21) and (24) to find the adjoint functions, which we use to characterize the optimal control. We rewrite the sensitivity PDEs as where and We rewrite the sensitivity system in Eqs. 27 and 28 as where and The adjoint operators L * i are related to sensitivity operators L i by L 2 inner products in our weak solution space, formally written as we use L 1 and L 2 to get expressions for L * 1 p and L * 2 q, the adjoint operators, where p and q are the adjoint variables. By integration by parts on the sensitivity operators, we obtain the adjoint operators, where Note that the δ terms will be used in the differentiation of J with respect to h and take into account the discount factor in J . Lastly, we have the transversality conditions, p(x, T ) = 0 and q(x, T ) = 0, and the boundary conditions, p(x, t) = 0 and q(x, t) = 0 on ∂ × (0, T ).

Optimal control characterization
We now use the sensitivity and adjoint PDEs to find the optimal control characterization. We differentiate our objective functional J at h * in the direction l. Since J is maximized at h * , we obtain = lim Then, using the weak form of the sensitivity and adjoint PDEs, we obtain Using this with appropriate variations l and the required bounds on h * , we obtain our desired optimal control characterization, Since the states u, k and the adjoints are 0 at the boundary of the spatial domain, we expect the fishing effort to be small (or vanish) near that boundary. Note that increasing σ , the vulnerability of the habitat to harvesting, would likely decrease the optimal fishing effort h * .

Numerical simulations
We now turn to numerical simulations to illustrate optimal solutions to the problem. We solved the optimality system (state and adjoint systems coupled with optimal control characterization) with a forward-backward sweep method (Hackbusch 1978;Lenhart and Workman 2007), using explicit finite differences to solve the PDEs. In addition to the optimal spatiotemporal distribution of effort, habitat quality, and stock density, we also computed the total effort the average stock density and the average habitat quality For our numerical solutions, we assume (for simplicity) that all of the parameters are constant in space and time. We also assume that habitat quality grows more slowly than the stock (r 2 < r 1 ) and that habitat quality diffuses more slowly than stock (a 2 < a 1 ). The parameter values we used in our simulations are listed in Table 1.
We considered two initial conditions. For most simulations, we assumed pristine initial conditions; i.e., we choose u 0 (x) and k 0 (x) as the equilibrium of an unexploited stock (Fig. 1, left panel). Because many stocks are far from pristine, we also computed harvesting strategies starting with an overexploited, open-access stock (Fig. 1, right panel). Under open access, we assume that effort increases until, at equilibrium, rent is dissipated at every location in space; i.e., Thus, the open access effort levelh(x) is given bȳ We substituteh for h in the partial differential equations (2) and (5), and, with corresponding boundary conditions, find the equilibrium solution.
In the remainder of this section, we explore the consequences of habitat quality dynamics, habitat vulnerability, habitat recovery rates, habitat existence values, and initial conditions on the optimal distribution of effort, and on the resulting distribution of stock and habitat quality. Given this, our model is relatively complex. Thus, we explore these consequences one at a time with a view towards improved understanding.
Habitat quality dynamics To assess the role that the dynamics of habitat degradation and recovery play in structuring the optimal harvesting strategy, we first computed the solution in a baseline case where there are no impacts of harvesting on habitat quality (σ = 0), there is no existence value associated with habitat quality (P K = 0), and in which habitat quality is constant in space and time (k = 1) (Fig. 2, top row). For the parameters, we chose (summarized in Table 1), the optimal distribution of harvesting begins with an initial period of intense harvest, which has the benefit of reducing negative effects of density dependence. Subsequently, a large central closed area gradually grows in  Table 1 Fig. 2 Stock dynamics and optimal harvesting strategies for three habitat dynamic scenarios: a constant, invulnerable habitat (k(x) ≡ 1, σ = 0, top row); a dynamic, invulnerable habitat (σ = 0, middle row); and a dynamic vulnerable habitat (σ = 0.5, bottom row). In each case the initial condition is the unexploited stock condition. This initial condition in the latter two cases is shown in Fig. 1 (left panel). Other parameter values as in Table 1 with P K = 0 Fig. 3 Optimal distributions of fishing effort (top row), stock density (middle row), and habitat quality (bottom row) for increasing levels of habitat vulnerability (σ = 0, 0.5, 1, 5, 10 from left to right). Initial conditions are an unexploited equilibrium (Fig. 1, left panel). Other parameters as in Table 1 with P K = 0 size. As the end of the time horizon approaches, the in situ value of the stock evaporates and as a result the central reserve gradually closes. 1 These results are consistent with the results in Joshi et al. (2009) andKelly et al. (2015).
In our simulations, the central reserve is also present, and is approximately the same size, when habitat quality is dynamic but invulnerable to fishing damage (i.e., when σ = 0; Fig. 2, middle row). There are small decreases in the present value of the fishery, average stock size, and total effort. These decreases are the result of reductions in habitat quality near the boundaries.
Habitat vulnerability Next, we assume that habitat quality is not only dynamic, but also is negatively impacted by harvesting (i.e., σ > 0; Fig. 2, bottom row; Fig. 3). The net present value and average stock size are both smaller in this scenario than when habitat is invulnerable (σ = 0). These reductions are primarily caused by the erosion of habitat quality towards the end of the time horizon. Interestingly, the total optimal effort is largest in this scenario and the increase comes largely from a reduction in the size and duration of the central closed area. This pattern does not continue; for larger habitat vulnerability, total effort decreases again (Fig. 4).
The optimal harvesting strategy in our model is particularly affected by the degree of habitat vulnerability (σ ). All else equal, more vulnerable habitats (or, equivalently, those subject to more destructive gear) are less valuable, and exhibit lower average stock densities, lower average habitat quality, and (usually) lower total effort (Fig. 4). In addition to affecting these aggregate quantities, the habitat vulnerability also affects the spatiotemporal distribution of fishing effort (Fig. 3). As habitat vulnerability increases, the central reserve becomes smaller and shorter in duration. For extreme sensitivities, a reserve is no longer optimal and the optimal strategy is to essentially extirpate the stock.
Habitat recovery rate An important timescale in our model is the time required for habitat quality to recover when harvesting ceases (1/r 2 ). This time ultimately affects the ability of the stock to grow and recover within the planning horizon. In effect, the economic timescale (T ) is longer if r 2 is large than it is when r 2 is small. (The discount rate will also play a role here.) As a result, increasing the recovery rate (r 2 ) increases the size and duration of reserves (Fig. 5). The net present value of the fishery, the stock density and the habitat quality all increase with habitat quality recovery rate (Online Resource 1). Conversely, effort remains essentially constant. Note that these results derive in part from the change in equilibrium initial conditions associated with the change in parameter.

Existence value of habitat quality
In the preceding results, we assumed that the only value that could be derived from the coupled stock-habitat system (1)-(8) was produced by the capture and sale of the stock. We now consider cases when habitat quality has additional intrinsic value, for example, because it provides for recreation or other ecosystem services. In these cases, P K > 0. Not surprisingly, our numerical results show a direct relationship between the intrinsic value of habitat quality and the extent and duration of no-take reserves (Fig. 6). As a direct consequence, the total effort varies inversely with habitat quality value. The net present value of the fishery, as well as the stock size and habitat quality, also increase monotonically with P K (Online Resource 1).
Initial conditions Because our model is dynamic, the optimal control will typically depend upon the initial data. In the previous examples, we took the stock and habitat quality distributions to be what they would be if they had been unharvested for a long time. Many marine fish stocks currently face a drastically different situation: they have been over-harvested for a long time. The classical example of over-harvesting is produced by the "tragedy of the commons" which results from unregulated open access (Clark 1990). In our model, open-access harvesting (see (54)) results in an extremely low stock Fig. 4 Effect of habitat vulnerability (σ ) on optimal levels of the net present value (J ), total harvesting effort (h tot ), average stock density (u avg ), and average habitat quality (k avg ). Parameters as in Table 1, with P K = 0. Initial conditions were for an unexploited fishery as in the left panel of Fig. 1  Fig. 5 Optimal distributions of fishing effort (top row), stock density (middle row), and habitat quality (bottom row) for increasing levels of habitat quality recovery rate; r 2 = [0.5, 1, 5, 10] from left to right. Other parameters as in Table 1, with P K = 0 population density as well as a lower habitat quality than would be found in an unexploited stock (Fig. 1, right  panel).
Open-access initial conditions also generate a very different pattern of optimal harvesting: there is no harvesting at the beginning of the time interval, allowing the  Table 1. Initial conditions were for an unexploited fishery as in the left panel of Fig. 1 stock to replenish. Eventually, harvesting (at maximal rates) opens up at the edges of the habitat and expands towards the center of the domain as the time horizon approaches (Fig. 7). A central reserve area persists almost until the end.

Discussion
Our bioeconomic analysis of the coupled stock-habitat quality model (1)-(11) has revealed several important conclusions. Some of these echo and reaffirm the conclusions of Neubert (2013, 2015); one would seem to stand in contrast to their results.
First, we, like Moeller and Neubert, found that the optimal spatiotemporal allocation of fishing effort can include no-take reserve areas. Our results add to a growing theoretical literature demonstrating that no-take reserves can theoretically have economic as well as conservation benefits (e.g., Neubert 2003;Sanchirico et al. 2006;Costello and Polasky 2008;White et al. 2008;Ding and Lenhart 2009;Hastings et al. 2017). Evidence for these conservation benefits in the real world seems to be clear (Costello 2014;Baskett and Barnett 2015;Sala and Giakoumi 2017). The  Table 1 with P K = 0 and T = 16. Initial conditions are illustrated in the right panel of Fig. 1 extent to which the theoretical economic benefits have materialized, or have the potential to materialize, in practice, however, is currently a topic of debate (Hilborn 2017; Pendleton et al. 2018). When there is a mismatch between potential and realized benefits, some of the difference can be attributed to imperfect implementation (Agardy 2018).
Also in agreement with Moeller and Neubert, we found that the optimal effort distribution can change dramatically when the effects of habitat damage are incorporated into the model. We extended the analysis of Moeller and Neubert by including habitat dynamics in our model to account for the fact that the cessation of harvesting does not result in an instantaneous restoration of habitat quality. These habitat dynamics introduce two new timescales in our model: the timescale over which harvest at a given intensity of effort h decreases habitat quality (1/σ h) and the characteristic timescale for habitat to recover when harvesting stops (1/r 2 ). As we highlighted in Figs. 3 and 5, the relationship between these two timescales (and the length of the planning horizon T ) plays an important role in structuring the optimal harvest, including the size and duration of no-take reserves.
If fishing is particularly destructive (i.e., habitat quality can be reduced quickly by harvesting, and the habitat recovery time is relatively long) then reserves are suboptimal for maximizing the net present value over a finite time horizon (Fig. 3, rightmost column). If, on the other hand, habitat quality is relatively resilient (i.e., habitat responds slowly to fishing and recovers relatively quickly, then reserves are necessary for optimal harvesting (Fig. 5, rightmost columns). This pattern would seem to contradict the findings of Neubert (2013, 2015) who found that the more damaging fishing was, the more likely it was that reserves would be part of a rent-maximizing effort distribution. This apparent disagreement is resolved, however, when one accounts for the fact that the objective in the papers by Moeller and Neubert was either the maximization of (undiscounted) sustainable equilibrium rent (in Moeller and Neubert 2013) or of the net present value over an infinite time horizon (in Moeller and Neubert 2015). Further, in both of these papers, habitat recovery was essentially instantaneous. Thus, habitat was always apparently relatively resilient in these previous analyses due to the effectively infinite time horizon. We conjecture that a significantly longer time horizon, possibly combined with a lower discount rate, would restore the optimality of reserves for vulnerable habitat in the model analyzed here.
Another conclusion (which probably does not need to be explained) is that optimal reserves are prominent when intrinsic habitat value is high. This result is consistent with the findings of Viana et al. (2017), who found the size of reserves should be larger with increasing tourism value, which was tied to increased stock density. In their analysis, however, tourism and fishing were essentially incompatible; that is, tourism only generated value when it occurred within a reserve. In the present model, habitat quality has value irrespective of fishing intensity. This would be the case if the value of habitat quality flowed from ecosystem services that were compatible with fishing. Incorporating incompatible services in our model would require the price of habitat quality to depend nonlinearly on harvest effort.
These conclusions and conjectures must be evaluated in light of the standard caveats that accompany any mathematical modeling. For example, our model has 18 parameters, so our exploration of the model's outcomes was necessarily limited to a very small parameter set (Table 1). In keeping our focus on the effects of habitat vulnerability and dynamics, we have ignored the effects of habitat size, advection, and of potentially fascinating geometric effects that might have been revealed by the analysis of a twodimensional habitat. (We note that it would not be difficult to extend the optimality system to two spatial dimensions (see, e.g., Joshi et al. 2009).) We have also only considered the case when congestion costs are small. As a result, the objective functional (9) is nearly linear in the control and the optimal effort distribution is approximately "bang-bang," with effort levels either at zero or at their maximum h max (cf. Figs. 3, 5, and 6). Solutions would doubtless change quantitatively if congestion costs were higher (as in Neubert and Herrera 2008) or if the effort constraint were relaxed. It remains for now an open question whether such adjustments would change our qualitative conclusions.
Another caveat is that our model is admittedly (and intentionally) stylized. Some extensions that have the potential to change or generalize our results, such as two spatial dimensions, could, in principle, be straightforward to implement. Others will require more thought. For example, our models for the dynamics of "habitat quality," and for the way habitat quality impacts the dynamics of the stock, are strictly phenomenological. We structured these models so that they would have certain properties that seemed reasonable and, to the extent that habitat quality modifies stock carrying capacity, matched with previous theoretical treatments and empirical findings (Fogarty 2005;Shephard et al. 2010;Foley et al. 2012). Habitat quality might also affect other vital rates (e.g., the density-independent birth and death rates that make up the intrinsic growth rate r 1 ) or the dispersal behavior of the fish-a topic that deserves further investigation (but see Langebrake et al. 2012). In our opinion, this phenomenological approach is the correct starting point. Nevertheless, a more concrete and mechanistic connection between habitat quality and the vital rates of the fish would be beneficial-particularly in an applied setting.
Finally, lest it escape the reader's attention, we note that all four of our indices of aggregate "bioeconomic health"-present value, total effort, average stock density, and average habitat quality-are significantly larger when fishing does not damage habitat (Fig. 4). This suggests that there may be substantial economic value, as well as ecological value, in improving and implementing technologies that reduce the impact of fishing gear.