Time delays and pollution in an open‐access fishery

We analyze the impacts of pollution on fishery sector using a dynamical system approach. The proposed model presupposes that the economic development causes emissions that either remediate or accumulate in the oceans. The model possesses a block structure where the solutions of the rate equations for the pollutant and the economic activity act as an input for the biomass and effort equation. We also account for distributed delay effects in both the pollution level and the economic activity level in our modeling framework. The weight functions in the delay terms are expressed in terms of exponentially decaying functions, which in turn enable us to convert the modeling framework to a higher‐order autonomous dynamical system by means of a linear chain trick. When both the typical delay time for the economic activity and the typical delay time for the pollution level are much smaller than the biomass time scale, the governing system is analyzed by means of the theory for singularly perturbed dynamical systems. Contrary to what is found for population dynamical systems with absolute delays, we readily find that the impact of the distributed time lags is negligible in the long‐run dynamics in this time‐scale separation regime.


| INTRODUCTION
Marine pollution is a widespread problem affecting the ocean health, the health of marine species, food security, and quality, and it contributes to climate change. As reported by Halpern et al. (2008), Islam and Tanaka (2004), Jones and Reynolds (1997), Todd et al. (2010), Vikas and Dwarakish (2015), and Lu et al. (2018), marine pollution is increasing and causes negative effects on commercial fisheries. A discussion of different forms of emission of pollution into the marine environment and their potential to be broken down can be found in Watson et al. (2016).
In a dynamic analysis of the interaction between pollution and fisheries, possible time delays are relevant. For the self-cleaning ability of waste both the present state as well as the past history are expected to influence the ongoing remediation capacity. Also for the growth of fish resources it is reasonable to expect that pollution history matters. Furthermore, when it comes to both fishermen and consumers, their behavior is affected by the present and the former economic conditions.
Economists have for a long time been engaged in modeling time delays in macro-and microeconomic models to reveal the consequences for the existence and stability of equilibrium states in economic systems. An early contribution dealing with this problem can be found in Frisch and Holme (1935). Examples on more recent papers focusing on time delays in economic models are Keller (2010), Bianca et al. (2013), and Gori et al. (2018). In bioeconomics, several authors have analyzed the effects of absolute time delays both in single-population models and in predator-prey models (see, e.g., Bretschger & Smulders, 2018;Chakraborty et al., 2011Chakraborty et al., , 2012Ferrara et al., 2019;Gazi & Bandyopadhyay, 2008;Kar & Chakraborty, 2010;Martin & Ruan, 2001;Yuan et al., 2016;Q. Zhang et al., 2011;Y. Zhang et al., 2014). These papers show that the incorporation of absolute time delay effect results in a more complicated dynamics than what one obtains in standard population dynamical systems without delay. The incorporation of time delays in bioeconomic models is also found in more detailed discrete time simulations (see, e.g., Dolder et al., 2020;Gomes et al., 2021). Dolder et al. (2020) assume a given rate of delay when they simulate the fish population movement in space and time, whereas Gomes et al. (2021) assume an absolute time delay in the way temperature affects the biological efficiency for fish species. Notice that such absolute time delays are often motivated by specific biological mechanisms (e.g., gestation period or reaction time related to for instance reproductive or consumption behavior). The problem with absolute time delays is that in economic and biological systems it is not possible to give a precise estimate of the time lags, however. On the basis of this observation, it is argued that distributed time delay effects provide more realistic descriptions in economic-growth models. See, for example, Guerrini et al. (2020; and some of the references therein).
This serves as a motivation for the present study: We analyze how a single species open-access commercial fishery is affected by marine pollution. Our model which is expressed in terms of a dynamical system is based on the simplifying assumption that production per capita in the economy grows with an exogenous given rate, and reaches an exogenous given saturation level in the long run. Furthermore, we assume that marine emissions are explained by means of this growth. Our objective is not to explain economic growth, but to analyze the possible consequences of economic growth for fisheries.
The modeling framework which we will present possesses a block structure: The solutions of the rate equations for the pollutant and the economic activity act as an input for the biomass and effort equation. This structure is reflected in the analysis of the existence and stability of the equilibrium states.
The proposed dynamical model deals with possible impacts of distributed time delays (instead of absolute time lags) on the fishery concerned. We incorporate these effects by expressing the pollutant level as a weighted decomposition of the instantaneous pollutant level and the corresponding time history of this level. The same type of decomposition applies to the economic activity state. We furthermore assume the weight functions in these delays to be modeled by means of exponentially decaying functions. This type of time delay modeling means that all earlier states matter, but the nearest history carries most weight. Similar ideas concerning the ways time delays may influence the dynamics are found in Guerrini et al. (2020) analyzing economic growth. The assumption of exponentially decaying weight functions, enables us to convert the governing system of equations to a higher-order autonomous dynamical system by using the linear chain trick (LCT; Cushing, 2013). The resulting dynamical system is investigated with respect to existence and stability of equilibrium states. Moreover, using this conceptual way of treating time delays makes it possible to consider separation of time scales related to the different variables. More specifically, we expect the typical timescale for the biomass evolution to be much greater than the typical delay time scales for the pollutant and economic activity level. In this time-scale separation regime we show that the modeling framework becomes a singularly perturbed dynamical system. This system is analyzed by means of the Fenichel theory for singularly perturbed systems (Fenichel, 1979) and its generalization to multiple time scales (Cardin & Teixeira, 2017). We explore and illuminate important features of the modeling framework by means of the outcome of this asymptotic analysis and numerical simulations and show that the impact of the time history is negligible in the long run. Also here the block structure of the model plays an important role. As far as we know, this type of analysis has not been carried out for similar models earlier.
The present paper is organized as follows: In Section 2 we discuss our modeling framework in detail which is a modified version of Bergland et al. (2019). We transform the original system to a dimensionless system by means of scaling in Section 3 and introduce the aforementioned timescale separation regimes from the view point of bioeconomics. Section 4 is devoted to the analysis of the subsystem describing the economic-growth and pollution part of the model, whereas we carry out the analysis of the full model in Section 5. The analysis shows that the results obtained in Section 4 can be used as input in the full model presented in Section 5. Section 6 contains concluding remarks and a list of problems for future investigations. Appendix A contains the proof of the persistency property of the model. In Appendix B we deal with the application of the well-known Fenichel theory for one-parameter singularly perturbed dynamical systems to the economic-growth and pollution block of the modeling framework, whereas the same problem is addressed for the full modeling framework in Appendix C.

