Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 815682, 11 pages

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

## Spatiotemporal Complexity of the Nutrient-Phytoplankton Model

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

Received 25 September 2014; Revised 3 January 2015; Accepted 13 January 2015

Academic Editor: Giovanni Garcea

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

#### Abstract

We consider the mathematical formulation, analysis, and numerical solution of a nonlinear system of nutrient-phytoplankton, which consists of a series of reaction-advection-diffusion equations. We derive the critical conditions for Turing instability without an advection term and define the range of Turing instability with the change of nutrient concentration. We show that horizontal movement of phytoplankton could influence the system and that it is unstable when the horizontal velocity exceeds a critical value. We also compare reaction-diffusion equations with reaction-advection-diffusion equations through simulations, with spotted, banded, and crenulate patterns produced from our model. We found that different spatial constructions could occur, impacted by the diffusion and sinking of nutrients and phytoplankton. The new model may help us better understand the dynamics of an aquatic community.

#### 1. Introduction

Environmental resources include clean air and water, forests, soil, climate, and many others. Ecological services exist within and are sustained by the ecosystem. They are generated by the continuous connection between organisms, populations, communities, and their physical and chemical environment [1]. Modeling the dynamics of phytoplankton is of great importance to many aspects of human interest, since phytoplankton provides the basis of the food chain in lakes, seas, and oceans [2]. In some cases lakes, reservoirs, and marine waters may experience plankton or algal blooms [3, 4]. Plankton blooms often lead to disorder of aquatic ecosystems, for example, the decrease of aquatic species, and then damage to the biological diversity. Therefore, research into plankton dynamics has become a topic of wide interest.

The historical work of Turing (1952), showed that complex spatial patterns can be produced from the reaction and diffusion of a number of chemicals [5], a phenomenon now known as diffusion-driven instability. A great deal of research has investigated spatial patterns in predator-prey systems. These arise through diffusion-driven instability and depend on significant differences between predator and prey diffusion coefficients [6, 7]. The effect of equal diffusion coefficients of spatial models has also been well studied with reaction-diffusion equations [8, 9].

Many researchers have investigated the predator-prey model among nutrients, phytoplankton, and zooplankton [10–20]. A series of studies simulated the dynamics of the plankton system. Zhang and Wang [21] studied a nutrient-phytoplankton-zooplankton model in an aquatic environment and showed that phytoplankton blooms may occur even if the nutrient input rate is low. Luo [22] investigated phytoplankton-zooplankton dynamics in periodic environments and showed that the eruption of algal blooms largely results from increasing of nutrient concentration. Sae-jie et al. [23] used a general forced dynamical model of nutrient-phytoplankton. They found that sinusoidal and periodic step function input of nutrients caused phytoplankton blooms with periodic behaviour and that changing the frequency of inputs produced blooms with a wide range of different dynamical behavior.

In recent years, numerous ecologists have paid increasing attention to spatial processes in practical contexts [24]. Klausmeier [25], for example, demonstrated a simple model of plant and water dynamics based on ecologically realistic assumptions, which showed that spatial patterns in ecological systems can result from both self-organization and amplification of an underlying heterogeneity. Other researchers have probed ecosystems, including vegetation [26–28] and phytoplankton systems. Asaeda et al. [29] developed a numerical model incorporating phytoplankton, submerged macrophytes,* Potamogeton pectinatus* L, which belongs to the Monocotyledon, their decomposition process, and the resulting nutrient dynamics in the overlying water, which investigated how external nutrient loading, water temperature, water depth, and retention time of water changed submerged macrophytes and phytoplankton biomass, phosphorous concentrations, and light penetration.

A large number of studies have simulated the processes of patchiness pattern formation [30, 31]. The research may be classified into two groups depending on what is considered as a major cause of these phenomena. One group concentrates on biological factors such as nonlinear growth and grazing interplay between phytoplankton and zooplankton [32]. The other emphasizes the influences of physical factors such as currents, eddies, and the turbulent stirring induced by them [33].

Medvinsky et al. [30] followed the first approach. They developed a model of phytoplankton and zooplankton with equal diffusivities and neglected turbulent stirring. This succeeded in reproducing patchiness-like patterns developed from regular spirals. In their simulation, the initial spiral pattern formation is followed by its collapse, which starts near the center, generating spatiotemporal chaos over the region. They also showed that the processes are almost independent of the initial distributions of the constituents. They found that biological factors play a very important role in the production of plankton patchiness [30].

Most researches following the second contention [31–33] adopt the three-component model comprising nutrients or their substitutes, phytoplankton, and zooplankton [34]. Abraham [31] developed a model consisting of the carrying capacities and population densities of phytoplankton and zooplankton. The maturation time of zooplankton enabled fine structures to emerge in the spatial distributions of the zooplankton population whereas the distribution of the phytoplankton population did not show such fine structures because of their rapid growth. He used the seeded-eddy model [35] as turbulent flows.

