# The Effect of Continuous and Pulse Input Nutrient on a Lake Model

Authors: Yongfeng Li; Dongliang Xie; Jing-An Cui
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

In an aquatic environment, a mathematical model consisting of nutrient, phytoplankton, and zooplankton has been considered, where excretion of zooplankton is considered as one of the sources of nutrient. We investigate the effect of the input rate of the limiting nutrient from outside on controlling algal bloom. First, we consider input limiting nutrient continuously, obtain the system has two boundary equilibria, and analyze the existence of the positive equilibrium by means of stability analysis, we get conditions for the stability of the equilibria. Then, we consider input limiting nutrient impulsively. We get the exact expression of the boundary periodic solution and obtain the condition for the stability of the periodic boundary solution. We also consider the effect of temperature on the system and give a model of Taihu Lake as an example. Finally, we give numerical simulation of our results and explain the effect of input limiting nutrient on controlling bloom of the lake system.


---

## Body

## 1. Introduction

In recent years, many mathematical models of ecosystems have been proposed and used for understanding complex aspects of marine ecosystems, especially eutrophication problem. Eutrophication is a serious “disease" of lakes around the world. It has badly damaged lake ecosystem health and has resulted in an imbalance between biological components, decreasing of biodiversity, and resilience. The adverse effects on human health, commercial fisheries, tourism, and environment are well established. A very important step towards protection and restoration of a lake is to limit, divert, or treat excessive nutrient, organic, and silt loads [1].Many models including nutrient concentration has been studied by [2–6]. Hallam [7] studied stability and persistence properties of a family of nutrient-controlled plankton models and obtained necessary and sufficient conditions for persistence. Gard [8] studied a nutrient-phytoplankton-zooplankton (NPZ) model with generalized functional response and obtained sharper criteria for persistence. A phytoplankton-zooplankton model was studied by Steffen et al. [9]; local and global behavior of the model were obtained. Busenberg et al. [10] demonstrated coexistence of plankton population in an orbitally stable oscillatory mode for a nutrient-plankton model.In this paper, we consider the model consisting of three ordinary differential equations for three dynamical variables, that is, the densities of limiting nutrient (N), phytoplankton (P), and zooplankton (Z). A NPZ model has been studied by [11, 12, 12–18] that modified the model of [11, 19] and investigated zooplankton mortality and dynamical behavior of the model. A top predator invasion into the NPZ model was studied by [15], the paper considered possible consequences of biological invasion in an epipelagic ecosystem. References [16, 18] used a general function to describe the nutrient uptake rates of phytoplankton and zooplankton and herbivore grazing.The model that we examined is based on NPZ model of [12–18]. The model differs from [12, 15] as a closed system is replaced by an open system; that is, we consider external inflow nutrient except for internal factors. In [16], nutrients absorbed by phytoplankton and zooplankton have been completely digested and absorbed and conversed phytoplankton growth; zooplankton-nutrient conversion and phytoplankton-nutrient conversion are total conversion. In [16, 18], the excretion of phytoplankton and zooplankton is nutrient. In fact, zooplankton-nutrient conversion and phytoplankton-nutrient conversion are part conversion, and part of phytoplankton and zooplankton excretion is nutrient. So, we consider modified phytoplankton nutrient conversion and zooplankton excretion.Water bloom takes place frequently in China, such as Taihu Lake, Chaohu Lake, and Dianchi Lake. The main reason for this phenomenon is that many kinds of pollutant have been discharged into lakes, therefore, it is very important to investigate nutrition input from outside. In the following, we take Taihu Lake as an example to explain the main reason of eutrophication.Through investigating the input and output balance relationship of a lake nitrogen (N), phosphorus (P) and chemical oxygen demand (COD), we show that sediment is the main source of nutrient salt. This result provide information for studying lake eutrophication. Table1 is statistical data of discharge amount for Taihu area in 2004 [20, page 25], where X, Y, and H
1 are industry pollution discharge, H
2 is life pollution discharge, H
3 is agricultural pollution discharge, H
4 is water and soil drained away, H
5 is poultry cultivation pollution, H
6 is aquatic cultivation pollution, H
7 is lake tourism pollution, H
8 is rain dropping, H
9 is dust falling, and H
10 is ships pollution.Table 1
System parameters and ranges of their values used in the model.
Parameters
Units
Meaning of the parameters
Literature values
λ
mg L−1 day−1
Input rate of limiting nutrient
0–1.5 [21]
d
day−1
Loss rate of limiting nutrient
0.1 [22]
r
day−1
Maximum growth rate of phytoplankton
0.3 [21]
μ
day−1
Maximum growth rate of zooplankton
0.6–1.4 [12]
k
1
mg L−1
Half saturation constant for nutrient uptake
0.02–0.15 [12, 19]
k
2
mg L−1
Half saturation constant for zooplankton grazing
0.02–0.1 [12, 19]
m
1
day−1
Death rate of phytoplankton
0.15 [22]
m
2
day−1
Death rate of zooplankton
0.1 [22]
d
1
day−1
Removal rate of phytoplankton from the system
0.1 [22]
d
2
day−1
Removal rate of zooplankton from the system
0.1 [22]
w
—
Phytoplankton efficiency coefficient
1 [12]
η
—
Zooplankton efficiency coefficient
0.2–0.5 [12]
α
—
Regeneration of nutrient from decomposition of phytoplankton
1 [12]
β
—
Regeneration of nutrient from decomposition of zooplankton
0.5–0.9 [12, 19]
γ
—
Zooplankton excretion coefficient
0.5–0.9 [12]From Table1, we can know that the main pollution source of Taihu Lake is the pollution from life and agriculture and poultry cultivation, so we can control water bloom by control these pollution inputs. So, in this paper, we consider the effect of nutrient coming from outside on the dynamics of the model and want to control bloom mainly by controlling the amount of outside nutrient. First, we consider input limiting nutrient continuously. Then we consider input limiting nutrient impulsively. System with impulsive effects describing evolution processes is characterized by the fact that at certain moments of time they abruptly experience a change of state. Processes of such character are studied in almost every domain of applied sciences. Theories of impulsive differential equations are found in the books [23, 24]. In recent years, their applications can be found in many domains of applied sciences [21, 22, 25–27].The paper is arranged like this. In Section2, we give a continuous input nutrient model, obtain the boundary equilibrium of the model, and analyze its stability; we also investigate the existence of the positive equilibrium and discuss its stability. In Section 3, pulse input nutrient is considered, nutrient flow into the lake every τ period. We obtain the exact boundary periodic solution of the impulsive input nutrient system. Using Floquet theory for the impulsive differential equation and small-amplitude perturbation skills and techniques of comparison, we prove that the boundary periodic solution is locally asymptotically stable if (
(
k
1
+
(
λ
-
k
1
)
exp
⁡
(
-
d
τ
)
)
/
(
k
1
+
(
λ
-
k
1
)
exp
⁡
(
-
d
τ
)
)
)
>
exp
⁡
(
-
d
τ
(
m
1
+
d
1
)
/
w
r
). But, in fact, the growth rate of algal is affected by temperature, so we consider the effect of temperature on the system in Section 4 and give an example of Taihu Lake. Finally, we give a brief discussion of our results and further numerical simulation; we also explain the effect of input limiting nutrient on controlling bloom of the lake system.


## 2. Continuous Input Nutrient

