Abstract

We investigated a nonlinear model of the interaction between nutrients and plankton, which was addressed using a pair of reaction-advection-diffusion equations. Based on numerical analysis, we studied a model without diffusion and sinking terms, and we found that the phytoplankton density (a stable state) increased with the increase of nutrient density. We analyzed the model using a linear analysis technique and found that the sinking of phytoplankton could affect the system. If the sinking velocity exceeded a certain critical value, the stable state became unstable and the wavelength of phytoplankton increased with the increase of sinking velocity. Furthermore, band patterns were also produced by our model, which was affected by the diffusion and sinking of phytoplankton. Thus, the change in the diffusion and sinking of phytoplankton led to different spatial distributions of phytoplankton. All of these results are expected to be useful in the study of plankton dynamics in aquatic ecosystems.

1. Introduction

Plankton play an important role in the ecology of the ocean and climate because of their participation in the global carbon cycle at the base of the food chain [1]. In certain environmental conditions, lakes, reservoir, and marine waters may experience plankton or algal blooms [2, 3]. However, the local and global impacts of plankton blooms on water quality, carbon cycling, and climate may be damaging. If nutrient source is abundant, and some conditions are satisfied, blooms may become long-term events that affect ecosystems. Plankton blooms can change the types of species present at the base of the aquatic food web and affect human health. Thus, the study of plankton dynamics is currently of major interest.

In the past years, there were many researches on the model between nutrient and phytoplankton and zooplankton [46]. A larger number of researchers have attempted to model the relationship between nutrient and phytoplankton and zooplankton, to investigate the dynamics in plankton model. Truscott and Brindley [7] presented a model for the evolution of phytoplankton and zooplankton populations which resembles models for the behavior of excitable media. Luo [8] investigated phytoplankton-zooplankton dynamics in periodic environments, where eutrophication was considered. El Saadi and Bah [9] modeled phytoplankton aggregation using numerical treatment and explored the asymptotic behavior of the model. Banerjee and Venturino [10] studied a phytoplankton-toxic phytoplankton-zooplankton model and found that the toxic phytoplankton does not drive the zooplankton population towards extinction under a certain mechanism. The result is very important for study on plankton. These works make contributions for the study on plankton.

In recent years, many ecologists have paid increasing attention to spatial processes in a wide variety of practical contexts [11]. For example, theoretical community ecologists have explored ecosystems, including vegetation systems [12] and phytoplankton systems [13]. In particular, the modeling of plankton systems is becoming increasingly important because of their roles in carbon cycling and temperature control, particularly their major impacts on global climate change [14]. These modeling strategies are focused on two areas: (i) studies of large and complex systems, which are eventually used to fit field data or to forecast future changes; and (ii) studies of skeleton models for various mechanisms, which can provide insights or stimulate new experiments [14].

The present study belongs to the latter area. We propose a model on phytoplankton using a pair of reaction-advection-diffusion equations, which allow spatial phenomena, such as sinking, and turbulence to be described directly, thereby enabling spatial structures to be studied. It is known that the sinking and mixing of phytoplankton have pronounced effects on the tendency of different phytoplankton to increase. Experimental studies indicated that most fresh water diatoms and other phytoplankton sink in undisturbed water [15]. Theoretical results also demonstrated the importance of sinking, mixing, and diffusion [16, 17]. Theoretical models predicted that a process with reduced vertical mixing may induce oscillations and chaos in the phytoplankton of the deep chlorophyll maxima, which leads to differences in the timescale between the sinking flux of phytoplankton and the upward flux of nutrients [18]. A remarkable finding was the survival of a sinking phytoplankton population even when the diffusivity in the deep layers could not prevent population washout [1921]. Mellard et al. demonstrated how externally imposed heterogeneity in the form of resource gradients and mixing interacted with internally generated heterogeneity in the form of competition, population dynamics, and movement to determine the spatial distribution of phytoplankton [22]. Ryabov et al. showed that the upper mixed layer was an important factor that had the potential to shape the spatial distribution and species composition of phytoplankton, but it also provided insights with general ecological importance [23]. van de Koppel et al. studied self-organized spatial patterning in an algae-mussel model, where regular spatial patterns were formed in young mussel beds on soft sediments in the Wadden Sea [24]. Self-organized spatial patterns are of considerable interest to theoretical biology [2530] because of the basic paper by Turing [31] on the role of nonequilibrium reaction-diffusion prepatterns in biomorphogenesis. Furthermore, recent modeling studies of plankton support the self-organized spatial patterns, such as patchiness [32, 33] and bands [25, 34, 35].

