Mathematical modeling and analysis of harmful algal blooms in flowing habitats

1.   Introduction

In the past decades, harmful algal blooms (HAB) have become important water quality issues in both coastal and inland waters, and the frequency and intensity of HAB are apparently increasing worldwide. Blooms of the haptophyte algae Prymnesium parvum have become more and more common in the world, and it is referred to as golden algae ^[1,2], which were documented to cause large fish kills ^[3,4]. Recent studies suggest possible potential techniques for managing and mitigating harmful algal blooms through flow manipulations in some river systems ^[5,6,7,8]. This motivates the theoretical exploration of harmful algal dynamics in riverine reservoirs. To understand longitudinal patterns arising along the axis of flow, advection-dispersion-reaction systems were employed to study the effects of spatial variations of harmful algae and its toxin production and decay in riverine reservoirs ^{[9,10,11,12,13]}. The models are one-dimensional systems with simple habitat geometry and transport processes. The flow reactor model with transport of nutrient and organisms by both advection and diffusion was first proposed in ^[14]. Recently, the flow reactor model in ^[14] has been further incorporated with a hydraulic storage zone for persistence and coexistence of competing populations ^[9]. Such systems with/without a hydraulic storage zone become more and more attractive since they can be regarded as a surrogate model for riverine reservoirs with strong advective flows ^[10].

It should be pointed out that the flow reactor system presented in ^[14] was used to model the influences of bacterial motility, fluid advection and other spatial variations on the competition between different strains of bacteria for the limiting nutrient in the large intestine of animals. Differently, our main purpose here is to use the flow reactor system to describe the dynamics/interactions of harmful algae and nutrient(s) in the river/stream. In ^[15], the author extended the model in ^[14] by considering two species competition for two essential/complementary nutrients, such as nitrogen and phosphorus. The complementary nutrient model is highly relevant since the limiting nutrient(s) in many ecosystems should be multiple, and hence, the single-nutrient model can be seen as a special case. Since the environment of the plug-flow reactor is the intestine or a river, it is much more realistic to assume that the input nutrient concentration is time-dependent. Thus, the periodically varying input nutrient concentration is incorporated into the model of ^[14] and the model of ^[15] in ^[16] and ^[17], respectively.

There is a persistence paradox in the river ecology, namely, rapid advective flow in such habitats can prevent persistence of one species for realistic parameters. This motivates us to incorporate the factor of hydraulic storage zones in flowing water habitats ^[9,11] since it might resolve this persistence paradox ^[18]. Introducing a hydraulic storage zone into the flow reactor model not only makes sense biologically but also causes mathematical challenges. Some equations in the flow reactor with a hydraulic storage zone have no diffusion terms, and hence, the associated solution maps are not compact. In order to obtain the existence of global attractor, we show that solution maps are asymptotically smooth by using the Kuratowski measure of noncompactness. Note that the existence of global attractor is assumed in the theory of uniform persistence and coexistence states (see, e.g., ^[19]).

Another problem is about the local stability of the trivial and semi-trivial solutions of the model, which are usually determined by the sign of the principal eigenvalue(s) of the associated linearized system at these states. Although the associated linearized system is cooperative, the "compactness" is required when one uses the classical Krein–Rutman theory to obtain the existence of the principal eigenvalue. Very recently, Wang and Zhao developed some sharp abstract results (see [20,Theorem 2.3] and [20,Remark 2.2]) on the existence of principal eigenvalues for an elliptic eigenvalue problem with some zero diffusion coefficients. Two closely related applications can be found in [12,Lemma 3.3] and [21,Theorem 2.1]. The authors in ^[22] further studied the existence of the principal eigenvalue for degenerate periodic reaction-diffusion systems (see [22,Theorem 2.16 and Remark 2.21]). An alternative approach is to directly utilize the generalized Krein–Rutman Theorem developed by Nussbaum ^[23], see, e.g., [12,Lemma 4.4].

The rest of the paper is organized as follows. In section 2, a flow reactor model proposed by Kung and Baltzis ^[14] and its extensions are briefly reviewed. Section 3 is devoted to the survey of a model of a flowing water habitat with a hydraulic storage zone in which no diffusive and advective flow occurs. The input nutrient concentration can be a constant ^[9] or time-dependent ^[11]. In section 4, we review a model of interactions of a single limiting nutrient, harmful algae, toxins, and zooplankton ^[10,12]. Coexistence of harmful algae and zooplankton was also investigated in ^[12]. We further discuss a model of algal toxins and nutrient recycling (see ^[10] and ^[13]), which is based on the fact that many cyanotoxins produced by cyanobacteria can get recycled back into the system as a nutrient resource after they decompose. Finally, a brief discussion section completes the paper.

2.   The flow reactor model

We first review a model of microbial competition for a single limited nutrient in a riverine reservoir occupying a simple channel of longitudinally invariant cross-section, which was formulated by Kung and Baltzis in ^[14] and analyzed by the authors in ^[24,25]. The channel is assumed to have constant cross-sectional area $A$ and length $L$, yielding volume $V$. A flow of water enters at the upstream end ($x = 0$), with discharge $F$ (dimensions length$^3$ / time). An equal flow exits at the downstream end ($x = L$), which is assumed to be a dam. Based on this flow, a dilution rate $D$ (dimensions time$^{-1}$) is defined as $F / V$. The advective flow within the channel is set to maintain water balance, by transporting water with a net velocity $\nu = DL$.

The reactor occupies the portion of the channel from $x = 0$ to $x = L$ in which the microbial populations $N_i$, $i = 1, 2$, compete for nutrient $R$. A flow of medium in the channel with velocity $\nu$ in the direction of increasing $x$ brings fresh nutrient at a constant concentration $R^{(0)}$ into the reactor at $x = 0$ and carries medium, unutilized nutrient and organisms out of the reactor at $x = L$. Nutrient and organisms are assumed to diffuse throughout the vessel with the same diffusivity $\delta$.

