Dynamics of a periodic stoichiometric model with application in predicting and controlling algal bloom in Bohai Sea off China

1.   Introduction

Reported from China's Xinhua news service in July 2011, "Almost 200 square miles of the Yellow Sea off China are covered by a massive bloom of green algae." The algal bloom was first seen in the Yellow Sea in 2007. However, researchers do not know what causes the massive algal blooms. One main culprit for algal bloom is the increase of nitrogen. Compared to the Yellow Sea, Bohai Sea is more severe in algal bloom, while East China Sea is in much better situation because nitrogen is difficult to accumulate in such an open sea ^[5].

The algal bloom is a sudden increase of algae in an aquatic ecosystem. It results in hypoxic conditions which greatly threatens marine life. The algal bloom can cause negative impacts to other living organisms via toxin production. The harmful algal bloom (HAB) can have large and varied impacts on marine ecosystems although only 2% of algal species are harmful ^[2]. HAB has been increasing and spreading along the coast of China in the past two decades ^[4,6,7]. To protect the Seas off China, it is urgent to take action now.

In this paper, we apply and modify a recently developed stoichiometric model ^[3] to study the algal bloom event. Our modeling framework suggests that blocking nutrient input from rivers may be a more effective way than the commonly used control method --- removing algae periodically. In addition, our mathematical model fits 2007 Bohai Sea data well, validating the model and its findings.

The remainder of this paper is organized as follows. In section 2, we develop a stoichiometric algal growth model incorporating a season-driven light intensity. In section 3, we characterize the model dynamics, via the method of qualitative analysis, coincidence degree theory and Liapunov direct method, such as the positive invariance, dissipativity, boundary and internal dynamics. Section 4 includes some numerical simulations to confirm the theoretical results and further explore the impacts of the seasonal light intensity and the nutrient availability on the algal dynamics. We discuss two control methods including removing algae (RA) periodically and blocking nutrient (BN) input from rivers constantly via our modeling approach in section 5. By comparison, the BN method is a more effective way to terminate algal bloom in Yellow Sea off China. In section 6, We fit our model to the Bohai Sea data for validation.

2.   Model formulation

Similar to the stoichiometrically derived algal growth model in ^[3], we define the variables: $x$ - the algal carbon concentration, $n$ - the algal nitrogen concentration. The model is provided by

where $g(T-n) = c(T-n)/[a+(T-n)]$. The parameters and their biological significance are listed in Table 1. The algal growth model captures key biological features of light- and nutrient- dependent algal growth and characterizes both carbon and phosphorus limitations.

Letting $r: = \mu-d$ and $\bar{q}: = \mu q/(\mu-d)>q$ leads to

The key mathematical findings are listed as follows ^[3].

● The open trapezoid domain $\Omega = \{(x, n)\in \mathbb{R}_+^2: 0<x<k, qx<n<T\}$, where $k = \min\{K, T/\bar{q}\}$, is positively invariant for the system (2.3)-(2.4).

● When $K\geq T/\bar{q}$ and $g(T)\leq d\bar{q}$, the origin $E_0 = (0, 0)$ is globally asymptotically stable.

● When $K\geq T/\bar{q}$ and $g(T)>d\bar{q}$, the internal equilibrium $E_1 = (\bar{x}, \bar{n})$ with $\bar{n} = \bar{q}\bar{x}$ and $g(T-\bar{n}) = d\bar{q}$ is globally asymptotically stable.

● When $n^*/\bar{q}<K\leq T/\bar{q}$ and $g(T)\leq d\bar{q}$ where $g(T-n^*)k = dn^*$, the origin $E_0$ is globally asymptotically stable.

● When $n^*/\bar{q}<K\leq T/\bar{q}$ and $g(T)>d\bar{q}$, the internal equilibrium $E_1$ is globally asymptotically stable.

● When $K\leq n^*/\bar{q}$, the boundary equilibrium $E_2 = (K, n^*)$ is globally asymptotically stable.

Given that $g(T-n) = c(T-n)/(a+(T-n))$, we have

Assume that light is sufficient (i.e., $K$ is large enough) and $g(T)>d\bar{q}$ (i.e., $T>ad\bar{q}/(c-d\bar{q})$), then

