Abstract

The stability and spatiotemporal dynamics of a diffusive nutrient-algae model are investigated mathematically and numerically. Mathematical theoretical studies have considered the positivity and boundedness of the solution and the existence, local stability, and global stability of equilibria. Turing instability has also been studied. Furthermore, a series of numerical simulations was performed and a complex Turing pattern found. These results indicate that the nutrient input rate has an important influence on the density and spatial distribution of algae populations. This may help us to obtain a better understanding of the interactions of nutrient and algae and to investigate plankton dynamics in aquatic ecosystems.

1. Introduction

Diffusive processes are often used to represent the formation of spatial patterns in biological systems [1]. Pattern formation in nonlinear complex systems is one of the central problems of the natural, social, and technological sciences. The occurrence of multiple steady states and transitions from one to another after critical fluctuations, the phenomena of excitability, oscillations, and waves, and the emergence of macroscopic order from microscopic interactions in various nonlinear nonequilibrium systems in nature and society have been the subject of many theoretical and experimental studies [2, 3]. Plankton plays an important role in ocean ecology and climate because of their participation in the global carbon cycle at the base of the food chain [4]. Now, using the reaction-diffusion equation and patterns to study the spatiotemporal dynamics of the planktonic ecosystem has aroused the interest of many researchers.

In recent years, the relationships among nutrients, phytoplankton, zooplankton, and fish have become the focus of researchers [58], including the control of phytoplankton [9, 10], and dynamic analysis of the relationships [11, 12]. Numerous theoretical ecologists have built many models to reveal the inner relationships among these populations and to investigate their dynamics. Drago et al. [13] presented a three-dimensional numerical model to analyze the dispersion of suspended solids and conservative pollutants released into ambient water and their effect on trophic behavior. Luo [14] derived and analyzed a mathematical model for interactions between phytoplankton and zooplankton in a periodic environment, where the growth rate and the intrinsic carrying capacity of phytoplankton are changing with respect to time and nutrient concentration. James et al. [15] built a model of the evolution of phytoplankton, zooplankton, and fish to investigate the mechanisms that trigger plankton blooms. Alvarez-Vázquez et al. [16] presented a mathematical model involving nutrients, phytoplankton, zooplankton, organic detritus, and dissolved oxygen to simulate eutrophication processes in aquatic media.

The Sanyang wetland is located in Wenzhou, in a subtropical region. It covers an area of 13 square kilometers and contains many streams. The ratio of water area to land in the Sanyang wetland is as high as 1.1 : 1. With rapid economic development in Wenzhou, most rivers and wetlands are in a state of eutrophication. Because of eutrophication, nuisance cyanobacterial blooms have taken place several times in the Sanyang wetland in recent years.

Recently, more and more researchers have investigated the dynamics of the planktonic ecosystem by means of spatial patterns [1724]. Self-organizing spatial patterns appeal to the theoretical biologist because of the pioneering work of Turing [25]. Many self-organizing spatial patterns have been found in their research, including patchiness, bands, and spiral waves. Wang et al. [26] used a reaction-diffusion-advection model of algae and mussels to demonstrate that young mussel beds on soft sediments can display large-scale regular spatial patterns and banded patterns. Their work is significant to study spatial patterns. Serizawa et al. [27] presented a minimal nutrient-phytoplankton model that can exhibit various types of spatial patterns, including patchiness. Although the model is simple, their work is still of interest. Liu et al. [28] studied a spatially extended nutrient-phytoplankton-zooplankton-fish reaction-diffusion system and found spiral waves and spatial chaos patterns. Their work is important to plankton ecological system. van de Koppel et al. [29] proposed a simple spatial model of algae and mussels and found simulated spatial patterns which are consistent with previous observations. Their works vigorously promote the study of emergent spatial patterning at larger spatial scales.

The nutrients in the Sanyang wetland comes mainly from sewage output by nearby residents and wastewater discharges from industries and businesses. Therefore, the input rate of nutrients flowing into the Sanyang wetland can be assumed constant. Many researchers have used to describe natural nutrient removal [27, 30, 31]. Therefore, the term has been used here to describe the natural removal of nutrients (mainly in the form of sedimentation) in the Sanyang wetland because of various complex long-term factors. However, this term has another important meaning in the mathematical model. Without this factor, the nutrient concentration will increase without limit when the algae are completely extinct. Wenzhou belongs to a subtropical monsoon climate zone and is perennially windy. Therefore, the effects of wind and waves on nutrients in the Sanyang wetland must be considered. Because of wind and wave action, a positive nutrient input is released into the water from bottom sediments. Ecologists usually approximate this input by the following sigmoid function [32]:

The Holling functional response function is suitable for algae, cells, and lower organisms. Many studies have shown that algae density in the Sanyang wetland is positively correlated with nutrient concentration within certain limits. Hence, in this paper, the Holling response function has been used to describe absorption of nutrients by algae. In the Sanyang wetland, nutrients are limited, and density restrictions exist in the algae population. Therefore, it is suitable to use a logistic growth function to describe algae population growth in the Sanyang wetland.

Based on this discussion, the following dynamic reaction-diffusion model has been proposed for nutrients and algae in the Sanyang wetland:where and denote the nutrient concentration and the algae density. Research has shown that phosphorus is limited and nitrogen is abundant in this wetland. Therefore, the main factor affecting algae population growth is phosphorus, and, in this paper, the nutrient has been assumed to be phosphorus. Let be the input rate of nutrients flowing into the water, the maximum growth rate of the algae population, the efficiency of conversion, the intrinsic growth rate, the carrying capacity of the algae population, the death rate of the algae population, and and the diffusion coefficients of nutrients and algae, respectively. is the usual Laplacian operator in two-dimensional space, and , , , , , , and are positive constants.

Model (2) will be studied under the following nonzero initial conditions,and the following zero-flux boundary conditions,where and are the size of the system in the - and -directions and is the outward unit normal vector of the boundary , which is assumed to be smooth.

For the reaction-diffusion nutrient-algae system (2), the reduced system is an ordinary differential equation of the form

Next, the ODE system (5) will be analyzed.

2. Dynamics Analysis of the Reduced ODE System

2.1. Positivity and Boundedness of the Solutions

Theorem 1. Suppose that holds; then all the solutions of system (5) with initial conditions are positive and bounded for all .

Proof. From the first equation of system (5), It follows that .
Hence, .
From the second equation of system (5), it can be determined thatleading to .
Then, If , then, for , there exists such that, for any , Moreover, is continuous when .
This means that there exists such that .
Therefore, when , and the proof is complete.

2.2. Existence of Equilibria

Theorem 2. System (5) has boundary equilibrium (extinction of algae population).

Proof. If is the boundary equilibrium of (5), thenAssume the following function:Obviously,It is easy to establish that the curve of function intersects the -axis, and, therefore, the existence of the boundary equilibrium is guaranteed, which completes the proof.

Theorem 3. System (5) has a unique interior equilibrium (coexistence of nutrients and algae) if and .

Proof. If is the interior equilibrium of (5), thenConsequently, is the positive root of the fourth-degree equation:Consider the following function:It is easy to establish thatIt is also easy to determine that if and .
Obviously,It can easily be shown that the curve of function intersects the -axis only once. Hence, the existence and uniqueness of the interior equilibrium is guaranteed, and this completes the proof.

2.3. Local Stability of Equilibria

Theorem 4. The boundary equilibrium of system (5) is locally asymptotically stable if The Jacobian matrix of system (5) is The Jacobian matrix of system (5) at isIt is clear that the Jacobian matrix has two eigenvalues: Hence, the equilibrium is locally asymptotically stable if

Theorem 5. The interior equilibrium of system (5) is locally asymptotically stable if .

Proof. The Jacobian matrix of system (5) at isMoreover, where , , and .
It can be easily determined that and if . This means that the two eigenvalues of the Jacobian matrix have a negative real part. Therefore, the interior equilibrium is locally asymptotically stable.

2.4. Global Stability of Equilibria

Theorem 6. Assuming thatthen the boundary equilibrium of system (5) is globally asymptotically stable.

Proof. Obviously, satisfies the equation .
To prove the above statement, let us consider the following Lyapunov function:It is easy to show that for all . Differentiating along the solutions of system (5), Assuming the following function,then, .
Let ; then, .
Obviously, , .
Therefore, .
Then Obviously, if and , then strictly for all except the boundary equilibrium , where .
Therefore, satisfies Lyapunov’s asymptotic stability theorem, and the boundary equilibrium of system (5) is globally asymptotically stable. This completes the proof.