| MODELING FRAMEWORK
Our model consists of two main parts. One part models economic activity (production per capita) Y t ( ) and stock of pollutant in the marine environment S t ( ), whereas the other part describes a modified fishery model that constitutes the impact of pollution and economic growth on the fish stock X t ( ) and dynamic effort variable E t ( ). Without any time lag effects involved, the actual model reads with nonnegative initial conditions X E S Y (0) > 0, (0) > 0, (0) > 0, (0) > 0. The last two equations in the model (1) account for the assumption that production and consumption activities are responsible for the emission of harmful substances in the marine environment, and model the pollution absorptive capacity of the environment. The instantaneous change in the production of pollution in the marine ecosystem is proportional to the production per capita Y with an emission rate ϱ 0 . The nonmonotonic function ρ represents the remediation capacity function which is given as The remediation capacity function ρ which models the pollution absorptive capacity of the environment, thus depends on the size of the pollutant density. It constitutes a two-parameter family of functions parameterized by means of the remediation parameters a and b. Here the nonnegative function ρ is a twice differential function of S and has a global maximum value at a unique point S S . Here we notice that ρ S a b S ( ; , ) = 1 2 max max ∕ . Furthermore, ρ is a strictly increasing function for S S 0 max ≤ ≤ and strictly decreasing for S S max ≥ . Notice that the function ρ also satisfies the homogeneity properties, that is, Hence, when the pollution level S is below the threshold value S max , the self-cleaning ability of the environment will increase with S. However, when exceeding this pollution threshold, the self-cleaning ability is expected to decrease and is negligible if pollution concentration in the marine ecosystem is high enough. A notable feature is that the remediation capacity function evaluated at S max is inverse proportional to S max . The last equation in the model (1) states that the change in production (and income) per capita is governed by the economic-growth function f Y Y Y Y ( ;¯) =¯− and the economic-growth rate g, where the parameter Ȳ is the unique zeros of the economic-growth function f . Moreover, Notice that the function f satisfies the homogeneity property, that is, where C is a positive parameter. It is worth mentioning that parameters Ȳ and C have the same dimension. 1 The assumption of modeling an exogenous growth rate seems adequate as the production from the fishery only is assumed to contribute marginally to the total economic activity in the society. At the same time, we assume that emissions from fishing activities are insignificant compared with emissions from overall production and consumption in the economy. For the other part of the model we assume a modified logistic growth equation for the development of the fish biomass with the negative term βS − incorporated to account for the impact of pollution on the fish biomass. The parameter r measures the intrinsic growth rate of fish biomass and K is the carrying capacity. The coupling parameter β is referred to as the growth-retardation parameter. The model next considers the fish harvesting H where the fish supply from the harvesting follows the Gordon-Schaefer production function H qXE = . The instantaneous change in the effort variable E is assumed to be the proportion of the timedependent sector profit π PH cE = − , where P is the time-dependent unit product fish price and c cost per unit price of fishing effort. The assumptions for the function P that represents the time-dependent marginal willingness to pay for fish products, include market demand-income-mechanisms and negative impact on fish demand due to pollution.
Inspired by the conjecture that the distributed delay framework appears to be more realistic and meaningful in the biological as well as in the economical context, we incorporate distributed time delays in the different modeling blocks of the model (1). First, we assume that the prehistory of the cleaning ability will affect the pollutant level S and should therefore be accounted for in the model extension. We take this property into account by making the replacement in the pollution equation of the model (1). Here S   is the temporal mean defined as Thus, the new pollution equation which accounts for the prehistory of cleaning ability reads Next we take into account the impact of the distributed time delay in the saturation term of the fish biomass equation in (1). We do this by making the extension of the biomass equation. Here the positive parameters r, K , and β have the same meaning as in (1). For a nonzero β the term β ωS ω S K + (1 −   in (8) thus measures a possible loss in fish biomass depending upon the accumulated pollution and its corresponding past history. Notice that this term is a convex combination of the instantaneous effect of the pollution (modeled by means of S and the time history of the pollution development, represented by the mean value S  ). The convexity parameter ω measures the degree of influence of the distributed time delay effect: For ω = 1, no such effect is present, whereas for the other extreme case ω = 0 the impact of the pollution effect on the wild fish population solely depends on the past history.
Regarding the price function P in the model (1), we assume that the price of the fish product is affected negatively by pollution and positively by economic activity (income per capita). The time-dependent marginal willingness to pay function P reads Here Y   is the temporal mean of the economic activity level Y defined as Moreover, we have assumed that the price function P depends on the pollutant density S and its corresponding time history S  , that is, Here F is defined as which ensures that the function P 1 is positive for all t. Thus we have assumed that distributed time lags are present in both the pollution level S and the income level Y in a way similar to the time lag effects in the biomass-pollution part of the modeling framework when it comes to the marginal willingness to pay function P. The case ν = 1 means that the effect of the past pollution history is negligible. To model P 1 we tacitly assume that accumulated pollution concentration affects the quality of the fish products and that it reduces the consumers beliefs in paying higher prices for the fish products. Therefore we assume that P 1 decreases for accumulating harmful substances in the marine ecosystem and it saturates to level A a B − > 0 0 1 0 . Here the positive constant A 0 measures a fixed market price for per capita income (Y ). The nonnegative parameters a 1 and B 0 are responsible for modeling the consumers beliefs in the fish products. It promotes the willingness of consumers to pay less prices for such harmful fish products (Chen et al., 2015;Fonner & Sylvia, 2014;Garzon et al., 2016;Wessells & Anderson, 1995;Whitehead, 2006). Here we also assume that the history of pollution might affect the demand for fish, by incorporating a distributed time effect ( ν 0 < 1 ≤ ). With the assumption that the higher the income per capita, the higher aggregated demand for fish is, we model per capita economy production or the general income level per capita by mean of state variable Y . We model the positive influence of income level per capita Y on consumers demands for fish products that also includes a possible time delay income effect and an ordinary income impact on fish demand. By inserting the expression (9) for P into the rate equation for E in the model (1), we end up with the distributed time delay equation for the effort variable E. Regarding the temporal evolution of the economical activity level (production per capita), we retain the same rate equation as in the model (1), that is, we assume that with an economic-growth rate g. By making the additional assumption that the weight functions in the mean values S   and Y   defined by (6) and (10) are exponentially decaying functions, that is, and the LCT produces the rate equations for the mean values S   and Y  . See Cushing (2013) for a detailed exposition on LCT.
To summarize, the modeling framework consists of the pollution dynamics equation (7), the fish biomass equation (8), the effort equation (13), and the production and income equation (14). The rate equation (17) for the mean value of the pollutant density and the rate equation (18) for the mean value of the economic activity level represent past history in the respective model variables. Thus our model comprises a 6D autonomous nonlinear dynamical system. The biological and/or economical interpretation of variables and the parameters of the proposed model are summarized in Table 1.

| Scaling and general properties of the model
We start out by observing that the exponentially decaying integral kernels α γ and α μ defined by (15) and (16) obey the scaling laws T A B L E 1 Interpretation of variables and parameters associated with the fishery-pollution models (8), (13), (7), and (14) and rate equations (17)  The scaled exponentially decaying kernels α ε 1 γ ∕ and α ε 1 μ ∕ are normalized functions. Thus they can be used as weight functions in integrals that define the temporal means θ   and ψ   of θ and ψ, respectively, that is, where the functions θ   and ψ   are defined in Table 2. Notice here that we can define θ   and ψ   by means of (20) for any normalized integral kernels α γ and α μ that satisfy the scaling laws (19). We interpret the dimensionless parameters ε γ and ε μ : T γ = 1 γ ∕ is the decay time of the weight function α γ in the temporal mean S   of the pollutant density S, whereas T r = 1 r ∕ is the logistic timescale of the biomass X . Thus the ratio ε γ is the ratio between these two time scales: Thus ε γ and ε μ parametrize the temporal delay effects in our modeling framework. These two parameters (and consequently also the three timescales T μ , T γ , and T r ) will play a crucial role in the forthcoming analysis. By means of the substitutions defined in Table 2 we readily derive the nondimensional version of the distributed delay system consisting of Equations (8), (7), (13), and (14). Here , , , and  are the functions defined as Here the notation ′ means differentation with respect to the time variable τ. Table 2 summarizes the definitions and the interpretations of all the dimensionless variables and parameters involved in the nondimensionalization process of our modeling framework. Imposing the homogeneity assumptions (3) and (4) on the economic-growth function f and the remediation capacity ρ we get R θ γ ρ θ γ aρ aθ a b ( , )˜( ; 1, ) = ( ; , )).