In a real open marine ecosystem, we consider the following continuous NPZ model:(1)
N
˙
(
t
)
=
λ
-
d
N
(
t
)
-
r
N
(
t
)
P
(
t
)
k
1
+
N
(
t
)
+
α
m
1
P
(
t
)
N
˙
(
t
)
=
+
β
m
2
Z
(
t
)
+
γ
μ
P
(
t
)
Z
(
t
)
k
2
+
P
(
t
)
,
P
˙
(
t
)
=
w
r
N
(
t
)
P
(
t
)
k
1
+
N
(
t
)
-
m
1
P
(
t
)
-
d
1
P
(
t
)
-
μ
P
(
t
)
Z
(
t
)
k
2
+
P
(
t
)
,
Z
˙
(
t
)
=
η
μ
P
(
t
)
Z
(
t
)
k
2
+
P
(
t
)
-
m
2
Z
(
t
)
-
d
2
Z
(
t
)
,
where N
(
t
) is the amount of limiting nutrient at time t, P
(
t
) is the biomass of phytoplankton at time t, and Z
(
t
) is the biomass of phytoplankton at time t. The parameter λ represents the input rate of the limiting nutrient from the environment, d is loss rate of the nutrient, r
(
μ
) is maximum phytoplankton (zooplankton) growth rate, m
1
(
m
2
) is phytoplankton (zooplankton) natural mortality and respiration rate, d
1
(
d
2
) is removal rate of phytoplankton (zooplankton), k
1 is half saturation constant for nutrient uptake, k
2 is half saturation constant for zooplankton grazing, w
(
η
) is phytoplankton (zooplankton) efficiency coefficient, γ is zooplankton excretion coefficient, and α
(
β
) is regeneration of nutrient from decomposition of phytoplankton (zooplankton).Here the uptake kinetics of the nutrient and consumption of phytoplankton by zooplankton are described by a Michaelis-Menten-Monod function which is nonnegative, increasing, and equal to zero in the absence of nutrient (phytoplankton); see [13]. Obviously, we have w
≤
1, α
≤
1, β
≤
1, and η
+
γ
≤
1. Natural mortality and respiration rate of the plankton are assumed by other authors to be linear losses.First, we show that solutions of system (1) always exist for all t
>
0. Since P
˙
|
P
=
0
=
0, Z
˙
|
Z
=
0
=
0, and N
˙
|
N
=
0
≥
λ
>
0, we have N
(
t
)
≥
0, P
(
t
)
≥
0, and Z
(
t
)
≥
0 for all t
>
0. Then the solutions of system (1) are nonnegative for t
>
0.We use realistic values of model parameters obtained from different sources and summarized in Table2. In the following, we investigate the existence and stability of equilibria of the system (1).Table 2
Pollutant discharge amount of Taihu Lake area in 1998.
Source of pollution
CODcr
TN
TP
X
/
(
t
/
a
)
Y
/
%
X
/
(
t
/
a
)
Y
/
%
X
/
(
t
/
a
)
Y
/
%
H
1
160 478.2
21.51
4041.8
3.71
264.74
1.70
H
2
197 470.2
26.47
32 911.7
30.22
7477.10
47.89
H
3
272 805.7
36.59
53 970.0
49.54
3990.41
25.56
H
4
8585.0
1.15
800.0
0.73
192.00
1.22
H
5
37 850.0
5.08
6710.0
6.16
2820.00
18.06
H
6
40 791.9
5.47
5703.9
5.24
634.92
4.07
H
7
150.2
0.02
81.1
0.07
5.60
0.04
H
8
23 595.0
3.16
2759.5
2.53
60.10
0.39
H
9
—
—
420.9
0.39
33.00
0.21
H
10
3883.7
0.05
1538.2
1.41
132.02
0.86
Total
745 609.9
100.0
108 937.1
100.0
15 609.89
100.0


### 2.1. Existence of Equilibria of the System (1)

It is easy to calculate the system (1) that always has a boundary equilibrium denoted by E
0
=
(
N
0
,
0,0
)
=
(
λ
/
d
,
0,0
). It also has a boundary equilibrium E
^
=
(
N
^
,
P
^
,
0
) if w
r
-
m
1
-
d
1
>
0 and λ
>
(
(
d
k
1
(
m
1
+
d
1
)
)
/
(
w
r
-
m
1
-
d
1
)
)
=
Δ
λ
0 hold, where
(2)
N
^
=
k
1
(
m
1
+
d
1
)
w
r
-
m
1
-
d
1
,
P
^
=
λ
-
d
N
^
d
1
=
w
(
λ
w
r
-
(
λ
+
d
k
1
)
(
m
1
+
d
1
)
)
(
w
r
-
m
1
-
d
1
)
(
m
1
(
1
-
w
α
)
+
d
1
)
.In the following, we consider the positive equilibrium denoted byE
*
=
(
N
*
,
P
*
,
Z
*
); then N
*, P
*, Z
* must satisfy
(3)
λ
-
d
N
-
r
N
P
k
1
+
N
+
α
m
1
P
+
β
m
2
Z
+
γ
μ
P
Z
k
2
+
P
=
0
,
(4)
w
r
N
k
1
+
N
-
m
1
-
d
1
-
μ
Z
k
2
+
P
=
0
,
(5)
η
μ
P
k
2
+
P
-
m
2
-
d
2
=
0
.
From (5) we have P
*
=
k
2
(
m
2
+
d
2
)
/
(
η
μ
-
m
2
-
d
2
), so P
* is positive if η
μ
-
m
2
-
d
2
>
0. In the following, we consider the existence of N
* and Z
* under the condition η
μ
-
m
2
-
d
2
>
0. Substituting P
* into (3) and (5), we get
(6)
λ
-
d
N
-
r
k
2
(
m
2
+
d
2
)
N
(
η
μ
-
m
2
-
d
2
)
(
k
1
+
N
)
+
α
m
1
k
2
(
m
2
+
d
2
)
η
μ
-
m
2
-
d
2
+
β
m
2
Z
+
γ
(
m
2
+
d
2
)
Z
η
=
0
,
w
r
N
k
1
+
N
-
m
1
-
d
1
-
(
η
μ
-
m
2
-
d
2
)
Z
k
2
η
=
0
.
That is,
(7)
-
d
N
2
+
a
1
N
+
k
1
a
2
Z
+
a
2
N
Z
+
a
3
=
0
,
(8)
Z
=
c
1
-
c
2
k
1
+
N
,
where a
1 = λ
-
d
k
1
+
(
k
2
(
m
2
+
d
2
)
(
α
m
1
-
r
)
/
(
η
μ
-
m
2
-
d
2
)
), a
2
=
(
β
m
2
+
(
γ
(
m
2
+
d
2
)
/
η
)
)
>
0, a
3
=
(
λ
k
1
+
(
α
m
1
k
1
k
2
(
m
2
+
d
2
)
/
(
η
μ
-
m
2
-
d
2
)
)
)
>
0, c
1
=
k
2
η
(
w
r
-
m
1
-
d
1
)
/
(
η
μ
-
m
2
-
d
2
), and c
2
=
w
r
k
1
k
2
η
/
(
η
μ
-
m
2
-
d
2
)
>
0. Equation (7) can also be rewritten as
(9)
Z
=
b
1
(
k
1
+
N
)
-
b
2
k
1
+
N
-
g
,
with
(10)
b
1
=
d
a
2
=
d
β
m
2
+
(
γ
(
m
2
+
d
2
)
/
η
)
>
0
,
b
2
=
a
3
-
k
1
(
d
k
1
+
a
1
)
a
2
b
2
=
k
1
k
2
(
m
2
+
d
2
)
r
c
c
c
c
c
×
(
(
η
μ
-
m
2
-
d
2
)
(
β
m
2
+
γ
(
m
2
+
d
2
)
η
)
)
-
1
>
0
,
g
=
2
d
k
1
+
a
1
a
2
g
=
(
(
λ
+
d
k
1
)
(
η
μ
-
m
2
-
d
2
)
c
c
c
c
c
c
c
+
k
2
(
m
2
+
d
2
)
(
α
m
1
-
r
)
)
c
c
c
c
×
(
(
η
μ
-
m
2
-
d
2
)
(
β
m
2
+
γ
(
m
2
+
d
2
)
η
)
)
-
1
.From (8), we have Z
′
|
N
=
c
2
/
(
k
1
+
N
)
2
>
0, so, Z is a monotonous increasing function about N, and the asymptotes are N
=
-
k
1 and Z
=
c
1. In the same way, for (9), Z
′
|
N
=
b
1
+
(
b
2
/
(
k
1
+
N
)
2
)
>
0, Z is a monotonous increasing function about N, and the asymptotes are N
=
-
k
1 and Z
=
b
1
(
k
1
+
N
)
-
g. Neglecting the left branches of figures of (8) and (9), the Figures of (8) and (9) are shown in Figures 1(a) and 1(b). If we want the curve of (8) and the curve of (9) to intersect in the first quadrant, we must have c
1
>
0; that is, w
r
-
m
1
-
d
1
>
0. If Z
=
0 in (8) and (9), we can calculate A
=
Δ
(
N
1
,
0
)
=
(
(
c
2
/
c
1
)
-
k
1
,
0
) and C
=
Δ
(
N
2
,
0
)
=
(
(
1
/
2
b
1
)
(
g
-
2
b
1
k
1
+
g
2
+
4
b
1
b
2
)
,
0
). While N
1
=
(
c
2
/
c
1
)
-
k
1
=
(
m
1
+
d
1
)
k
1
/
(
w
r
-
m
1
-
d
1
)
>
0 if w
r
-
m
1
-
d
1
>
0; that is, N
1
=
N
^.Figures of (8) and (9), where (a) is the figure of (9) and (b) is the figure of (8).
(a)
(b)IfN
1
<
N
2, then the figures of (8) and (9) only have a positive intersection point; that is,
(11)
k
2
(
m
2
+
d
2
)
(
w
r
-
m
1
-
d
1
)
(
m
1
(
1
-
w
α
)
+
d
1
)
+
w
(
m
2
+
d
2
-
η
μ
)
(
w
r
λ
-
(
λ
+
d
k
1
)
(
m
1
+
d
1
)
)
<
0
.
That is, λ
>
(
k
2
(
m
2
+
d
2
)
(
w
r
-
m
1
-
d
1
)
(
m
1
(
1
-
w
α
)
+
d
1
)
+
w
(
η
μ
-
m
2
-
d
2
)
d
k
1
(
m
1
+
d
1
)
) / w
(
η
μ
-
m
2
-
d
2
)
(
w
r
-
m
1
-
d
1
)
=
Δ
λ
*. Obviously, λ
*
>
λ
0. If N
1
=
N
2, then (8) and (9) only intersect on the positive N axis; that is, Z
=
0 and λ
=
λ
*; then P
*
=
P
^. In this condition, the equilibrium only has E
^
=
(
N
^
,
P
^
,
0
).According to the above statement, we know that the system (1) only has a positive equilibrium if η
μ
-
m
2
-
d
2
>
0, w
r
-
m
1
-
d
1
>
0, and λ
>
λ
*.The existence of equilibria is summarized as follows:(i)
E
0
=
(
N
0
,
0,0
) always exists;
(ii)
E
^
=
(
N
^
,
P
^
,
0
) exists if w
r
-
m
1
-
d
1
>
0 and λ
>
λ
0 hold;
(iii)
E
*
=
(
N
*
,
P
*
,
Z
*
) exists if w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
>
λ
*.


