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

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.


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.
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][18][19][20][21][22][23][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-phytoplanktonzooplankton-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 reactiondiffusion model has been proposed for nutrients and algae in the Sanyang wetland: where  and V 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 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.
Proof.If  1 = (, 0) is the boundary equilibrium of ( 5), then Assume 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 𝐸
Proof.If  * = ( * , V * ) is the interior equilibrium of (5), then Consequently,  * is the positive root of the fourth-degree equation: Consider the following function: It is easy to establish that It is also easy to determine that ℎ  () > 0 ( > 0) if  2  −  −  > 0 and ( − ) +  > 0. 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.

Local Stability of Equilibria
Theorem 4. The boundary equilibrium  1 = (, 0) of system ( 5) is locally asymptotically stable if The Jacobian matrix of system ( 5) is The Jacobian matrix of system ( 5) at  1 = (, 0) is It is clear that the Jacobian matrix  ( 1 ) has two eigenvalues: Hence, the equilibrium  1 = (, 0) is locally asymptotically stable if Theorem 5.The interior equilibrium  * = ( * , V * ) of system ( 5) is locally asymptotically stable if Proof.The Jacobian matrix of system (5) at Moreover, where tr( It can be easily determined that tr( ( * ) ) < 0 and det( ( * ) ) > 0 if  11 < 0. This means that the two eigenvalues of the Jacobian matrix  ( * ) have a negative real part.Therefore, the interior equilibrium  * = ( * , V * ) is locally asymptotically stable.

Global Stability of Equilibria
Theorem 6.Assuming that then the boundary equilibrium  1 = (, 0) of system ( 5) is globally asymptotically stable.
Remark 7. It is easy to verify that if condition (25) holds, then condition ( 18) is true.
By simple computation, This shows that the positive equilibrium ( * , V * ) is the only extremum of the function  (,V) in the first quadrant.Obviously, Equations ( 32) and (33) show that the positive equilibrium ( * , V * ) is the global minimum in the first quadrant.Therefore,  (,V) ≥  ( * ,V * ) = 0 for all , V > 0.

Stability Analysis of Equilibria in the Presence of Diffusion
Theorem 9.If  +  −  > 0, the steady state  1 = (, 0) of the diffusive system (2) is unstable.
Proof.To investigate the stability of the equilibrium  1 , let us consider the corresponding eigenvalue problem of the linearized operator around  1 .The linearization of system (2) at steady state  1 = (, 0) can be expressed as where  = (, V)  ,  = diag( 1 ,  2 ), and From the above, the linearized result of system (2) around  1 is The corresponding characteristic equation of (37) can be obtained as where  is an eigenvalue of (38) and the corresponding eigenvector is (, V).If V ̸ = 0, then  is an eigenvalue of the operator  2 Δ + ( +  − ) with a homogeneous Neumann boundary condition.Therefore,  must be real.In the same way,  is also real if  ̸ = 0. Hence, all eigenvalues of (38) are real.Let  max be the largest eigenvalue.Consider the principal eigenvalue λ of the following equation: From the above equation, λ > 0 and the associated eigenvector V > 0 if  +  −  > 0. It is claimed here that λ is also an eigenvalue of (38).In fact, û > 0 is taken to be a solution of then (, V) = (û, V) satisfies (38) with  = λ.Therefore, λ > 0 is an eigenvalue of (38).Hence,  max ≥ λ > 0, and the steady state  1 of the diffusive system (2) is unstable.This completes the proof.
To prove Theorem 12, let us introduce the following lemma [33].
Applying Lemma 11, the first conclusion of Theorem 12 follows.
On the other hand, if 2) is unstable by Lemma 11.This completes the proof.
Consider the following function: Then, Thus,  1 () satisfies Lyapunov's asymptotic stability theorem, and the positive equilibrium ( * , V * ) 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.

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,  1 and  2 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 tr 0 =  11 +  22 < 0 and Δ 0 =  11  22 −  12  21 > 0. It is clear that tr  < tr 0 < 0. Therefore, the stability of equilibrium  * of system (2) changes with the sign of Δ  .It is easy to establish that Δ  < 0 for  1 2 < If  1,2 2 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 In Figure 1, the green line corresponds to the critical Turing value,  * = 0.006195.When  = 0.0055 <  * (the yellow line in Figure 1), Turing instability occurs, whereas when  = 0.008 >  * (the red line in Figure 1), the Turing instability decays.

Numerical Simulations
On the basis of actual monitoring data from the Sanyang wetland from 2012 to 2013, the following parameter set is considered:  = 0.008, initial value ((0), V(0)) = (0.1, 0.1),  = 40, 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  * = (0.02462459569, 0.5574801727) 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   =   = 200.The spatial step size is ℎ = 1/3, and the temporal step size is  = 0.01.The initial value of system ( 2) is placed at the interior equilibrium  * = ( * , V * ), 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 V when the value of  increases from 0.005 to 0.006.
Figure 4 shows the process of pattern formation in system (2) with  = 0.0055 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 ( * , V * ) = (0.7722029183, 0.7931146159) (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  = 0.005, the homogeneous state is ( * , V * ) = (0.7498330425, 0.7860636982), and the spot pattern is made up of "hot spots" around the edge (Figure 5(a)).However, with  = 0.006, ( * , V * ) = (0.7860324719, 0.7974736497), and the spot-stripe pattern consists of a certain number of "hot spots" and stripes (Figure 5(b)).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.

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 − + 2/(1 +  2 ) 2 < 0 and  +  −  < 0, then the nonspatial system and the diffusive system are stable at  1 .When − − V * + 2 * /(1 +  * 2 ) 2 < 0, the interior equilibrium  * = ( * , V * ) 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  * = ( * , V * ) of system (2).However, when 0 < − − V * + 2 * /(1 +  * 2 ) 2 < (/) 2  * , 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.
. The other parameters come from the actual monitoring data from the Sanyang wetland from 2012 to 2013 and from the literature:

Figure 1 :Figure 2 :
Figure 1: Variation of the dispersion relation of model (2) around  * .The yellow line corresponds to  = 0.0055, the green to  * = 0.006195, and the red to  = 0.008.

Figure 3 :
Figure 3: Phase trajectories of nutrients and algae beginning at different initial levels.
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  1 and  2 the diffusion coefficients of nutrients and algae, respectively.Δ =  2 / 2 +  2 / 2 is the usual Laplacian operator in twodimensional space, and , , , , , , and  are positive constants.