provides the algal density (measured by carbon biomass). Plot $\bar{x}$ versus $T$ (nutrient level) to indicate how nutrient richness enhances the algal bloom. The slope, indicating how fast algal density increases with respect to increasing nutrient density, is $1/\bar{q} = (\mu-d)/\mu q = 1/q-d/\mu q$, which depends on algal species. It is less likely for algal bloom to occur when the algal species has smaller $\mu$ (maximum growth rate), larger $d$ (recycling rate), or larger $q$ (minimum structural N:C ratio). Since there is no officially defined threshold level for algal bloom, we assume the algal bloom level $x_b$ (some threshold of algal density) in this simulation, then we can compute the threshold of nutrient availability

over which algal bloom occurs. The threshold is an increasing function of $\bar{q}$, thus the larger $\bar{q}$ (i.e., smaller $1/\bar{q}$) increases $T_b$ such that algal bloom is less likely to occur. Figure 1 shows that the $T_b$ can have a large variation for different algal species. This explains the biological phenomenon that typically only one or a small number of algal species are involved in algal bloom. The algal species on the left-up corner of Figure 1 is difficult to have outbreak since the higher threshold the lower probability of algal bloom. This result may explain algal bloom rarely occurs in some coasts or some lakes. The nutrient density in algae is given by

When algal density reaches the carrying capacity in the case $K\leq n^*/\bar{q}$, its internal nutrient density is given by

The environmental variation is time-dependent and is one of the processes that the autonomous model ignores. Natural environment is physically highly variable, and in response, some vital rates of populations, are usually subject to seasonal fluctuations and vary greatly in time, which is more biologically practical. Therefore, realistic models should take into account the seasonal effect. A desired model taking into account such environmental fluctuations must be nonautonomous in mathematics. In order to study the resulting dynamics, one must ascribe some properties to the time dependence of parameters in the models, for example, one might assume they are periodic.

The more realistic model incorporating the season-driven light intensity is represented by

where the carrying capacity $K(t)$ is a periodic function with period one year and peak at summer time. Assume $g(x)$ is a smooth function that satisfies

3.   Dynamics of system (2.6)-(2.7)

In this section, we explore the dynamics of (2.6)-(2.7) and present some results on the positive invariance, the boundary and internal dynamics. For convenience, we introduce the following notations

3.1.   Positive invariance and dissipativity

Theorem 3.1. The set defined by

is positively invariant with respect to (2.6)-(2.7).

Proof. To prove the positive invariance of $\Gamma$, we only need to examine the direction fields on the boundary of $\Gamma$. On the upper boundary of $\Gamma$, one has $n = T$, $x \in [0, k]$, $\frac{dn}{dt} = -dT<0$, where $k = \min\{K_{\max}, T/\bar{q}\}$. On the right boundary of $\Gamma$, $x = k$, $n \in (0, T)$, $\frac{dx}{dt} = rk(1-\frac{k}{\min\{K(t), n/\bar{q}\}})\leq0$ hold. Therefore, all orbits starting from $\Gamma$ cannot escape from these two boundaries.

On the bottom boundary $\{(x, n) : 0 < x < k, n = qx\}$ of $\Gamma$. We claim if $n(0)/x(0) > q$, then $n(t)/x(t) > q$ for all $t > 0$. Otherwise, there exists a $t_2 > 0$ such that $Q(t_2) = n(t_2)/x(t_2) = q$ and $Q(t) = n(t)/x(t) > q$ for $t \in [0, t_2)$. It is obvious that $\frac{dQ(t_2)}{dt}\leq0$. When $n(t)/\bar{q}\geq K_{\max}$, $\min\{K(t), n(t)/\bar{q}\} = K(t)$, one has

It follows that

which contradicts $\frac{dQ(t_2)}{dt}\leq0$.

It remains to show that all orbits starting from $\Gamma$ cannot leave from the left boundary of $\Gamma$, we assume $x(t_1) = 0$ and then

where $\alpha$ is a constant. Thus, $x(t_1)\geq x(0)\exp\{\alpha t_1\}>0$, which contradicts $x(t_1) = 0$. Therefore, no trajectory can touch the left boundary of $\Gamma$. This completes the proof.

Theorem 3.2. (2.6)-(2.7) is ultimately bounded with respect to $\Gamma\cup\partial\Gamma$, i.e. (2.6)-(2.7) is dissipative.