Variable/parameter definition Interpretation
Relative unit cost of effort in fishery remains in r for all τ τ* ≥ , where τ* 0 ≤ , that is, Σ is invariant under the solution flow of the systems (21) and (22). For the sake of completeness we have relegated the details of this proof to Appendix A (Theorem 1).
In the forthcoming analysis we refer to the parameter ι γ γ 5 2 ≡ as the emission-remediation ratio. By restoring to the original dimensional parameters we readily find that Notice that the inequality (12) which is a sufficient condition for positive willingness to pay translates into the condition γ γ.
We observe that the functions  and  defined by means of (25) and (26), respectively, are independent of ξ and η. This means that we can view the solutions of the subsystem as input functions to the subsystem We further notice that the component functions , , and  in (35) are independent of the mean value ψ  . This means that the solution of the subsystem (35) can be determined by first obtaining the solution to the 3D system and thereafter finding the temporal mean ψ   by solving the linear rate equation using ψ as an input function. Figure 1 summarizes schematically the relationship between the different blocks in the present model. (21)-(30) is of particular interest from the viewpoint of bioeconomics: According to the previous discussion of the logistic timescale T r and the delay time scales T μ and T μ , we observe that the parameter regime This regime can be analyzed directly by means of the standard theory for singularly perturbed problems (Fenichel, 1979). If we impose the additional scale separation ε ε γ μ (which can be motivated from viewpoint of bioeconomics), we end up with the two-parameter singularly perturbed system. Multiscale singularly perturbed dynamical systems can be dealt with using a generalization of Fenichel's theory (Cardin & Teixeira, 2017). In subsequent sections we will deal with these two parameter regimes.