The rest of this paper is organized as follows. In the next section, we present a model based on the theoretical ecology and partial differential equations, which is addressed using a pair of reaction-advection-diffusion equations. In Section 3, we analyze stable behavior of the nonspatial system firstly. What is more, the stable behavior of the spatial system is analyzed. And we obtain the condition under which the steady state becomes unstable. Finally, a series of simulations are given. Using simulation, we investigate the effect of critical factor on the system. In Section 4, discussion and conclusion are presented.

2. The Model

Natural ecosystems of plankton exhibit great variability in space and time. The growth of phytoplankton is dependent mainly on nutrients and light. After the mortality of phytoplankton, nutrients are returned to the system over short time scales with minimal losses [39] through microbial decomposition. In addition, biological factors such as higher predation and physical factors such as the sinking of phytoplankton into the water column also affect the ecosystem, which has been examined previously [40]. Turbulence also affects these systems [3, 19, 41]. Vertical mixing brings nutrients from the lower layers of the ocean into the mixed layer. Based on the previous analysis, the following general structure is obtained:where is the nutrient density and is the phytoplankton population density.

Dugdale proposed the use of Michaelis-Menten enzyme kinetics to describe nutrient-phytoplankton interactions [42]. The Michaelis-Menten equations have the same form as the well-known Monod equations [43], which are used in the Droop equations, and they have formed the basis of a number of modeling studies that aimed to simulate phytoplankton blooms [44]. Thus, we employed Michaelis-Menten kinetics in terms of “uptake.” Furthermore, a Holling type II functional response has been used widely to describe zooplankton predation in various theoretical studies [45, 46]. It has also been reported that the Holling type II functional response shows good concordance with experimental data [33, 47, 48]. Hence, in the present paper, Holling II functional response is adopted to describe zooplankton grazing on phytoplankton. Therefore, a pair of specific models is defined as follows:where a vertical water column is considered. Let indicate the depth in the water column; is the width in the water column. For vertical mixing, we assume that is a constant concentration, which includes the nutrient input flowing into the system and the nutrient from below the mixed layer, is the rate of exchange between the lower and upper layers, is the nutrient content of phytoplankton, denotes the maximum growth rate of phytoplankton, is the half-saturation constant for nutrients, is the half-saturation constant for phytoplankton, denotes the maximum predation rate of zooplankton on phytoplankton, is the mortality of phytoplankton, is the proportion of nutrients in dead phytoplankton that is recycled, is the sinking velocity of phytoplankton, and and are the diffusion rates of nutrients and phytoplankton, respectively, which are caused by mixing and turbulence. In addition, is the Laplacian operator. Table 1 provides the parameter values used and their units, which is obtained from published studies [19, 3638].

3. Results

3.1. Stable Behavior of the Nonspatial System

In the nonspatial system (i.e., system (2a), (2b) without spatial derivatives), according to and , vertical isocline and horizontal isocline can be obtained, respectively, as follows. Vertical isocline : . Horizontal isocline :. Obviously, the line, , is asymptote of vertical isocline , and the line, , is asymptote of horizontal isocline . In the following discussion, it is assumed that the condition always holds; otherwise phytoplankton become extinct eventually. For vertical isocline , , which is derivative. There are two roots in when the condition holds, root is . Then, it is obvious that the asymptote, , is on the left of line . Therefore, the vertical isocline is continuous when . And holds when ; holds when . From the horizontal isocline , if , then there is a positive equilibrium in the nonspatial system at least; if , then there is a positive equilibrium in the nonspatial system at least when the condition holds.

When the condition holds, there is no root in if . From , vertical isocline is monotone decreasing when , and holds when ; holds when . According to horizontal isocline , if , then there is a positive equilibrium in the nonspatial system at least when the condition holds; if , then there is no positive equilibrium in the nonspatial system when the condition holds.