### 2.2. Stability of Equilibria of the System (1)

First, we know that all solutions of (1) are uniformly ultimately bounded, that is,Theorem 1.
For each positive solutionx
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) of system (1) there exists a constant M
>
0, such that N
(
t
)
≤
M, P
(
t
)
≤
M, and Z
(
t
)
≤
M with t large enough.Proof.
Define a functionL
(
t
,
x
) such that L
(
t
,
x
)
=
N
(
t
)
+
P
(
t
)
+
Z
(
t
). Then L
(
t
,
x
)
∈
V
0 and the upper right derivative of L
(
t
,
x
) along a solution of system (1) is described as
(12)
D
+
L
(
t
,
x
)
|
(
1
)
=
λ
-
d
N
+
(
α
-
1
)
m
1
P
+
(
β
-
1
)
m
2
Z
-
d
1
P
-
d
2
Z
+
(
w
-
1
)
r
N
P
k
1
+
N
+
(
γ
+
η
-
1
)
μ
P
Z
k
2
+
P
≤
λ
-
d
N
-
(
d
1
+
(
1
-
α
)
m
1
)
P
-
(
d
2
+
(
1
-
β
)
m
2
)
Z
≤
λ
-
d
^
L
(
t
,
x
)
,
where d
^
=
min
⁡
{
d
,
d
1
+
(
1
-
α
)
m
1
,
d
2
+
(
1
-
β
)
m
2
}. Hence D
+
L
(
t
,
x
)
|
(
1
)
≤
-
d
^
L
(
t
,
x
)
+
λ. So, we have
(13)
L
(
t
,
x
)
≤
λ
d
^
+
(
L
(
0
)
-
λ
d
^
)
e
-
d
^
t
⟶
λ
d
^
,
t
⟶
∞
.
Therefore, L
(
t
,
x
) is ultimately bounded by a constant and there exists a constant M
>
0, such that N
(
t
)
≤
M, P
(
t
)
≤
M, and Z
(
t
)
≤
M, for each solution x
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) of system (1) with t large enough. The proof is complete.In the following, we letE
0
=
(
N
0
,
P
0
,
Z
0
) be any equilibrium of system (1); then the Jacobian matrix about E
0 is given by(
*
)
J
|
E
0
=
(
-
d
-
r
k
1
P
0
(
k
1
+
N
0
)
2
-
r
N
0
k
1
+
N
0
+
α
m
1
+
γ
μ
k
2
Z
0
(
k
2
+
P
0
)
2
β
m
2
+
γ
μ
P
0
k
2
+
P
0
w
r
k
1
P
0
(
k
1
+
N
0
)
2
w
r
N
0
k
1
+
N
0
-
m
1
-
d
1
-
μ
k
2
Z
0
(
k
2
+
P
0
)
2
-
μ
P
0
k
2
+
P
0
0
η
μ
k
2
Z
0
(
k
2
+
P
0
)
2
η
μ
P
0
k
2
+
P
0
-
m
2
-
d
2
)
,denoted by
(14)
J
|
E
0
=
(
a
11
a
12
a
13
a
21
a
22
a
23
0
a
32
a
33
)
;
then the characteristic equation about E
0 is: λ
3
+
p
1
λ
2
+
p
2
λ
+
p
3
=
0, where
(15)
p
1
=
-
(
a
11
+
a
22
+
a
33
)
,
p
2
=
a
11
a
22
+
a
11
a
33
+
a
22
a
33
-
a
21
a
12
-
a
23
a
32
,
p
3
=
a
12
a
21
a
33
-
a
11
a
22
a
33
-
a
21
a
32
a
13
+
a
11
a
23
a
32
.About the stability of the equilibriumE
0, E
^, and E
*, we have the following results.Theorem 2.
The steady stateE
0 of system (1) is locally asymptotically stable if m
1
+
d
1
>
(
w
r
N
0
/
(
k
1
+
N
0
)
); that is, w
r
-
m
1
-
d
1
>
0 and λ
<
λ
0.Proof.
For equilibriumE
0
=
(
N
0
,
0,0
)
=
(
λ
/
d
,
0,0
), the Jacobian matrix evaluated at E
0 is
(16)
J
|
E
0
=
(
-
d
-
r
N
0
k
1
+
N
0
+
α
m
1
β
m
2
0
w
r
N
0
k
1
+
N
0
-
m
1
-
d
1
0
0
0
-
(
m
2
+
d
2
)
)
;
obviously, two of eigenvalues are λ
1
=
-
d
<
0 and λ
2
=
-
(
m
2
+
d
2
)
<
0; the other eigenvalue is λ
3
=
(
w
r
N
0
/
(
k
1
+
N
0
)
)
-
m
1
-
d
1. The steady state E
0 is asymptotically stable if and only if the eigenvalues of the Jacobian matrix at E
0 have negative real parts; the proof is complete.Theorem2 implies that if the nutrient-phytoplankton conversion rate is less than phytoplankton loss rate, then both phytoplankton and zooplankton population will become extinct.Theorem 3.
The steady stateE
^ of system (1) is locally asymptotically stable if m
2
+
d
2
-
η
A
2
>
0; that is, w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
<
λ
*.Proof.
For equilibriumE
^
=
(
N
^
,
P
^
,
0
), from the second equation of the system (1), we have (
w
r
N
^
/
(
k
1
+
N
^
)
)
-
m
1
-
d
1
=
0; that is, (
r
N
^
/
(
k
1
+
N
^
)
)
=
(
(
m
1
+
d
1
)
/
w
), so -
(
r
N
^
/
(
k
1
+
N
^
)
)
+
α
m
1
+
(
γ
μ
k
2
Z
^
/
(
k
2
+
P
^
)
2
)
=
-
(
(
m
1
+
d
1
)
/
w
)
+
α
m
1
+
(
γ
μ
k
2
Z
^
/
(
k
2
+
P
^
)
2
)
=
-
(
(
1
-
w
α
)
m
1
+
d
1
)
/
w
<
0. Hence (
*
) reduces to(17)
J
|
E
^
=
(
-
d
-
r
k
1
P
^
(
k
1
+
N
^
)
2
-
(
1
-
w
α
)
m
1
+
d
1
w
β
m
2
+
γ
μ
P
^
k
2
+
P
^
w
r
k
1
P
^
(
k
1
+
N
^
)
2
0
-
μ
P
^
k
2
+
P
^
0
0
η
μ
P
^
k
2
+
P
^
-
m
2
-
d
2
)
≜
(
-
d
-
A
1
-
d
~
β
m
2
+
γ
A
2
w
A
1
0
-
A
2
0
0
η
A
2
-
m
2
-
d
2
)
,with A
1
=
r
k
1
P
^
/
(
k
1
+
N
^
)
2
>
0, A
2
=
μ
P
^
/
(
k
2
+
P
^
)
>
0, and d
~
=
(
(
1
-
w
α
)
m
1
+
d
1
)
/
w
>
0. Then λ
3
+
p
1
λ
2
+
p
2
λ
+
p
3
=
0, where
(18)
p
1
=
(
d
+
A
1
)
+
(
m
2
+
d
2
-
η
A
2
)
,
p
2
=
(
d
+
A
1
)
(
m
2
+
d
2
-
η
A
2
)
+
d
~
w
A
1
,
p
3
=
d
~
w
A
1
(
m
2
+
d
2
-
η
A
2
)
;
then p
1
p
2
-
p
3
=
(
d
+
A
1
)
2
(
m
2
+
d
2
-
η
A
2
)
+
d
1
w
A
1
(
d
+
A
1
)
+
(
d
+
A
1
)
(
m
2
+
d
2
-
η
A
2
)
2. Obviously, p
i
>
0 (i
=
1,2
,
3) and p
1
p
2
-
p
3
>
0 if m
2
+
d
2
-
η
A
2
>
0. According to Routh-Hurwitz criteria, E
^ is locally stable if m
2
+
d
2
-
η
A
2
>
0. The proof is complete.Theorem3 implies that if the growth rate of zooplankton is less than its loss rate, then the zooplankton population will die out and the phytoplankton population will survive on the nutrient.Theorem 4.
The steady stateE
* of system (1) is locally asymptotically stable if h
1
>
0 and h
1
h
2
-
h
3
>
0.Proof.
For equilibriumE
*
=
(
N
*
,
P
*
,
Z
*
), (
*
) reduces to(19)
J
|
E
*
=
(
-
d
-
r
k
1
P
*
(
k
1
+
N
*
)
2
-
r
N
*
k
1
+
N
*
+
α
m
1
+
γ
μ
k
2
Z
*
(
k
2
+
P
*
)
2
β
m
2
+
γ
μ
P
*
k
2
+
P
*
w
r
k
1
P
*
(
k
1
+
N
*
)
2
w
r
N
*
k
1
+
N
*
-
m
1
-
d
1
-
μ
k
2
Z
*
(
k
2
+
P
*
)
2
-
μ
P
*
k
2
+
P
*
0
η
μ
k
2
Z
*
(
k
2
+
P
*
)
2
0
)
≜
(
-
d
-
B
1
-
B
2
+
α
m
1
+
γ
B
3
β
m
2
+
γ
B
4
w
B
1
w
B
2
-
m
1
-
d
1
-
B
3
-
B
4
0
η
B
3
0
)
,with B
1
=
r
k
1
P
*
/
(
k
1
+
N
*
)
2, B
2
=
r
N
*
/
(
k
1
+
N
*
), B
3
=
μ
k
2
Z
*
/
(
k
2
+
P
*
)
2, and B
4
=
μ
P
*
/
(
k
2
+
P
*
). Then λ
3
+
h
1
λ
2
+
h
2
λ
+
h
3
=
0, where
(20)
h
1
=
(
d
+
B
1
)
+
(
m
1
+
d
1
+
B
3
-
w
B
2
)
,
h
2
=
(
d
+
B
1
)
(
m
1
+
d
1
+
B
3
-
w
B
2
)
ccccc
+
w
B
1
(
B
2
-
α
m
1
-
γ
B
3
)
+
B
4
η
B
3
,
h
3
=
-
w
B
1
η
B
3
(
β
m
2
+
γ
B
4
)
+
B
4
(
d
+
B
1
)
η
B
3
;
then
(21)
h
1
h
2
-
h
3
=
(
d
+
B
1
)
2
(
m
1
+
d
1
+
B
3
-
w
B
2
)
+
(
d
+
B
1
)
w
B
1
(
B
2
-
α
m
1
-
γ
B
3
)
+
(
d
+
B
1
)
(
m
1
+
d
1
+
B
3
-
w
B
2
)
2
+
w
B
1
(
B
2
-
α
m
1
-
γ
B
3
)
·
(
m
1
+
d
1
+
B
3
-
w
B
2
)
+
B
4
η
B
3
(
m
1
+
d
1
+
B
3
-
w
B
2
)
+
w
B
1
η
B
3
(
β
m
2
+
γ
B
4
)
.
From (5), we have w
B
2
-
m
1
-
d
1
-
(
μ
Z
*
/
(
k
2
+
P
*
)
)
=
0, so w
B
2
-
m
1
-
d
1
-
B
3
=
(
μ
Z
*
/
(
k
2
+
P
*
)
)
-
B
3
=
(
μ
Z
*
/
(
k
2
+
P
*
)
)
-
(
μ
k
2
Z
*
/
(
k
2
+
P
*
)
2
)
=
μ
P
*
Z
*
/
(
k
2
+
P
*
)
2
>
0. According to Routh-Hurwitz criteria, we have Theorem 4.
The existence and local stability of equilibria are summarized as follows:(i)
E
0
=
(
N
0
,
0,0
) always exists; E
0 is local stable if w
r
-
m
1
-
d
1
>
0 and λ
<
λ
0; E
0 is unstable if w
r
-
m
1
-
d
1
>
0 and λ
>
λ
0;
(ii)
E
^
=
(
N
^
,
P
^
,
0
) exists if w
r
-
m
1
-
d
1
>
0 and λ
>
λ
0 hold; E
^ is local stable if w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
0
<
λ
<
λ
*; E
^ is unstable if w
r
-
m
1
-
d
1
>
0
,
η
μ
-
m
2
-
d
2
>
0 and λ
>
λ
*;
(iii)
E
*
=
(
N
*
,
P
*
,
Z
*
) exists if w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
>
λ
*; E
* is local stable if h
1
>
0 and h
1
h
2
-
h
3
>
0.
From analysis, we know that the number of equilibrium increases with the increase ofλ: from one equilibrium to three equilibria. The stability of equilibria changes with changing λ; when equilibrium E
* becomes unstable, algal bloom will break out. So, controlling the amount of input nutrient is very important in controlling water bloom; for more detail we can see the Discussion section.


