# Nonlinear Dynamics of a Nutrient-Plankton Model

Authors: Yapei Wang; Min Zhao; Chuanjun Dai; Xinhong Pan
Journal: Abstract and Applied Analysis
 (2014)

---

## 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.


---

## Body

## 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 [4–6]. 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 [19–21]. 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 [25–30] 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 Section3, 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:(1a)∂N∂t=input-uptake+recycling+mixing,(1b)∂P∂t=growth-mortality-predation-sinking+mixing,where N is the nutrient density and P 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:(2a)∂N∂t=f(N,P)+dNΔN=k(N0-N)-αβNHN+NP+εmP+dNΔN,(2b)∂P∂t=g(N,P)-v∂P∂z+dPΔP=βNHN+NP-mP-fPPHP+P-v∂P∂z+dPΔP,where a vertical water column is considered. Let z indicate the depth in the water column; x is the width in the water column. For vertical mixing, we assume that N0 is a constant concentration, which includes the nutrient input flowing into the system and the nutrient from below the mixed layer, k is the rate of exchange between the lower and upper layers, α is the nutrient content of phytoplankton, β denotes the maximum growth rate of phytoplankton, HN is the half-saturation constant for nutrients, HP is the half-saturation constant for phytoplankton, fP denotes the maximum  predation rate of zooplankton on phytoplankton, m is the mortality of phytoplankton, ε is the proportion of nutrients in dead phytoplankton that is recycled, v is the sinking velocity of phytoplankton, and dN and dP 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, 36–38].Table 1
Parameter values used.
Symbol
Value
Unit
N
0
0.5
g·m−3
k
0.08
day−1
α
0.02
dimensionless
β
0.5
day−1
ε
0.01
dimensionless
m
0.24
day−1
f
P
2
g·m−2·day−1
H
N
0.005
g·m−2
H
P
4
g·m−2
v
1.008
m·day−1
d
N
1.038
m2·day−1
d
P
1.038
m2·day−1
Note: parameter valueε was estimated, parameter value HP was estimated based on de Angelis et al. [36], and the other parameter values were derived from previous studies [19, 36–38].


## 3. Results

### 3.1. Stable Behavior of the Nonspatial System