When the condition holds, there are two roots in if . From , vertical isocline is monotone increasing, when . And holds, when ; holds, when . According to horizontal isocline , if , then there is a positive equilibrium state in the nonspatial system at least when the condition holds; if , then there is no positive equilibrium state in the nonspatial system when the condition holds.

It is noted that these conditions only confirm that there exists positive equilibrium state in the nonspatial system when these conditions are satisfied. It does not mean that there must be no positive equilibrium state in the nonspatial system when these conditions are not satisfied. In addition, it is not difficult to find that there is always a trivial steady state, , consisting of bare nutrients without phytoplankton in the nonspatial model. The Jacobian matrix of nonspatial system at the equilibrium is

It is obvious that the index of equilibrium is , when the condition holds, which is stable. In particular, when the conditions and hold, if or and , then there is no positive equilibrium in nonspatial system. Under these conditions, equilibrium is locally asymptotically stable. Furthermore, equilibrium is globally asymptotically stable in . Because first quadrant is a positive invariant set according to and , there is no limit cycle because there is no equilibrium in first quadrant.

Based on previous discussion, there exists positive equilibrium in the nonspatial model under some conditions, which is defined by . The Jacobian matrix of nonspatial model at the equilibrium is

From , it is easy to find that the equilibrium is unstable when , which is saddle. When , the index of equilibrium is when the condition holds, and it is locally asymptotically stable using Roth-Hurwitz criterion when the condition holds.

Although the expression of equilibrium can hardly be obtained, the stable behavior of the nonspatial system is determined when some parameters are given in Table 1. In the present paper, our interest is how some factors, such as the nutrient concentration and nutrient cycling effect , affect the system. Hence, the stable behavior of the nonspatial system is analyzed using the graph (see Figure 1). In Figure 1, when the nutrient concentration or nutrient cycling effect increases, there is always a trivial equilibrium consisting of bare nutrients without phytoplankton. In Figures 1(a) and 1(b), when the nutrient concentration is , there is only a trivial steady state in the nonspatial system. When the nutrient concentration , there were two other steady states: one is always unstable (green dashed, saddle), while the other is stable (red solid line, focus). When the nutrient concentration is , the focus disappears, and a node emerges (blue line, node). In Figures 1(c) and 1(d), a similar analysis is obtained, and the difference among Figures 1(a), 1(b), 1(c), and 1(d) is indicated by the grey zone. The trivial steady state and saddle coexist in the grey zone. In the following discussion, the nontrivial homogeneous steady state (focus or node) is defined by .

3.2. Stable Behavior of the Spatial System

In this section, we consider the sensitivity of the system (2a), (2b) to change in the parameter values. A linear analysis technique is employed to focus on the parameters essential for the system behavior [49]. Our interest is how the nutrient concentration , nutrient cycling effect , sinking velocity , and the diffusion rate of phytoplankton affect the system. Symmetry breaking occurred when the stationary homogeneous solution, , is linearly unstable to small spatial perturbations in the presence of diffusion and advection, but that is linearly stable to perturbations in the absence of the diffusion and advection terms. To analyze the spatial system and determine how a small heterogeneous perturbation of the homogeneous steady state developed within a large time period, the following perturbation is considered [41]: where is the perturbation growth rate, is the wave-number, and is an imaginary unit (). Substituting expression (5) into (2a), (2b) and neglecting all nonlinear terms in and , the following characteristic equation is obtained for the eigenvalues : where the elements of the Jacobian determinant of the nonspatial system are taken at the stationary homogeneous solution , as follows:

By , it is not difficult to find that , , and , when is positive equilibrium. When , . The characteristic equation (8) can be described as where and . In the previous analysis, the parameters , , and are allowed to vary, but the other parameters are fixed in Table 1. In Figure 2(a), the value of is given when the parameters and are changed. In addition, in Figure 2(b), in zone III, , so for ; in zone II, , but for , so ; in zone I, , and for . In zone I, to determine the sign of for different values of and , we need to analyze further because of the for .

From expression (8), we can obtain where and .

To analyze the spatial system, the real and imaginary parts of the eigenvalues must be obtained, which are described as follows:where . The solution is stable when the real parts of all eigenvalues are less than zero; that is, . The solution is unstable when one of the real parts with a finite wave number is greater than zero at least. The critical point is got when . However, the analytical expression for the critical point is difficult to be obtained. Indeed, we only need to consider the maximum value of . Thus, the critical condition can be obtained using , as follows:

In expression (11), the solution is unstable if the right-hand side of the equal sign is always less than zero. Otherwise, a necessary condition for expression (11) to hold is that . We consider the following case: the nutrient concentration is allowed to vary, but the values of other parameters are in Table 1. Then, the sinking velocity is a function related to the nutrient concentration and the wave number . By Figure 2(b), if the diffusion rate of phytoplankton, , is larger than 0.2 cm2·s−1, then . Thus, the right-hand side of expression (11) is positive within , and the sinking velocity has a minimum at the point . Figure 3(a) confirms expression (11). In Figure 3(a), the neutral curve is convex with a unique minimum in the range .

The effects of sinking velocity and nutrient concentration on the behavior of system (2a), (2b) are shown in Figure 3(b), which shows the transition from a no phytoplankton state through a banded phytoplankton state to a homogeneous phytoplankton state when the nutrient concentration increases and the sinking velocity is fixed. However, when the nutrient concentration is fixed, a homogeneous phytoplankton state becomes a banded phytoplankton state with the increase of sinking velocity , that is, a banded self-organized spatial pattern emerges because of the sinking velocity .

In zone I of Figure 3(b), there is only a trivial steady state, that is, no phytoplankton. In zone II, the steady state is stable, while the steady state becomes unstable in zone III because of the effect of sinking velocity of phytoplankton, which is further confirmed by Figure 3(c), where  (m·day−1) and . It is obvious that the maximal real part of is larger than zero when the sinking velocity of phytoplankton is larger than the critical value of the sinking velocity of phytoplankton; that is, the instability of the steady state will occur. The imaginary value of is not equal to zero. In Figure 3(b), the red line represents the critical value of the sinking velocity of phytoplankton. The critical value of the sinking velocity of phytoplankton increases with the increase of nutrient concentration ; that is, when the nutrient concentration increases, the sinking velocity cannot affect the stability of the stable state.

3.3. Effects of the Parameters on the Wavelength of the Banded Pattern

Section 3.2 has described how a banded phytoplankton state emerges when the sinking velocity has reached a certain critical value. The determination of the wavelength of the banded pattern is a key issue. In particular, how do the parameters, such as the nutrient concentration , and the sinking velocity , affect the change of wavelength? The relationships among the wavelength, the nutrient concentration , and the sinking velocity are shown in Figure 4, which shows that the wavelength increases when the sinking velocity exceeds the critical value , but the wavelength decreases when the nutrient concentration increases.

3.4. The Simulation

In the previous sections, we discussed the effects of parameters, including the nutrient concentration and the sinking velocity , on the system (2a), (2b). In this section, we discuss the numerical solution of the system (2a), (2b) in one-dimensional and two-dimensional spaces. In a one-dimensional space, a periodic boundary condition is employed, and system (2a), (2b) is solved on a rectangular spatial grid of 1 × 200 points. In two-dimensional space, system (2a), (2b) is studied in a horizontal (x, z)-plane with zero-flux boundary conditions (left and right) and periodic boundary conditions (top and bottom), which is solved on a rectangular spatial grid of 100 × 300 points. The initial conditions comprise a homogeneous state which is randomly perturbed. Furthermore, we assume that the diffusion rate of phytoplankton is larger or smaller than that of nutrients because of the viscosity and living of phytoplankton. Of course, it is also feasible that the diffusion rate of phytoplankton is equal to that of nutrients.

Firstly, the one-dimensional solution of system (2a), (2b) is shown in Figure 5. In Figure 5, we consider a vertical water column, where the depth of the water column is 120 m and the time is 600 days. We found that oscillations did occur; that is, the stable state became unstable because of spatial effects. Figure 6 shows the analysis of relationship between nutrients and phytoplankton further. In Figure 6(a), the relationship between the spatial distributions of nutrients and phytoplankton is given at the 600th day, which shows that nutrient concentration reaches the minimal value when the density of phytoplankton reaches the maximal value. The nutrient concentration affects the density of phytoplankton, and the density of phytoplankton increases with the increases of nutrient concentration. Thus, eutrophication may explain phytoplankton blooms. Furthermore, the effects of phytoplankton sinking on the relative maxima for nutrients and phytoplankton are shown in Figure 6(b). The relative maxima of phytoplankton increase with the increase of sinking velocity, whereas the relative maxima of nutrients decrease with the increase of sinking velocity. Therefore, the sinking flux has an important role in the increase of the density of phytoplankton.