| THE ECONOMIC-GROWTH AND THE POLLUTION PART OF THE MODEL
Here we discuss the properties of (35) and examine the impact of Group 1 parameters and delay weight parameter p on the dynamics of the subsystem (35). We start out by studying the existence of equilibrium points. An equilibrium point   without any restriction. Hence the equilibrium problem boils down to a study of the equilibrium states θ γ ( , ) e 5 of the 2D system We make use of the properties of this system when dealing with the study of the slow-fast regime (35) We denote the possible solutions of this equation as θ e . Let us define a function Γ : and introduce the subsets , , non eq,1 eq,2    , and  defined as Here the curve  represents the nontransversality condition ( ) Notice that the region non  produces no equilibrium points whereas the parameter region eq,2  produces exactly two equilibrium points. eq,1  is the region which gives rise to one negative equilibrium point. Now, let us assume that ( ) γ ι , 3 2 eq,1 eq,1    ∈ ∪ ∪ which means that we are guaranteed the existence of at least one equilibrium point (and maximum two equilibrium points) of the type θ γ ( , ) e 5 of the 2D system (38). The stability problem of (38) can be resolved with the help of the Jacobian matrix The eigenvalues of the above stability matrix are computed as γ R θ γ − ′( ; )  Figure 2 shows a numerical example with the phase portrait of the 2D system (38)  , where θ e satisfies (39). The Jacobian matrix for the 3D subsystem (36) at 3 is given as is the 2 × 2-block matrix given as Here and in the sequel we will often make use of the abbreviation R . The eigenvalues of (43) are given as Thus the spectral properties of p 3 ( ) show that equilibrium 3 is asymptotically stable if to a saddle-node bifurcation. This means that the weight parameter p does not change the stability properties as compared with the case p = 1. Thus the stability properties of the 2D subsystem (38) carry over to the 3D subsystem (36). We notice that = p 2 ( =1) 2 , where 2 is given by (42). The system (36) becomes a decoupled system in the case p = 1, in which the dynamical evolution of the state variables θ and ψ is governed by the 2D system (38), whereas the mean value θ   can be computed by solving the linear rate equation ε θ θ θ ′ = − γ     by direct integration treating θ as an input function obtained by solving (38).
Similarly, we linearize the 4D subsystem (35) around the equilibrium 4 to resolve the stability problem. One can easily compute the eigenvalues of the following Jacobian matrix to examine the stability of subsystem (35)  Here we have tacitly assumed that the ( ) . We also find that one eigenvalue becomes zero if R θ γ ′( ; ) = 0 e 3 . Therefore the equilibrium 4 is locally asymptotically stable if R θ γ ′( ; ) > 0 e 3 , whereas it will be a saddle point for the case R θ γ ′( ; ) < 0 e 3 . Simple computation reveals that points located on the curve  correspond to a saddle-node bifurcation. It is worth noting that also in this case qualitative properties of the subsystem (35) follow the dynamical properties of the 2D subsystem (38).
In Figure 3 we present numerical results associated with a subsystem (35)  Here we see that the dynamical variable ψ and its temporal mean ψ   converge to the equilibrium value ψ = 1.1 e . However, the variable θ and its temporal mean θ   increase as time progresses (see Figure 3b). Obviously these findings are consistent with our theoretical outcomes. Here a moderate initial production per capita promotes stability. On the other hand, the accumulated pollution goes unbounded when the initial per capita production is high enough. A similar situation is depicted in Figure 3c when initial production per capita exceeds a certain threshold.