Proof. Let $(x(t), n(t))$ be a solution of (2.6)-(2.7) through $x(t_0)>0$ and $n(t_0)>0$. It is not difficult to verify that $\liminf _{t\rightarrow\infty}x(t)\geq0, ~\liminf _{t\rightarrow\infty}n(t)\geq0, \limsup _{t\rightarrow\infty}x(t)\leq k, ~\limsup _{t\rightarrow\infty}n(t)\leq T$.

From the proof of Theorem 3.1, we have

It follows that $\liminf_{{t\rightarrow\infty}} Q(t) = \liminf_{{t\rightarrow\infty}} \frac{n(t)}{x(t)}\geq q$. Therefore, (2.6)-(2.7) is ultimately bounded with respect to $\Gamma\cup\partial\Gamma$, i.e. (2.6)-(2.7) is dissipative. The proof is complete.

3.2.   Boundary dynamics

In this subsection, we explore the boundary dynamics of (2.6)-(2.7) and prove that all orbits tend to the origin $E_0 = (0, 0)$, i.e., algae go extinct. Actually, the origin is not a steady state but plays a similar role as a steady state.

Theorem 3.3. Assume $g(T)\leq d\bar{q}$ and $T<\bar{q}K_{\min}$, all solutions in $\Omega$ tend to the origin $E_0$, i.e., $E_0$ is globally asymptotically stable.

Proof. we introduce the transformation

which converts the system (2.6)-(2.7) in $\Omega$ into the new system

Here,

This system has two equilibria $(0, 0)$ and $(u_0, 0)$ with $u_0 = \frac{d+r}{\bar{q}r+g(T)}$. Note that $g(T)\leq d\bar{q}$, then $u_0<1/\bar{q}<1/q$, and the equilibrium $(u_0, 0)$ lies on the left boundary of $\Phi(\Omega)$.

The $u$-nullcline is $l_1$ : $u = 0$ and $u = u_1(n) = \frac{d+r}{\bar{q}r+g(T-n)}, ~0 \leq n < T$. The p-nullcline is $l_2$ : $n = 0$ and $u = u_2(n) = \frac{d}{g(T-n)}, 0\leq n<T .$ Since $1/q >\frac{d}{g(T-n)}>\frac{d+r}{\bar{q}r+g(T-n)}$, $l_2$ is above $l_1$.

Define the regions

In $D_1$, $\frac{du}{dt}>0, ~\frac{dn}{dt}<0$. In $D_2$, $\frac{du}{dt}<0, ~\frac{dn}{dt}<0$. In $D_3$, $\frac{du}{dt}<0, ~\frac{dn}{dt}>0$. Thus, any solution starting from the region $D_1$ tends to the equilibrium $(u_0, 0)$. Any solution starting from the region $D_3$ first enters the region $D_2$, and then either directly tends to the equilibrium $(u_0, 0)$ or passes through the region $D_1$ and then tends to the equilibrium $(u_0, 0)$. Thus, any solution $(u(t), n(t))$ of the system (3.2)-(3.3) tends to the equilibrium $(u_0, 0)$, i.e.,

Thus,

which implies the origin of the system (2.6)-(2.7) is globally asymptotically stable.

3.3.   Internal dynamics

In this subsection, we investigate the internal dynamics of (2.6)-(2.7) by establishing some sufficient criteria for the existence and GAS of internal equilibrium and positive periodic solution of (2.6)-(2.7).

Theorem 3.4. When $g(T)>d\bar{q}$ and $T<\bar{q}K_{\min}$, there exists a unique internal equilibrium $E_1 = (\bar{x}, \bar{n})$ satisfying $\bar{n} = \bar{q}\bar{x}$ and $g(T-\bar{n}) = d\bar{q}$. Furthermore, $E_1$ is globally asymptotically stable.

The proof of this theorem is similar to the proof of Theorem 3. in ^[3], so we omit it here.

Next we explore the existence and GAS of positive periodic solution of (2.6)-(2.7). For the reader's convenience, we shall summarize below a few concepts and results from ^[8] that will be essential for the following discussion.