In the  nonspatial system (i.e., system (2a), (2b) without spatial derivatives), according to f(x,y)=0 and g(x,y)=0, vertical isocline and horizontal isocline can be obtained, respectively, as follows. Vertical isocline l1: P=(k(N0-N)(HN+N))/((αβ-εm)N-εmHN). Horizontal isocline l2:P=(fP(HN+N)/((β-m)N-mHN))-HP. Obviously, the line, N=εmHN/(αβ-εm), is asymptote of vertical isocline l1, and the line, N=(mHN)/(β-m), is asymptote of horizontal isocline l2. In the following discussion, it is assumed that the condition β>m always holds; otherwise phytoplankton become extinct eventually. For vertical isocline l1, P′=((εm-αβ)N2+2εmHNN+(εmHN-αβN0)HN)/((εm-αβ)N+εmHN)2, which is derivative. There are two roots in P′=0 when the condition εm>αβ holds, root is   N=(εmHN±αβHN(εmHN+(εm-αβ)N0))/(αβ-εm). Then, it is  obvious that the asymptote, N=εmHN/(αβ-εm), is on the left of line N=0. Therefore, the vertical isocline l1 is continuous when N>0. And P>0 holds when N>N0; P<0 holds when 0<N<N0. From the horizontal isocline l2, if fP>HP(β-m), then there is a positive equilibrium in the nonspatial system at least; if fP<HP(β-m), then there is a positive equilibrium in the nonspatial system at least when the condition mHN/(β-m)>N0 holds.When the conditionεm<αβ holds, there is no root in P′=0 if N0>(εmHN/(αβ-εm)). From P′, vertical isocline l1 is monotone decreasing when N>(εmHN/(αβ-εm)), and P<0 holds when N∈(0,(εmHN/(αβ-εm))); P>0 holds when N∈((εmHN/(αβ-εm)),N0). According to horizontal isocline l2, if fP>HP(β-m), then there is a positive equilibrium in the nonspatial system at least when the condition (m/β)ε<α<ε holds; if fP<HP(β-m), then there is no positive equilibrium in the nonspatial system when the condition (mHN/(β-m))>N0 holds.When the conditionεm<αβ holds, there are two roots in P′=0 if N0<(εmHN/(αβ-εm)). From P′, vertical isocline l1 is monotone increasing, when N∈(N0,(εmHN/(αβ-εm))). And P<0 holds, when N∈(0,N0)∪((εmHN/(αβ-εm)),+∞); P>0 holds, when N∈(N0,(εmHN/(αβ-εm))). According to horizontal isocline l2, if fP>HP(β-m), then there is a positive equilibrium state in the nonspatial system at least when the condition (m/β)ε<α<ε holds; if fP<HP(β-m), 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,E0=(N0,0), consisting of bare nutrients without phytoplankton in the nonspatial model. The Jacobian matrix of nonspatial system at the equilibrium E0=(N0,0) is
(3)A=(-kεm-αβN0HN+N00(β-m)N0-mHNHN+N0-fPHP).It is obvious that the index of equilibriumE0 is +1, when the condition βN0HP<(mHP+fP)(N0+HN) holds, which is stable. In particular, when the conditions εm<αβ and fP<HP(β-m) hold, if (mHN/(β-m))>N0>(εmHN/(αβ-εm)) or N0<(εmHN/(αβ-εm)) and α>ε, then there is no positive equilibrium in nonspatial system. Under these conditions, equilibrium E0 is locally asymptotically stable. Furthermore, equilibrium E0 is globally asymptotically stable in Ω=[0,+∞)×[0,+∞]. Because first quadrant is a positive invariant set according to f(x,y) and g(x,y), 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 byE*=(N*,P*). The Jacobian matrix of nonspatial model at the equilibrium E*=(N*,P*) is
(4)T=(-k-αβHNP*(HN+N*)2k(N*-N0)P*βHNP*(HN+N*)2fPP*(HP+P*)2).FromT, it is easy to find that the equilibrium E*=(N*,P*) is unstable when N*>N0, which is saddle. When N*<N0, the index of equilibrium E* is +1 when the condition kβHN(N0-N*)(HP+P*)2>(k(HN+N*)2+αβHNP*)fPP* holds, and it is locally asymptotically stable using Roth-Hurwitz criterion when the condition (fPP*/(HP+P*)2)<k+(αβHNP*/(HN+N*)2) holds.Although the expression of equilibriumE* 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 N0 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 N0 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 0≤N0<0.489, there is only a trivial steady state in the nonspatial system. When the nutrient concentration 0.489<N0<0.609, 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 0.609<N0<2, 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 E*=(N*,P*).The existence and stability of an equilibrium versus the parameterN0 (a, b) and parameter ε (c, d). The other parameters are given in Table 1.
(a)
(b)
(c)
(d)