To further analyze the dynamic behavior of system (2a), (2b), we consider the solution of system (2a), (2b) in two-dimensional space. The band patterns are observed in the field, as shown in Figure 7. Figures 7(a), 7(b), and 7(c) show the patterns of phytoplankton at the 1000th day in the two-dimensional space. As discussed in Section 3.2, our numerical results confirm the predictions of the linear analysis that a band pattern of phytoplankton occurs if the nutrient concentration and the sinking velocity of phytoplankton satisfy some conditions. Figure 7(a) shows the emergence of parallel and crossed patterns, which indicate that band patterns with different speeds coexist in the system (2a), (2b) where the wavelength of patterns is different. By contrast, the patterns in Figure 7(b) are much more regular and almost parallel. In the real world, the sinking velocity of phytoplankton varies at different spatial points. Thus, to add more realism to the system, we forced the model to undergo periodic changes in the sinking velocity of phytoplankton, that is, , where denotes the width of water column. The results are shown in Figure 7(c), where the values of the parameters are the same as those used in Figures 7(a) and 7(b), except the sinking velocity . However, it is obvious that the patterns in Figure 7(c) are very different from those in Figures 7(a) and 7(b). Thus, the sinking flux of phytoplankton has an important role in the system.

4. Discussion and Conclusion

Banded patterns have been described in several resource-limited ecosystems around the world. In the real world, numerous population patterns have been observed, including banded vegetation, patches, and spiral waves, which can be regular or irregular. Physical factors may cause these types of pattern, such as wind, water flow, and turbulence. Internal factors in populations also force these patterns to occur.

In the present study, we used a nutrient-plankton model with both diffusion and advection to investigate the interaction between nutrient and plankton. Our model was simple because it was only an abstraction of real-world phenomena but the model reproduced many features of real-world phenomena. Our explanation focuses on a predator-prey interaction between phytoplankton and their nutrient source. In particular, how do the sinking of phytoplankton and the input of nutrients affect the interaction? Our analytical results showed that the homogeneous steady state became unstable because of the sinking of phytoplankton. The critical value of the sinking phytoplankton led to an instability in the homogeneous steady state, which depended on the input of nutrients. Our numerical results showed that the homogeneous steady state was unstable against small spatially heterogeneous perturbations.

Figure 1(b) shows that when the nutrient concentration increased beyond a critical value, the increase in the concentration of phytoplankton was stable; that is, the concentration of phytoplankton tended toward a certain stable state. Spatial effects did not influence the stable state when the sinking flux was below a critical value, as shown in Figure 3(b). Thus, an abundance of nutrient inputs flowed into the system, which led to the high-level reproduction of phytoplankton, which may trigger phytoplankton blooms.

Figures 5 and 7 show that oscillation could occur because of the sinking flux. In particular, Figure 5 shows that both spatial and temporal oscillations arose in the nutrients and the phytoplankton. The sinking of phytoplankton can also lead to the increase in the phytoplankton density and wavelength when the input of nutrients is fixed. It is possible that the phytoplankton sinks from the surface of water until it reaches a depth where the nutrient conditions are suitable for growth. Figure 6(a) shows that the relationship between phytoplankton and nutrients is mutually constrained. Thus, abundant nutrition leads to the mass propagation of phytoplankton, which consumes large amounts of nutrient, thereby depleting the nutrient levels. Thus, the sinking of phytoplankton and the input of nutrients can change the spatial distribution of phytoplankton under these conditions, which may promote the increase of phytoplankton density. In particular, eutrophication may promote phytoplankton blooms.

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 Key Program of Zhejiang Provincial Natural Science Foundation of China (Grant no. LZ12C03001), the National Natural Science Foundation of China (Grant no. 31170338), and the National Key Basic Research Program of China (973 Program, Grant no. 2012CB426510).