Analysis of Mathematics and Sustainability in an Impulsive Eutrophication Controlling System

Abstract

Eutrophication removal problems have captured the attention of biologists, mathematicians, and environmental scientists. Within this framework, an impulsive eutrophication controlling system is studied analytically and numerically. A key advantage of the eutrophication system is that it can be quite accurate to describe the interaction effect of some critical factors (fishermen catch and releasing small fry, etc.), which enables a systematic and logical procedure for fitting eutrophication mathematical system to real monitoring data and experiment data. Mathematical theoretical works have been pursuing the investigation of two threshold functions of some critical parameters under the condition of all species persistence, which can in turn provide a theoretical basis for the numerical simulation. Using numerical simulation works, we mainly focus on how to choose the best value of some critical parameters to ensure the sustainability of the eutrophication system so that the eutrophication removal process can be well developed with maximizing economic benefit. These results may be further extended to provide a basis for simulating the algal bloom in the laboratory and understanding the application of some impulsive controlling models about eutrophication removal problems.

1. Introduction

The eutrophication of lakes and reservoirs has been widely and intensively studied for several decades, which is a degradation process originating from introduction of nutrients from agricultural runoff and untreated industrial and urban discharges [1]. The eutrophication has been identified as a major and serious water quality management issue [2], not only in the developing countries [3–5], but also in the developed countries [6–8]. Impairment of water quality due to eutrophication can lead to a series of problems and result in losses of ecological integrity, sustainability, and safe use of aquatic ecosystems [9]. It is characterized by frequently recurring algal blooms, zooplankton aggregation, and reduction in species diversity in water bodies at all trophic levels [2].

In 1932, Bertalanffy firstly proposed a method of using mathematical model to study the biological system [10]. He pointed out that this approach is a combination of coordination, order, and purpose, which in turn formed three basic ideas of studying the biological system, namely, trophic-level analysis, system perspective, and dynamic view. The main aim of modeling population dynamics is to improve the understanding of the interactions between populations and their dependence on internal and external conditions [11]. Hence, the mathematical models of biological population dynamics have to account not only for the growth and interactions between populations but also for the forecast of the disasters which can be caused by the plankton like blue-green algae, Moina and Brachionus. Early attempts began with the logistic growth, exponential growth, and Cui growth (this growth can explain the relationship between population increment and limiting resources) for biological population. With the increasing effects of the eutrophication, zooplankton aggregation, and algal blooms, the history of mathematical modeling plankton or zooplankton dynamics and biological eutrophication removal process is already quite long and has been initiated by biological science. These approaches have been refined to more realistic descriptions of the development of biological natural populations, especially algae population. In recent decades, the main researches focused on how to control the eutrophication, how to predict the outbreak of algae, how to simulate the spread tendency of algae, and how to develop eutrophication control strategy based on actual data. Hence, more and more researches have given an attempt to combine the forces of ecologist and mathematical engineers in an evaluation of possible control strategies for the management of eutrophication of water bodies, which is becoming an urgent water-management problem of drinking water reservoirs [12]. Malmaeus and Håkanson have depicted the development of a lake eutrophication model; they indicated that all changes that provided significant differences and improved predictions have been adopted [13]. Chen et al. said that their mathematical description provides a different insight into how to use nitrogen and phosphorus more effectively in the eutrophication control and elimination of blue-green blooms [14].

Zeya reservoir of Wenzhou is located in subtropical regions. Because of eutrophication, the nuisance algal blooms and the serious zooplankton aggregation frequently come forth, which caused clogging and blocking of the filtration system and resulted in millions of people without drinking water. Eutrophication removal is achieved by two major processes, physicochemical and biological treatment techniques; specifically biological eutrophication removal from drinking water in the lake and reservoir is usually considered to accomplish optimal and economic eutrophication treatment [15]. In order to understand and improve the performance of biological eutrophication removal treatment, it is necessary to briefly look into the mechanisms of eutrophication removal. The removal process of eutrophication in the Zeya reservoir is complex and dynamic with many variables. Since successful management of eutrophication in Zeya reservoir is often dependent on proper identification of the most likely growth limiting amount for algae, special emphasis is given to the limiting amount concept. In order to apply the biological principle to control eutrophication in the Zeya reservoir, we give good biological control agents for algae. First, we bring up two species of filter-feeding fish, silver carp (Hypophthalmichthys molitrix), and bighead carp (Aristichthys nobilis) [16], which are thought to be suited to control algae and zooplankton biomass directly in freshwater reservoir, especially cyanobacteria. Second, because the decision-maker must give cost-benefit and maintenance considerations, fish should be carried out and controlled at an appropriate time. At the same time, in order to prevent the rapid growth of plankton, especially algae, and maintain a good circulating development of biological process, we must release a certain amount of filter-feeding fish at a fixed time. These processes are subject to short-term perturbations whose durations are too small compared to the duration of the process. Consequently, it is natural to assume that these perturbations act instantaneously in the form of impulses. Thus impulsive differential equations and differential equations involving impulsive effects appear as a natural mathematical description (model) of observed evolution phenomena of several real world problems [17–21]. These works have introduced the theory of impulsive differential equations and its applications and provided a theoretical basis of the application of impulsive differential equations. In recent years, Fan et al. [22] have proposed the dynamics of a class of mutualism-competition-predator interaction models with Beddington-DeAngelis functional responses and impulsive perturbations; this work is important for solving the global stability and the globally exponential stability of system. Yu and Lu [23] mainly focused on the permanence and global attractivity of a discrete almost periodic ratio-dependent Leslie system with time delays and feedback controls; these works vigorously promote the application of feedback control. Shen et al. [24] mainly have studied the permanence of a single species system with distributed time delay and feedback controls; it is a very important fact that the feedback control is harmless to the permanence of species; these works are significant for the application of feedback control in biological control. Furthermore, the main aim of this investigation is to answer that what is the fundamental mathematical controlling mechanism, which can support the biological principle to control eutrophication in the Zeya reservoir and to provide a different insight into how to make the effects of environmental factors in connection with the ability of controlling eutrophication.