## 2.1. Existence of Equilibria of the System (1)

It is easy to calculate the system (1) that always has a boundary equilibrium denoted by E
0
=
(
N
0
,
0,0
)
=
(
λ
/
d
,
0,0
). It also has a boundary equilibrium E
^
=
(
N
^
,
P
^
,
0
) if w
r
-
m
1
-
d
1
>
0 and λ
>
(
(
d
k
1
(
m
1
+
d
1
)
)
/
(
w
r
-
m
1
-
d
1
)
)
=
Δ
λ
0 hold, where
(2)
N
^
=
k
1
(
m
1
+
d
1
)
w
r
-
m
1
-
d
1
,
P
^
=
λ
-
d
N
^
d
1
=
w
(
λ
w
r
-
(
λ
+
d
k
1
)
(
m
1
+
d
1
)
)
(
w
r
-
m
1
-
d
1
)
(
m
1
(
1
-
w
α
)
+
d
1
)
.In the following, we consider the positive equilibrium denoted byE
*
=
(
N
*
,
P
*
,
Z
*
); then N
*, P
*, Z
* must satisfy
(3)
λ
-
d
N
-
r
N
P
k
1
+
N
+
α
m
1
P
+
β
m
2
Z
+
γ
μ
P
Z
k
2
+
P
=
0
,
(4)
w
r
N
k
1
+
N
-
m
1
-
d
1
-
μ
Z
k
2
+
P
=
0
,
(5)
η
μ
P
k
2
+
P
-
m
2
-
d
2
=
0
.
From (5) we have P
*
=
k
2
(
m
2
+
d
2
)
/
(
η
μ
-
m
2
-
d
2
), so P
* is positive if η
μ
-
m
2
-
d
2
>
0. In the following, we consider the existence of N
* and Z
* under the condition η
μ
-
m
2
-
d
2
>
0. Substituting P
* into (3) and (5), we get
(6)
λ
-
d
N
-
r
k
2
(
m
2
+
d
2
)
N
(
η
μ
-
m
2
-
d
2
)
(
k
1
+
N
)
+
α
m
1
k
2
(
m
2
+
d
2
)
η
μ
-
m
2
-
d
2
+
β
m
2
Z
+
γ
(
m
2
+
d
2
)
Z
η
=
0
,
w
r
N
k
1
+
N
-
m
1
-
d
1
-
(
η
μ
-
m
2
-
d
2
)
Z
k
2
η
=
0
.
That is,
(7)
-
d
N
2
+
a
1
N
+
k
1
a
2
Z
+
a
2
N
Z
+
a
3
=
0
,
(8)
Z
=
c
1
-
c
2
k
1
+
N
,
where a
1 = λ
-
d
k
1
+
(
k
2
(
m
2
+
d
2
)
(
α
m
1
-
r
)
/
(
η
μ
-
m
2
-
d
2
)
), a
2
=
(
β
m
2
+
(
γ
(
m
2
+
d
2
)
/
η
)
)
>
0, a
3
=
(
λ
k
1
+
(
α
m
1
k
1
k
2
(
m
2
+
d
2
)
/
(
η
μ
-
m
2
-
d
2
)
)
)
>
0, c
1
=
k
2
η
(
w
r
-
m
1
-
d
1
)
/
(
η
μ
-
m
2
-
d
2
), and c
2
=
w
r
k
1
k
2
η
/
(
η
μ
-
m
2
-
d
2
)
>
0. Equation (7) can also be rewritten as
(9)
Z
=
b
1
(
k
1
+
N
)
-
b
2
k
1
+
N
-
g
,
with
(10)
b
1
=
d
a
2
=
d
β
m
2
+
(
γ
(
m
2
+
d
2
)
/
η
)
>
0
,
b
2
=
a
3
-
k
1
(
d
k
1
+
a
1
)
a
2
b
2
=
k
1
k
2
(
m
2
+
d
2
)
r
c
c
c
c
c
×
(
(
η
μ
-
m
2
-
d
2
)
(
β
m
2
+
γ
(
m
2
+
d
2
)
η
)
)
-
1
>
0
,
g
=
2
d
k
1
+
a
1
a
2
g
=
(
(
λ
+
d
k
1
)
(
η
μ
-
m
2
-
d
2
)
c
c
c
c
c
c
c
+
k
2
(
m
2
+
d
2
)
(
α
m
1
-
r
)
)
c
c
c
c
×
(
(
η
μ
-
m
2
-
d
2
)
(
β
m
2
+
γ
(
m
2
+
d
2
)
η
)
)
-
1
.From (8), we have Z
′
|
N
=
c
2
/
(
k
1
+
N
)
2
>
0, so, Z is a monotonous increasing function about N, and the asymptotes are N
=
-
k
1 and Z
=
c
1. In the same way, for (9), Z
′
|
N
=
b
1
+
(
b
2
/
(
k
1
+
N
)
2
)
>
0, Z is a monotonous increasing function about N, and the asymptotes are N
=
-
k
1 and Z
=
b
1
(
k
1
+
N
)
-
g. Neglecting the left branches of figures of (8) and (9), the Figures of (8) and (9) are shown in Figures 1(a) and 1(b). If we want the curve of (8) and the curve of (9) to intersect in the first quadrant, we must have c
1
>
0; that is, w
r
-
m
1
-
d
1
>
0. If Z
=
0 in (8) and (9), we can calculate A
=
Δ
(
N
1
,
0
)
=
(
(
c
2
/
c
1
)
-
k
1
,
0
) and C
=
Δ
(
N
2
,
0
)
=
(
(
1
/
2
b
1
)
(
g
-
2
b
1
k
1
+
g
2
+
4
b
1
b
2
)
,
0
). While N
1
=
(
c
2
/
c
1
)
-
k
1
=
(
m
1
+
d
1
)
k
1
/
(
w
r
-
m
1
-
d
1
)
>
0 if w
r
-
m
1
-
d
1
>
0; that is, N
1
=
N
^.Figures of (8) and (9), where (a) is the figure of (9) and (b) is the figure of (8).
(a)
(b)IfN
1
<
N
2, then the figures of (8) and (9) only have a positive intersection point; that is,
(11)
k
2
(
m
2
+
d
2
)
(
w
r
-
m
1
-
d
1
)
(
m
1
(
1
-
w
α
)
+
d
1
)
+
w
(
m
2
+
d
2
-
η
μ
)
(
w
r
λ
-
(
λ
+
d
k
1
)
(
m
1
+
d
1
)
)
<
0
.
That is, λ
>
(
k
2
(
m
2
+
d
2
)
(
w
r
-
m
1
-
d
1
)
(
m
1
(
1
-
w
α
)
+
d
1
)
+
w
(
η
μ
-
m
2
-
d
2
)
d
k
1
(
m
1
+
d
1
)
) / w
(
η
μ
-
m
2
-
d
2
)
(
w
r
-
m
1
-
d
1
)
=
Δ
λ
*. Obviously, λ
*
>
λ
0. If N
1
=
N
2, then (8) and (9) only intersect on the positive N axis; that is, Z
=
0 and λ
=
λ
*; then P
*
=
P
^. In this condition, the equilibrium only has E
^
=
(
N
^
,
P
^
,
0
).According to the above statement, we know that the system (1) only has a positive equilibrium if η
μ
-
m
2
-
d
2
>
0, w
r
-
m
1
-
d
1
>
0, and λ
>
λ
*.The existence of equilibria is summarized as follows:(i)
E
0
=
(
N
0
,
0,0
) always exists;
(ii)
E
^
=
(
N
^
,
P
^
,
0
) exists if w
r
-
m
1
-
d
1
>
0 and λ
>
λ
0 hold;
(iii)
E
*
=
(
N
*
,
P
*
,
Z
*
) exists if w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
>
λ
*.