Both advective and diffusive transport occur at the upstream boundary ($x = 0$). The inflow contains dissolved nutrient $R(x, t)$ at a concentration $R^{(0)}$. The downstream boundary is assumed to be a dam, over which there is advective flow but through which no diffusion can take place. These assumptions lead to boundary conditions for nutrient:

Boundary conditions for the nutrient are given by equations (2.1), while those for population densities are

These boundary conditions mean that no inflow of the populations occurs, and there is no dispersive transport over the dam at the downstream end. Under these assumptions, we have the following reaction-diffusion system describing the densities $R(x, t)$, $N_1(x, t)$ and $N_2(x, t)$:

with boundary conditions

and initial conditions

where $q_i$ is the constant nutrient quota for species $i$. The nonlinear functions $f_{i}(R)$ describes the nutrient uptake rate and the growth rate of the organisms $N_{i}$ at nutrient concentration $R$. We assume that these functions satisfy

The usual example is the Monod function

where $\mu_\text{max, i}$ (resp. $K_{i}$) represents the maximal growth rate (resp. the half saturation constant) of species $i$. In [25,Chapter 8], the author showed that both species in (2.2)-(2.4) can coexist under suitable conditions. The condition (2.3) is called as the Danckwerts' boundary condition by Aris ^[26]. For a detailed derivation of it, we refer to a review paper ^[27].

The authors in ^[24] extended (2.2)-(2.4) to the following system

with boundary conditions (2.3), and initial conditions (2.4). Here $\delta_0$ and $\delta_i$ stand for the random motility coefficients of nutrient and species $i$, respectively; $m_i$ is the death rate of species $i$. The effects of random motility on the extinction/persistence of a single population model, and the influences of random motility on the competition outcomes between two species were investigated in ^[24]. Note that if we assume $\delta_0 = \delta_1 = \delta_2 = \delta$ and $m_1 = m_2 = 0$ in (2.5), then it becomes (2.2). The authors in ^[16] further incorporated a periodically varying input nutrient concentration into system (2.5) with boundary conditions (2.3) and initial conditions (2.4), where the input concentration $R^{(0)}$ is replaced by a $\tau$-periodic function $R^{(0)}(t)$. Then they used the theories of monotone dynamical systems and uniform persistence to obtain some analytic results about the extinction/persistence of a single population model and coexistence of two species system in terms of the principal eigenvalue(s) of the associated periodic-parabolic eigenvalue problem(s).

To address the multiple nutrients in ecosystems, the author in ^[15] also generalized model (2.2)-(2.4) to an evolution system of two species competition for two essential nutrients with constant input concentrations. Later on, a model of two species competition for two essential nutrients with periodically varying input concentrations was studied in ^[17].

3.   Models with storage zones

This section is devoted to the survey of models with hydraulic storage zones, which partition the cross-section of the channel into a flowing zone of area $A$, and a static zone of area $A_S$. Exchange of nutrient and populations between the flowing and storage zones occurs by Fickian diffusion with rate $\alpha$. Nutrient concentration and population densities vary with location in both the flowing channel and the storage zone, however, we assume that advective and diffusive transport occur only in the flowing zone, not the storage zone. Nutrient concentration and population densities in the flowing channel (resp. the storage zone) are denoted by $R(x, t)$ and $N_i(x, t)$ (resp. $R_{S}(x, t)$, $N_{S, i}(x, t)$). Then we generalize system (2.2)-(2.4), by adding a hydraulic storage zone and including the seasonality of $R^{(0)}$, to the following form ^[9,11]:

with boundary conditions

and initial conditions

Here $R^{(0)}(t)$ satisfies

We should point out that the authors in ^[9] considered system (3.1)-(3.3) under the assumption that $R^{(0)}(t)\equiv R^{(0)}$ is a positive constant.

By [33,Theorem 1 and Remark 1.1], we have the following results.

Lemma 3.1. ([11,LEMMA 2.1]) System (3.1)-(3.3) has a unique noncontinuable solution and solutions of (3.1)-(3.3) remain non-negative on their interval of existence if they are non-negative initially.

In the following, we demonstrate that system (3.1)-(3.3) have a mass conservation in the flow and storage zones. Let

Then $W(x, t)$ and $W_{S}(x, t)$ satisfy the following system

The following result is concerned with the global dynamics of system (3.6).

Lemma 3.2. ([11,LEMMA 2.3]) System (3.6) admits a unique positive $\tau$-periodic solution $(W^{*}(x, t), W^{*}_{S}(x, t))$ and for any $(W_{0}, W_{0}^S)\in C([0, L], \mathbb{R}^2)$, the unique mild solution $(W(x, t), W_{S}(x, t))$ of (3.6) with $(W(x, 0), W_{S}(x, 0)) = (W_{0}(x), W_{0}^S(x))$ satisfies $\lim\limits _{t \rightarrow \infty}((W(x, t), W_{S}(x, t))-(W^{*}(x, t), W^{*}_{S}(x, t))) = (0, 0)$ uniformly for $x \in [0, L]$.

We should point out that [11,Lemma 2.2] is based on the assumption that the associated eigenvalue problem admits a principal eigenvalue, and there is a gap in the arguments for the existence of the principal eigenvalue in the paragraph above [11,Lemma 2.2]. However, this gap can be easily filled by using [21,Theorem 2.1] or the arguments similar to those in [12,Lemma 3.3] combined with [20,Theorem 2.3 and Remark 2.2].

By Lemma 3.1, the relation (3.5) and Lemma 3.2, we have the following result.

Lemma 3.3. ([11,LEMMA 2.4]) Any solution of the system (3.1)-(3.3) exists globally on $[0, \infty)$. Moreover, solutions are ultimately bounded and uniformly bounded.

