Mathematical Problems in Engineering

Volume 2015, Article ID 251937, 12 pages

http://dx.doi.org/10.1155/2015/251937

## The Dynamics of a Diffusive Nutrient-Algae Model Based upon the Sanyang Wetland

^{1}School of Mathematics and Information Science, Wenzhou University, Wenzhou, Zhejiang 325035, China^{2}Zhejiang Provincial Key Laboratory for Water Environment and Marine Biological Resources Protection, Wenzhou University, Wenzhou, Zhejiang 325035, China^{3}School of Life and Environmental Science, Wenzhou University, Wenzhou, Zhejiang 325035, China

Received 25 September 2014; Accepted 21 January 2015

Academic Editor: P. Balasubramaniam

Copyright © 2015 Yi Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 [5–8], 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 [17–24]. 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 is**It 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*

*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:*