## 2.2. Stability of Equilibria of the System (1)

First, we know that all solutions of (1) are uniformly ultimately bounded, that is,Theorem 1.
For each positive solutionx
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) of system (1) there exists a constant M
>
0, such that N
(
t
)
≤
M, P
(
t
)
≤
M, and Z
(
t
)
≤
M with t large enough.Proof.
Define a functionL
(
t
,
x
) such that L
(
t
,
x
)
=
N
(
t
)
+
P
(
t
)
+
Z
(
t
). Then L
(
t
,
x
)
∈
V
0 and the upper right derivative of L
(
t
,
x
) along a solution of system (1) is described as
(12)
D
+
L
(
t
,
x
)
|
(
1
)
=
λ
-
d
N
+
(
α
-
1
)
m
1
P
+
(
β
-
1
)
m
2
Z
-
d
1
P
-
d
2
Z
+
(
w
-
1
)
r
N
P
k
1
+
N
+
(
γ
+
η
-
1
)
μ
P
Z
k
2
+
P
≤
λ
-
d
N
-
(
d
1
+
(
1
-
α
)
m
1
)
P
-
(
d
2
+
(
1
-
β
)
m
2
)
Z
≤
λ
-
d
^
L
(
t
,
x
)
,
where d
^
=
min
⁡
{
d
,
d
1
+
(
1
-
α
)
m
1
,
d
2
+
(
1
-
β
)
m
2
}. Hence D
+
L
(
t
,
x
)
|
(
1
)
≤
-
d
^
L
(
t
,
x
)
+
λ. So, we have
(13)
L
(
t
,
x
)
≤
λ
d
^
+
(
L
(
0
)
-
λ
d
^
)
e
-
d
^
t
⟶
λ
d
^
,
t
⟶
∞
.
Therefore, L
(
t
,
x
) is ultimately bounded by a constant and there exists a constant M
>
0, such that N
(
t
)
≤
M, P
(
t
)
≤
M, and Z
(
t
)
≤
M, for each solution x
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) of system (1) with t large enough. The proof is complete.In the following, we letE
0
=
(
N
0
,
P
0
,
Z
0
) be any equilibrium of system (1); then the Jacobian matrix about E
0 is given by(
*
)
J
|
E
0
=
(
-
d
-
r
k
1
P
0
(
k
1
+
N
0
)
2
-
r
N
0
k
1
+
N
0
+
α
m
1
+
γ
μ
k
2
Z
0
(
k
2
+
P
0
)
2
β
m
2
+
γ
μ
P
0
k
2
+
P
0
w
r
k
1
P
0
(
k
1
+
N
0
)
2
w
r
N
0
k
1
+
N
0
-
m
1
-
d
1
-
μ
k
2
Z
0
(
k
2
+
P
0
)
2
-
μ
P
0
k
2
+
P
0
0
η
μ
k
2
Z
0
(
k
2
+
P
0
)
2
η
μ
P
0
k
2
+
P
0
-
m
2
-
d
2
)
,denoted by
(14)
J
|
E
0
=
(
a
11
a
12
a
13
a
21
a
22
a
23
0
a
32
a
33
)
;
then the characteristic equation about E
0 is: λ
3
+
p
1
λ
2
+
p
2
λ
+
p
3
=
0, where
(15)
p
1
=
-
(
a
11
+
a
22
+
a
33
)
,
p
2
=
a
11
a
22
+
a
11
a
33
+
a
22
a
33
-
a
21
a
12
-
a
23
a
32
,
p
3
=
a
12
a
21
a
33
-
a
11
a
22
a
33
-
a
21
a
32
a
13
+
a
11
a
23
a
32
.About the stability of the equilibriumE
0, E
^, and E
*, we have the following results.Theorem 2.
The steady stateE
0 of system (1) is locally asymptotically stable if m
1
+
d
1
>
(
w
r
N
0
/
(
k
1
+
N
0
)
); that is, w
r
-
m
1
-
d
1
>
0 and λ
<
λ
0.Proof.
For equilibriumE
0
=
(
N
0
,
0,0
)
=
(
λ
/
d
,
0,0
), the Jacobian matrix evaluated at E
0 is
(16)
J
|
E
0
=
(
-
d
-
r
N
0
k
1
+
N
0
+
α
m
1
β
m
2
0
w
r
N
0
k
1
+
N
0
-
m
1
-
d
1
0
0
0
-
(
m
2
+
d
2
)
)
;
obviously, two of eigenvalues are λ
1
=
-
d
<
0 and λ
2
=
-
(
m
2
+
d
2
)
<
0; the other eigenvalue is λ
3
=
(
w
r
N
0
/
(
k
1
+
N
0
)
)
-
m
1
-
d
1. The steady state E
0 is asymptotically stable if and only if the eigenvalues of the Jacobian matrix at E
0 have negative real parts; the proof is complete.Theorem2 implies that if the nutrient-phytoplankton conversion rate is less than phytoplankton loss rate, then both phytoplankton and zooplankton population will become extinct.Theorem 3.
The steady stateE
^ of system (1) is locally asymptotically stable if m
2
+
d
2
-
η
A
2
>
0; that is, w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
<
λ
*.Proof.
For equilibriumE
^
=
(
N
^
,
P
^
,
0
), from the second equation of the system (1), we have (
w
r
N
^
/
(
k
1
+
N
^
)
)
-
m
1
-
d
1
=
0; that is, (
r
N
^
/
(
k
1
+
N
^
)
)
=
(
(
m
1
+
d
1
)
/
w
), so -
(
r
N
^
/
(
k
1
+
N
^
)
)
+
α
m
1
+
(
γ
μ
k
2
Z
^
/
(
k
2
+
P
^
)
2
)
=
-
(
(
m
1
+
d
1
)
/
w
)
+
α
m
1
+
(
γ
μ
k
2
Z
^
/
(
k
2
+
P
^
)
2
)
=
-
(
(
1
-
w
α
)
m
1
+
d
1
)
/
w
<
0. Hence (
*
) reduces to(17)
J
|
E
^
=
(
-
d
-
r
k
1
P
^
(
k
1
+
N
^
)
2
-
(
1
-
w
α
)
m
1
+
d
1
w
β
m
2
+
γ
μ
P
^
k
2
+
P
^
w
r
k
1
P
^
(
k
1
+
N
^
)
2
0
-
μ
P
^
k
2
+
P
^
0
0
η
μ
P
^
k
2
+
P
^
-
m
2
-
d
2
)
≜
(
-
d
-
A
1
-
d
~
β
m
2
+
γ
A
2
w
A
1
0
-
A
2
0
0
η
A
2
-
m
2
-
d
2
)
,with A
1
=
r
k
1
P
^
/
(
k
1
+
N
^
)
2
>
0, A
2
=
μ
P
^
/
(
k
2
+
P
^
)
>
0, and d
~
=
(
(
1
-
w
α
)
m
1
+
d
1
)
/
w
>
0. Then λ
3
+
p
1
λ
2
+
p
2
λ
+
p
3
=
0, where
(18)
p
1
=
(
d
+
A
1
)
+
(
m
2
+
d
2
-
η
A
2
)
,
p
2
=
(
d
+
A
1
)
(
m
2
+
d
2
-
η
A
2
)
+
d
~
w
A
1
,
p
3
=
d
~
w
A
1
(
m
2
+
d
2
-
η
A
2
)
;
then p
1
p
2
-
p
3
=
(
d
+
A
1
)
2
(
m
2
+
d
2
-
η
A
2
)
+
d
1
w
A
1
(
d
+
A
1
)
+
(
d
+
A
1
)
(
m
2
+
d
2
-
η
A
2
)
2. Obviously, p
i
>
0 (i
=
1,2
,
3) and p
1
p
2
-
p
3
>
0 if m
2
+
d
2
-
η
A
2
>
0. According to Routh-Hurwitz criteria, E
^ is locally stable if m
2
+
d
2
-
η
A
2
>
0. The proof is complete.Theorem3 implies that if the growth rate of zooplankton is less than its loss rate, then the zooplankton population will die out and the phytoplankton population will survive on the nutrient.Theorem 4.
The steady stateE
* of system (1) is locally asymptotically stable if h
1
>
0 and h
1
h
2
-
h
3
>
0.Proof.
For equilibriumE
*
=
(
N
*
,
P
*
,
Z
*
), (
*
) reduces to(19)
J
|
E
*
=
(
-
d
-
r
k
1
P
*
(
k
1
+
N
*
)
2
-
r
N
*
k
1
+
N
*
+
α
m
1
+
γ
μ
k
2
Z
*
(
k
2
+
P
*
)
2
β
m
2
+
γ
μ
P
*
k
2
+
P
*
w
r
k
1
P
*
(
k
1
+
N
*
)
2
w
r
N
*
k
1
+
N
*
-
m
1
-
d
1
-
μ
k
2
Z
*
(
k
2
+
P
*
)
2
-
μ
P
*
k
2
+
P
*
0
η
μ
k
2
Z
*
(
k
2
+
P
*
)
2
0
)
≜
(
-
d
-
B
1
-
B
2
+
α
m
1
+
γ
B
3
β
m
2
+
γ
B
4
w
B
1
w
B
2
-
m
1
-
d
1
-
B
3
-
B
4
0
η
B
3
0
)
,with B
1
=
r
k
1
P
*
/
(
k
1
+
N
*
)
2, B
2
=
r
N
*
/
(
k
1
+
N
*
), B
3
=
μ
k
2
Z
*
/
(
k
2
+
P
*
)
2, and B
4
=
μ
P
*
/
(
k
2
+
P
*
). Then λ
3
+
h
1
λ
2
+
h
2
λ
+
h
3
=
0, where
(20)
h
1
=
(
d
+
B
1
)
+
(
m
1
+
d
1
+
B
3
-
w
B
2
)
,
h
2
=
(
d
+
B
1
)
(
m
1
+
d
1
+
B
3
-
w
B
2
)
ccccc
+
w
B
1
(
B
2
-
α
m
1
-
γ
B
3
)
+
B
4
η
B
3
,
h
3
=
-
w
B
1
η
B
3
(
β
m
2
+
γ
B
4
)
+
B
4
(
d
+
B
1
)
η
B
3
;
then
(21)
h
1
h
2
-
h
3
=
(
d
+
B
1
)
2
(
m
1
+
d
1
+
B
3
-
w
B
2
)
+
(
d
+
B
1
)
w
B
1
(
B
2
-
α
m
1
-
γ
B
3
)
+
(
d
+
B
1
)
(
m
1
+
d
1
+
B
3
-
w
B
2
)
2
+
w
B
1
(
B
2
-
α
m
1
-
γ
B
3
)
·
(
m
1
+
d
1
+
B
3
-
w
B
2
)
+
B
4
η
B
3
(
m
1
+
d
1
+
B
3
-
w
B
2
)
+
w
B
1
η
B
3
(
β
m
2
+
γ
B
4
)
.
From (5), we have w
B
2
-
m
1
-
d
1
-
(
μ
Z
*
/
(
k
2
+
P
*
)
)
=
0, so w
B
2
-
m
1
-
d
1
-
B
3
=
(
μ
Z
*
/
(
k
2
+
P
*
)
)
-
B
3
=
(
μ
Z
*
/
(
k
2
+
P
*
)
)
-
(
μ
k
2
Z
*
/
(
k
2
+
P
*
)
2
)
=
μ
P
*
Z
*
/
(
k
2
+
P
*
)
2
>
0. According to Routh-Hurwitz criteria, we have Theorem 4.
The existence and local stability of equilibria are summarized as follows:(i)
E
0
=
(
N
0
,
0,0
) always exists; E
0 is local stable if w
r
-
m
1
-
d
1
>
0 and λ
<
λ
0; E
0 is unstable if w
r
-
m
1
-
d
1
>
0 and λ
>
λ
0;
(ii)
E
^
=
(
N
^
,
P
^
,
0
) exists if w
r
-
m
1
-
d
1
>
0 and λ
>
λ
0 hold; E
^ is local stable if w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
0
<
λ
<
λ
*; E
^ is unstable if w
r
-
m
1
-
d
1
>
0
,
η
μ
-
m
2
-
d
2
>
0 and λ
>
λ
*;
(iii)
E
*
=
(
N
*
,
P
*
,
Z
*
) exists if w
r
-
m
1
-
d
1
>
0, η
μ
-
m
2
-
d
2
>
0, and λ
>
λ
*; E
* is local stable if h
1
>
0 and h
1
h
2
-
h
3
>
0.
From analysis, we know that the number of equilibrium increases with the increase ofλ: from one equilibrium to three equilibria. The stability of equilibria changes with changing λ; when equilibrium E
* becomes unstable, algal bloom will break out. So, controlling the amount of input nutrient is very important in controlling water bloom; for more detail we can see the Discussion section.