3.1.   Single species growth

This subsection is devoted to the investigation of the single population model. Mathematically, it means that we set $(N_1, N_{S, 1}) = (0, 0)$ or $(N_2, N_{S, 2}) = (0, 0)$ in the model system (3.1)-(3.3). In order to simplify notation, we drop all subscripts in the remaining equations and then consider

with boundary conditions

and initial conditions

where $R^{(0)}(t)$ satisfies (3.4).

Let

Then $\mathcal{W}(x, t)$ and $\mathcal{W}_{S}(x, t)$ satisfy (3.6). By Lemma 3.2, we see that the limiting system of (3.7)-(3.9) takes the following form:

with boundary conditions

and initial conditions

From the biological view of point, the feasible domain $\Lambda(t)$ for (3.10)-(3.12) should be

We further have the following basic properties of the set $\Lambda(t)$.

Lemma 3.4. ([11,LEMMA 3.1]) For any $\phi: = (\phi_1, \phi_2)\in \Lambda(0)$, system (3.10)-(3.12) has a unique mild solution $(N(\cdot, t), N_{S}(\cdot, t))$ with $(N(\cdot, 0), N_{S}(\cdot, 0)) = \phi$ and $(N(\cdot, t), N_{S}(\cdot, t))\in \Lambda(t)$, for all $t\geq 0$.

By Lemma 3.4, we can define solution maps $\Psi_{t}:\Lambda(0)\rightarrow \Lambda(t)$ associated with (3.10)-(3.12) by

Note that $\Psi_{\tau}:\Lambda(0)\rightarrow \Lambda(\tau) = \Lambda(0)$ is the Poincaré map associated with (3.10)-(3.12).

For convenience, we let

Since one equation in (3.10)-(3.12) has no diffusion term, its solution map $\Psi_{t}$ is not compact. Due to the lack of compactness, we need to impose the following condition:

Recall that the Kuratowski measure of noncompactness (see, e.g., ^[28]), $\kappa$, is defined by

for any bounded set $B$. We set $\kappa(B) = \infty$ whenever $B$ is unbounded. Note that $B$ is precompact(i.e., $\bar{B}$ is compact) if and only if $\kappa(B) = 0$.

Lemma 3.5. ([11,LEMMA 3.2]) Let (3.13) hold. Then $\Psi_{\tau}$ is $\kappa$-contracting in the sense that $\lim_{n\to \infty}\kappa(\Psi_{\tau}^nB) = 0$ for any bounded set $B \subset Y^{+}$.

By Lemma 3.3, Lemma 3.5 and [29,Theorem 2.6], we have the following result.

Theorem 3.1. ([11,THEOREM 3.1]) $\Psi_{\tau}$ admits a global attractor on $Y^{+}$ that attracts each bounded set in $Y^{+}$ provided that (3.13) holds.

Note that $(0, 0)$ is the trivial solution of (3.10)-(3.12). Linearizing system (3.10)-(3.12) at $(0, 0)$, we have

Substituting $N(x, t) = e^{-\mu t}\phi_1(x, t)$ and $N_{S}(x, t) = e^{-\mu t}\phi_2(x, t)$, we obtain the associated eigenvalue problem

As in the proof of [11,Lemma 3.3], we let $\Pi_t:C([0, L], \mathbb{R}^2)\rightarrow C([0, L], \mathbb{R}^2)$ be the solution maps associated with (3.15). Then $\mathcal{P}: = \Pi_\tau$ is the Poincaré map associated with system (3.15). Let $r(\mathcal{P})$ be the spectral radius of $\mathcal{P}$. By the proof of [11,Lemma 3.3], we further see that $\Pi_\tau$ is an $\kappa$-contraction on $C([0, L], \mathbb{R}^2)$ in the sense that

for any bounded set $B$ in $C([0, L], \mathbb{R}^2)$, where $r_0$ is a positive number such that

By (3.17) and the arguments similar to those in [12,Lemma 4.4]) or [30,Lemma 3.1], we can use the generalized Krein-Rutman Theorem ^[23] to obtain the following result, which is a corrected version of [11,Lemma 3.3].

Lemma 3.6. Define $\mu^{*}: = -\frac{1}{\tau}\ln r(\mathcal{P})$ and let (3.13) hold. If $r(\mathcal{P})\geq1$, then $\mu^{*}$ is the principal eigenvalue of the eigenvalue problem (3.16) with a strongly positive eigenfunction $\phi^{*} = (\phi^{*}_1, \phi^{*}_2)\gg 0$.

The following result is concerned with the global dynamics of system (3.10)-(3.12).

Theorem 3.2. ([11,THEOREM 3.2]) Assume that (3.13) holds. Let $(N(x, t), N_{S}(x, t))$ be the solution of (3.10)-(3.12) with initial data $(N^0(\cdot), \ N_{S}^0(\cdot))\in Y^{+}$. Then the following statements are valid:

$\rm{(1)}$If $\mu^{*} > 0$, then $\lim\limits _{t \rightarrow \infty}|(N(x, t), N_{S}(x, t))| = 0$ uniformly for $x \in [0, L]$;

$\rm{(2)}$If $\mu^{*} < 0$, then (3.10)-(3.12) admit a unique positive $\tau$-periodic solution

$(N^{*}(x, t), N^{*}_{S}(x, t))$ and for any $(N^0(\cdot), N_{S}^0(\cdot))\in Y_{0}$, we have

By Theorem 3.2 and the theories of chain transitive sets (see, e.g., ^[19,31]), one can obtain a threshold type result on the global dynamics of the single population model (3.7)-(3.9).

3.2.   Two species competition

This subsection focuses on the investigation of the possibility of coexistence for system (3.1)-(3.3). In view of the relation (3.5) and Lemma 3.2, we see that the limiting systems of (3.1)-(3.3) take the forms:

in $(0, L)\times (0, \infty)$, with boundary conditions

and initial conditions

From the biological view of point, the feasible domain $D(t)$ for (3.18)-(3.20) should be

The following result indicates that $D(t)$ is positively invariant for the solution maps associated with (3.18)-(3.20).

Lemma 3.7. ([11,LEMMA 2.5]) For any $\phi: = (\phi_1, \phi_2, \phi_3, \phi_4)\in D(0)$, system (3.18)-(3.20) has a unique mild solution $(N_1(\cdot, t), N_{S, 1}(\cdot, t), N_2(\cdot, t), N_{S, 2}(\cdot, t))\in D(t)$, for all $t\geq 0$, whenever $(N_1(\cdot, 0), N_{S, 1}(\cdot, 0), N_2(\cdot, 0), N_{S, 2}(\cdot, 0)) = \phi$.

Since two equations in (3.18)-(3.20) have no diffusion terms, its solution maps are not compact. So we require the following conditions in this subsection:

Fix $i\in\{1, 2\}$, we consider the following linear system

Then the associated eigenvalue problem takes the form

Let $\mathcal{P}_i$, $i = 1, 2$, be the Poincaré map associated with system (3.22), and $r(\mathcal{P}_i)$ be the spectral radius of $\mathcal{P}_i$. By the same arguments as in Lemma 3.6, we have the following result.

Lemma 3.8. Define $\mu_i^{*}: = -\frac{1}{\tau}\ln r(\mathcal{P}_i)$ and let (3.21) hold. If $r(\mathcal{P}_i)\geq1$, then $\mu_i^{*}$ is the principal eigenvalue of the eigenvalue problem (3.23) with a strongly positive eigenfunction.

Note that Then system (3.18)-(3.20) admits the following possible trivial/semi-trivial solutions:

$\rm{(i)}$ Trivial solution $\hat{0}: = (0, 0, 0, 0)$ always exists;

$\rm{(ii)}$ Semi-trivial solution $(N_{1}^{*}(x, t), N_{S, 1}^{*}(x, t), 0, 0)$ exists provided that $\mu_{1}^{*} < 0$;

$\rm{(iii)}$ Semi-trivial solution $(0, 0, N_{2}^{*}(x, t), N_{S, 2}^{*}(x, t))$ exists provided that $\mu_{2}^{*} < 0$;

$\rm{(iv)}$ There may be additional $\tau$-periodic solutions as well and these must be positive.

Here $(N_{i}^{*}(x, t), N_{S, i}^{*}(x, t))$ denotes the unique positive $\tau$-periodic solution of (3.10)-(3.12) resulting from putting $f = f_{i}$ and $q = q_{i}$. The two organisms can coexist if a positive $\tau$-periodic solution exists.

In view of Lemma 3.7, we let $\Phi_{t}:D(0)\rightarrow D(t)$ be the solution map of system (3.18)-(3.20). Let $K = C([0, L], \mathbb{R}_{+}^2)\times (-C([0, L], \mathbb{R}_{+}^2))$ and denote its induced order by $\leq_K$. Thus, the solution map $\Phi_t$ is monotone ^[25] with respect to the partial order $\leq_K$. Note that $\Phi_{\tau}:D(0)\rightarrow D(\tau) = D(0)$ and for the Poincaré map $S: = \Phi_{\tau}$, we have $S^{n}(P) = \Phi_{n \tau}(P)$, for all $n\in \mathbb{Z}$. Set $\mathbb{Y}^{+} = D(0)$,

and $\partial \mathbb{Y}_{0}: = \mathbb{Y}^{+}\backslash \mathbb{Y}_{0}$. For convenience, we further set

and

where ${\bf{u}}: = (u_1, u_2)\in \mathbb{R}_{+}^2$ and ${\bf{v}}: = (v_1, v_2)\in \mathbb{R}_{+}^2$. With the assumption (3.21), it follows whenever $\alpha\frac{A}{A_{S}}$ is sufficiently large, there exists a constant $r > 0$ such that

where ${\bf{g}}(t, x, {\bf{u}}, {\bf{v}}): = (g_{1}(t, x, u_1, u_2, v_1, v_2), g_{2}(t, x, u_1, u_2, v_1, v_2))$.

Lemma 3.9. ([11,LEMMA 4.1]) Let (3.21) and (3.24) hold. Then the map $\Phi_{\tau}$ is $\kappa$-contracting in the sense that $\lim_{n\to \infty}\kappa (\Phi_{\tau}^n(B)) = 0$ for any bounded set $B \subset \mathbb{Y}^{+}$, where $\kappa$ is the Kuratowski measure of noncompactness as defined in (3.14).

By Lemma 3.3, Lemma 3.9, and [29,Theorem 2.6], we have the following result.

Theorem 3.3. ([11,THEOREM 4.1]) Let (3.21) and (3.24) hold. Then $\Phi_{\tau}$ admits a global attractor on $\mathbb{Y}^{+}$ that attracts each bounded set in $\mathbb{Y}^{+}$.

Fix $i\in\{1, 2\}$, let $\hat{\mathcal{P}}_i$ be the Poincaré map associated with system (3.22) when $(f_i(W^{*}(x, t)), f_i(W^{*}_{S}(x, t)))$ in (3.22) is replaced by

Let $r(\hat{\mathcal{P}}_i)$ be the spectral radius of $\hat{\mathcal{P}}_i$. By the same arguments as in Lemma 3.6, we have the following observation.

Lemma 3.10. Define $\eta_{i}^{*}: = -\frac{1}{\tau}\ln r(\hat{\mathcal{P}}_i)$ and let (3.21) hold. If $r(\hat{\mathcal{P}}_i)\geq1$, then $\eta_{i}^{*}$ is the principal eigenvalue of (3.23) with $(f_i(W^{*}(x, t)), f_i(W^{*}_{S}(x, t)))$ replaced by the one in (3.25).