| Slow-fast analysis
The subsystem (35) shows that the economic activity ψ and the pollution level θ do not depend on the past history of the economic activities. Therefore we will not include the equation of temporal mean ψ   when studying the slow-fast progression analysis of the subsystem associated with the economic activity and the pollution part. A detailed mathematical analysis of a slow-fast version of the subsystems (35) and (36) is appended in Appendix B.
Here we show that numerical simulations of the dynamics of the slow-fast version of the subsystem (36) (i.e., the scaled subsystem B2) confirm the mathematical findings in Appendix B. This is illustrated in of 3 . During this initial phase the variables θ and ψ stay almost constant. In the second phase they approach very slowly the asymptotically stable equilibrium on the slow manifold a  as predicted by the reduced problem. In Figure 4b one can see that ψ (red curve) approaches its equilibrium value ψ = 1.1 e slowly. Moreover, the mean value θ   (yellow curve) rapidly converges to θ (blue curve) following the dynamics of the fast subsystem on the slow-manifold a  . However, θ grows unbounded following the dynamics of the differential algebraic equation  Uncontrolled pollution in case of initially high economic activity. Moreover, the economic activity ψ stabilizes as time progresses.
In this section we first discuss the existence and stability of equilibrium points of the full systems (21)-(30). If any such equilibrium points exist, then they must be in the form where θ e is a solution of Equation ( (21)-(32) to have equilibrium states. We keep this condition in mind during the upcoming analysis of the proposed model.
In Section 5.1 we study the possibility of having equilibrium points of the system located in the hyperplanes ξ = 0 and η = 0 with the corresponding stability properties, whereas we in Section 5.2 deal with conditions ensuring the existence of equilibrium points in the first orthant. In the latter subsection we address the stability of the interior equilibrium points. The stability assessment of the equilibrium points (49) is inferred from the spectral properties of the Jacobian matrix. We also illustrate the theoretical findings by means of numerical simulations in these subsections. where the functions  and  are defined by means of (23) and (24), respectively. Assume that ξ = 0 e . We readily find that η θ γ θ γ γ (0, , , , , ) = − < 0

| Existence and stability of equilibrium points in the hyperplanes
For the purpose of investigating the stability of Q 0 , we first derive the Jacobian 6 evaluated at Q 0 . Simple computation reveals that where 0 m n × denotes the zero m n × matrix. The eigenvalues of the Jacobian matrix (51) are given by means of (47) and represent bifurcation points for which the stability assessment cannot be based on the linearization procedure.
In Figure 5, we numerically illustrate the stability of the equilibrium point Q 0 and try to examine how the pollution and the economic activity part affect the fishery dynamics. Here we observe that for a certain set of parameter values and initially low pollution density, the economic activities and the pollution level are sustained and remain in the vicinity of their respective equilibrium values (θ ψ , e e ). We also observe the eradication of the fish biomass together with a vanishing of the harvesting effort consistent with the stability properties of the equilibrium point Q 0 . In the case the pollution goes uncontrolled we notice a rapid extinction of the fish population and the vanishing of the harvesting effort (see Figure 5b).
Let us next consider the hyperplane η = 0. In this case we either get ξ γ θ = 1 − , respectively. Notice that if γ θ 1 h 1 ≤ , we will have γ θ 1 l 1 ≤ Thus we will have the coexistence of two equilibrium points of the type Q 1 for the case γ θ 1 h 1 ≤ and θ θ < l h . We next compute the Jacobian 6 evaluated at Q 1 : always is fulfilled. Figure 6 demonstrates the impact of the pollution on the stability properties of the equilibrium point Q 1 . In Figure 6a we observe that a small perturbation around Q 1 does not harm its stability. Notice that this long-run evolution of the pollution level makes the fishery unprofitable, thus resulting in a vanishing of the effort even though the fish stock has not gone extinct. However, Figure 6b shows that an initially high pollution level also causes the extinction of the fish biomass. This means that the solution is attracted towards the equilibrium Q 0 .

| Existence and stability of positive equilibrium points
Here we discuss the existence and stability of positive equilibrium points, that is, equilibrium points of the type (49) for which ξ η θ , , > 0 e e e . By using the equilibrium condition . (a) Initial condition (0.4, 0.4, 0.6, 0.6, 0.6, 0.6) belongs to the attractor basin of Q 0 . The fishery sector is wiped out and the effort vanishes. (b) Initial condition (0.4, 0.4, 1.2, 0.6, 0.6, 0.6) does not belong to the attractor basin of Q 0 . The fish stock is depleted very quickly and the pollution grows uncontrolled due to high initial pollution production. Here ξ represents fish biomass, η is the effort variable, θ measures the concentration of pollution, and ψ is the economic activity.
will always be fulfilled. This means that we are guaranteed the existence of positive equilibrium points Q e provided the subsystem (36) has at least one positive equilibrium point.
Remark 1. We notice that the first inequality in (57) is equivalent to the boundedness condition for the equilibrium biomass X e and the equilibrium pollutant density S e (corresponding to ξ e and θ e , respectively) when restoring to the dimensional parameters. The result (58) shows that X e is bounded by the carrying capacity K whereas S e is bounded by K β ∕ . The second inequality in (57) is the nondimensional version of the positive marginal willingness to pay, that is, the condition, By restoring to the dimensional parameters, we find that the equilibrium effort E e and the equilibrium population density X e (corresponding to Equation 56) are given as Here P e is the marginal willingness to pay given as (59), also interpreted as the market value per unit catch. We notice that the equilibrium biomass X e and the equilibrium fishery effort E e are determined by fish value, cost, harvest efficiency, and growth parameters in a similar manner as in the standard Gordon-Schaefer model (e.g., Clark, 2010;Flaaten, 2018;Hannesson, 1978). For instance, it follows that the special case with no biomass growth damage (β = 0) our equilibrium (corresponding to the effort and stock level) is referred to as the bioeconomic equilibrium of the unregulated open-access fishery in the standard Gordon-Schaefer model. The cost-value ratio ( c P e ) is crucial when determining the profitability in the fishery, and thereby the equilibrium biomass and the equilibrium fishery effort. If X X > e , fishing is profitable and the effort increases causing a reduction of the biomass. As soon as X falls below X e , fishing becomes unprofitable and the effort will be withdrawn. Now to determine the stability of the positive equilibrium point Q e we first derive the Jacobian matrix 6 evaluated at this point. We readily obtain The eigenvalues of stability matrix (61) are given by and those eigenvalues which can be extracted from submatrix p 4 ( ) , that is, the expressions (47).
From Remark 1 we know that Q ( ) > 0 ξ e  (see condition 57). We observe that eigenvalues λ e ± remain negative if the condition ξ Q > 4 ( ) e ξ e  is satisfied. Combining this condition with (57) we get the following restriction: which ensures negative values of λ e ± . Here the structure of the function F θ γ ( ; ) e 8 is defined in (12). Remark 2 gives a detailed bioeconomic interpretation of the condition (64). It is noteworthy that the condition (64) is interlinked with the stability dynamics of the fisheries and economic effort part which are automatically satisfied if the condition (57) is satisfied. Assuming that the condition (64) is fulfilled, the stability assessment of the positive equilibrium point Q e simplifies the study of the slope of the remediation capacity R evaluated at θ e : Q e is a locally asymptotically stable equilibrium if the positive slope condition R θ γ ′( ; ) > 0 e 3 holds true. A saddle-node bifurcation occurs when R θ γ ′( ; ) = 0 e 3 . Q e becomes unstable if R θ γ ′( ; ) < 0 e 3 . We illustrate this result numerically. To this end, we run simulations for a set of parameters satisfying the restrictions stated in (56) which ensure the positivity of the equilibrium point Q e . In Figure 7, one can easily see how an initially high pollutant level may alter the dynamics of a stable equilibrium point. In Figure 7a, we observe that a small disturbance of the equilibrium point Q e does not harm the stability of the system. However, an initially high pollutant density promotes uncontrolled pollution which may harm the stability of the system. It is worth mentioning that a high initial pollution level may give rise to the bifurcation from the interior equilibrium point Q e to the fisheries-effort eradication equilibrium point Q 0 (see Figure 7a).
We also signify the effect of Group 3 parameters on the stability of the interior equilibrium point Q e . We numerically observe that Group 3 parameters have no impact on the equilibrium value of the state variables θ ψ , . Therefore we only capture the time evolution of state variables ξ and η in the upcoming simulations. Remark 2. By restoring to the dimensional parameters, we find that the inequality (64) translates into the boundedness condition P qK Here P e is the marginal willingness to pay defined by means of (59). We notice the appearance of the market value per unit, in addition to the harvest efficiency (q), the growth parameters (r and K ), and the friction in effort adjustment (k) in this expression. We conclude from the above analysis that the equilibrium Q e is asymptotically stable as long as the market value P e in the fishery is below this threshold level and for small and moderate values of the equilibrium pollutant level.

| Slow-fast analysis
In this subsection we apply geometrical singular perturbation theory to investigate the dynamics of slow-fast progression in the original system. We discuss different cases depending on the values of parameters μ γ and . ≪ . This regime can be investigated by means of Fenichel's theory (Fenichel, 1979) and its generalization to multitime-scale problems (Cardin & Teixeira, 2017) for slow-fast systems. The application of this theory to our model is relegated to Appendix C.

| Case I: T T T
Here we focus on how the numerical simulation matches the theoretical findings presented in Appendix C. To this end we consider the flow of the slow part of the system (35) (i.e., the subsystem C1) on the 4D critical manifold 0  defined by This evolution is governed by the system Here we see that the economic activity variable ψ evolves independently of the remaining state variables as time progresses. Moreover it is always attracted to its equilibrium ψ γ = e 5 . We recall the results obtained in previous sections to determine the stability of equilibrium points of the system (67) is unstable within the framework of the model (67). Notice that this stability result is independent of ε γ . The mathematical analysis carried out in Appendix C suggests that the impact of distributed delay over the long period of time is negligible in the time-scale regime T T T μ γ r ≪ . To illustrate the mathematical observations appended in Appendix C we run numerical simulations for the same set of parameter values and initial data underlying the computations leading to Figure 7. We observe a rapid attraction of fast variables θ   and ψ   towards their equilibrium points. The slow variables evolve very slowly and almost look like constant in this phase (see Figure 8).
In Figure 9, we present numerical results for the regime ε ε 0 < = 1 μ γ ≪ showing the temporal evolution of the biomass ξ and the effort variable η by assuming the ) and initial data are the same as in Figure 7a,b. Here ξ represents fish biomass, η is the effort variable, θ measures the concentration of pollution, and ψ is the economic activity. . The initial condition is the same as in Figure 7a. Here ξ represents fish biomass, η is the effort variable, θ measures the concentration of pollution, and ψ is the economic activity.
pollution-economic activity block as input information to the fishery part. In Figure 9a we observe the well-known oscillations where the profitability in the fishery determines the effort, and the following harvest volume causes biomass variations. Figure 9 displays the temporal evolution of the slow variable which is evolved on slow manifold ε γ  . These simulations demonstrate how a general trajectory of the system (67)   ∈       of 6 acts as a normally attracting hyperbolic manifold. The dynamical evolution in the next phase takes place on 1  and is governed by the system which is a singularly perturbed system with ε γ as a perturbation parameter. As ε 0 γ → , the asymptotic approximation of the dynamical evolution prescribed by (68)  Here ⋅ denotes the time derivative with respect to τ where time variable τ τ ε =˜μ ∕ and other state variables are transformed similarly as in (B1). Notice that 1  and 2  are the critical and two-critical manifolds on which the intermediate and the reduced problems of systems (21)-(30) are defined, respectively. Notice that 1  and 2  are also the equilibrium points of the layered and the intermediate problem of the systems (21)-(30). By appealing to Theorems A, B, C, and Corollary 1 in Cardin and Teixeira (2017) we conclude that for small parameters ε μ and ε γ , there is an invariant manifold under the flow of the systems (21)-(30) which is the perturbed counterparts of the manifold 2  with Hausdorff distance ε ε ( + ) μ γ  and, it is diffeomorphic to 2  . Moreover, the stability property of ε ε 2 ( , ) μ γ  is inherited from the attracting critical manifold 2  . We omit the proof of multitime-scale convergence results for the slow-fast systems as these results follow directly from Theorems A, B, and C in Cardin and Teixeira (2017).
We finally run simulations for the singularly perturbed original systems (21)-(30) in the parameter regime ε ε 0 < < 1 μ γ ≪ . The outcome of these simulations is summarized in Figure 10. Here we observe that the mean variables θ   and ψ   converge to the trajectories of θ and ψ very quickly, respectively. However, the slow variables evolve very slowly on the same time interval in comparison. These variables remain almost constant in the initial phase (see magnified view in Figure 10a). However, in the long run the slow variables converge to their respective equilibrium values. The numerical findings thus agree surprisingly well with the predictions deduced by means of the theory for singularly perturbed systems, even though ε ε = 5 γ μ in the numerical simulations. We conclude that neither the past history of economic activities nor the past history of pollution growth will affect the fishery dynamics in the long run in this time separation regime. Moreover, high initial pollution level gives rise to uncontrolled pollution growth which in turn promotes unprofitable fisheries and could be the reason for the extinction of the fish biomass. ) and the initial data are the same as in Figure 7a,b. Here ξ represents the fish biomass, η is the effort variable, θ measures the concentration of pollution, ψ is the economic activity, and τ is the time scale for the equivalent fast system of the original models (21)-(30). 6 | CONCLUDING REMARKS AND EXTENSIONS