## 3. Impulsive Input Nutrient

Pulse input nutrient can be defined as the repeated application of input nutrient. Some factory discharge pollution periodicity, and rains are also seasonally. So in this section, we assume that the pulse scheme proposes to input nutrientλ in a single pulse, applied every τ period (
τ
>
0
). So the expression λ in system (1) is presented in pulses; we have the following impulsive differential equations:(22)
N
˙
=
-
d
N
-
r
N
P
k
1
+
N
+
α
m
1
P
+
β
m
2
Z
+
γ
μ
P
Z
k
2
+
P
,
P
˙
=
w
r
N
P
k
1
+
N
-
m
1
P
-
d
1
P
-
μ
P
Z
k
2
+
P
,
Z
˙
=
η
μ
P
Z
k
2
+
P
-
m
2
Z
-
d
2
Z
,
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
t
≠
n
τ
,
N
(
t
+
)
=
N
(
t
-
)
+
λ
,
P
(
t
+
)
=
P
(
t
-
)
,
Z
(
t
+
)
=
Z
(
t
-
)
,
c
c
c
c
c
c
c
c
t
=
n
τ
,
where n
=
0,1
,
2
,
…
,
τ is the period of input nutrient, n
τ
+ is the time at which we apply the nth pulse, and n
τ
- is the time just before applying the nth pulse. Other denotations are the same as system (1).Before we consider the stability of the boundary periodic solution, we need the following lemmas.Lemma 5.
Supposex
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) is a solution of system (22) with initial values x
(
0
+
)
≥
0; then x
(
t
)
≥
0; that is, N
(
t
)
≥
0
,
P
(
t
)
≥
0 and Z
(
t
)
≥
0; further x
(
t
)
>
0 for all t
≥
0 if x
(
0
+
)
>
0.The result can be easily obtained from system (22), so we omit it.Theorem 6.
For each positive solutionx
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) of system (22) there exists a constant M
>
0, such that N
(
t
)
≤
M, P
(
t
)
≤
M, and Z
(
t
)
≤
M with t large enough.Proof.
Define a functionL
(
t
,
x
) such that L
(
t
,
x
)
=
N
(
t
)
+
P
(
t
)
+
Z
(
t
). Then L
(
t
,
x
)
∈
V
0 and the upper right derivative of L
(
t
,
x
) along a solution of (22) is described as
(23)
D
+
L
(
t
,
x
)
|
(
22
)
=
-
d
N
+
(
α
-
1
)
m
1
P
+
(
β
-
1
)
m
2
Z
-
d
1
P
-
d
2
Z
+
(
w
-
1
)
r
N
P
k
1
+
N
+
(
γ
+
η
-
1
)
μ
P
Z
k
2
+
P
≤
-
d
N
-
d
1
P
-
d
2
Z
≤
-
d
^
L
(
t
,
x
)
,
where d
^
=
min
⁡
{
d
,
d
1
,
d
2
}. Hence D
+
L
(
t
,
x
)
|
(
22
)
≤
-
d
^
L
(
t
,
x
). When t
=
n
T, we obtain L
(
n
τ
+
)
≤
L
(
n
τ
)
+
λ. So, we have
(24)
L
(
t
)
=
L
(
0
+
)
e
-
d
^
t
+
λ
e
d
^
τ
e
d
^
τ
-
1
(
e
-
d
^
(
t
-
n
τ
)
-
e
-
d
^
t
)
⟶
λ
e
d
^
τ
e
d
^
τ
-
1
,
t
⟶
∞
.
Therefore L
(
t
,
x
) is ultimately bounded by a constant and there exists a constant M
>
0, such that N
(
t
)
≤
M, P
(
t
)
≤
M, and Z
(
t
)
≤
M, for each solution x
(
t
)
=
(
N
(
t
)
,
P
(
t
)
,
Z
(
t
)
) of system (22) with t large enough. The proof is complete.For convenience, we give some basic properties of the following system:(25)
N
˙
(
t
)
=
-
d
N
(
t
)
,
t
≠
n
τ
,
Δ
N
(
t
)
=
λ
,
t
=
n
τ
,
N
(
0
+
)
=
N
0
≥
0
.Like [27], we have the following Lemma.Lemma 7.
System (25) has a positive periodic solution N
~
(
t
) and, for every solution N
(
t
) of system (25) with initial value N
(
0
+
)
=
N
0
≥
0, we have N
(
t
)
→
N
~
(
t
) as t
→
∞; moreover, N
(
t
)
≥
N
~
(
t
) if N
0
≥
λ
/
(
1
-
exp
⁡
(
-
d
τ
)
) and N
(
t
)
<
N
~
(
t
) if N
0
<
λ
/
(
1
-
exp
⁡
(
-
d
τ
)
), where N
~
(
t
)
=
λ
exp
⁡
(
-
d
(
t
-
n
τ
)
)
/
(
1
-
exp
⁡
(
-
d
τ
)
), t
∈
(
n
τ
,
(
n
+
1
)
τ
], and n
∈
Z
+, N
~
(
0
+
)
=
λ
/
(
1
-
exp
⁡
(
-
d
τ
)
).Therefore, the system (22) has a boundary eradication periodic solution
(26)
(
N
~
(
t
)
,
0,0
)
=
(
λ
exp
⁡
(
-
d
(
t
-
n
τ
)
)
1
-
exp
⁡
(
-
d
τ
)
,
0,0
)
,
cccccccccccccccccc
for
t
∈
(
n
τ
,
(
n
+
1
)
τ
]
.
About the stability of the boundary periodic solution of the system (22), we can use the similar method of [27] and get the following theorem.Theorem 8.
If(
(
k
1
+
(
λ
-
k
1
)
exp
⁡
(
-
d
τ
)
)
/
(
k
1
+
λ
-
k
1
exp
⁡
(
-
d
τ
)
)
)
>
exp
⁡
(
-
(
d
τ
(
m
1
+
d
1
)
)
/
w
r
), then the periodic boundary solution (
N
~
(
t
)
,
0,0
) of system (22) is locally asymptotically stable.We denoteλ
^ = (
k
1
(
1
-
exp
⁡
(
-
d
τ
)
)
(
-
1
+
exp
⁡
(
-
(
m
1
+
d
1
)
d
τ
/
w
r
)
)
)
/
(
exp
⁡
(
-
d
τ
)
-
exp
⁡
(
-
(
m
1
+
d
1
)
d
τ
/
w
r
)
), τ
^ = 2
r
w
λ
(
w
r
-
m
1
-
d
1
)
/
d
(
λ
w
2
r
2
-
λ
(
m
1
+
d
1
)
2
+
2
w
r
k
1
(
m
1
+
d
1
)
). According to Theorem 8 we can easily obtain the following results.Corollary 9.
(i) Ifw
r
-
m
1
-
d
1
>
0, then (
N
~
(
t
)
,
0,0
) is locally asymptotically stable provided that τ
>
τ
^. (ii) If w
r
-
m
1
-
d
1
>
0, then (
N
~
(
t
)
,
0,0
) is locally asymptotically stable provided that λ
<
λ
^.Corollary9 implies that the phytoplankton and zooplankton will disappear if the length of impulsive period is large enough or the impulsive input rate of the nutrient is small enough.