The following result is concerned with the coexistence of system (3.18)-(3.20),

Theorem 3.4. Let (3.21) and (3.24) hold, and assume that $\mu_{i}^{*} < 0$ and $\eta_{i}^{*} < 0, \ i = 1, 2$. Then system (3.18)-(3.20) admits at least one (componentwise) positive $\tau$-periodic solution and there exists a positive constant $\zeta > 0$ such that for any solution $(N_1(x, t), N_{S, 1}(x, t), N_2(x, t), N_{S, 2}(x, t))$ of system (3.18)-(3.20) with the initial data in $\mathbb{Y}_{0}$ satisfies $\liminf_{t\rightarrow \infty}\min_{x\in [0, L]}N_i(x, t)\geq \zeta$ and $\liminf_{t\rightarrow \infty}\min_{x\in [0, L]}N_{S, i}(x, t)\geq \zeta$, for all $i = 1, 2$.

We remark that Theorem 3.4 follows from [11,Theorem 4.2], where the theory of monotone dynamical systems have been used. Instead, we can also obtain Theorem 3.4 by using the theory of uniform persistence. In [11,Section 5], the authors further lifted the dynamics of the limiting system (3.18)-(3.20) to the full system (3.1)-(3.3) by the theory of chain transitive sets (see, e.g., ^[19,31]).

4.   Interactions of harmful algae and zooplankton

In this section, we survey systems modeling the interactions of nutrient, harmful algae, toxins, and zooplankton, in which the input concentration $R^{(0)}$ is always a constant.

4.1.   A model of harmful algae and their toxins

This subsection is devoted to the study of the influences of spatial variations on the growth of harmful algae and the production/decay of their toxins in riverine reservoirs. Suppose $R(x, t)$, $N(x, t)$ and $C(x, t)$ (resp. $R_{S}(x, t)$, $N_{S}(x, t)$ and $C_{S}(x, t)$) denote dissolved nutrient concentration, algal abundance and dissolved toxin concentration at location $x$ and time $t$ in the flowing channel, respectively (resp. in the storage zone). The authors in ^[10] propose the following advection-dispersion-reaction system:

in $(x, t)\in (0, L)\times (0, \infty)$ with boundary conditions

and initial conditions

in $x\in (0, L)$. Here the mortality of algae is assumed to be a constant rate $m$; $q_N$ represents the constant quota of algae. For simplicity, we have assumed that toxin degradation follows first order kinetics with a decay coefficient $k$. We point out that system (4.1)-(4.3) applies to many flagellate toxins ^[32].

There are two types of productions for dissolved toxins ^[10]. The first assumes that the algae produce toxin more rapidly when there is little nutrient in the system,

where $\epsilon$ is a constant coefficient and $\mu_\text{max}$ represents the maximal growth rate. It has been observed that toxins produced by Prymnesium parvum (toxic flagellates) are proportional to the degree of algal nutrient limitation. The second type of toxin production assumes that the toxin is produced proportional to the algal productivity,

This case assumes that toxin is produced in proportion to other cellular products and released into the water at a constant rate. We refer to this as the case of cylindrospermopsin, which is a cyanotoxin produced by a variety of freshwater cyanobacteria.

By [33,Theorem 1 and Remark 1.1], we have the following result.

Lemma 4.1. ([12,LEMMA 3.1]) System (4.1)-(4.3) has a unique noncontinuable solution and solutions of (4.1)-(4.3) remain non-negative on their interval of existence if they are non-negative initially.

In the following, we will demonstrate that mass conservation is satisfied in the flow and storage zones for the equations given by (4.1)-(4.3). Let

Then $W(x, t)$ and $W_{S}(x, t)$ satisfy the following system

Then one can show that (e.g., ^[9] and [11,Lemma 2.3]) system (4.4) admits a unique positive steady-state solution $(R^{(0)}, R^{(0)})$ and

It is not hard to see that $(R^{(0)}, 0, 0, R^{(0)}, 0, 0)$ is the trivial steady-state solution of (4.1)-(4.3). Linearizing system (4.1)-(4.3) around $(R^{(0)}, 0, 0, R^{(0)}, 0, 0)$, we get the following cooperative system for the algae population:

Substituting $N(x, t) = e^{\lambda t}\phi(x)$ and $N_{S}(x, t) = e^{\lambda t}\phi_S(x)$ into (4.5), we obtain the associated eigenvalue problem

Due to the noncompactness of the model system, we impose the following condition