### 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 N0, nutrient cycling effect ε, sinking velocity v, and the diffusion rate of phytoplankton dP affect the system.  Symmetry breaking occurred when the stationary homogeneous solution, E*=(N*,P*), 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]:
(5)(NP)=(N*P*)+δ(N0P0)exp⁡(λt+ikz)+c.c.+O(ε2),
where λ is the perturbation growth rate, k is the wave-number, and i is an imaginary unit (i2=-1). Substituting expression (5) into (2a), (2b) and neglecting all nonlinear terms in N and P, the following characteristic equation is obtained for the eigenvalues λ:
(6)|a11-k2dN-λa12a21a22-k2dP-ivk-λ|=0,
where the elements of the Jacobian determinant of the nonspatial system are taken at the stationary homogeneous solution E*=(N*,P*), as follows:
(7)J=(-k-αβHNP*(HN+N*)2k(N*-N0)P*βHNP*(HN+N*)2fPP*(HP+P*)2)=(a11a12a21a22).ByJ, it is not difficult to find that a11<0, a21>0, and a22>0, when E* is positive equilibrium. When N*<N0, a12<0. The characteristic equation (8) can be described as
(8)λ2-(trk-ivk)λ+Δk+ivk(k2dN-a11)=0,
where trk=(a11+a22)-(dN+dP)k2 and Δk=a11a22-a21a12-k2(a11dP+a22dN)+k4dNdP. In the previous analysis, the parameters N0,ε, v, and dP are allowed to vary, but the other parameters are fixed in Table 1. In Figure 2(a), the value of a22 is given when the   parameters N0 and ε are  changed. In addition, in Figure 2(b), in zone III, a11dP+a22dN<0, so Δk>0 for k>0; in zone II, a11dP+a22dN>0, but min⁡(Δk)>0 for k>0, so Δk>0; in zone I, a11dP+a22dN>0, and min⁡(Δk)<0 for k>0. In zone I, to determine the sign of Δk for different values of N0 and dP, we need to analyze Δk further because of the min⁡(Δk)<0 for k>0.(a) An illustration ofa22>0 versus ε and N0. The other parameters are given in Table 1. (b) An illustration of the sign of Δk versus dP and N0, v=0.6 (m·day−1).
(a)
(b)From expression (8), we can obtain
(9)λ=12[trk-ivk±Φ+iΘ],
where Φ=trk2-v2k2-4Δk and Θ=-2(trk+2(k2dN-a11))vk.To analyze the spatial system, the real and imaginary parts of the eigenvalues must be obtained, which are described as follows:(10a)Re(λ)=12[trk+j12(Φ2+Θ2+Φ)],(10b)Im⁡(λ)=12[-vk+jsign⁡(Φ)12(Φ2+Θ2-Φ)],where j=±1. The solution is stable when the real parts of all eigenvalues are less than zero; that is, Re(λ)<0. The solution is unstable when one of the real parts with a finite wave number k>0 is greater than zero at least. The critical point is got when Re(λ)=0. However, the analytical expression for the critical point is difficult to be obtained. Indeed, we only need to consider the maximum value of Re(λ). Thus, the critical condition can be obtained using Re(λ)=0, as follows:
(11)v2=trk2Δk(k2dN-a11)(a22-dPk2)k2.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 Δk(a22-dPk2)>0. We consider the following case: the nutrient concentration N0 is allowed to vary, but the values of other parameters are in Table 1. Then, the sinking velocity v is a function related to the nutrient concentration N0 and the wave number k. By Figure 2(b), if the diffusion rate of phytoplankton, dP, is larger than 0.2 cm2·s−1, then Δk>0. Thus, the right-hand side of expression (11) is positive within 0<k<a22/dP, and the sinking velocity v has a minimum vc at the point k=kc. Figure 3(a) confirms expression (11). In Figure 3(a), the neutral curve is convex with a unique minimum in the range 0<k<a22/dP.(a) A typical neutral curvev, defined using expression (11) for different values of N0. (b) Numerical calculation of the stability on (v,N0) space. (c) An illustration of the dispersion relation (Re(λ),Im⁡(λ) versus the wave number k), where dP=0.3 (m2·day−1) and the other parameters are given in Table 1. Blue line: v=0.6>vc≈0.358 (m·day−1); red line: v=vc≈0.358 (m·day−1); green line: v=0.1<vc≈0.358 (m·day−1).
(a)
(b)
(c)The effects of sinking velocityv and nutrient  concentrationN0 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 N0 increases and the sinking velocity v is fixed. However, when the nutrient concentration N0 is fixed, a homogeneous phytoplankton state becomes a banded phytoplankton state with the increase of sinking velocity v, that is, a banded self-organized spatial pattern emerges because of the sinking velocity v.In zone I of Figure3(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 vc≈0.358 (m·day−1) and kc≈0.283. 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 N0; that is, when the nutrient concentration N0 increases, the sinking velocity v cannot affect the stability of the stable state.


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

Section3.2 has described how a banded phytoplankton state emerges when the sinking velocity v 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 N0, and the sinking velocity v, affect the change of wavelength? The relationships among the wavelength, the nutrient concentration N0, and the sinking velocity v are shown in Figure 4, which shows that the wavelength increases when the sinking velocity v exceeds the critical value vc, but the wavelength decreases when the nutrient concentration N0 increases.Figure 4
An illustration of the variation in the pattern wavelength withv and N0, where the symbol (■) represents the critical value of the sinking velocity v, dP=0.3 (m2·day−1), and the other parameters are given in Table 1.


### 3.4. The Simulation

In the previous sections, we discussed the effects of parameters, including the nutrient concentrationN0 and the sinking velocity v, 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.Model simulations based on depth and time: (a) nutrient density (unit g·m−3); (b) phytoplankton density (unit g·m−3), where dP=1.2(m2·day−1), v=1.2(m·day−1), and the other parameters are given in Table 1.
(a)
(b)(a) Spatial distributions of the nutrient density and phytoplankton density on the 600th day, where the left vertical axis denotes the nutrient density, the right vertical axis denotes the phytoplankton density,dP=0.3 (m2·day−1), v=1.2 (m·day−1), and the other parameters are given in Table 1. (b) Relationship between nutrient density, phytoplankton density, and the sinking velocity v, where the right vertical axis denotes the maximal nutrient density, and the left vertical axis denotes the maximal phytoplankton density. The purple field shows that the equilibrium state is stable, where N0=0.5 (g·m−3) and dP=0.3 (m2·day−1).
(a)
(b)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 N0 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, v=v*+A×sin(4π×x/L), where L 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 v. 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.Simulation in two-dimensional space. The figures show the density levels of the phytoplankton on the 1000th day, where the width is 100 m and the depth is 300 m. (a)dP=0.3 (m2·day−1) and v=0.9 (m·day−1). (b) dP=1.038 (m2·day−1) and v=1.2 (m·day−1). (c) dP=1.038 (m2·day−1) and v=0.9+0.25*sin(4*π*x/100). The other parameters are given in Table 1.
(a)
(b)
(c)


## 3.1. Stable Behavior of the Nonspatial System

In the  nonspatial system (i.e., system (2a), (2b) without spatial derivatives), according to f(x,y)=0 and g(x,y)=0, vertical isocline and horizontal isocline can be obtained, respectively, as follows. Vertical isocline l1: P=(k(N0-N)(HN+N))/((αβ-εm)N-εmHN). Horizontal isocline l2:P=(fP(HN+N)/((β-m)N-mHN))-HP. Obviously, the line, N=εmHN/(αβ-εm), is asymptote of vertical isocline l1, and the line, N=(mHN)/(β-m), is asymptote of horizontal isocline l2. In the following discussion, it is assumed that the condition β>m always holds; otherwise phytoplankton become extinct eventually. For vertical isocline l1, P′=((εm-αβ)N2+2εmHNN+(εmHN-αβN0)HN)/((εm-αβ)N+εmHN)2, which is derivative. There are two roots in P′=0 when the condition εm>αβ holds, root is   N=(εmHN±αβHN(εmHN+(εm-αβ)N0))/(αβ-εm). Then, it is  obvious that the asymptote, N=εmHN/(αβ-εm), is on the left of line N=0. Therefore, the vertical isocline l1 is continuous when N>0. And P>0 holds when N>N0; P<0 holds when 0<N<N0. From the horizontal isocline l2, if fP>HP(β-m), then there is a positive equilibrium in the nonspatial system at least; if fP<HP(β-m), then there is a positive equilibrium in the nonspatial system at least when the condition mHN/(β-m)>N0 holds.When the conditionεm<αβ holds, there is no root in P′=0 if N0>(εmHN/(αβ-εm)). From P′, vertical isocline l1 is monotone decreasing when N>(εmHN/(αβ-εm)), and P<0 holds when N∈(0,(εmHN/(αβ-εm))); P>0 holds when N∈((εmHN/(αβ-εm)),N0). According to horizontal isocline l2, if fP>HP(β-m), then there is a positive equilibrium in the nonspatial system at least when the condition (m/β)ε<α<ε holds; if fP<HP(β-m), then there is no positive equilibrium in the nonspatial system when the condition (mHN/(β-m))>N0 holds.When the conditionεm<αβ holds, there are two roots in P′=0 if N0<(εmHN/(αβ-εm)). From P′, vertical isocline l1 is monotone increasing, when N∈(N0,(εmHN/(αβ-εm))). And P<0 holds, when N∈(0,N0)∪((εmHN/(αβ-εm)),+∞); P>0 holds, when N∈(N0,(εmHN/(αβ-εm))). According to horizontal isocline l2, if fP>HP(β-m), then there is a positive equilibrium state in the nonspatial system at least when the condition (m/β)ε<α<ε holds; if fP<HP(β-m), 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,E0=(N0,0), consisting of bare nutrients without phytoplankton in the nonspatial model. The Jacobian matrix of nonspatial system at the equilibrium E0=(N0,0) is
(3)A=(-kεm-αβN0HN+N00(β-m)N0-mHNHN+N0-fPHP).It is obvious that the index of equilibriumE0 is +1, when the condition βN0HP<(mHP+fP)(N0+HN) holds, which is stable. In particular, when the conditions εm<αβ and fP<HP(β-m) hold, if (mHN/(β-m))>N0>(εmHN/(αβ-εm)) or N0<(εmHN/(αβ-εm)) and α>ε, then there is no positive equilibrium in nonspatial system. Under these conditions, equilibrium E0 is locally asymptotically stable. Furthermore, equilibrium E0 is globally asymptotically stable in Ω=[0,+∞)×[0,+∞]. Because first quadrant is a positive invariant set according to f(x,y) and g(x,y), 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 byE*=(N*,P*). The Jacobian matrix of nonspatial model at the equilibrium E*=(N*,P*) is
(4)T=(-k-αβHNP*(HN+N*)2k(N*-N0)P*βHNP*(HN+N*)2fPP*(HP+P*)2).FromT, it is easy to find that the equilibrium E*=(N*,P*) is unstable when N*>N0, which is saddle. When N*<N0, the index of equilibrium E* is +1 when the condition kβHN(N0-N*)(HP+P*)2>(k(HN+N*)2+αβHNP*)fPP* holds, and it is locally asymptotically stable using Roth-Hurwitz criterion when the condition (fPP*/(HP+P*)2)<k+(αβHNP*/(HN+N*)2) holds.Although the expression of equilibriumE* 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 N0 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 N0 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 0≤N0<0.489, there is only a trivial steady state in the nonspatial system. When the nutrient concentration 0.489<N0<0.609, 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 0.609<N0<2, 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 E*=(N*,P*).The existence and stability of an equilibrium versus the parameterN0 (a, b) and parameter ε (c, d). The other parameters are given in Table 1.
(a)
(b)
(c)
(d)


## 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 N0, nutrient cycling effect ε, sinking velocity v, and the diffusion rate of phytoplankton dP affect the system.  Symmetry breaking occurred when the stationary homogeneous solution, E*=(N*,P*), 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]:
(5)(NP)=(N*P*)+δ(N0P0)exp⁡(λt+ikz)+c.c.+O(ε2),
where λ is the perturbation growth rate, k is the wave-number, and i is an imaginary unit (i2=-1). Substituting expression (5) into (2a), (2b) and neglecting all nonlinear terms in N and P, the following characteristic equation is obtained for the eigenvalues λ:
(6)|a11-k2dN-λa12a21a22-k2dP-ivk-λ|=0,
where the elements of the Jacobian determinant of the nonspatial system are taken at the stationary homogeneous solution E*=(N*,P*), as follows:
(7)J=(-k-αβHNP*(HN+N*)2k(N*-N0)P*βHNP*(HN+N*)2fPP*(HP+P*)2)=(a11a12a21a22).ByJ, it is not difficult to find that a11<0, a21>0, and a22>0, when E* is positive equilibrium. When N*<N0, a12<0. The characteristic equation (8) can be described as
(8)λ2-(trk-ivk)λ+Δk+ivk(k2dN-a11)=0,
where trk=(a11+a22)-(dN+dP)k2 and Δk=a11a22-a21a12-k2(a11dP+a22dN)+k4dNdP. In the previous analysis, the parameters N0,ε, v, and dP are allowed to vary, but the other parameters are fixed in Table 1. In Figure 2(a), the value of a22 is given when the   parameters N0 and ε are  changed. In addition, in Figure 2(b), in zone III, a11dP+a22dN<0, so Δk>0 for k>0; in zone II, a11dP+a22dN>0, but min⁡(Δk)>0 for k>0, so Δk>0; in zone I, a11dP+a22dN>0, and min⁡(Δk)<0 for k>0. In zone I, to determine the sign of Δk for different values of N0 and dP, we need to analyze Δk further because of the min⁡(Δk)<0 for k>0.(a) An illustration ofa22>0 versus ε and N0. The other parameters are given in Table 1. (b) An illustration of the sign of Δk versus dP and N0, v=0.6 (m·day−1).
(a)
(b)From expression (8), we can obtain
(9)λ=12[trk-ivk±Φ+iΘ],
where Φ=trk2-v2k2-4Δk and Θ=-2(trk+2(k2dN-a11))vk.To analyze the spatial system, the real and imaginary parts of the eigenvalues must be obtained, which are described as follows:(10a)Re(λ)=12[trk+j12(Φ2+Θ2+Φ)],(10b)Im⁡(λ)=12[-vk+jsign⁡(Φ)12(Φ2+Θ2-Φ)],where j=±1. The solution is stable when the real parts of all eigenvalues are less than zero; that is, Re(λ)<0. The solution is unstable when one of the real parts with a finite wave number k>0 is greater than zero at least. The critical point is got when Re(λ)=0. However, the analytical expression for the critical point is difficult to be obtained. Indeed, we only need to consider the maximum value of Re(λ). Thus, the critical condition can be obtained using Re(λ)=0, as follows:
(11)v2=trk2Δk(k2dN-a11)(a22-dPk2)k2.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 Δk(a22-dPk2)>0. We consider the following case: the nutrient concentration N0 is allowed to vary, but the values of other parameters are in Table 1. Then, the sinking velocity v is a function related to the nutrient concentration N0 and the wave number k. By Figure 2(b), if the diffusion rate of phytoplankton, dP, is larger than 0.2 cm2·s−1, then Δk>0. Thus, the right-hand side of expression (11) is positive within 0<k<a22/dP, and the sinking velocity v has a minimum vc at the point k=kc. Figure 3(a) confirms expression (11). In Figure 3(a), the neutral curve is convex with a unique minimum in the range 0<k<a22/dP.(a) A typical neutral curvev, defined using expression (11) for different values of N0. (b) Numerical calculation of the stability on (v,N0) space. (c) An illustration of the dispersion relation (Re(λ),Im⁡(λ) versus the wave number k), where dP=0.3 (m2·day−1) and the other parameters are given in Table 1. Blue line: v=0.6>vc≈0.358 (m·day−1); red line: v=vc≈0.358 (m·day−1); green line: v=0.1<vc≈0.358 (m·day−1).
(a)
(b)
(c)The effects of sinking velocityv and nutrient  concentrationN0 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 N0 increases and the sinking velocity v is fixed. However, when the nutrient concentration N0 is fixed, a homogeneous phytoplankton state becomes a banded phytoplankton state with the increase of sinking velocity v, that is, a banded self-organized spatial pattern emerges because of the sinking velocity v.In zone I of Figure3(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 vc≈0.358 (m·day−1) and kc≈0.283. 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 N0; that is, when the nutrient concentration N0 increases, the sinking velocity v cannot affect the stability of the stable state.


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

Section3.2 has described how a banded phytoplankton state emerges when the sinking velocity v 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 N0, and the sinking velocity v, affect the change of wavelength? The relationships among the wavelength, the nutrient concentration N0, and the sinking velocity v are shown in Figure 4, which shows that the wavelength increases when the sinking velocity v exceeds the critical value vc, but the wavelength decreases when the nutrient concentration N0 increases.Figure 4
An illustration of the variation in the pattern wavelength withv and N0, where the symbol (■) represents the critical value of the sinking velocity v, dP=0.3 (m2·day−1), and the other parameters are given in Table 1.


## 3.4. The Simulation

In the previous sections, we discussed the effects of parameters, including the nutrient concentrationN0 and the sinking velocity v, 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.Model simulations based on depth and time: (a) nutrient density (unit g·m−3); (b) phytoplankton density (unit g·m−3), where dP=1.2(m2·day−1), v=1.2(m·day−1), and the other parameters are given in Table 1.
(a)
(b)(a) Spatial distributions of the nutrient density and phytoplankton density on the 600th day, where the left vertical axis denotes the nutrient density, the right vertical axis denotes the phytoplankton density,dP=0.3 (m2·day−1), v=1.2 (m·day−1), and the other parameters are given in Table 1. (b) Relationship between nutrient density, phytoplankton density, and the sinking velocity v, where the right vertical axis denotes the maximal nutrient density, and the left vertical axis denotes the maximal phytoplankton density. The purple field shows that the equilibrium state is stable, where N0=0.5 (g·m−3) and dP=0.3 (m2·day−1).
(a)
(b)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 N0 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, v=v*+A×sin(4π×x/L), where L 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 v. 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.Simulation in two-dimensional space. The figures show the density levels of the phytoplankton on the 1000th day, where the width is 100 m and the depth is 300 m. (a)dP=0.3 (m2·day−1) and v=0.9 (m·day−1). (b) dP=1.038 (m2·day−1) and v=1.2 (m·day−1). (c) dP=1.038 (m2·day−1) and v=0.9+0.25*sin(4*π*x/100). The other parameters are given in Table 1.
(a)
(b)
(c)


## 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.Figure1(b) shows that when the nutrient concentration N0 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.Figures5 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.



---

*Source: 451757-2014-01-16.xml*