Speaking on eutrophication of Zeya reservoir, as the algae population grows to a certain biomass level, the dominant zooplankton species (Moina, Cyclops, and Brachionus are mainly grazing algae) are bound to sire and grow, so that these results can lead to damage or even crash of aquatic ecosystems. Furthermore, based on the above discussion, modeling is playing an increasing role in helping to identify causes and effective control measures through what if studies, and strategies are being identified for the management and control of eutrophication [25, 26]; the paper considers an eutrophication mathematical model, which can be described by the following differential equations: 𝑑𝑥(𝑡)𝑑𝑡=𝑟₁𝑥(𝑡)𝑥(𝑡)1−𝑘₁−𝑎₁𝑥(𝑡)𝑦(𝑡)𝑏₁(𝑡)−𝑎+𝑥²𝑥(𝑡)𝑧(𝑡)𝑏₂(𝑡),(𝑡)+𝑥𝑑𝑦𝑑𝑡=𝑟₂𝑦(𝑡𝑦(𝑡))1−𝑘₂+𝑒₁𝑎₁𝑥(𝑡)𝑦(𝑡)𝑏₁(𝑡)−𝑎+𝑥³𝑦(𝑡)𝑧(𝑡)𝑏₃(𝑡),(𝑡)+𝑦𝑑𝑧=𝑎𝑑𝑡²𝑒₂𝑥(𝑡)𝑧(𝑡)𝑏₂(𝑡)+𝑎+𝑥³𝑒₃𝑦(𝑡)𝑧(𝑡)𝑏₃(𝑡)(𝑡),()(𝑡)+𝑦−𝑚𝑧𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,Δ𝑥=−𝛿₁𝑥(𝑡),Δ𝑦(𝑡)=−𝛿₂𝑦(𝑡),Δ𝑧(𝑡)=−𝛿₃𝑧(𝑡),()(𝑡)(𝑡)(𝑡)𝑡=𝑛+𝑙−1𝑇,Δ𝑥=0,Δ𝑦=0,Δ𝑧=𝑝,𝑡=𝑛𝑇, (1) where𝑥(𝑡),𝑦(𝑡), and𝑧(𝑡)are the biomass of the algae population (cyanobacteria) (10⁶ cell/L), the dominant zooplankton (Moina, Cyclops, and Brachionus) (10² piece/L), and the filter-feeding fish (silver carp and bighead carp) (1 fish/Cubic meter) at time𝑡, respectively,Δ𝑥(𝑡)=𝑥(𝑡⁺)−𝑥(𝑡),Δ𝑦(𝑡)=𝑦(𝑡⁺)−𝑦(𝑡), andΔ𝑧(𝑡)=𝑧(𝑡⁺)−𝑧(𝑡).0<𝑙<1is used to describe the intervals of time between the pulsed use of controls.𝑟_𝑖(𝑖=1,2) are the intrinsic growth rate,𝑎_𝑖(𝑖=1,2,3) and𝑏_𝑖(𝑖=1,2,3) measure the efficience of the prey in evading a predator attack.𝑒_𝑖  (𝑖=1,2,3) denote the efficiency with which resources are converted to new consumers,𝑘_𝑖(𝑖=1,2) are carrying capacity in the absence of the filter-feeding fish, and𝑚is the mortality rates for filter-feeding fish. The decision-maker of biological treatment technology must implement an impulsive control strategy that can better maintain biological treatment process, so we not only release a certain amount of the filter-feeding fishr𝑧(𝑝>0is the release amount of the filter-feeding fishr𝑧) at𝑡=𝑛𝑇to graze the algae and the dominant zooplankton, but also consider the algae and the dominant zooplankton losses due to other reasons as well as the filter-feeding fish harvest at𝑡=(𝑛+𝑙−1)𝑇(𝛿_𝑖(𝑖=1,2,3) (0<𝛿_𝑖<1) represent the fraction of three species), where𝑇is the period of the impulsive effect,𝑛∈𝑁, and𝑁is the set of all nonnegative integers.

The main purpose of this paper is to enhance a theoretical and constructive framework by making the attempt to deal with eutrophication controlling problem for some drinking reservoirs using an impulsive eutrophication controlling system. Meanwhile, the theory discussions on the global asymptotical stability and permanent would help to design some efficient parameter constraint representations, which can in turn provide a theoretical basis for the numerical simulation. Further, in the context of population growth dynamics, how to implement impulsive control strategy to prevent and control the algal bloom is studied by simulating the dynamics of the impulsive controlling eutrophication system.

2. Mathematical Analysis

Now we will pursue the investigation of the constraint expression of the release amount𝑝and the average harvest𝛿₃of the filter-feeding fish. Let𝑅₊=[0,∞), and𝑅³₊={𝑋∈𝑅³∣𝑋≥0}. Denote𝑓=(𝑓₁,𝑓₂,𝑓₃)as the map defined by the right hand of the first, second, and third equations of system (1). Let𝑉:𝑅₊×𝑅³₊→𝑅₊; then𝑉is said to belong to class𝑉₀if(1)𝑉is continuous in(𝑛𝑇,(𝑛+𝑙−1)𝑇)∪((𝑛+𝑙−1)𝑇,(𝑛+1)𝑇]×𝑅³₊, and for each𝑋∈𝑅³₊,𝑛∈𝑁,lim_{(𝑡,𝑦)→(𝑛𝑇}⁺_,𝑋)𝑉(𝑡,𝑦)=𝑉(𝑛𝑇⁺,𝑋)exists;(2)𝑉is locally Lipschitzian in𝑋.

Definition 1. Let𝑉∈𝑉₀; then for(𝑡,𝑥)∈(𝑛𝑇,(𝑛+𝑙−1)𝑇)∪((𝑛+𝑙−1)𝑇,(𝑛+1)𝑇]×𝑅³₊, the upper right derivative of𝑉(𝑡,𝑋)with respect to the impulsive differential system (1) is defined as 𝐷⁺𝑉(𝑡,𝑋)=lim_ℎ→0+1supℎ𝑉𝑡+ℎ,𝑋+ℎ𝑓(𝑡,𝑥)−𝑉(𝑡,𝑋). (2)
The solution of system (1) is a piecewise continuous function𝑋∶𝑅₊→𝑅³₊,𝑋(𝑡)is continuous on(𝑛𝑇,(𝑛+𝑙−1)𝑇)∪((𝑛+𝑙−1)𝑇,(𝑛+1)𝑇],𝑛∈𝑁, and𝑋(𝑛𝑇⁺)=lim_{𝑡→𝑛𝑇}⁺𝑋(𝑡)exists. The smoothness properties of𝑓guarantee the global existence and uniqueness of solution of system (1); for details see [8, 9].

Definition 2. System (1) is said to be permanent if there exists a compactΩ⊂int𝑅³₊such that every solution(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡))of system (1) will eventually enter and remain in the regionΩ.

The following lemma is obvious.

Lemma 3. Let𝑋(𝑡)be a solution of system (1) with𝑋(0⁺)≥0; then𝑋(𝑡)≥0for all𝑡≥0, and further𝑋(𝑡)>0,𝑡>0if𝑋(0⁺)>0.

We will use an important comparison theorem on impulsive differential equation.