Let $X, Z$ be normed vector spaces, $L:{\rm Dom}~ L\subset X\rightarrow Z$ be a linear mapping, and $N:X\rightarrow Z$ be a continuous mapping. The mapping $L$ will be called a Fredholm mapping of index zero if $\dim~{\rm Ker}~L = {\rm codim}~{\rm Im}~L<+\infty$ and ${\rm Im}~L$ is closed in $Z$. If $L$ is a Fredholm mapping of index zero, then there exist continuous projectors $P:X\rightarrow X$ and $Q:Z\rightarrow Z$ such that ${\rm Im}~P = {\rm Ker}~L$, ${\rm Im}~L = {\rm Ker}~Q = {\rm Im}(I-Q)$. It follows that $L|{\rm Dom} L\cap~{\rm Ker}~P:(I-P)X\rightarrow~{\rm Im}~L$ is invertible. We denote the inverse of that map by $K_P$. If $\Omega$ is an open bounded subset of $X$, the mapping $N$ will be called $L$-compact on $\overline{\Omega}$ if $QN(\overline{\Omega})$ is bounded and $K_P(I-Q)N:\overline{\Omega}\rightarrow X$ is compact. Since ${\rm Im}~{Q}$ is isomorphic to ${\rm Ker}~L$, there exist isomorphisms $J:{\rm Im}~Q\rightarrow~{\rm Ker}~L.$

Lemma 3.5 (Continuation theorem ^[8]). Let $L$ be a Fredholm mapping of index zero and let $N$ be L-compact on $\overline{\Omega}$. Suppose

(a) for each $\lambda\in(0, 1)$, every solution $x$ of $Lx = \lambda Nx$ is such that $x\notin \partial\Omega$;

(b) $QNx\neq 0$ for each $x\in\partial\Omega\cap \text{Ker} L$ and $\deg \{JQN, \Omega\cap \text{Ker} L, 0 \} \neq 0.$ Then the equation $Lx = Nx$ has at least one solution lying in ${\rm Dom} L\cap\overline{\Omega}.$

Theorem 3.6. Assume that $g(T)>2d\bar{q}$, then (2.6)-(2.7) has at least one strictly positive $\omega$-periodic solution.

Proof. Considering the biological significance of (2.6)-(2.7), we specify $x(0)>0$ and $n(0)>0$. Let $x(t) = \exp\{y_1(t)\}$ and $n(t) = \exp\{y_2(t)\}$, then (2.6)-(2.7) turns into

Define

and $\|y\| = \|(y_1, y_2)^T\| = \max_{t\in[0, \omega]}|y_1(t)|+\max_{t\in[0, \omega]}|y_2(t)|$ for any $y\in X(\text{or}~Z)$. Then $X$ and $Z$ are Banach spaces endowed with the norm $\|\cdot\|$. Let

Then, it is not difficult to show that ${\rm Ker}~L = R^2, ~{\rm Im}~L = \{z\in Z:\int_0^{\omega}z(t)dt = 0\}$ is closed in $Z$, $\dim~{\rm Ker}~L = 2 = {\rm codim}~{\rm Im}~L$, and $P$, $Q$ are continuous projectors such that ${\rm Im}~P = {\rm Ker}~L$, ${\rm Ker}~Q = {\rm Im}~L = {\rm Im}(I-Q)$. Therefore, $L$ is a Fredholm mapping of index zero. Furthermore, the generalized inverse (to $L$)$K_P:{\rm Im}~L\rightarrow~{\rm Ker}~P\cap~{\rm Dom} L$ is given by $K_P(z) = \int_0^t{z(s)}ds-1/\omega\int_0^\omega\int_0^tz(s)dsdt.$ By the Arzela-Ascoli theorem, it is not difficult to show that $N$ is $L$-compact on $\overline{\Omega}$ with any open bounded set $\Omega\in X$.

Now we are right here to search for an appropriate open bounded subset $\Omega$ for the application of the continuation theorem. Corresponding to the operator equation $Ly = \lambda Ny$, $\lambda\in(0, 1)$, one has

Let $y = y(t)\in X$ be a solution of (3.6) for some certain $\lambda\in(0, 1)$. Integrating (3.6) from $0$ to $\omega$ produces

then

From (3.6)-(3.10), it follows that

and

Since $y(t) = (y_1(t), y_2(t))^T\in X$, there exist $\xi_i, \eta_i\in[0, \omega]$ such that

By (3.9), one has

that is, $y_1(\xi_1)\leq\ln[K_{\max}],$ then

Let $F(x) = \frac{g(T-x)}{x}$, it is easy to verify that $F(x)$ is decreasing when $x>0$. By (3.10), one has

then

Since $F(0) = \infty>d\exp\{-K_1\}, ~F(T) = 0<d\exp\{-K_1\}$, therefore $F^{-1}(d\exp\{-K_1\})$ exists within $(0, T)$. Then one can verify that $y_2(\xi_2)\leq F^{-1}(d\exp\{-K_1\}),$ then

