Abstract
A reaction-diffusion-advection model is proposed for the Zeya Reservoir to study interactions between algae and zooplankton, including the diffusive spread of algae and zooplankton and the sinking of algae. The model is investigated both with and without sinking. Conditions of Hopf and Turing bifurcation in the spatial domain are obtained, and conditions for differential-flow instability that gives rise to the formation of spatial patterns are derived. Using numerical simulation, the authors examine the impacts on algae of different nutrient concentrations, different sinking rates, and various diffusive spreading patterns. Finally, the models with and without sinking are compared, revealing that the sinking of algae plays an important role in the oscillations of algae and zooplankton. All these results may help to achieve a better understanding of the impact of algae in the Zeya Reservoir.
1. Introduction
With the economic development of human society, waters in lakes and reservoirs are experiencing more and more serious eutrophication, which can cause sustained algal growth. With a high degree of eutrophication, the rapid growth of algae may form algal blooms, bringing about ecological failure and even causing harm to humans. For example, due to the eutrophication, algal blooms frequently appear in the Zeya Reservoir in Wenzhou, which is located in a subtropical region; this can cause deterioration in water quality, which can deprive the drinking water of millions of people.
Because of the local and global impacts of algae blooms on water quality, the studies of dynamic model of plankton are of current interest, and numerous studies have addressed the dynamics of phytoplankton, zooplankton, and algal growth [1–3]. In a number of these studies, the researchers have shown great interest in identifying patterns (or spatiotemporal complexity) of phytoplankton and zooplankton growth [4–6].
A number of empirical and theoretical studies on pattern formation in reaction-diffusion systems [7–12] have been carried out since the pioneering work of Turing [13]. Nowadays, many applications use Turing’s scenario to attempt to explain patterns in fish skin, vegetation, plankton, and so on [14–16]. However, the reaction-diffusion system cannot explain well certain phenomena in vegetation, plankton, and other biological systems because these systems include flow. Therefore, it is necessary and interesting to investigate reaction-diffusion systems involving flow, which are called reaction-diffusion-advection systems. Systems exhibiting flow as a part of pattern formation, such as phytoplankton-nutrient systems with sinking, have been recently studied [17]. Therefore, it is practical and significant to study reaction-diffusion-advection systems.
The organization of this paper is as follows. In Section 2, we present the nutrient-algae-zooplankton model, explaining the biological model. In Section 3, the model is analyzed, including the one with diffusion and the one with diffusion and sinking, and the critical conditions of unstability are obtained. In addition, a series of simulations are given in Section 4. Using simulation, we investigate the effect of critical factor on the model. Finally, the paper ends with a brief discussion and conclusions in Section 5.
2. Model
The reservoir ecology is a relatively ideal ecological environment which includes a number of species, mainly phytoplankton (with algae composing the majority of the phytoplankton), zooplankton, and fish. In a reservoir, two food chains typically exist: (I) the grazing food chain: nutrients-phytoplankton-zooplankton-fish and (II) the detritus food chain: dead organisms-bacteria or benthonic organisms-zooplankton-fish [18]. These two food chains constitute the food web, which maintains the ecological balance of the reservoir. In this food web, once the trophic level is abnormal, the balance may be disturbed, for example, by an algal bloom.
In recent years, more and more researchers have been interested in studying the nutrients-phytoplankton-zooplankton system because of the serious impact of algal blooms on the ecological environment and on human health. A number of researchers have simulated the relations among phytoplankton, zooplankton, and fish using reaction-diffusion models [19–22]. Recently, some studies have included the sinking of phytoplankton in their models [23–25]. Sedimentation losses of phytoplankton seem to be an important factor that may limit phytoplankton biomass at low mixing depths [26–28]. Because most phytoplankton taxa have a higher specific mass than water, all nonmotile and negatively buoyant taxa tend to sink out of the epilimnion with time [23]. The specific sedimentation loss rate of a given taxon depends strongly on its sinking velocity, and therefore sinking velocity is an important determinant of sedimentation losses; the actual sedimentation loss rate also depends on mixing depth [23].
Based on the above analysis, consider a vertical water column. Let indicate the depth of the water column. Let denote the algae population density and the zooplankton population density. The algae-zooplankton model can be described by a reaction-diffusion system with advection: where is the nutrient level of the system, is the specific growth rate, is the grazing rate of zooplankton on algae, is the prey assimilation efficiency of zooplankton, is the nutrient availability constraint, is the carrying capacity, is the handling time, is the mortality and respiration rate of zooplankton, is the algae sinking velocity, and are the respective diffusion coefficients, . In model (1), all parameters are positive constants.
The nutritional growth function of algae, as well as the dependence of the zooplankton grazing rate on algae density, is in agreement with the Monod type [29]. Hence, in the absence of zooplankton, algal growth will saturate at , where is the competition coefficient of the algae population. Growth constraints imposed by different nutrients are not considered separately, but an overall carrying capacity, which depends on the total nutrient level, is assumed [30]. Algae predation by zooplankton follows a Michaelis-Menten (or Michaelis-Menten-Holling) type model; in fact, in recent years, many experiments and observations have shown that the functional response should be described by a ratio-dependent function [31, 32]. In addition, many works in the literature have examined the so-called Michaelis-Menten-type ratio-dependent predator-prey model [33–38].
3. Analysis of the Model
To perform the analysis, model (1) is rewritten as follows: where −.
3.1. Model without Diffusion and Sinking
The nonspatial model has two positive equilibriums under the conditions and , which correspond to spatially homogeneous equilibrium of the full model (1). The equilibriums of the system are(i) (extinction of the zooplankton);(ii), where , .
The Jacobian matrix, , of nonspatial system at the equilibrium is
From , when the condition, , holds, the equilibrium is locally stable, whose index is +1, and then there is no positive equilibrium in the nonspatial model. When the condition, , holds, the equilibrium is saddle and then there is a positive equilibrium, , in the nonspatial model when the condition, , holds.
The Jacobian matrix, , of nonspatial system at the equilibrium is
By , the index of is +1, and the equilibrium is locally stable when the condition, , holds.
3.2. Model with Diffusion
Now let us investigate the most interesting equilibrium point . To perform a linear stability analysis, model (2) is rewritten as follows:
Model (5) can be linearized around the spatially homogeneous fixed point for small space- and time-dependent fluctuations and expanded in Fourier space: The following characteristic equation, then, can be obtained: where
Note that (7) can be solved using the characteristic polynomial of the original problem: where
The roots of (9) can be expressed in the following form:
At the bifurcation point, two equilibrial states of the model intersect and exchange their stability. The space-independent Hopf bifurcation breaks the temporal symmetry of the system, giving rise to oscillations that are uniform in space and periodic in time. In addition, the Turing bifurcation breaks the spatial symmetry, forming patterns that are stationary in time and oscillatory in space.
The Hopf bifurcation occurs when Then the critical value of the transition, the critical parameter on Hopf bifurcation, can be obtained as The Turing bifurcation occurs when The critical value of the bifurcation parameter is where , .
From previous analysis, we can know that the diffusion affects the stability of equilibrium. Next, we will discuss the effect of sinking on equilibrium.
3.3. Model with Diffusion and Sinking
In this subsection, the model with sinking will be discussed. To analyze this model, following Murray [15], zero-flux boundary conditions will be assumed and initial conditions specified. By substituting and into model (1) and assuming , a linear approximation of model (1) can be obtained as follows:
The initial conditions are assumed as follows: and , where the functions and decay rapidly for . Following the standard approach, Laplace transformation of the linearized equations over the two independent variables and is performed. For , the so-called two-sided version of the transformation is used [39]. The relations for the forward and backward transforms are where and are complex variables. In (17) for the backward transformation, the integration contour in the -plane is the imaginary axis. In the -plane, the contour is parallel to the imaginary axis and located to the right of all singularities of the integrand.
After this transformation, the kinetic equations read as follows: where and are the transforms of and . To reveal the presence of instability and disclose its nature, it is sufficient to consider one variable. Further calculations yield where the denominator is Hence, the dispersion equation is For the sake of convenience, the one-dimensional case is considered first. From the quadratic equation (22), the roots can be obtained, where Straightforward manipulation of (23) yields where and denote the real and imaginary parts of , respectively.
The condition for a spatial mode to be unstable is that model (1) grows into a pattern; that is, . If for any values of , then pattern formation will occur. Therefore, it is necessary to find the critical value of for . In fact, what is needed is the maximum value of . Thus, the critical condition can be obtained using , as follows: where
Generically, (27) has a strictly positive minimum . When (27) holds, the neutral curve has a minimum at a nonzero value of .
From previous analysis, we can know that the sinking affects the stability of equilibrium.
4. Numerical Results
In the previous, we analyzed the model, including nonspatial model, the model with diffusion, and the model with diffusion and sinking. In this section, in order to carry out numerical simulations, it is necessary to discretize the space and time variables in the problem. So-called “discretization” means that the continuous problem is transformed into a discrete problem; that is, the continuous space becomes a discrete domain with lattice sites. The spacing between the lattice points is defined by the lattice constant . The following discussion uses an explicit Euler method with time step for time integration and a finite-difference method for the diffusion and advection terms, in which the difference approaches the derivatives when . The numerical simulation is performed over mesh points with spatial resolution , and the zero-flux boundary condition is used. Because the spatiotemporal dynamics of the system depend on the choice of initial conditions, random spatial distributions of the species around the nontrivial stationary state were chosen, which seems to be a fairly general assumption from the biological point of view. In the rest of this paper, the model parameters are set to
Firstly, we consider the Hopf and Turing bifurcations of model (1) without sinking in the parameter space spanned by the parameters and , as illustrated in Figure 1(a). The Hopf bifurcation line and the Turing bifurcation line intersect at a codimension-2 bifurcation point, called the Turing-Hopf bifurcation point. The parametric space can be separated into four domains. Domain I, located above both bifurcation lines, corresponds to systems with homogeneous equilibrium which is unconditionally stable. Domain II is the region of pure Hopf instability, and domain III is that of pure Turing instability. In Domain IV, both Hopf and Turing instabilities occur. Figure 1(b) shows the dispersion relationships of the unstable Hopf mode and the Turing mode as they move from stable to unstable states, using model (2).
(a)
(b)
By (23), the critical wavelength () can be obtained such that , where and are chosen as control parameters, and their corresponding values are and , as shown in Figure 2. Figure 2 depicts the range of values of for certain parameter values, and it is apparent that spatial patterns can occur when . At the same time, it is easy to see that instability occurs when (Figure 2(a)), while instability also occurs when (Figure 2(b)).
(a)
(b)
According to (27), the neutral curve has a minimum at a nonzero value of , as shown in Figure 3(a). From Figure 3(a), it is easy to see that increases with increasing nutrient concentration . The reason for this is that the algae may pursue nutrients that have been deposited on the bottom. In this way, nutrients can influence the sinking of algae. Moreover, using (27), a parameter space ( space) can be obtained in which the stable and unstable domains of model (1) can be determined, as shown in Figure 3(b). Figure 3(b) also shows that increases with increasing nutrient concentration . Moreover, instability will occur for any when is fixed, as is shown in Figure 3(b) as domain II (identified as the “spatial patterns” domain).
(a)
(b)
It is obvious that detailed numerical investigation of the model in the whole parameter space is virtually impossible because the model has a relatively large number of parameters. However, the numerical simulations show that the type of dynamics in the system is determined by the values of, , and , which represent the effect of nutrients, the diffusion of zooplankton, and the sinking of algae.
Figure 4 shows the stripelike patterns that spontaneously form at . From Figure 4, it is clear the stripelike spatial patterns arise from random initial conditions. After some time, the stripelike patterns form from the initial state (Figure 4(b)) and grow steadily over time until they reach a certain width. Finally, the stripelike spatial patterns prevail over the whole domain (Figure 4(c)).
(a)
(b)
(c)
Figure 5 shows snapshots of algae spatial patterns with at 0, 50000, and 100000 iterations. From Figure 5, it is interesting to note that the spotted spatial patterns arise from the random initial conditions. The stripelike pattern forms from the initial state (Figure 5(a)) after some time has passed. Surprisingly, however, the final pattern becomes a spotted spatial pattern. Compared with Figure 4, only the amount of nutrients has been changed, but this is clearly an essential difference. Therefore, nutrient concentration evidently plays an important role in the spatial distribution of algae.
(a)
(b)
(c)
Figure 6 shows snapshots of algae spatial patterns with at 0, 50000, and 100000 iterations. Comparing Figures 5 and 6, the final pattern in Figure 6 is similar to that in Figure 5 (they are both spotted), but their formative processes are different. Moreover, it is obvious that the density of spots is sparser in Figure 6(c) than in Figure 5(c). In fact, the difference between Figures 5 and 6 is the change in diffusion coefficient. These results show that the diffusion of zooplankton is a key factor in the formative processes of algae patterns.
(a)
(b)
(c)
Now including sinking in model (1), that is, , the effects of sinking are simulated, and the results are shown in Figures 7 and 9. Figure 7 shows snapshots of algae spatial patterns with at 0, 50000, and 100000 iterations. In Figure 7, vertical stripes move from top to bottom and then exit from the bottom. Compared to Figure 4, it is apparent that sinking exists in Figure 7. To see the effect of sinking, is replaced by , and the resulting patterns are shown in Figure 8. If Figure 9 is compared with Figure 6, the spotted pattern becomes vertical stripes. Note that Figures 7, 8, and 9 are quite similar at the end.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Based on the above analysis, it is apparent that parameters , and are critical factors for controlling the amount and the spatial distribution of algae. Moreover, the sinking of algae can play an important role in the oscillations of algae and zooplankton. Such evolutionary developments of the reaction-diffusion-advection system may well be taking place in the real world.
5. Discussion and Conclusions
In the present paper, we used a reaction-diffusion-advection model to investigate the interaction between algae and zooplankton, where the sinking of algae was considered which was described by the advection term. Using linear analysis technique, the model was analyzed, and the effects of critical factors, including , , and , on the model were discussed. Then the instability conditions of Hopf and Turing were obtained. In addition, from the analysis described in Section 3 and Figures 1 and 9, it was not difficult to find that the theoretical results were consistent with the numerical simulations.
The model was relatively simple because it was only an abstraction of real-world phenomena, but it reproduced many features of real-world phenomena. In the real world, many patterns have been observed, such as banded vegetation, patches, and spiral waves, which could be regular or irregular. The patterns may be forced to occur by physical factors or internal factors. While our explanation focused on a predator-prey interaction between algae and zooplankton, especially, how the sinking of algae affected the interaction and patterns, our analytical results showed that the homogeneous steady state became unstable due to the sinking of algae. In addition, diffusion also led to instability of the homogeneous steady state. But, the sinking of algae forced band patterns to occur. That is, when the sinking flux of algae was beyond a certain critical value which could be obtained from our analyzed results, the spatial distribution of population was not homogeneous because the sinking of algae existed.
From the numerical results, it was obvious that different , , or can lead to different patterns at the same time. Many field studies indicated that the oscillation of population density occurs in reality, while Figures 4 and 9 implied that nutrient, diffusion, and the sinking can result in oscillation of population density. It means that these factors may play an important role in the change of algae population density. These results may provide better understanding on the study of algal dynamics in water environment.
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. 31170338), by the Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001), and by the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510).