Remark 7. It is easy to verify that if condition (25) holds, then condition (18) is true.

Theorem 8. The positive equilibrium of system (5) is globally asymptotically stable if .

Proof. Obviously, satisfies the equationsNow let us construct the following Lyapunov function:Obviously, is continuous for all .
By simple computation,This shows that the positive equilibrium is the only extremum of the function in the first quadrant. Obviously, Equations  (32) and (33) show that the positive equilibrium is the global minimum in the first quadrant.
Therefore, for all .
Differentiating along the solutions of model (5), Obviously, . Furthermore, if , , and, therefore, strictly for all except the positive equilibrium , where .
Therefore, satisfies Lyapunov’s asymptotic stability theorem, and the positive equilibrium of system (5) is globally asymptotically stable. This completes the proof.

3. Stability Analysis in the Presence of Diffusion and Turing Instability

3.1. Stability Analysis of Equilibria in the Presence of Diffusion

Theorem 9. If  , the steady state of the diffusive system (2) is unstable.

Proof. To investigate the stability of the equilibrium , let us consider the corresponding eigenvalue problem of the linearized operator around .
The linearization of system (2) at steady state can be expressed as where , , and From the above, the linearized result of system (2) around isThe corresponding characteristic equation of (37) can be obtained aswhere is an eigenvalue of (38) and the corresponding eigenvector is . If , then is an eigenvalue of the operator with a homogeneous Neumann boundary condition. Therefore, must be real. In the same way, is also real if . Hence, all eigenvalues of (38) are real. Let be the largest eigenvalue. Consider the principal eigenvalue of the following equation:From the above equation, and the associated eigenvector if . It is claimed here that is also an eigenvalue of (38). In fact, is taken to be a solution ofthen satisfies (38) with . Therefore, is an eigenvalue of (38). Hence, , and the steady state of the diffusive system (2) is unstable. This completes the proof.

Theorem 10. If and , the diffusive system (2) is stable at .

Proof. In the case and , let be the principal eigenvector of (38) corresponding to the largest eigenvalue .
If , then is also an eigenvalue of (39). Therefore, can be obtained if holds.
If , then . Hence, is an eigenvalue of Obviously, the largest eigenvalue of (41) is . Therefore, if . Therefore, the diffusive system (2) is stable at if and . This completes the proof.

To prove Theorem 12, let us introduce the following lemma [33].

Lemma 11. Consider the following equation:Let be a steady state of (42); that is, .
If , , , , wherethen is uniformly asymptotically stable. Furthermore, if , then is unstable.

Theorem 12. The following conclusions hold.If then the positive constant solution of system (2) is uniformly asymptotically stable.If then the positive constant solution of system (2) is unstable.

Proof. ConsiderObviously, , , , if .
Applying Lemma 11, the first conclusion of Theorem 12 follows.
On the other hand, if , it can easily be established that . Therefore, of system (2) is unstable by Lemma 11. This completes the proof.

Theorem 13. If then the positive equilibrium of system (2) is globally asymptotically stable.

Proof. To prove Theorem 13, it is necessary to construct a Lyapunov function. Define Obviously, , if and only if . Calculating the derivative of along the solutions of model (2), Obviously, , , and .
Consider the following function:Then, . Because , thenTherefore, if , , which is equivalent to .
Thus, satisfies Lyapunov’s asymptotic stability theorem, and the positive equilibrium of system (2) is globally asymptotically stable. This completes the proof.

Remark 14. It is easy to verify that if condition (47) holds, then condition (44) is true.

3.2. Turing Instability

To analyze the spatial system and how a small heterogeneous perturbation of a homogeneous steady state develops over a long time period, small space- and time-dependent perturbations for of system (2) will be considered:where , are small enough and is the wave number. Substituting (52) into (2) yields its characteristic matrix:Correspondingly, and are the roots of the following characteristic equation:, .

The roots of (54) yield the dispersion relation:

As is well known, Turing instability means that the stable equilibrium point is driven to become unstable by the local dynamics and diffusion of species. The stability conditions for of system (5) are and . It is clear that . Therefore, the stability of equilibrium of system (2) changes with the sign of . It is easy to establish that for , where