Moreover, when $\exp\{y_2(t)\}\geq2\bar{q}\exp\{y_1(t)\}$, it follows from (3.9) that

which implies

When $\exp\{y_2(t)\}<2\bar{q}\exp\{y_1(t)\}$, it follows from (3.10) that

which implies

Thus,

Whence,

Therefore, $\max_{t\in[0, \omega]}|y_1(t)|\leq\max \{|K_1|, |K_3|\}: = K_5.$

On the other hand, by (3.10),

It is not difficult to show $F^{-1}(d\exp\{-K_3\})$ exists within $(0, T)$. Thus,

then,

Therefore, $\max_{t\in[0, \omega]}|y_2(t)|\leq\max \{|K_2|, |K_4|\}: = K_6.$ Clearly, $K_i(i = 1, ..., 6)$ are independent of $\lambda$.

Considering the following algebraic equations

Under the assumption of $g(T)>2d\bar{q}>d\bar{q}$, we will divide the discussion into the following three cases:

● If $T\leq\bar{q}K_{\min}$, there exists a unique solution $(v_1^*, v_2^*)$ satisfying $0<v_1^*\leq K_{\min}, ~0<v_2^*\leq \bar{q}K_{\min}.$

● If $T>\bar{q}K_{\min}$ and $F(\bar{q}K_{\max})<\tilde{K}$, where $\tilde{K} = \frac{1}{\omega}\int_0^\omega\frac{1}{K(t)}dt$, there exist at least one solution. We denote the largest one as $(v_1', v_2')$, which satisfies $K_{\min}\leq v_1'\leq 1/\tilde{K}, ~\bar{q}K_{\min}\leq v_2'< T$.

● If $T>\bar{q}K_{\min}$ and $F(\bar{q}K_{\max})>\tilde{K}$, there exists a largest solution $(v_1^{**}, v_2^{**})$ satisfying $v_1^{**} = \frac{1}{\tilde{K}}, ~\bar{q}K_{\max}\leq v_2^{**}< T$.

Hence, the solutions of the system of algebraic equations $(\tilde{v}_1, \tilde{v}_2)$, satisfy $v_1^*\leq\tilde{v}_1\leq v_1^{**}$, $v_2^*\leq\tilde{v}_2\leq v_2^{**}$, thus there exists a sufficiently large number $K_7$ such that $\|(\ln(\tilde{v}_1), \ln(\tilde{v}_2))^T\| = |\ln(\tilde{v}_1)|+|\ln(\tilde{v}_2)|\leq \max\{|\ln v_1^*|, |\ln v_1^{**}|\}+\max\{|\ln v_2^*|, |\ln v_2^{**}|\}\leq K_7.$ Denote $K = K_5+K_6+K_7$ and take $\Omega = \{ y(t) = (y_1(t), y_2(t))^T\in X:\|y\|<K \},$ then $\Omega$ verifies the requirement (a) in Lemma 3.5. When $x\in\partial\Omega\cap~{\rm Ker}~L = \partial\Omega\cap R^2$, $x$ is a constant vector in $R^2$ with $\|y\| = K$. Then

Define $H_\nu(y_1, y_2) = \nu QN(y_1, y_2)+(1-\nu)G(y_1, y_2),$ where $\nu\in[0, 1]$ and $G(y_1, y_2) = (1-\bar{q}\exp\{y_1\}/\exp\{y_2\}, ~d\exp\{y_2\}-g(T-\exp\{y_2\})\exp\{y_1\})^T$. Obviously, $H_\nu(y_1, y_2)$ is homotopic mapping, here $H_0(y_1, y_2) = QN(y_1, y_2), ~H_1(y_1, y_2) = G(y_1, y_2)$ and $0\notin H_\nu(\partial\Omega\cap~{\rm Ker}~L).$ Furthermore, in view of the assumptions in Theorem 3.6, direct calculations show that $\deg\{JQNy, \Omega\cap~{\rm Ker}~L, 0\} = \deg\{G, \Omega\cap~{\rm Ker}~L, 0\}\neq0$ when $g(T)>d\bar{q}$. Here $J$ can be the identity mapping since ${\rm Im}~P = {\rm Ker}~L$. By now we have proved that $\Omega$ verifies all the requirements in Lemma 3.5. Hence, (3.4) has at least one solution $y^*(t) = (y_1^*(t), y_2^*(t))^T$ in ${\rm Dom} L\cap\overline{\Omega}$. Set $x^*(t) = \exp\{y_1^*(t)\}$, $n^*(t) = \exp\{y_2^*(t)\}$, then $(x^*(t), n^*(t))^T$ is an $\omega$-periodic solution of system (2.6)-(2.7) with strictly positive components. Thus, there exists an $\omega$-periodic solution of system (2.6)-(2.7). The proof of Theorem 3.6 is complete.