## 4. The Effect of Temperature on a Lake Model

Temperature is an important factor in the growth of algae. So, in order to accurately describe the growth of algae, the effect of temperature on the algal growth was considered in [20, 28–30]. In [28], f
(
T
)
=
e
-
(
2.3
/
15
)
|
T
-
T
opt
| was considered as the coefficient of temperature, where T
opt is the optimal temperature for phytoplankton growth; f
(
T
)
=
(
T
/
T
opt
)
e
1
-
(
T
/
T
opt
) and f
(
T
)
=
e
-
2.3
|
(
T
-
T
opt
)
/
(
T
opt
-
T
min
⁡
)
| were considered in [29, 30], respectively, where T
min
⁡ is minimum temperature for phytoplankton growth. In fact, f
(
T
)
=
e
-
(
2.3
/
15
)
|
T
-
T
opt
| and f
(
T
)
=
e
-
2.3
|
(
T
-
T
opt
)
/
(
T
opt
-
T
min
⁡
)
| are the same form and the latter is more generality, so, we only need to compare f
(
T
)
=
e
-
2.3
|
(
T
-
T
opt
)
/
(
T
opt
-
T
min
⁡
)
| with f
(
T
)
=
(
T
/
T
opt
)
e
1
-
(
T
/
T
opt
). If we let T
opt
=
28, T
min
⁡
=
5, we can get Figure 2. From this figure, we can know that f
(
T
)
=
e
-
(
2.3
/
23
)
|
T
-
28
| changes more rapidly than f
(
T
)
=
(
T
/
28
)
e
1
-
(
T
/
28
). According to interrelated literature, we take the growth of the algae as f
(
T
)
=
e
-
2.3
|
(
T
-
T
opt
)
/
(
T
opt
-
T
min
⁡
)
|.Figure 2
f
1
(
T
)
=
(
T
/
28
)
e
1
-
(
T
/
28
), f
2
(
T
)
=
e
-
(
2.3
/
23
)
|
T
-
28
|.In fact, the effect of temperature on the algal growth can be described as follows: algae can not grow if the temperature is too low; the metabolism rate will become large if the temperature is too high; organic matter of photosynthesis synthesis is decomposed by metabolism rapidly, so it also limits the increasing of algae biomass. So, the effect of temperature on the growth of algae in general can be expressed as the following form:f
(
T
)
=
e
-
K
T
|
T
-
T
opt
| [20, page 349], where K
T is the temperature impact coe±cient of the growth of algae. Taking into account thatMicrocystis aeruginosa population is the dominant algae species in Taihu Lake, generally the optimal temperature forMicrocystis aeruginosa growth is thought of as 2
8
°C [20, page 213], so we take T
opt
=
2
8
°C.The effect of temperature on the growth of zooplankton is similar to the effect of temperature on the growth of phytoplankton [30].Taking into account the effect of temperature on the growth of plankton and incorporating the paper of Taihu Lake [20], we consider the following model of Taihu Lake:
(27)
N
˙
(
t
)
=
λ
-
d
N
(
t
)
-
u
(
t
)
P
(
t
)
+
α
m
1
P
(
t
)
+
β
m
2
Z
(
t
)
+
γ
v
(
t
)
Z
(
t
)
,
P
˙
(
t
)
=
w
u
(
t
)
P
(
t
)
-
m
1
P
(
t
)
-
d
1
P
(
t
)
-
v
(
t
)
Z
(
t
)
,
Z
˙
(
t
)
=
η
v
(
t
)
Z
(
t
)
-
m
2
Z
(
t
)
-
d
2
Z
(
t
)
,
where u
(
t
)
=
u
max
⁡
f
1
(
T
)
f
2
(
N
), v
(
t
)
=
v
max
⁡
f
3
(
T
)
f
4
(
P
), m
1
=
m
1
max
⁡
f
1
(
T
), m
2
=
m
2
max
⁡
f
3
(
T
), f
1
(
T
)
=
e
-
2.3
|
(
T
-
T
1
opt
)
/
(
T
1
opt
-
T
1
min
⁡
)
|, f
2
(
N
)
=
N
/
(
k
1
+
N
), f
3
(
T
)
=
e
-
2.3
|
(
T
-
T
2
opt
)
/
(
T
2
opt
-
T
2
min
⁡
)
|, and f
4
(
P
)
=
P
/
(
k
2
+
P
). The meaning of parameters is as follows: u
max
⁡
(
v
max
⁡
) is the maximum growth rate of phytoplankton (zooplankton), T is temperature, m
1
max
⁡
(
m
2
max
⁡
) is the maximum death and respiration rate of phytoplankton (zooplankton), T
1
opt
(
T
2
opt
) is the optimal temperature of phytoplankton (zooplankton) growth, and T
1
min
⁡
(
T
2
min
⁡
) is the minimum temperature of phytoplankton (zooplankton) growth. From [27, 28, 31, 32], we have T
1
opt
=
2
8
°C, T
2
opt
=
2
8
°C, T
1
min
⁡
=
5
°C, T
2
min
⁡
=
5
°C, u
max
⁡
=
1.27
d
-
1, v
max
⁡
=
0.35
d
-
1, m
1
max
⁡
=
1.02
d
-
1, and m
2
max
⁡
=
0.145
d
-
1. We take d
=
d
1
=
d
2, where d
=
Q
/
V, Q is the outflow of the lake, V is the volume of the lake. For Taihu Lake [27], k
1
=
0.35
mg
L
-
1, k
2
=
0.5
mg
L
-
1, Q
=
1.56
×
1
0
7
m
3
d
-
1, and V
=
4.43
×
1
0
9
m
3, so we have d
=
0.323
d
-
1. If we let w
=
1, α
=
1, β
=
0.5, η
=
0.8, γ
=
0.1, and λ
=
0.08, then the positive equilibrium is unstable and we can get Figures 3(a), 3(b), and 3(c) for T
=
10,18,27, respectively. From comparing these figures, it consumes different lengths of time for the amount of algae to get the maximum. The time at which algae get the peak is the shortest for T
=
27, so comparatively suitable temperature is easier to break out water bloom. This agrees with the reality and provides a theoretical basis for understanding water bloom.(a), (b), and (c) are time series of nutrient, phytoplankton, and zooplankton forT
=
10,18,27, respectively, where the thick line is the phytoplankton curve, the middle line is the nutrient curve, and the thin line is zooplankton curve.
(a)
(b)
(c)Water bloom often breaks out in Taihu Lake in recent years. Last year, water bloom broke out in Taihu Lake on May; it caused Wuxi citizen drinking water crisis and brought to the lives of people a lot of inconvenience. In order to control (or reduce) water bloom of Taihu Lake immediately, one of measures which the government took is to make Yangtze River water flow into Taihu Lake. This measure reflected in model (27) is that parameter d is changing, that is, increasing coefficient d.If we increase the parameterd and let d
=
0.01, T
=
18, then we can know that the positive equilibrium is unstable for d
=
0.0035 (Figure 4(a)) and stable for d
=
0.01 (Figure 4(b)). Again, water bloom corresponding to the unstable positive equilibrium, normal state corresponding to the stable equilibrium, we know that water bloom disappears if we increase d from 0.0035 to 0.01.(a), (b), and (c) are time series of nutrient, phytoplankton, and zooplankton ford
=
0.0035,0.01,0.014, respectively, where the thick line is the phytoplankton curve, the middle line is the nutrient curve, and the thin line is zooplankton curve.
(a)
(b)
(c)And the maximum biomass of algae is about8.5 for d
=
0.0035 and is about 3.3 for d
=
0.01 nearby t
=
400; we also see that the biomass of algae is decreasing also. If we increase d to 0.014 again (Figure 4(c)), then the positive equilibrium disappears and the boundary equilibrium E
^ is stable. The biomass of algae gets the maximum value 1.824 and gets stable finally. We can see that the biomass of algae is decreasing with the increase of parameter d. So, from above analysis, we give theoretical explanation for the government measure to Taihu Lake, and this indicates that our model is reasonable for Taihu Lake at the same time.Controlling eutrophication of Taihu Lake fundamentally is a long term project. In the long term, we need to control the nutrient from outside and make a series of strict standards of discharging sewage step by step, that is, decreasing parameterλ; this will take long time to realize it. In order to control water bloom in a short term, reduce the negative impact on people’s lives immediately; we can increase the flow of the water, that is increasing parameter d; this agrees with the government measure.