Lemma 4 (Bainov and Simeonov, 1993 [17]). Suppose𝑉∈𝑉₀. Assume that 𝐷⁺𝑉()((()𝑉𝑡𝑡,𝑋≤𝑔𝑡,𝑉𝑡,𝑋)),𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,𝑡,𝑋⁺≤𝜓¹_𝑛(𝑉(()𝑉𝑡𝑡,𝑋)),𝑡=𝑛+𝑙−1𝑇,𝑡,𝑋⁺≤𝜓²_𝑛(𝑉(𝑡,𝑋)),𝑡=𝑛𝑇, (3) where𝑔∶𝑅₊×𝑅²₊→𝑅²₊is continuous in((𝑛−1)𝑇,(𝑛+𝑙−1)𝑇)and((𝑛+𝑙−1)𝑇,𝑛𝑇]for𝑢∈𝑅²₊,𝑛∈𝑁,lim_{(𝑡,𝑣)→(𝑛𝑇}⁺_,𝑢)𝑔(𝑡,𝑣)=𝑔(𝑛𝑇⁺,𝑢), andlim_{(𝑡,𝑣)→((𝑛+𝑙−1)𝑇}⁺_,𝑢)𝑔(𝑡,𝑣)=𝑔((𝑛+𝑙−1)𝑇⁺,𝑢)exists;𝜓^𝑖_𝑛(𝑖=1,2)∶𝑅₊→𝑅²₊is nondecreasing. Let𝑟(𝑡)be maximal solution of the scalar impulsive differential equation 𝑑𝑢(𝑡)((𝑡()𝑢𝑡𝑑𝑡=𝑔𝑡,𝑢)),𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₊=𝜓¹_𝑛(𝑢(𝑡()𝑢𝑡)),𝑡=𝑛+𝑙−1𝑇,₊=𝜓²_𝑛(𝑢(𝑡𝑢0)),𝑡=𝑛𝑇,₊=𝑢₀, (4) existing on[0,∞). Then𝑉(0⁺,𝑋₀)≤𝑢₀implies that𝑉(𝑡,𝑋(𝑡))≤𝑟(𝑡),𝑡≥0, where𝑋(𝑡)is any solution of system (1).

Finally, we give some basic properties about the following subsystems of system (1): 𝑑𝑧(𝑡)(𝑡)()𝑧𝑡𝑑𝑡=−𝑚𝑧,𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₊=1−𝛿₃𝑧(𝑡)()𝑧𝑡,𝑡=𝑛+𝑙−1𝑇,₊(𝑡)0=𝑧+𝑝,𝑡=𝑛𝑇,₊=𝑧₀. (5)

Clearly 𝑧^∗(𝑡)=𝑝exp(−𝑚(𝑡−(𝑛−1)𝑇))1−1−𝛿₃()()()𝑝exp−𝑚𝑇𝑛−1𝑇<𝑡≤𝑛+𝑙−1𝑇,1−𝛿₃exp(−𝑚(𝑡−(𝑛−1)𝑇))1−1−𝛿₃()()exp−𝑚𝑇𝑛+𝑙−1𝑇<𝑡≤𝑛𝑇, (6)𝑛∈𝑁,and 𝑧^∗0⁺=𝑧^∗𝑛𝑇⁺=𝑝1−1−𝛿₃(),𝑧exp−𝑚𝑇_∗𝑙𝑇⁺=𝑝1−𝛿₃exp(−𝑚𝑙𝑇)1−1−𝛿₃exp(−𝑚𝑇) (7) is a positive periodic solution of system (5). Since 𝑧(𝑡)=1−𝛿₃^𝑛−1𝑧0⁺−𝑝1−1−𝛿₃()×()exp−𝑚𝑇exp−𝑚𝑡+𝑧^∗(𝑡()()),𝑛−1𝑇<𝑡≤𝑛+𝑙−1𝑇,1−𝛿₃^𝑛𝑧0⁺−𝑝1−1−𝛿₃()()exp−𝑚𝑇×exp−𝑚𝑡+𝑧^∗(𝑡),(𝑛+𝑙−1)𝑇<𝑡≤𝑛𝑇 (8) is the solution of system (5) with initial value𝑧₀=𝑧(0⁺)≥0, where𝑡∈((𝑛−1)𝑇,(𝑛+𝑙−1)𝑇)∪((𝑛+𝑙−1)𝑇,𝑛𝑇],𝑛∈𝑁, we get the following.

Lemma 5. For a positive periodic solution𝑧^∗(𝑡)of system (5) and every solution𝑧(𝑡)of system (5) with𝑧₀≥0, one has|𝑧(𝑡)−𝑧^∗(𝑡)|→0,𝑡→∞.

Lemma 6. There exists a constant𝑀such that𝑥(𝑡)≤𝑀,𝑦(𝑡)≤𝑀,𝑧(𝑡)≤𝑀, and𝑝(𝑡)≤𝑀for each solution(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡),𝑝(𝑡))of system (1) with all𝑡largest enough.

Therefore, we obtain the complete expression for a periodic solution(0,0,𝑧^∗(𝑡))of system (1). Now, we study the stability of the periodic solution, which explains the threshold expression of the release amount𝑝and the average harvest𝛿₃of the filter-feeding fish under the condition of some species extinction.

Theorem 7. Let(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡))be any solution of system (1); then(0,0,𝑧^∗(𝑡))is said to be globally asymptotically stable if 𝑟₁𝑎𝑇−²𝑏₂𝑚𝑝1−1−𝛿₃()×exp−𝑚𝑇1−𝛿₃exp(−𝑚𝑙𝑇)−1−𝛿₃()1exp−𝑚𝑇<ln1−𝛿₁,𝑟₂𝑎𝑇−³𝑏₃𝑚𝑝1−1−𝛿₃()×exp−𝑚𝑇1−𝛿₃exp(−𝑚𝑙𝑇)−1−𝛿₃()1exp−𝑚𝑇<ln1−𝛿₂. (9)