Our present research builds on the second contention. We propose a model with a series of reaction-advection-diffusion equations in regard to nutrients and phyto-plankton, which allow spatial phenomena, such as sinking and turbulence, to be described directly, therefore enabling research on spatial structures. The sinking, turbulence, and mixing of phytoplankton have an important influence on the dynamic behavior of the system. Wang et al. [36] developed a relatively simple reaction-diffusion-advection model of the interaction between algae and mussels, which included the diffusive spread of mussels and the advection of algae. Theoretical results also showed the significance of sinking, turbulence, and diffusion [37, 38].

In Section 2 we present the mathematical formulation of the model incorporating the biology in the partial differential equations. In Section 3, we analyze the positive equilibrium state, the necessary conditions of Turing instability are deduced, and we analyze the stable behavior of the nonspatial and spatial systems. We derive the condition where the stationary homogeneous solution is linearly unstable. In Section 4, a number of simulations are discussed. From this we show the key factors and the critical conditions in the system. In Section 5, we discuss the outcomes and draw our final conclusion.

#### 2. Mathematical Formulation of the Model

Fresh water wetland ecosystems depend on the available nutrients. Most often, man can contribute to the nutrient through sewage and agricultural fertilizer runoffs. Phytoplankton growth, in particular, is dependent on light and nutrients [1]. Our model is described aswhere the variables and represent the concentration of nutrients and phytoplankton, respectively, which are functions of the actual time, ; all the parameters are positive constants, with regard to the parameters; is the (constant) input of nutrients from the environment; is the rate of the nutrient concentration loss; is the consumption rate of nutrient by phytoplankton; is the maximum nutrient intake rate for phytoplankton; is the phytoplankton mortality; is the maximum predation rate of zooplankton on phytoplankton; is the half-saturation constant for phytoplankton; is the sinking velocity of phytoplankton, which is caused by turbulent stirring; is the width in the water column; is the depth of water, and are the diffusion rates of nutrients and phytoplankton, respectively, which arise from turbulent flow; and is the usual Laplacian operator in two-dimensional space.

The consumption of nutrients by phytoplankton is represented by the Holling I Functional response function. The implicitly assumed zooplankton grazing on phytoplankton is combined as a nondynamic term, consisting of a Holling II Functional response function.

The term addressing removal of nutrients, , is added to estimate the losses except for the intake by phytoplankton, such as the sedimentation to a deep bottom or a runoff from the system. This term has another meaning in the case of mathematical modelling; if the phytoplankton is finally extinct, the nutrient concentration continues to increase infinitely. To avoid this situation, we add the removal term. However, if the term was omitted, the simulation results for spatial pattern formation are unchanged.

The term of represents the feedback of the bottom sediments [1]. When the bottom phytoplankton is decreased, the bottom sediments will be more easily affected by winds, waves, and bottom zooplankton. Finally, nutrients will be more easily released into the water by turbulent stirring and turbulent flow. It is obvious that the term is nondecreasing in . Furthermore, when increases, it will reach an upper bound due to when all green plants in the bottom disappearing, increasing will not release any more nutrients from the bottom sediments. Generally, ecologists make use of the sigmoid function, and hence is the feedback.

As a response to the interaction between nutrients and phytoplankton, we utilize a Holling type I Functional response function in our model. Holling explained the functional response of predators to prey density and the impact on mimicry and population regulation [39].

It is also important to use a Holling type II Functional response to describe zooplankton grazing on phytoplankton, because this is a nondynamic term. This type of functional response has been broadly applied to represent zooplankton predation in many theoretical studies [30, 39–41] and shows good agreement with experimental data [42, 43].

#### 3. Stability Analysis of the System

##### 3.1. Analysis of the Spatially Uniform System

The system of (1a) and (1b), without spatial derivatives, can be described as

The vertical isocline, where , is , which is continuous when . When is close to zero, is close to the positive infinite. On the other side, when is close to positive infinite, is close to . Because is continuous when , there must exist a point at least that satisfies . If , we can obtain , then it is monotonically decreasing when the . Since when , it follows that when .

The horizontal isocline, where , is , and . When , the line, , is the asymptote of . The horizontal isocline is continuous when . Also, is satisfied when ; is satisfied when .

When the system consists of bare nutrients without phytoplankton in the nonspatial model, , and . From the discussion above, if , there is always a boundary equilibrium state, which can be denoted as .

On the other hand, where and , if , there exists a positive equilibrium state in the nonspatial system when ; whereas there is no positive equilibrium in the nonspatial system when . The positive equilibrium can be written as .

Thus, we confirm that positive equilibrium states exist in the nonspatial system. Of course, it does not mean that there cannot be positive equilibrium states in the nonspatial system when these conditions are not satisfied.

The Jacobian matrix of the nonspatial system at the equilibrium is

The index of equilibrium is when the condition, , , is satisfied, which is stable. According to the eigenvalues of , it is locally asymptotically stable when the above conditions are satisfied. In particular, the equilibrium is unstable, when or .

From the above discussion, we find a positive equilibrium in the nonspatial system under some preconditions, defined by . The Jacobian matrix of the nonspatial system at the equilibrium, is

When , the index of equilibrium is if . It is locally asymptotically stable using the Routh-Hurwitz criteria when the above conditions are satisfied. In particular, the equilibrium is unstable, when , which is a saddle. Figure 1 shows the phase diagram of system which is locally asymptotically stable.