Lemma 3.7. Let $h$ be a real number and $f$ be a nonnegative function defined on $[h, +\infty)$ such that $f$ is integrable on $[h, +\infty)$ and is uniformly continuous on $[h, +\infty)$. Then $\lim_{t\rightarrow +\infty}f(t) = 0$.

Theorem 3.8. Let $(x^*(t), n^*(t))$ be the positive periodic solution of (2.6)-(2.7). If

then $(x^*(t), n^*(t))$ is globally asymptotically stable.

Proof. Let $(x(t), n(t))$ be any solution of (2.6)-(2.7) with positive initial value. Since $\Gamma\cup\partial\Gamma$ is an ultimately bounded region of (2.6)-(2.7), there exists $T_1>0$ such that $(x(t), n(t))\in \Gamma\cup\partial\Gamma$ and $(x^*(t), n^*(t)) \in \Gamma\cup\partial\Gamma$ for $t \geq t_0+T_1$. Consider the Liapunov function defined by

A direct calculation of the right derivative of $V(t)$ along the solution of (2.6)-(2.7) leads to

where $\theta$ is a positive constant. Integrating both sides of above inequality from $t_0+T_1$ to $t$ produces

It follows from (3.11)

Therefore, $|x(s)-x^*(s)|+|n(s)-n^*(s)|\in L^1([t_0+T_1, +\infty))$.

The boundedness of $x^*(t)$ and $n^*(t)$ and the ultimate boundedness of $x(t)$ and $n(t)$ imply that $x(t), n(t), x^*(t), n^*(t)$ all have bounded derivatives for $t>t_0+T_1$ (from the equations satisfied by them). Then it follows that $|x(t)-x^*(t)|+|n(t)-n^*(t)|$ is uniformly continuous on $[t_0+T_1, +\infty)$. By Lemma 3.7, one has

The proof is complete.

In summary, we have established some sufficient criteria for the dynamics of (2.6)-(2.7). Figure 2 shows such scenarios, where the parameter values are deliberately specified such that Theorems 3.3, 3.4, 3.6 and 3.8 are satisfied. Figure 2(a) and (b) show that the algal density tends to 0 and a positive steady state under the conditions of Theorems 3.3 and 3.4, respectively. In Figure 2(c), there exists a positive periodic solution that is globally asymptotically stable when Theorems 3.6 and 3.8 hold. The solution curves in Figure 2(c) converges rapidly, and these solutions in the early time window are exhibited in Figure 2(d). Of course, these results can be improved. For instance, Figure 2(e) and (f) illustrates that (2.6)-(2.7) still admits a positive periodic solution that is globally asymptotically stable, even when the conditions in Theorem 3.8 are not satisfied.

4.   Numerical analysis and applications

In this section, based on the above theoretical analysis, numerical simulations show that the ecological impact and consequence of the seasonal light intensity and the nutrient availability on the model variables. Here we choose

Figure 3(a) shows the season-dependent function $K(t)$ that reaches maximum in summer time and minimum in winter time. Figure 3(b) illustrates the solution of the nonautonomous system (2.6)-(2.7) with $T = 0.03, ~0.08, ~0.11, ~0.3$ and other parameters chosen in Table 1. The seasonal carrying capacity leads to periodic fluctuations of algae.