By [20,Theorem 2.3] or [21,Theorem 2.1] (see also the arguments in [12,Lemma 3.3], it follows that the eigenvalue problem (4.6) has a principal eigenvalue, denoted by $\lambda^{0}$.

We are in a position to adopt the results developed in ^[20] to define the basic reproduction ratio for algae. Let $S(t):C([0, L], \mathbb{R}^2)\rightarrow C([0, L], \mathbb{R}^2)$ be the $C_0$-semigroup generated by the following system

Note that $S(t)$ is a positive $C_0$-semigroup on $C([0, L], \mathbb{R}^2)$. We further assume that both algae individuals in the flow and storage zones are near the trivial steady-state solution of (4.5), and introduce fertile individuals at time $t = 0$, where the distribution of initial algae individuals in the flow and storage zones is described by $\varphi: = (\varphi_2, \varphi_5)\in C(\bar{\Omega}, \mathbb{R}^2)$. Thus, $S(t)\varphi$ represents the distribution of fertile algae individuals at time $t\geq 0$.

Let ${\bf{L}}:C([0, L], \mathbb{R}^2)\rightarrow C([0, L], \mathbb{R}^2)$ be defined by

It then follows that ${\bf{L}}(\varphi)(\cdot)$ represents the distribution of the total new population generated by initial fertile algae individuals $\varphi: = (\varphi_2, \varphi_5)$, and hence, ${\bf{L}}$ is the next generation operator. We define the spectral radius of ${\bf{L}}$ as the basic reproduction ratio for algae, that is,

By ^[34] or [20,Theorem 3.1 (ⅰ) and Remark 3.1], we have the following observation.

Lemma 4.2. $\mathcal{R}_0-1$ and $\lambda^{0}$ have the same sign.

We first consider the following auxiliary system:

in $(x, t)\in (0, L)\times (0, \infty)$ with boundary conditions

and initial conditions

The biologically relevant domain for the system (4.8)-(4.10) is given by

For convenience, we let $\mathcal{Y}_{0} = \mathcal{Y}^{+}\backslash \{(0, 0)\}, \ \partial \mathcal{Y}_{0}: = \mathcal{Y}^{+}\backslash \mathcal{Y}_{0} = \{(0, 0)\}$. By Lemma 4.2 and the arguments similar to those in [11,Lemma 3.2,Theorems 3.1 and 3.2], we have the following result.

Lemma 4.3. ([12,LEMMA 3.6])Assume that (4.7) holds. For any $(N^0(\cdot), N_{S}^0(\cdot))\in \mathcal{Y}^{+}$, let $(N(\cdot, t), N_{S}(\cdot, t))$ be the solution of (4.8)-(4.10). Then the following statements are valid:

$\rm{(i)}$ If $\mathcal{R}_0\leq 1$, then $\lim\limits _{t \rightarrow \infty}(N(x, t), N_{S}(x, t)) = (0, 0)$ uniformly for $x \in [0, L]$;

$\rm{(ii)}$ If $\mathcal{R}_0 > 1$, then (4.8)-(4.10) admit a unique positive steady-state solution $(N^{*}(x), N^{*}_{S}(x))$ and for any $(N^0(\cdot), \ N_{S}^0(\cdot))\in \mathcal{Y}_{0}$, we have

Recall that $X^{+} = C([0, L], \mathbb{R}_{+}^6)$ is the biologically relevant domain for the system (4.1)-(4.3). For convenience, we set $X_{0}: = X^{+}\backslash \{(R^{(0)}, 0, 0, R^{(0)}, 0, 0)\}, \ \partial X_{0}: = X^{+}\backslash X_{0} = \{(R^{(0)}, 0, 0, R^{(0)}, 0, 0)\}$. By Lemma 4.3 and the theory of chain transitive sets (see, e.g., ^[19,31]), one can lift the threshold type result of (4.8)-(4.10) to the full system (4.1)-(4.3).

Theorem 4.1. ([12,THEOREM 3.2]) Assume that (4.7) holds. Let

be the solution of (4.1)-(4.3) with initial data in $X^{+}$. Then the following statements are valid:

$\rm{(i)}$ If $\mathcal{R}_0\leq 1$, then

uniformly for $x \in [0, L]$.

$\rm{(ii)}$ If $\mathcal{R}_0 > 1$, then (4.1)-(4.3) admit a unique positive steady-state solution $(R^{*}(x), N^{*}(x), C^{*}(x), R_{S}^{*}(x), N_{S}^{*}(x), C_{S}^{*}(x)),$ and for any

we have

Next, we consider a model incorporating nutrient recycling. Cyanobacteria excrete some toxins that contain nitrogen, a potential limiting nutrient for algae. Hence, chemical decomposition of the toxin results in nutrient recycling ^[10]. We assume that $\epsilon$ represents a dimensionless coefficient that specifies the allocation to toxin production ^[10]. Accordingly, the authors in ^[10] proposed another reaction-diffusion-advection system:

for $(x, t)\in (0, L)\times (0, \infty)$ with boundary conditions

and initial conditions

for $x\in (0, L)$, where $q_N$ ($q_C$) represents the nutrient quota of algae (toxin). The terms $kq_{C}C$ and $kq_{C}C_S$ in (4.11) reflect that the toxin can get recycled back into the system as available nutrient. From the second and fifth equations of (4.11), we realize that only a part, $(1-\epsilon)$, of the nutrient consumed is used for algal growth, which is discounted by the cost of toxin production.

In ^[13], the extinction/persistence of system (4.11)-(4.13) is investigated in terms of a reproduction number by the comparison arguments and the theory of uniform persistence. Due to the introduction of nutrient recycling, the mathematics becomes more challenging. For example, the uniqueness and global attractivity of the positive steady state of system (4.11)-(4.13) are unclear in general. With an additional assumption, we can establish the uniqueness and global attractivity of the positive steady state, see [13,Section 4].

4.2.   A model of harmful algae and zooplankton

In [12,Section 4], the zooplankton is further incorporated into system (4.1)-(4.3). Suppose $Z$ and $Z_{S}$ represent the densities of zooplankton in the flow and storage zones, respectively; $q_{Z}$ is the constant nutrient quota for zooplankton; $m_{Z}$ is the mortality of zooplankton. Then the governing equations take the following form:

in $(x, t)\in (0, L)\times (0, \infty)$ with boundary conditions

and initial conditions

in $x\in (0, L)$. Here $\eta > 0$ is a constant and represents the effect of the inhibitor on zooplankton, the term $e^{-\eta C}$ represents the degree of inhibition of $C$ on the growth rate of zooplankton, and the function $g(N)$ has the following form:

Let $\mathbb{X}^{+}: = C([0, L], \mathbb{R}_{+}^8)$. By comparison arguments, one can show that solutions of system (4.14)-(4.16) exist globally on $[0, \infty)$, and ultimately bounded and uniformly bounded in $\mathbb{X}^{+}$ (see Lemma 4.1 and Lemma 4.2 in ^[12]). Then we define the solution semiflow $\Theta(t):\mathbb{X}^{+}\rightarrow\mathbb{X}^{+}$ of (4.14)-(4.16) by

where $u(x, t, \phi)$ is the solution of (4.14)-(4.16) with $u(\cdot, 0, \phi) = \phi\in \mathbb{X}^{+}$. We can further find a bounded set $\mathbb{D}$ in $\mathbb{X}^{+}$ and a $t_0 > 0$ such that

and $\mathbb{D}$ is positively invariant for $\Theta(t)$ in the sense that

In view of the assumption (4.7), it follows whenever $\alpha\frac{A}{A_{S}}$ is sufficiently large, there exists a constant $r > 0$ such that

where $\mathcal{M}(R, N, C, Z, R_S, N_S, C_S, Z_S) =$

and

We note that the last four equations in system (4.14)-(4.16) have no diffusion terms, and hence, its solution map $\Theta(t)$ is not compact. By arguments similar to those in [11,Lemma 4.1], we have the following result.

Lemma 4.4. ([12,LEMMA 4.3]) Let (4.7) and (4.17) hold. Then the solution semiflow $\Theta(t)$ is $\kappa$-contracting in the sense that $\lim_{t\to \infty}\kappa (\Theta(t)(B)) = 0$ for any bounded set $B \subset \mathbb{X}^{+}$, where $\kappa$ is the Kuratowski measure of noncompactness.

By [29,Theorem 2.6], we have the following result.

Theorem 4.2. ([12,THEOREM 4.1]) Let (4.7) and (4.17) hold. Then $\Theta(t)$ admits a global attractor on $\mathbb{X}^{+}$ that attracts each bounded set in $\mathbb{X}^{+}$.

We note that the system (4.14)-(4.16) admits the following trivial/semitrivial steady states: $E_{0}: = (R^{(0)}, 0, 0, 0, R^{(0)}, 0, 0, 0)$ and

where $\mathcal{R}_0$ is the algal reproduction ratio for system (4.1)-(4.3), and

is the unique positive steady-state solution of (4.1)-(4.3). Linearizing system (4.14)-(4.16) around the state $E_{1}$, we get the following system for the zooplankton compartments ($Z, Z_S$):

The eigenvalue problem associated with (4.18) takes the form:

Due to the loss of compactness, we need to impose the following condition:

The following result is a straightforward consequence of [21,Theorem 2.1].

Lemma 4.5. Assume that condition (4.20) holds. Then the eigenvalue problem (4.19) admits the principal eigenvalue, denoted by $\Lambda^{*}$.

We remark that in [12,Lemma 4.4], the authors used a generalized Krein-Rutman Theorem (see, e.g., ^[23]) to show that (4.19) admits the principal eigenvalue if (4.20) holds and one additional condition is satisfied. Combining [21,Lemmas 2.1-2.3] with [20,Theorem 2.3], one can obtain [21,Theorem 2.1] and hence Lemma 4.5. Thus, Lemma 4.5 is an improved version of [12,Lemma 4.4] since that additional condition is removed. Here we emphasize that Lemmas 2.1-2.3 in ^[21] hold true only for the autonomous system. So the arguments in [12,Lemma 4.4] are still useful for us to establish the existence of the principal eigenvalue for degenerate periodic reaction-diffusion systems.

In the following, we shall adopt the theory developed in ^[20] to define the basic reproduction ratio for zooplankton. Let $\mathbb{S}(t):C([0, L], \mathbb{R}^2)\rightarrow C([0, L], \mathbb{R}^2)$ be the $C_0$-semigroup generated by the following system

Note that $\mathbb{S}(t)$ is a positive $C_0$-semigroup on $C([0, L], \mathbb{R}^2)$. Assume that both zooplankton individuals in the flow and storage zones are near the trivial steady-state solution $(0, 0)$ for (4.18), and introduce fertile individuals at time $t = 0$, where the distribution of initial zooplankton individuals in the flow and storage zones is described by $\varphi: = (\varphi_4, \varphi_8)\in C([0, L], \mathbb{R}^2)$. Thus, $\mathbb{S}(t)\varphi$ represents the distribution of fertile zooplankton individuals at time $t\geq 0$.

Let $\mathbb{L}:C([0, L], \mathbb{R}^2)\rightarrow C([0, L], \mathbb{R}^2)$ be defined by

It then follows that $\mathbb{L}(\varphi)(\cdot)$ represents the distribution of the total new population generated by initial fertile zooplankton individuals $\varphi: = (\varphi_4, \varphi_8)$, and hence, $\mathbb{L}$ is the next generation operator. We define the spectral radius of $\mathbb{L}$ the basic reproduction ratio of zooplankton compartments for system (4.14)-(4.16), that is,

By ^[34] or [20,Theorem 3.1 (ⅰ) and Remark 3.1], we have the following observation.

Lemma 4.6. $\mathcal{R}_0^z-1$ and $\Lambda^{*}$ have the same sign.

Recall that $\mathbb{X}^{+}: = C([0, L], \mathbb{R}_{+}^8)$. Let

and

Now we are in a position to state the main result of this subsection.

Theorem 4.3. ([12,THEOREM 4.2]) Assume (4.7), (4.17) and (4.20) hold. Then the following statements are valid:

$\rm{(i)}$ If $\mathcal{R}_0 < 1$, then the trivial solution $E_{0}$ is globally attractive in $\mathbb{X}^{+}$ for (4.14)-(4.16).

$\rm{(ii)}$ If $\mathcal{R}_0 > 1$ and $\mathcal{R}_0^z > 1$, then system (4.14)-(4.16) admits at least one (componentwise) positive equilibrium

and there is a positive constant $\zeta > 0$ such that every solution

of (4.14)-(4.16) with

satisfies $\liminf_{t\rightarrow \infty}\min_{x\in [0, L]}Z(x, t)\geq \zeta \ \mathit{\mbox{and}}\ \liminf_{t\rightarrow \infty}\min_{x\in [0, L]}Z_S(x, t)\geq \zeta$.

It is still an open problem whether $E_{1}$ is globally attractive in $\mathbb{X}^{+}$ for system (4.14)-(4.16) in the case where $\mathcal{R}_0 > 1$ and $\mathcal{R}_0^z < 1$.

5.   Discussion

This paper surveys mathematical models describing the spatial variation of population dynamics of harmful algae and toxin production and decay in flowing-water habitats ^[9,10,14]. Previous mathematical models have been somewhat simplified, and raise many paradoxes ^[14,18]. One of the paradoxes is the persistence of harmful algae in the river/stream. Intuitively, phytoplankton populations in riverine reservoirs should be washed out by the strong flow, however, we did observe the occurrences of harmful algal blooms. This persistence paradox may be resolved by the complexity of the channel. In fact, the shoreline features and the bed of the channel can retard flow, producing slow-flowing regions. These slow-flowing regions constitute a hydraulic storage zone that may promote algal persistence ^[9,10]. The authors in ^[9] proposed and analyzed system (3.1)-(3.3) under the case where $R^{(0)}(t)\equiv R^{(0)}$ is a positive constant. The analytical and numerical results in ^[9] confirm that the system with a storage zone can enhance the persistence of phytoplankton populations. More precisely, numerical work in ^[9] shows that persistence is possible at higher advective flows for biologically reasonable parameters in the system with a storage zone. The authors in ^[10] also proposed two-vessel gradostat models of algal dynamics, in which one compartment is a small cove connected to a larger lake. Incorporating seasonal temperature variations into two-vessel gradostat models, rigorous analysis of the time-periodic two-vessel gradostat models are given in ^[35], and their numerical simulations on the basic reproduction number also indicate that seasonality can play a central role in the extinction/persistence of harmful algae.

Some previous mathematical models closely related to this survey paper, using ordinary or partial differential equations and integro-differential or integro-difference equations, can be found in ^{[36,37,38,39,40,41,42]}. Those works focus on the investigations of spatial spread and persistence of populations in the river/stream. Recently, the authors in ^[40,43] also studied reaction-diffusion-advection systems describing the growth of a single species where the species lives in both flowing water and river benthos, respectively. The next generation operator mapping the population from one generation to its next generation offsprings was also used to define three different measures that can determine the extinction/persistence of population in a river. The global dynamics and spreading properties were also investigated in ^[22,30] for time-periodic benthic-drift population models. In ^[44], the authors further studied a reaction-diffusion-advection system of two species competing in a river environment where the populations grow and compete in the benthic zone and disperse in the drifting water zone. The influences of advection rates, diffusion rates, river length, competition rates, transfer rates, and spatial heterogeneity on the persistence/coexistence of species were also numerically investigated. Comparing with the models reviewed in this paper, those in ^{[40,41,43,44,45]} neglect the classes of nutrient(s) and toxin(s). Mathematically, these models are similar to our limiting systems (3.10) and (3.18). Recently, the authors in ^{[46,47,48,49,50,51,52]} (and the related references therein) also considered two-species competition models in a one-dimensional advective environment, where the governing equations are restricted to Lotka-Volterra type reaction-diffusion-advection systems. Assuming that the two species share the same resources, these authors focused on the study of different evolution strategies reflected by their different random dispersal rates and/or advection rates.

In a real ecosystem, the interactions of nutrients, harmful algae, toxins and zooplankton can be very complex. For example, in a real reservoir, P. parvum competes for nitrogen and phosphorus with cyanobacteria, which also excrete allelopathic cyanotoxins that inhibit the growth of P. parvum. A small-bodied zooplankton population consume both types of algae for growth, but the dissolved toxins produced by P. parvum also inhibits zooplankton ingestion, growth and reproduction. In order to understand such complex interactions and reactions in an ecosystem, the authors in ^[53] proposed a well-mixed chemostat system to explore the dynamics of nutrients, P. parvum, toxin(s) produced by P. parvum, cyanobacteria, cyanotoxin(s) produced by cyanobacteria, and zooplankton.

In ^[54], the authors further modify the model in ^[53] to an unstirred chemostat model of the dynamics of P. parvum, cyanobacteria, and a zooplankton population, in which spatial variations are included, but the compartments of algal toxins produced by P. parvum and cyanobacteria are neglected. The strength of inhibition/allelopathy is directly determined by the densities of P. parvum and cyanobacteria, respectively, which reduces the numbers of the modeling equations. It turns out that this model system admits a coexistence steady state and is uniformly persistent provided that the trivial steady state, two semi-trivial steady states and a global attractor on the boundary are all weak repellers.

The factors of seasonal temperature, salinity and vertical variations (due to light limitation in deeper riverine systems) also have been known to have crucial influences on the evolution dynamics of harmful algae. It will be challenging and interesting projects if the aforementioned mechanisms are added into the models reviewed in this paper.

Acknowledgments

We are grateful to Professor James P. Grover for the continued discussions on the mathematical models of harmful algal blooms. Research of S.-B. Hsu is supported in part by Ministry of Science and Technology, Taiwan. Research of F.-B. Wang is supported in part by Ministry of Science and Technology, Taiwan, and National Center for Theoretical Sciences (NCTS), National Taiwan University and Chang Gung Memorial Hospital (CRRPD3H0011, BMRPD18, NMRPD5F0543 and CLRPG2H0041). Research of X.-Q. Zhao is supported in part by the NSERC of Canada. We would like to express our thanks to the anonymous reviewers for their careful reading and helpful suggestions which led to improvements of our original manuscript.

Conflict of interest

The authors declared that they have no conflicts of interest to this work.