Proof. The local stability of periodic solution(0,0,𝑧^∗(𝑡))may be determined by considering the behavior of small amplitude perturbation of the solution. Define 𝑥(𝑡)=𝑢(𝑡),𝑦(𝑡)=𝑣(𝑡),𝑧(𝑡)=𝑤(𝑡)+𝑧^∗(𝑡). (10)
We will put (10) into (1). The linearization of the system becomes 𝑑𝑢(𝑡)=𝑟𝑑𝑡¹−𝑎₂𝑧^∗(𝑡)𝑏₂𝑢(𝑡),𝑑𝑣(𝑡)=𝑟𝑑𝑡²−𝑎₃𝑧^∗(𝑡)𝑏₃𝑣(𝑡)𝑑𝑤(𝑡)(𝑡),()(𝑡)𝑑𝑡=−𝑚𝑤𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,Δ𝑢=−𝛿₁𝑢(𝑡),Δ𝑣(𝑡)=−𝛿₂𝑣(𝑡),Δ𝑤(𝑡)=−𝛿₃𝑤(𝑡),()(𝑡)(𝑡)(𝑡)𝑡=𝑛+𝑙−1𝑇,Δ𝑢=0,Δ𝑣=0,Δ𝑤=0,𝑡=𝑛𝑇. (11)
Therefore, we have 𝑢(𝑡)𝑣(𝑡)𝑤(𝑡)(𝑡)𝑢(0)𝑣(0)𝑤(0)=Φ,0≤𝑡<𝑇, (12) whereΦ(𝑡)satisfies 𝑑Φ=𝑟𝑑𝑡¹−𝑎₂𝑧^∗(𝑡)𝑏₂000𝑟₂−𝑎₃𝑧^∗(𝑡)𝑏₃0Φ(𝑡)00−𝑚 (13) andΦ(0)=𝐼, the identity matrix, and 𝑢(𝑛+𝑙−1)𝑇⁺𝑣(𝑛+𝑙−1)𝑇⁺𝑤(𝑛+𝑙−1)𝑇⁺=1−𝛿₁0001−𝛿₂0001−𝛿₃𝑢)𝑇)𝑣)𝑇)𝑤)𝑇),𝑢((𝑛+𝑙−1((𝑛+𝑙−1((𝑛+𝑙−1𝑛𝑇⁺𝑣𝑛𝑇⁺𝑤𝑛𝑇⁺=𝑢()𝑣()𝑤().100010001𝑛𝑇𝑛𝑇𝑛𝑇 (14)
The stability of the periodic solution(0,0,𝑧^∗(𝑡))is determined by the eigenvalues of Θ=1−𝛿₁0001−𝛿₂0001−𝛿₃Φ(𝑇)100010001. (15)
If each of these eigenvalues is less than one in magnitude, then the periodic solution(0,0,𝑧^∗(𝑡))is locally stable. Therefore all eigenvalues ofΘ, namely, 𝑢₁=1−𝛿₁exp^𝑇₀𝑟₁−𝑎₂𝑧^∗(𝑡)𝑏₂,𝑢𝑑𝑡₂=1−𝛿₂exp^𝑇₀𝑟₂−𝑎₃𝑧^∗(𝑡)𝑏₃,𝑢𝑑𝑡₃=1−𝛿₃exp(−𝑚𝑇)<1. (16)
According to Floquet theory of impulsive differential equation,(0,0,𝑧^∗(𝑡))is locally asymptotically stable if|𝑢₁|<1and|𝑢₂|<1; that is to say, 𝑟₁𝑎𝑇−²𝑏₂𝑚𝑝1−1−𝛿₃()×exp−𝑚𝑇1−𝛿₃exp(−𝑚𝑙𝑇)−1−𝛿₃()1exp−𝑚𝑇<ln1−𝛿₁,𝑟₂𝑎𝑇−³𝑏₃𝑚𝑝1−1−𝛿₃()×exp−𝑚𝑇1−𝛿₃exp(−𝑚𝑙𝑇)−1−𝛿₃()1exp−𝑚𝑇<ln1−𝛿₂. (17)
In the following, we will prove the global attraction. Choose a𝜀>0such that 𝜉₁≜1−𝛿₁exp^𝑇₀𝑟₁−𝑎₂𝑧^∗(𝑡)−𝜖𝑏₂+𝑀𝑑𝑡<1. (18)
Note that𝑑𝑧(𝑡)/𝑑𝑡≥−𝑚𝑧(𝑡); from Lemma 5 and comparison theorem of impulsive equation, we have 𝑧(𝑡)>𝑧^∗(𝑡)−𝜖 (19) for all𝑡large enough. For simplification, we may assume that (19) holds for all𝑡≥0. From (1) and (19), we can obtain 𝑑𝑥(𝑡)≤𝑟𝑑𝑡¹−𝑎₂𝑧^∗(𝑡)−𝜖𝑏₂𝑥(𝑡)()𝑥𝑡+𝑀,𝑡≠𝑛+𝑙−1𝑇,₊=1−𝛿₁𝑥(𝑡),𝑡=(𝑛+𝑙−1)𝑇, (20) which leads to 𝑥((𝑛+𝑙)𝑇)≤𝑥(𝑛+𝑙−1)𝑇⁺×exp_()𝑇^{(𝑛+𝑙)𝑇}_{𝑛+𝑙−1}𝑟₁−𝑎₂𝑧^∗(𝑡)−𝜖𝑏₂)𝑇)+𝑀𝑑𝑡=𝑥((𝑛+𝑙−11−𝛿₁×exp_()𝑇^{(𝑛+𝑙)𝑇}_{𝑛+𝑙−1}𝑟₁−𝑎₂𝑧^∗(𝑡)−𝜖𝑏₂)𝑇)𝜉+𝑀𝑑𝑡=𝑥((𝑛+𝑙−1₁. (21)
Hence,𝑥((𝑛+𝑙)𝑇)≤𝑥(𝑙𝑇)𝜉^𝑛₁and𝑥((𝑛+𝑙)𝑇)→0as𝑛→∞. Therefore,𝑥(𝑡)→0as𝑛→∞since0<𝑥(𝑡)<𝑥((𝑛+𝑙−1)𝑇)(1−𝛿₁)exp(𝑟₁𝑇)for(𝑛+𝑙−1)𝑇<𝑡≤(𝑛+𝑙)𝑇. Similarly,𝑦(𝑡)→0as𝑛→∞.
Next, we prove that𝑧(𝑡)→𝑧^∗(𝑡)as𝑡→∞. For0<𝜖sufficiently small, there must exist a𝑇^′>0such that0<𝑥(𝑡)<𝜖and0<𝑦(𝑡)<𝜖,𝑡≥𝑇^′. Without loss of generality, we may assume that0<𝑥(𝑡)<𝜖and0<𝑦(𝑡)<𝜖for all𝑡≥0; then from system (1) we obtain (𝑡)≤(𝑡)−𝑚𝑧𝑑𝑧≤𝑧(𝑡)𝑒𝑑𝑡−𝑚+𝛾𝜖,𝛾=²𝑎₂𝑏₂+𝑒₃𝑎₃𝑏₃. (22) From Lemmas 4 and 5, we have𝑣₁(𝑡)≤𝑧(𝑡)≤𝑣₂(𝑡)and𝑣₁(𝑡)→𝑧^∗(𝑡),𝑣₂(𝑡)→𝑣^∗(𝑡)as𝑡→∞, where𝑣₁(𝑡)and𝑣₂(𝑡)are solutions of 𝑑𝑣₁(𝑡)𝑑𝑡=−𝑚𝑣₁(𝑡)()𝑣,𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₁𝑡⁺=1−𝛿₃𝑣₁(𝑡)()𝑣,𝑡=𝑛+𝑙−1𝑇,₁𝑡⁺=𝑣₁(𝑡)𝑣+𝑝,𝑡=𝑛𝑇,₁0⁺=𝑧0⁺,𝑑𝑣₂(𝑡)=𝑣𝑑𝑡−𝑚+𝛾𝜖₂(𝑡)()𝑣,𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₂𝑡⁺=1−𝛿₃𝑣₂(𝑡)()𝑣,𝑡=𝑛+𝑙−1𝑇,₂𝑡⁺=𝑣₂(𝑡)𝑣+𝑝,𝑡=𝑛𝑇,₂0⁺=𝑧0⁺, (23)(24) respectively, 𝑣^∗(𝑡)=𝑝exp−𝑚+𝛾𝜖(𝑡−(𝑛−1)𝑇)1−1−𝛿₃𝑇,()()𝑝exp−𝑚+𝛾𝜖𝑛−1𝑇<𝑡≤𝑛+𝑙−1𝑇,1−𝛿₃exp−𝑚+𝛾𝜖(𝑡−(𝑛−1)𝑇)1−1−𝛿₃𝑇,()exp−𝑚+𝛾𝜖𝑛+𝑙−1𝑇<𝑡≤𝑛𝑇. (25)
Then, for any𝜖₁>0there exists a𝑇₁>0such that 𝑧^∗(𝑡)−𝜖₁<𝑧(𝑡)<𝑧^∗(𝑡)+𝜖,𝑡>𝑇₁. (26)
Let𝜖→0; we have 𝑧^∗(𝑡)−𝜖₁<𝑧(𝑡)<𝑧^∗(𝑡)+𝜖₁ (27) for𝑡large enough, which implies𝑧(𝑡)→𝑧^∗(𝑡)as𝑡→∞. This completes the proof.

Now we study the stability of the periodic solution, which explains the threshold expression of the release amount𝑝and the average harvest𝛿₃of the filter-feeding fish under the condition of all species persistence.

Theorem 8. The system (1) is permanent if 𝑟₁𝑎𝑇−²𝑏₂𝑚𝑝1−1−𝛿₃()×exp−𝑚𝑇1−𝛿₃exp(−𝑚𝑙𝑇)−1−𝛿₃()1exp−𝑚𝑇>ln1−𝛿₁,𝑟₂𝑎𝑇−³𝑏₃𝑚𝑝1−1−𝛿₃()×exp−𝑚𝑇1−𝛿₃exp(−𝑚𝑙𝑇)−1−𝛿₃()1exp−𝑚𝑇>ln1−𝛿₂. (28)(29)

Proof. Suppose𝑋(𝑡)=(𝑥(𝑡),𝑦(𝑡),𝑧(𝑡),𝑝(𝑡))is any solution of the system (1) with𝑋(0)>0. From Lemma 6, we assume that𝑥(𝑡)≤𝑀,𝑦(𝑡)≤𝑀, and𝑧(𝑡)≤𝑀with𝑡≥0. From (19), we have𝑧(𝑡)>𝑧^∗(𝑡)−𝜖for all𝑡large enough and some𝑧(𝑡)≥𝑝(1−𝛿₃)exp(−𝑚𝑇)/(1−(1−𝛿₃)exp(−𝑚𝑇))−𝜖≜𝜁₁for𝑡large enough. Thus we only need to find𝜁₂and𝜁₈such that𝑥(𝑡)>𝜁₂and𝑦(𝑡)>𝜁₈for𝑡large enough. Now we will find a𝜁₂.
We will prove this in the following two steps. First, let0<𝜁₄,𝜖₁>0be small enough such that 𝜂₁≜1−𝛿₁exp_()𝑇^{(𝑛+𝑙)𝑇}_{𝑛+𝑙−1}𝑟₁𝜁1−⁴𝑘₁−𝑎₁𝑀𝑏₁−𝑎₂𝑏₂𝑣^∗₃(𝑡)+𝜖₁𝑑𝑡>1. (30)
Easy to prove that𝑥(𝑡)<𝜁₄cannot hold all𝑡≥0. Otherwise, 𝑑𝑧(𝑡)e𝑑𝑡≤−𝑚−²𝑎₂𝑥(𝑡)𝑏₂−𝑒₃𝑎₃𝑦(𝑡)𝑏₃𝑧(𝑡)𝑒≤−𝑚−²𝑎₂𝜁₄𝑏₂−𝑒₃𝑎₃𝑀𝑏₃𝑧(𝑡),()𝑧𝑡𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₊=1−𝛿₃𝑧(𝑡)()𝑧𝑡,𝑡=𝑛+𝑙−1𝑇,₊(𝑡)𝑧0=𝑧+𝑝,𝑡=𝑛𝑇,₊=𝑧₀. (31) Then we have𝑧(𝑡)≤𝑣₃(𝑡)and𝑣₃(𝑡)→𝑣^∗₃(𝑡)(𝑡→∞), where𝑣₃(𝑡)is the solution of 𝑑𝑣₃(𝑡)𝑒𝑑𝑡=−𝑚−²𝑎₂𝜁₄𝑏₂−𝑒₃𝑎₃𝑀𝑏₃𝑣₃(𝑡),()𝑣𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₃𝑡⁺=1−𝛿₃𝑣₃(𝑡)()𝑣,𝑡=𝑛+𝑙−1𝑇,₃𝑡⁺=𝑣₃(𝑡)𝑣+𝑝,𝑡=𝑛𝑇,₃0⁺=𝑧₀, (32)

𝑣^∗₃(𝑡)=𝑝exp−𝑚−𝑒₂𝑎₂𝜁₄/𝑏₂−𝑒₃𝑎₃𝑀/𝑏₃(𝑡−(𝑛−1)𝑇)1−1−𝛿₃exp−𝑚−𝑒₂𝑎₂𝜁₄/𝑏₂−𝑒₃𝑎₃𝑀/𝑏₃𝑇,()()𝑇𝑝𝑛−1𝑇<𝑡≤𝑛+𝑙−11−𝛿₃exp−𝑚−𝑒₂𝑎₂𝜁₄/𝑏₂−𝑒₃𝑎₃𝑀/𝑏₃(𝑡−(𝑛−1)𝑇)1−1−𝛿₃exp−𝑚−𝑒₂𝑎₂𝜁₄/𝑏₂−𝑒₃𝑎₃𝑀/𝑏₃𝑇,(𝑛+𝑙−1)𝑇<𝑡≤𝑛𝑇. (33)

Therefore, there exists a𝑇₁>0such that 𝑧(𝑡)≤𝑣₃(𝑡)<𝑣^∗₃(𝑡)+𝜖₁,𝑑𝑥(𝑡)≥𝑟𝑑𝑡¹𝜁1−⁴𝑘₁−𝑎₁𝑀𝑏₁−𝑎₂𝑏₂𝑣^∗₃(𝑡)+𝜖₁𝑥(𝑡),()𝑥𝑡𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₊=1−𝛿₁𝑥(𝑡),𝑡=(𝑛+𝑙−1)𝑇 (34)(35) for𝑡>𝑇₁. Let𝑁₁∈𝑁and(𝑁₁+𝑙−1)𝑇≥𝑇₁. Integrating (35) on((𝑛+𝑙−1)𝑇,(𝑛+𝑙)𝑇],𝑛>𝑁₁, we can get 𝑥((𝑛+𝑙)𝑇)≥𝑥(𝑛+𝑙−1)𝑇⁺×exp_()𝑇^{(𝑛+𝑙)𝑇}_{𝑛+𝑙−1}𝑟₁𝜁1−⁴𝑘₁−𝑎₁𝑀𝑏₁−𝑎₂𝑏₂𝑣^∗₃(𝑡)+𝜖₁)𝑇)𝑑𝑡=𝑥((𝑛+𝑙−11−𝛿₁×exp_()𝑇^{(𝑛+𝑙)𝑇}_{𝑛+𝑙−1}𝑟₁𝜁1−⁴𝑘₁−𝑎₁𝑀𝑏₁−𝑎₂𝑏₂𝑣^∗₃(𝑡)+𝜖₁)𝑇)𝜂𝑑𝑡=𝑥((𝑛+𝑙−1₁. (36)

Then𝑥((𝑁₁+𝑛+𝑙))𝑇≥𝑋((𝑁₁+𝑙)𝑇)𝜂^𝑛₁→∞as𝑛→∞, which is a contradiction to the boundedness of𝑥(𝑡). Hence there exists a𝑡₁>0such that𝑥(𝑡₁)≥𝜁₄. So, there exists a𝑡₃>𝑡₁such that𝑥(𝑡₃)≥𝜁₄.

Second, if𝑥(𝑡)≥𝜁₄for all𝑡≥𝑡₃, then our aim is obtained. Hence we only need to consider those solutions which leave the region𝑅={𝑋(𝑡)∈𝑅³₊∶𝑥(𝑡)≤𝜁₂,𝑦(𝑡)≤𝜁₈,𝑧(𝑡)≤𝜁₁}and reenter it again. Let𝑡^∗=inf_𝑡≥𝑡3{𝑥(𝑡)<𝜁₄}, there are two possible cases for𝑡^∗.

Case I. Consider𝑡^∗=(𝑛₁+𝑙−1)𝑇,𝑛₁∈𝑁. Then𝑥(𝑡)≥𝜁₄for𝑡∈[𝑡₃,𝑡^∗]and(1−𝛿₁)𝜁₄≤𝑥(𝑡^∗+)=(1−𝛿₁)𝑥(𝑡^∗)<𝜁₄. Select𝑛₂,𝑛₃∈𝑁such that 𝑛₂𝜖−1𝑇>ln₁/𝑀+𝑝−𝑚+𝑒₂𝑎₂𝜁₄/𝑏₂+𝑒₃𝑎₃𝑀/𝑏₃,1−𝛿₁^𝑛2exp𝑛₂𝜂₂𝑇𝜂^𝑛3₁>1−𝛿₁^𝑛2exp𝑛₂+1𝜂₂𝑇𝜂^𝑛3₁>1, (37) where𝜂₂=(𝑟₁(1−(𝜁₄/𝑘₁))−(𝑎₁𝑀/𝑏₁)−(𝑎₂𝑀/𝑏₂))<0. Let𝑇=𝑛₂𝑇+𝑛₃𝑇. We claim that there must be a𝑡₄∈(𝑡^∗,𝑡^∗+𝑇]such that𝑥(𝑡)>𝜁₄. Otherwise, consider (32) with𝑣₃(𝑡^∗+)=𝑧(𝑡^∗+), and we have𝑣₃(𝑡)=1−𝛿₃^{𝑛−(𝑛1+1)}𝑣₃𝑛₁𝑇⁺−𝑝1−1−𝛿₃exp−𝑚+𝑒₂𝑎₂𝜁₄/𝑏₂+𝑒₃𝑎₃𝑀/𝑏₃𝑇𝑒×exp−𝑚+²𝑎₂𝜁₄𝑏₂+𝑒₃𝑎₃𝑀𝑏₃𝑡−𝑛₁𝑇+𝑣^∗₃(𝑡),(𝑛−1)𝑇<𝑡≤(𝑛+𝑙−1)𝑇,1−𝛿₃^(𝑛−1)𝑣₃𝑛₁𝑇⁺−𝑝1−1−𝛿₃exp−𝑚+𝑒₂𝑎₂𝜁₄/𝑏₂+𝑒₃𝑎₃𝑀/𝑏₃𝑇𝑒×exp−𝑚+²𝑎₂𝜁₄𝑏₂+𝑒₃𝑎₃𝑀𝑏₃𝑡−𝑛₁𝑇+𝑣^∗₃(𝑡),(𝑛+𝑙−1)𝑇<𝑡≤𝑛𝑇 (38)and𝑛₁+1≤𝑛≤𝑛₁+𝑛₂+𝑛₃. Then 𝑣₃(𝑡)−𝑣^∗₃(𝑡<)𝑀+𝑝exp−𝑚₊𝑒₂𝑎₂𝜁₄𝑏₂+𝑒₃𝑎₃𝑀𝑏₃𝑡−𝑛₁𝑇<𝜖₁,𝑧(𝑡)≤𝑣₃(𝑡)≤𝑉^∗₃(𝑡)+𝜖₁, (39)(𝑛₁+𝑛₂+𝑛₃)𝑇≤𝑡≤𝑡^∗+𝑇,which implies that (35) holds for𝑡^∗+𝑛₂𝑇≤𝑡≤𝑡^∗+𝑇. As in the first step, we have 𝑥𝑡^∗+𝑛₂𝑇≥𝑥𝑡^∗+𝑛₂𝑇𝜂^𝑛3₁. (40)

The first equation of the system (1) gives 𝑑𝑥(𝑡)≥𝑟𝑑𝑡¹𝜁1−⁴𝑘−𝑎₁𝑀𝑏₁−𝑎₂𝑀𝑏₂𝑥(𝑡),()𝑥𝑡𝑡≠𝑛𝑇,𝑡≠𝑛+𝑙−1𝑇,₊=1−𝛿₁𝑥(𝑡),𝑡=(𝑛+𝑙−1)𝑇. (41) Integrating (41) on[𝑡^∗,𝑡^∗+𝑛₂𝑇], we have 𝑥𝑡^∗+𝑛₂𝑇≥𝜁₄1−𝛿₁^𝑛2exp𝑛₂𝜂₂𝑇. (42) Thus𝑥(𝑡^∗+𝑇)≥𝜁₄(1−𝛿₁)^𝑛2exp(𝑛₂𝜂₂𝑇)𝜂^𝑛3₁>𝜁₄, and this is a contradiction. Let𝑡=inf_{𝑡≥𝑡∗}{𝑥(𝑡)>𝜁₄}, and then𝑥(𝑡)≤𝜁₄,𝑥(𝑡)=𝜁₄, for𝑡∈(𝑡^∗,𝑡); we have 𝑥(𝑡)≥𝜁₄1−𝛿₁^𝑛2+𝑛3exp𝑛₂+𝑛₃𝜂₂𝑇≜𝜁₅. (43)

For𝑡>𝑡, the same arguments can be continued since𝑥(𝑡)≥𝜁₄.

Case II. Consider𝑡^∗≠(𝑛₁+𝑙−1)𝑇,𝑛₁∈𝑁. Then𝑥(𝑡)≥𝜁₄for𝑡∈[𝑡₃,𝑡^∗]and𝑥(𝑡^∗)=𝜁₄, and suppose𝑡^∗∈((𝑛^′₁+𝑙−1)𝑇,(𝑛^′₁+𝑙)𝑇),𝑛^′₁∈𝑁.There are two possible cases for𝑡∈(𝑡^∗,(𝑛^′₁+𝑙)𝑇).

Case IIa.(𝑥(𝑡)≤𝜁₄for all𝑡∈(𝑡^∗,(𝑛^′₁+𝑙)𝑇)).Similar to Case I, we can prove that there must be a𝑡^′₂∈[(𝑛^′₁+𝑙)𝑇,(𝑛^′₁+𝑙)𝑇+𝑇]such that𝑥(𝑡^′₁)>𝜁₄.Here we omit it. Let𝑡̃=inf_{𝑡≥𝑡∗}{𝑥(𝑡)>𝜁₄}; then𝑥(𝑡)̃≤𝜁₄. For𝑡∈(𝑡^∗,𝑡)̃, we have 𝑥(𝑡)≥𝜁₄1−𝛿₁^𝑛2𝑛3exp𝑛₂+𝑛₃+1𝜂₂𝑇. (44)

Let𝜁₆≜𝜁₄(1−𝛿₁)^𝑛2𝑛3exp((𝑛₂+𝑛₃+1)𝜂₂𝑇)<𝜁₅, so𝑥(𝑡)≥𝜁₆for𝑡∈(𝑡^∗,𝑡)̃.For𝑡>𝑡,̃the same arguments can be continued since𝑥(𝑡)̃≥𝜁₄.

Case IIb. There exists a𝑡∈(𝑡^∗,(𝑛^′₁+𝑙)𝑇)such that𝑥(𝑡)>𝜁₄.Let𝑡̃=inf_{𝑡≥𝑡∗}{𝑥(𝑡)>𝜁₄}; then𝑥(𝑡)̃≤𝜁₄for𝑡∈(𝑡^∗,𝑡)̃and𝑥(𝑡)̃=𝜁₄.For𝑡∈(𝑡^∗,𝑡)̃,(41) holds true. Integrating (41) on(𝑡^∗,𝑡)̃, we have 𝑥(𝑡)≥𝑥𝑡^∗exp𝜂₂𝑡−𝑡^∗≥𝜁₄exp𝜂₂𝑇>𝜁₆. (45) Since𝑥(𝑡)≥𝜁₄for𝑡>𝑡,̃the same arguments can be continued. Hence𝑥(𝑡)≥𝜁₆for all𝑡>𝑡.̃The proved method for𝑦(𝑡)is similar to the proved method for𝑥(𝑡). SetΩ={(𝑥,𝑦,𝑥)∶𝑥≥𝜁₂,𝑦≥𝜁₈,𝑧≥𝜁₁,𝑥+𝑦+𝑧≤3𝑀}.Obviously, we know that the setΩ∈int𝑅³₊is global attractors. every solution of system (1) will eventually enter and remain in regionΩ. Therefore, system (1) is permanent. The proof is completed.

Using the theoretical analysis, we obtain the threshold expression of the release amount𝑝and the average harvest𝛿₃of the filter-feeding fish under the condition of some species extinction and all species persistence, which in turn provides a theoretical basis for the numerical simulation.

3. Numerical Analysis

The effective methods used for eutrophication removal, algal bloom control process, and associated public perception are a crucial consideration in developing effective systems of eutrophication removal management for Zeya reservoir; hence the sustainability analysis of the impulsive eutrophication controlling system must be investigated in detail. In order to study how to choose the value of the parameters𝑝and𝛿₃so that they can satisfy the actual monitoring data and enable Theorems 7 and 8 to be feasible, it is necessary to present a more in-depth analysis of the relationship among𝑝and𝛿₃. Set𝑓₁(𝑝,𝛿₃)=𝑟₁𝑇−(𝑎₂/𝑏₂𝑚)(𝑝/(1−(1−𝛿₃)exp(−𝑚𝑇)))(1−𝛿₃exp(−𝑚𝑙𝑇)−(1−𝛿₃)exp(−𝑚𝑇))−ln(1/(1−𝛿₁)), and𝑓₂(𝑝,𝛿₃)=𝑟₂𝑇−(𝑎₃/𝑏₃𝑚)(𝑝/(1−(1−𝛿₃)exp(−𝑚𝑇)))(1−𝛿₃exp(−𝑚𝑙𝑇)−(1−𝛿₃)exp(−𝑚𝑇))−ln(1/(1−𝛿₂)). On the basis of the actual monitoring data of the Zeya reservoir from year 2005 to year 2011 and laboratory controlling experiment data, the values of the parameters can be taken as follows:𝑟₁=0.75,𝑟₂=0.8,𝑘₁=20,𝑘₂=10,𝑎₁=0.3,𝑎₂=0.35,𝑎₃=0.5,𝑒₁=0.85,𝑒₂=0.5,𝑒₃=0.45,𝑏₁=2,𝑏₂=2.2,𝑏₃=2.5,𝑚=0.15,𝛿₁=0.15,𝛿₂=0.25,𝑙=0.003, and𝑇=30. Hence we have𝑓₁(𝑝,𝛿₃)=−22.33748107+(1.060606061(0.9888910035−0.9754817198𝛿₃)𝑝)/(0.9888910035+0.01110899654𝛿₃)and𝑓₂(𝑝,𝛿₃)=−23.71231793+(1.333333333(0.9888910035−0.9754817198𝛿₃)𝑝)/(0.9888910035+0.01110899654𝛿₃). It is obvious to find from Figure 1(a) that there exist some values of𝑝and𝛿₃so that𝑓₁(𝑝,𝛿₃)<0and𝑓₂(𝑝,𝛿₃)<0or𝑓₁(𝑝,𝛿₃)>0and𝑓₂(𝑝,𝛿₃)>0, which in turn confirms the existence of Theorems 7 and 8. From the actual stocking amount of the little filter-feeding fish in Zeya reservoir to control the algal bloom, we can choose biologically feasible parametric numerical value of𝑝as𝑝=35. Further, it can be seen from Figure 1(b) that if the value of𝛿₃is less than0.4, there exist𝑓₁(𝑝,𝛿₃)>0and𝑓₂(𝑝,𝛿₃)>0, but if the value of𝛿₃is greater than or equal to 0.4, then the system (1) remains locally asymptotically stable. It is to be noted from Figure 2 that the time series solution of the filter-feeding fish population oscillates with stable limit cycle with increasing time. In contrast, both the algae population and the dominant zooplankton population rapidly decrease to their extinction with the increasing time. Meanwhile, these results further demonstrate that Theorems 7 and 8 are correct and feasible.

(a)

(b)

Figure 1 

(a) The relationship diagram of the function𝑓1(𝑝,𝛿3)and the function𝑓2(𝑝,𝛿3)with respect to the parameters𝑝and𝛿3; (b) the relationship diagram of the function𝑓1(𝑝,𝛿3)and the function𝑓2(𝑝,𝛿3)with respect to the parameter𝛿3with𝑝=35.

(a)

(b)

(c)

Figure 2 

The dynamics of the system (1).

Considering the great uncertainty that still surrounds biological treatment mechanism, it is most prudent to interpret the simulation output on the basis of information and data obtained from observation and experiment. Simulation design offers tremendous flexibility to study the behavior dynamic of biological treatment process, so it is necessary to investigate that the release amount𝑝is how to affect the dynamics of biological treatment process with the parameter𝛿₃=0.4and0<𝑝≤35. It is worth pointing out from Figure 3 that the system (1) with0<𝑝<35exhibits relatively complex dynamical behaviors: chaos, period-halving bifurcation, and periodic attractor, but the system (1) with𝑝=35remains locally asymptotically stable. Furthermore, it should be stressed that the dominant zooplankton population rapidly decreases to their extinction when the value of𝑝is greater than 4.4268; this phenomenon shows that the dominant zooplankton population will advance toward extinction with the increase of𝑝. But it is regretful for this phenomenon because the dominant zooplankton population is one of the main food for the filter-feeding fish and can graze the algae population, and the idea of using biological eutrophication removal treatment to control the algal bloom not only can control the algal bloom, but also consider the sustainable development of eutrophication removal process with the maximizing economic benefits. Thus, controlling the desirable biomass of the dominant zooplankton population is a primary aim of operation and maintenance, and how to well carry out the eutrophication removal process is an extremely vital problem of Zeya reservoir.

(a)

(b)

(c)

Figure 3 

Bifurcation diagram of the system (1) with0<𝑝≤35.

In real ecosystem, the effect of some critical factors can be depicted by the properties and the dynamic behaviors of species [27, 28]. To more carefully investigate how to legitimately implement eutrophication removal treatment, Figure 4 depicts the long-run dynamics of system (1) for0≤𝛿₃≤1with𝑝=5. It is easy to find from Figure 4 that the system exhibits relatively complex dynamical behaviors: chaos, period-halving bifurcation, and periodic attractor, but it is worth pointing out that, when the value of𝛿₃is beyond 0.6, the algae population will rapidly go to extinction and the dominant zooplankton population will go toward a reasonable state, and the biomass of the dominant zooplankton population is so small that it does not cause the aggregation of the dominant zooplankton. Meanwhile, the biomass of the filter-feeding fish is feasible to maintain the sustainable development of eutrophication removal process with the maximizing economic benefits because of𝑝=5. This result suggests that we can choose the value of𝑝and𝛿₃as 5 and 0.6 for legitimately implementing eutrophication removal treatment. For further analysis of the feasibility of this result, Figure 5 gives the dynamics of the system (1) with the increase of𝑡as well as𝑝=5and𝛿₃=0.6. It is gratifying to see from Figure 5 that the algae population quickly goes to the extinct status with the increase of𝑡, and the dominant zooplankton population is close to a reasonable number with the increase of𝑡. It is worth encouraging that the filter-feeding fish can arrive to a viable amount so that the eutrophication removal process can be developed with the maximizing economic benefits. Therefore, it should be stressed that these results can perfectly confirm that the impulsive eutrophication controlling system can be well applied in the process of controlling the algal bloom in the theoretical work of eutrophication removal and mathematical model.

(a)

(b)

(c)

Figure 4 

Bifurcation diagram of the system (1) with0≤𝛿3≤1.

(a)

(b)

(c)

Figure 5 

The dynamics of the system (1) with𝑝=5and𝛿3=0.6.

Based on the foregoing analysis, it is successful for impulsive eutrophication controlling system to implement the eutrophication removal process. The impulsive eutrophication controlling system not only can promote all wanted species persistence, but also is convenient for the manager to maintain the normal development of biological treatment process. Moreover, simulation analysis provides an approximation of the real biological controlling system behaviors; hence it is very easy to obtain the best value of𝑝and𝛿₃, which is in accordance with observations and in turn shows that the selection of the best value about some critical parameters in the view of population dynamics is feasible.

4. Conclusions and Remarks

The theoretical ecologists and applied mathematicians are now working on controlling the algal bloom in ecosystem, but the optimal way to eradicate the algae population to control the algal bloom is yet to know, but in the context of population growth dynamics, how to implement impulsive control strategy to prevent and control the algal bloom with the maximizing economic benefits is worthy of studying. In the paper, we take an approach to study how to better implement impulsive eutrophication controlling system to control the algal bloom and maintain the normal development of biological treatment process with maximizing economic benefits by using mathematical analysis and simulating the dynamics. On the basis of eutrophication ecology and differential equation, an impulsive eutrophication controlling system is studied analytically and numerically. A key advantage of the impulsive eutrophication controlling system is that it can be quite accurate to describe the interaction effect of some critical factors (fishermen catch and releasing small fry, etc.), which enables a systematic and logical procedure for fitting eutrophication mathematical system to real monitoring data and experiment data. mathematical theoretical works have been pursuing the investigation of two threshold functions𝑓₁(𝑝,𝛿₃)and𝑓₂(𝑝,𝛿₃)of two critical parameters𝑝and𝛿₃, if𝑓₁(𝑝,𝛿₃)<0and𝑓₂(𝑝,𝛿₃)<0, the system (1) remain locally asymptotically stable, and if𝑓₁(𝑝,𝛿₃)>0and𝑓₂(𝑝,𝛿₃)>0, the system (1) can ensure all species persistence. These results in turn provide a theoretical basis for the numerical simulation.

From the viewpoint of population dynamics, numerical simulation works have been pursuing the investigation of confirming the existence and feasibility of Theorems 7 and 8 and how to choose the best value of two critical parameters𝑝and𝛿₃. Within this framework, it can be assumed that the sustainability of the impulsive eutrophication controlling system mainly depends on two critical parameters𝑝and𝛿₃. Then the direct and indirect effects caused by critical parameters𝑝and𝛿₃are investigated by means of bifurcation analysis. These results indicate that if𝑝=5and𝛿₃=0.6, the algae population will go to extinction, but the dominant zooplankton and the filter-feeding fish can survive in a feasible range, which can assure that the eutrophication removal process can be circulatingly developed with maximizing economic benefits. In addition, it is worthwhile to remark that mathematical theoretical results and numerical simulation results are complementary. It is hopeful that these results have an important implication for eutrophication removal management; hence in the follow-up works, we will use the system to character the growth rule and the correlation between the algae population, the dominant zooplankton, and the filter-feeding fish.

Acknowledgments

This work was supported by the Zhejiang Provincial Natural Science Foundation of China (Grant no. LY13A010010), by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001), by the National Natural Science Foundation of China (Grant no. 31370381), by the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510) and by the Nonprofit Technology Research Projects of Zhejiang Province (Grant no. 2011C23119).