In summer time, algal density is over the algal bloom level, and during this period, algal nutrient content reaches the maximum while N:C ratio of algae reaches the minimum, due to the fact that the speed of algal growth is faster than that of the algal nutrient content. Algal density reaches the minimum in winter time, during which algal nutrient content is lowest and N:C ratio of algae is highest. When there exists an cyclic oscillation ($T = 0.08, ~0.11, ~0.3$), Figure 3(b) also shows that the maximum of the algal density is positively related to the nutrient availability ($T$), while the minimum remains unchanged with varying $T$. In other words, $T$ may influence the peak value of the algal density, thus decreasing $T$ in summer time can be an effective measure to reduce the algal density under the algal bloom level. When $T = 0.03$, the solutions system (2.6)-(2.7) tend to a positive steady state with a low level. In order to further understand the continuous variation of the algal density with respect to $T$, we draw a bifurcation diagram (Figure 4(a)). One can see that there exists a unique internal equilibrium at the beginning, and then a positive periodic solution appears. The minimum of the cyclic oscillation keeps constant with respect to $T$, while the maximum increases as $T$ increases until a certain value of $T_1$ and stays unchanged beyond this value. It reveals that algal growth is restricted not only by light intensity and nutrient availability but also its own stoichiometry when the nutrient availability beyond $T_1$. The amplitude of the cyclic oscillation of algae increases first and stay constant with the increase in the strength of cyclic forcing of $T$ as shown in Figure 4(b).

There is another question: what is the threshold of $T$ for the occurrence of algae bloom? We consider the autonomous system (2.3)-(2.4) and the nonautonomous system (2.6)-(2.7) to give the thresholds of $T$, respectively. As mentioned above, one can compute $T_b^1 = 0.1003$ from (2.5) for system (2.3)-(2.4). For system (2.6)-(2.7), we find the threshold $T_b^2$ through

By simple calculations, $T_b^2 = 0.1010>T_b^1$. That is to say, if we consider the seasonal changing light intensity, there is a higher threshold than the autonomous system (2.3)-(2.4), but the difference is relatively small.

There is another question: what is the threshold of $T$ for the occurrence of algae bloom? We consider the autonomous system (2.3)-(2.4) and the nonautonomous system (2.6)-(2.7) to give the thresholds of $T$, respectively. As mentioned above, one can compute $T_b^1 = 0.1003$ from (2.5) for system (2.3)-(2.4). For system (2.6)-(2.7), we find the threshold $T_b^2$ through

By simple calculations, $T_b^2 = 0.1010>T_b^1$. That is to say, if we consider the seasonal changing light intensity, there is a higher threshold than the autonomous system (2.3)-(2.4), but the difference is relatively small.

5.   Control of algal bloom

The control method of the Chinese government is to remove algae (RA) periodically. Here we propose the other control method: to block nutrient (BN) input from rivers constantly. For the RA method, the model becomes

where $\delta(t)$ is a Haar (step) function, which is a positive constant in summer time (July and August) and zero in other seasons. Figure 5 shows that the algal bloom can never be controlled by removing algae periodically. Every year algal bloom occurs no matter how high summer removal rate is. Hence, the RA method is not an effective way to inhibit algal bloom in Yellow Sea off China. For the BN method, the model becomes

where $\sigma$ is the percentage of nutrient block rate from input rivers. Figure 6 shows that the algal bloom can potentially be controlled by reducing N concentration. The reduction of N concentration leads to lower peaks in summer time. Algal bloom clearly disappears when N is reduced by 70% or a higher percentage. Sufficiently N reduction excludes fluctuations, and algal density will almost stay in a very low constant (steady state) level, which is good for the control of algal bloom. Obviously the BN method is an effective way to terminate algal bloom in Yellow Sea off China.

6.   Application to Bohai Sea

Primary production in Bohai Sea exhibited an annual cycle with a large peak occurring in May or June. In the recent few years the cycle only had one peak, while during the years 1998-2006 the annual cycle had two clear peaks with a large one occurring in May or June and a small one occurring in August, September or October ^[5]. The small peak almost disappeared in the recent years, thus it is not as robust as the large peak. For this reason, we fit 2007 primary production data of Bohai Sea that clearly has one peak (see the data in Table 2)

We assume primary production is proportional to algal density and justify the light-dependent carrying capacity accordingly for Bohai Sea. We apply our model to fit the Bohai Sea data in Figure 7 using the least square method. The simulation time series fits the real data well with tiny mean relative error about 0.094.

Because the nitrogen concentration in Bohai Sea is usually increased by industrial pollution, we suggest the Chinese government to control pollution from industry to rivers, such as Yellow River, merging to Bohai Sea, and to control direct pollution in Bohai Sea such as production of oil and gas.

Acknowledgment

This research is partially supported by NSFC-11671072, NSFC-11801052, and NSERC grant RGPIN-2015-04581.

Conflict of interest

All authors declare no conflicts of interest in this paper.