| Main results
In the present paper we have investigated a conceptual causal bioeconomic model describing the impact of pollution on a fishery using a dynamical system approach. We have accounted for distributed time histories in both the pollution level and the economic activity level in this modeling framework. This means that the fundamental model assumes the form of a distributed time delay system. The model decomposes into two subsystems. The first subsystem which describes the interaction between the pollution and the economic activity plays the role as an input to the second subsystem which is defined as rate equations for the fish stock density and the effort variables.
The model is converted to a 6D autonomous dynamical system by means of the LCT when assuming the weight functions in the time delay terms to be exponentially decaying functions (Cushing, 2013). Existence and stability of equilibrium states are investigated. We also identify the characteristic time scales for these time lags and investigate the dynamical evolution as a function of these time scales. When the former time scales are much smaller than the logistic time scale, the modeling framework becomes a singularly perturbed dynamical system. We have analyzed this system using Fenichel's theory (Fenichel, 1979) and its generalization to multiscale problems (Cardin & Teixeira, 2017) together with numerical simulations. On the basis of this analysis and the simulation results we conclude that the prominent feature to be observed is that the impact of the distributed time lags is negligible in the long-run dynamical evolution, contrary to what is found for population dynamical systems with absolute time lags (Bretschger & Smulders, 2018;Chakraborty et al., 2011Chakraborty et al., , 2012Ferrara et al., 2019;Gazi & Bandyopadhyay, 2008;Kar & Chakraborty, 2010;Martin & Ruan, 2001;Yuan et al., 2016;Q. Zhang et al., 2011;Y. Zhang et al., 2014). Hence, our model analysis adds to the existing literature knowledge of how distributed time delays affect short-and long-run development of stock and effort in an open-access harvesting model.