If has positive values, then the range of instability for can be obtained, which is called the Turing space. To illustrate the Turing space, the dispersion relations corresponding to several values of the bifurcation parameter are plotted in Figure 1. The other parameters come from the actual monitoring data from the Sanyang wetland from 2012 to 2013 and from the literature:

In Figure 1, the green line corresponds to the critical Turing value, . When (the yellow line in Figure 1), Turing instability occurs, whereas when (the red line in Figure 1), the Turing instability decays.

4. Numerical Simulations

On the basis of actual monitoring data from the Sanyang wetland from 2012 to 2013, the following parameter set is considered: , initial value , , and other parameter values as given in (57). To verify the feasibility and correctness of the theoretical results, numerical simulations were performed. It is obvious from Figures 2 and 3 that the interior equilibrium of system (5) is locally asymptotically stable.

To investigate more thoroughly how the input rate of nutrients flowing into the water affects the spatiotemporal dynamics of system (2), spatial distribution diagrams of the system are generated for various values of . All these numerical simulations employ the zero-flux boundary conditions and a discrete two-dimensional domain with . The spatial step size is , and the temporal step size is . The initial value of system (2) is placed at the interior equilibrium , and the initial perturbation is 0.0005 space units per time unit. As is well known, in such numerical simulations, the spatial distributions of predator and prey are always of the same type. Therefore, it is necessary to analyze the spatial distribution of only one of these. Here, the spatial distribution of the algae population is considered. Some snapshots have been extracted and are shown in red (blue) corresponding to the high (low) value of algae density when the value of increases from 0.005 to 0.006.

Figure 4 shows the process of pattern formation in system (2) with and other parameters fixed as in (57).

From Figure 4, it is clear that the algae population pattern takes a long time to settle down, beginning with a uniform state (Figure 4(a)). After 40000 iterations, the spatial distribution of the algae population consists mainly of a number of spots and interconnected stripes (Figure 4(b)). The remaining four images consist of blue/red spots on a red/blue background. These are referred to here as “hot spots” and “cold spots” [34]. The hot spots are isolated zones with high algae density and low nutrient concentration, whereas the cold spots represent the opposite case. Figures 4(c) and 4(d) show some isolated “cold spots,” whereas Figures 4(e) and 4(f) contain some isolated “hot spots” and stripes.

Figure 5 shows two patterns obtained from system (2) at 300000 iterations. When , the homogeneous state is , and the spot pattern is made up of “hot spots” around the edge (Figure 5(a)). However, with , , and the spot-stripe pattern consists of a certain number of “hot spots” and stripes (Figure 5(b)).

By comparing Figures 5(a), 4(f), and 5(b), it is apparent that the algae population density increases with within certain limits. Moreover, as the value of increases, the spatial distribution of the algae population becomes relatively intensive and the distribution relatively wide.

5. Discussion and Conclusions

In this paper, a reaction-diffusion nutrient-algae model has been used to investigate the interaction between nutrients and algae mathematically and numerically. Mathematical theoretical work has examined the positivity and boundedness of solutions and the existence and stability of equilibria. Local and global stability analyses of equilibria in the presence of diffusion have also been performed. Mathematical analysis indicated that all solutions of the reduced ODE system are positive and bounded under some certain conditions and that the ODE system always has boundary equilibrium. If and , then the nonspatial system and the diffusive system are stable at . When , the interior equilibrium is always stable, regardless of whether the system is an ODE or a PDE system. Theorem 12 gives the stability condition of the positive constant solution of system (2). However, when , the stability of may vary. In this paper, only the positive constant solution of system (2) is discussed. Other questions remain to be answered: whether system (2) has other nonnegative constant solutions or nonconstant positive steady states; if there are positive solutions, what is their exact multiplicity and how stable are any positive solutions.

Numerical simulations have indicated that the input rate of nutrients flowing into the water has an important influence on the density and spatial distribution of the algae population. The spatial distribution of the algae population becomes relatively intensive and the distribution relatively wide as increases. These results may help to provide a better understanding of the interactions of nutrients and algae and the variations in the spatial distribution of algae over time. This research is also expected to contribute to exploring the mechanisms of eutrophication.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 31370381) and by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001).