## 5. Discussion

In the present research, the phenomenon of phytoplankton bloom has received much attention among experimental ecologists and mathematical ecologists. There are many papers where phytoplankton blooms have been modelled through different aspects [2–8]. In this paper, we study the bloom phenomenon by controlling outside input nutrient in an aquatic environment. Motivated by [12, 13], we have also incorporated the phenomena of nutrient recycling which have established importance in the context of a real open marine ecosystem; that is, we consider decomposition of phytoplankton, decomposition of zooplankton, and excretion of zooplankton in the nutrient equation. The difference between [13] and this paper is that [13] only had numerical results and this paper gives not only theoretical analysis but also numerical simulations. Different recycling effects make the situation more complicated than a simple food chain system, considering the amount of input nutrient to be a control parameter which is seldom.First, we consider a continuous input nutrient model, we obtain that the system has two boundary equilibria, and we analyze their stability. If the nutrient-phytoplankton conversion rate is less than phytoplankton loss rate, that is, the amount of input nutrient is less than some critical value, then both phytoplankton and zooplankton population will become extinct; see Theorem2. If the growth rate of zooplankton is less than its loss rate, that is, the amount of input nutrient is less than some critical value and larger than another critical value, then the zooplankton population will die out and the phytoplankton population will survive on the nutrient; see Theorem 3. We also investigate the existence of the positive equilibrium and discuss the stability of positive equilibrium; see Theorem 4.According to Table1, if we let d
=
0.1, d
1
=
0.1, d
2
=
0.1, r
=
0.3, k
1
=
0.1, k
2
=
0.05, m
1
=
0.15, m
2
=
0.1, β
=
0.5, γ
=
0.5, μ
=
0.8, η
=
0.3, w
=
1, and α
=
1, then we can calculate λ
0
=
0.05, λ
*
=
0.075. If we let λ
=
0.01, then the condition of Theorem 2 is satisfied (
0.01
<
0.05
), and then both phytoplankton and zooplankton population will become extinct; that is, E
0
=
(
N
0
,
0,0
)
=
(
0.1,0
,
0
) is locally asymptotically stable. If we increase λ to 0.06, then the condition of Theorem 3 is satisfied (
0.06
<
0.075
), and then the zooplankton population will die out and the phytoplankton population will survive on the nutrient; that is, E
^
=
(
N
^
,
P
^
,
0
)
=
(
0.5,0.1,0
) is locally asymptotically stable. If we continue increasing λ to 0.08, then the condition of Theorem 4 is satisfied (
h
1
=
0.11498
>
0
,
h
1
h
2
-
h
3
=
0.00017
>
0
), and the positive equilibrium E
*
=
(
N
*
,
P
*
,
Z
*
)
=
(
0.5462,0.25,0.00134
) is local asymptotically stability. Once again, if we continue increasing λ to 0.2, then the condition of Theorem 4 is not satisfied (
h
1
=
0.07438
>
0
,
h
1
h
2
-
h
3
=
-
0.00022
<
0
), and then the positive equilibrium E
*
=
(
N
*
,
P
*
,
Z
*
)
=
(
1.7144,0.25,0.01255
) will become unstable, species densities will oscillate vast scale, and then algae bloom will break out; this corresponds to eutrophic. So, one way controlling algae bloom is to control the amount of the limiting nutrient input from the environment, that is, parameter λ.Then, we consider the model of pulse input nutrient, that is, nutrient flow into the lake everyτ period. We obtain the exact boundary periodic solution of the impulsive input nutrient system. Using Floquet theory for the impulsive equation, small-amplitude perturbation skills, and techniques of comparison, we prove that the boundary periodic solution is locally asymptotically stable if some conditions are satisfied, that is, Theorem 8. If we let λ
=
0.1
,
τ
=
10, other parameters are the same as above; then λ
^
=
0.5316
(
0.1
<
0.5316
); that is, Corollary 9 (ii) is satisfied. From numerical simulation, with the increase of λ, we guess there is nutrient-phytoplankton coexistence and zooplankton extinct periodic solution denoted by (
N
¯
(
t
)
,
P
¯
(
t
)
,
0
) (see Figure 5). And the behavior of the system changes with the increase of λ again (see Figure 6). We leave these for future study.Existence of the periodic boundary solution(
N
¯
(
t
)
,
P
¯
(
t
)
,
0
) of system (22). (a), (b), and (c) are time series of nutrient, phytoplankton, and zooplankton, respectively, where λ
=
1.5.
(a)
(b)
(c)Another phenomenon of system (22). (a), (b), and (c) are time series of nutrient, phytoplankton, and zooplankton, respectively, where λ
=
1.7.
(a)
(b)
(c)Finally, the effect of temperature on the algal growth was considered; it is more accurate to describe the growth of algae with consideration of temperature. From the example of Taihu Lake, we analyze the effect of temperature on the system and give some numerical simulation. We also give theoretical analysis of measures which the government took during water bloom of Taihu Lake. Our results suggest that our lake model is correct and controlling the amount of input rate is an important measure in controlling water bloom.



---

*Source: 462946-2014-05-06.xml*