| Possible extensions
Here we will point to possible extensions and modifications of the present modeling framework which could be topics for future research.
First, in the present paper we have assumed that the emission rate ϱ 0 in the pollution equation is constant, compare Equation (7). A natural extension of this assumption consists of modeling the emission rate by means of a strictly decreasing, which is bounded from above and below. The motivation for assuming a decrease in the emission rate as a function of the pollutant density and the economic activity goes as follows: First, when the production increases (which implies a higher income in the economy), it seems reasonable to expect that the individuals' preferences for reducing the pollution become relatively more important. If such reasoning holds, the environmental policies will become stricter. Second, an increase in emissions may also lead to an increased willingness to implement stricter environmental standards to avoid the negative impacts from pollution. However, it is reasonable to believe that the emission rate remains positive and saturates when assuming such decreasing effects.
Second, we have assumed the time histories to be modeled by means of exponentially decaying functions. We have argued that this assumption mimics a realistic time history when it comes to economic activity and pollution emission. We do not exclude the possibility that other time histories can be relevant, however. One possible future research topic consists of assuming that the time histories can be modeled as a Taylorpolynomial multiplied with an exponentially decaying function. This type of time history can model situations where certain times in the past carry more weights than others when it comes to the influence of the time history on the present situation. The LCT is applicable in this case. See Cushing (2013). It will produce a higher hierarchy of ordinary differential equations (ODEs). Here a detailed numerical examination of the change in the temporal development of the state variables would be of interest. One should, however, be aware of that the change to this type of weight function will not alter the fundamental stability properties of the detected equilibrium states. Notice also that the time-scale separation regimes which also in this case can be analyzed by means of Fenichel's theory for singularly perturbed dynamical systems (and its extension to multiscale problems) will also lead to the conclusion that the impact of the delay effects on the dynamical evolution will be negligible in the long run.
The argument for this proceeds in the same way as in Sections 4 and 5. Notice that it is possible to extend the LCT method to generally distributed time lags as suggested in Ponosov et al. (2004). However, in this case LCT may produce an infinite dimensional system of ODEs, so that it may be a mathematical challenge to perform stability analysis, and there is no guarantee that the properties of the present model would be preserved under such a farreaching generalization. However, if the lag distribution is polynomial-based, the LCT gives a finite-dimensional system, and the stability analysis can be performed in the standard way, as explained above.
Third, a natural and well-motivated extension consists of incorporating spatial effects in the modeling framework. This means to consider effects like local as well as nonlocal advection-diffusion effects in the description of the biomass evolution as well as in the pollutant movement. We believe that such spatial effects are not present in the rate equations for the economic activity and the effort variables. Here we could follow the same line of thought as in Marciniak-Czochra et al. (2017) and Heilmann et al. (2018). This means that the modeling framework assumes the form of an ODE-reaction diffusion coupled model with advection effects incorporated. The key problem here will be to detect linear instabilities of equilibrium states as a route to pattern formation. The outcome of this process is expected to be stationary spatial patterns, standing, and traveling waves, and possibly also chaotic spatiotemporal patterns. It will thus be of interest to find out if the resulting system permits coherent structures like stable traveling wave solutions. The role of such structures in the pattern-forming process should be examined. Here we should notice that nonlocal diffusion effects incorporated in the biomass and the pollutant equation in a way analogous to Cheng and Yuan (2017) could result in the existence of traveling waves. Another interesting problem for further investigation is to study the proposed model under biomass diffusion as an optimal fisheries management problem. The main idea is to investigate the effects of optimal control of biomass diffusion processes in the emergence of spatial heterogeneity through diffusion-driven instability. Optimal control in harvesting and diffusion in the emergence of spatial heterogeneity can be studied based on the footprints of William and Xepapadeas (2010).
Fourth, we will expect that the determination of the maximal remediation capacity as well as other parameters will be subject to uncertainties, and thus the stochastic processes should be included in the modeling framework. Hence the modeling framework becomes a system of stochastic differential equations. Following the thoughts as in Evans (2012) and Øksendal (2003), the model can be rewritten as an autonomous dynamical system of first-order equations with the stochastic effects accounting for additive noise terms in the system.
In addition to more sophisticated and realistic assumptions regarding the pollution and the cleaning abilities of the ecosystem, alternative and more reasonable assumptions regarding the economic behavior and the functioning of the fish market also represent possible model extensions. For instance, modifications of Verhulst's population growth model and the Gordon-Schaefer harvest production function which are discussed in Hannesson (1978), Clark (2010), Flaaten (2018), or Eide (2018, could be incorporated in the model formulation.
Seeking for an overall welfare optimal resource allocation, environmental, and fishery policies based on the predictions derived from a model like ours would also be an interesting task for future research.
To sum up, our comments regarding possible extensions may serve as a relevant background for future studies concerning marine pollution and fishery resources within the framework of bioeconomic models.
2. If for some τ 0 the components ξ and η of the solution of the systems (21) and (22) are positive, then they remain positive for all τ τ 0 ≥ even if the other components may be negative.
The first statement follows from the assumption that f ψ γ( ; ) > 0 slow manifold or critical manifold (Fenichel, 1979). Notice the role of 0  : It plays the role as the equilibrium state of the subsystem Θ = Θ − Θ , Ψ = Ψ − Ψ         of (C4). On the fastest time scale, the state variables ξ , η, θ, and ψ are constant, in accordance with the system (C4). The slow manifold 0  is a normally attracting hyperbolic manifold. Therefore the flow on the perturbed slow manifold is diffeomorphic to the associated slow flow on the critical manifold for sufficiently small ε γ (Fenichel, 1979).