Abstract

A wetland ecosystem is studied theoretically and numerically to reveal the rules of dynamics which can be quite accurate to better describe the observed spatial regularity of tussock vegetation. Mathematical theoretical works mainly investigate the stability of constant steady states, the existence of nonconstant steady states, and bifurcation, which can deduce a standard parameter control relation and in return can provide a theoretical basis for the numerical simulation. Numerical analysis indicates that the theoretical works are correct and the wetland ecosystem can show rich dynamical behaviors not only regular spatial patterns. Our results further deepen and expand the study of dynamics in the wetland ecosystem. In addition, it is successful to display tussock formation in the wetland ecosystem may have important consequences for aquatic community structure, especially for species interactions and biodiversity. All these results are expected to be useful in the study of the dynamic complexity of wetland ecosystems.

1. Introduction

Ecosystems consist of organisms, the environment, and the interactions with each other [1]. Many ecosystems exhibit remarkable complexity in both structure and dynamics, which has been used to predict the response of ecosystems to human interference [2, 3]. Complexity theory has successfully predicted that localized ecological interactions may strongly affect the organization of ecosystems on larger spatial scales, which may induce spatial self-organization. There is a growing body of literature emphasizing the importance of self-organization in determining the spatial complexity of ecosystems. Localized ecological interactions can generate striking large-scale spatial patterns in ecosystems through self-organization [4ā€“7]. In recent decades, regular spatial patterning, which is common in ecosystems, has been a dominant topic in the study of spatial self-organization [8, 9]. In addition, a number of theoretical ecologists have proposed various self-organization mechanisms to explain regular spatial patterns in real ecosystems [10ā€“14]. Many regular spatial patterns can be observed in various ecosystems, including arid ecosystems [12], savanna ecosystems [15], mussel beds [9], coral reefs [16], ribbon forests [17], intertidal mudflats [18], and marsh tussocks [8, 19]. In the study of regular spatial patterns, reaction-diffusion equations [20ā€“22] have proved to be an effective tool and are extensively used.

Tussock formation is a global phenomenon that enhances microtopography and increases biodiversity by adding structure to ecological communities [19]. Recently, many ecologists have paid increasing attention to experimental investigation of Carex stricta [8, 19, 23, 24]. Carex stricta, the tussock sedge, is a species with widespread distribution in fresh water marshes of North America. Lawrence and Zedler [19] studied Carex stricta tussock size in relation to elevation at a range of sites in southern Wisconsin, USA, and tested the effect of five hydroperiods and N+P addition on tussock formation during a three-year mesocosm experiment. Their work is significant to the study of tussock development in relation to environmental factors. van de Koppel and Crain [8] presented an experimental investigation of regular spatial patterning in Carex stricta. When van de Koppel and Crain studied a marsh tussocks ecosystem, they proposed a model to investigate alternative theoretical mechanisms driving regular spatial patterning in general. To obtain the best theoretical explanation of the observed spatial regularity of tussock vegetation, van de Koppel and Crain analyzed and compared these alternative mechanisms, using three simple, spatially explicit models to describe the growth of Carex stricta. These models, which describe the interaction between plants and wrack, have the following general structure [8]:where is plant biomass, is wrack biomass, is a function describing the positive effect of plant biomass on its own growth, is the specific rate of plant senescence, is a function describing the inhibiting effect of wrack on plant growth as a function of plant and wrack biomass, is the decay rate of wrack, and are diffusion constants describing lateral movement of plants and wrack, and is the usual Laplacian operator in two-dimensional space.

van de Koppel and Crain [8] compared three alternative models to examine potential mechanisms driving regularity and found that scale-dependent inhibition provides the best explanation for the observed regular tussock spacing. van de Koppel and Crain emphasized the importance of scale-dependent feedback in the formation of regular spatial patterns in ecosystems. Their work is very important in the study of the spatial distribution of Carex stricta. However, if and vary, what are the dynamics? And what are the distributions of plants and wrack? Furthermore, if other parameters vary, what are the resulting changes in the dynamics and spatial distribution of vegetation? These problems are what we are interested in. One model proposed by van de Koppel and Crain [8] was investigated in detail using analysis theory. Although the model is relatively simple, it is important and meaningful because it can reveal the theoretical mechanism that explains tussock formation in wetland ecosystems. The model can be expressed as follows:where is added to the inhibition term to lower inhibition as increases, is the level of plant biomass where inhibition is reduced by half, and is an inhibition coefficient [8].

Model (2) was analyzed under the following nonzero initial conditions:and the zero-flux boundary conditions:where is the outward unit normal vector of the boundary .

The rest of the paper is organized as follows. In Section 2, the properties of constant steady-state solutions of (2) are analyzed, including dissipation, persistence properties, and local and global stability. In Section 3, we give a priori estimates of the positive solutions and establish the existence and nonexistence of nonconstant positive steady-state solutions of this reaction-diffusion system. In Section 4, we perform bifurcation analysis, including Hopf bifurcation and Turing bifurcation. In Section 5, some numerical simulations are presented to verify the feasibility of the theoretical results. The paper ends with some conclusions and a brief discussion.

2. Constant Steady States

Obviously, is the trivial constant steady state of system (2). If , there exists a unique positive constant steady state , andwhere . Furthermore, if , there are no positive constant steady states of system (2).

2.1. Dissipation

Theorem 1. If , any solution of system (2) satisfies Therefore, for any , the rectangle is a global attractor of (2) in .

Proof. From the first equation of system (2),and then from the comparison principle, . This implies that, for any , there exists such that for all and , where . Then it follows that satisfies:Let be a solution of the following ordinary differential equation:Then . It follows from the comparison principle that and . This completes the proof.

2.2. Persistence

Theorem 2. If and , then system (2) is persistent.

Proof. From the proof of Theorem 1, it can be concluded that there exists such that for any , , where . From the first equation of system (2), for all and ,Because , then from the comparison principle, for any , there exists such that, for any ,where .
Then satisfiesLet be a solution of the following ordinary differential equation:Then . It follows from the comparison principle that there exists such that, for any ,where .
From (11) and (14), it can be easily obtained that and . Therefore, together with Theorem 1, we can conclude that system (2) is persistent.

2.3. Local and Global Stability of Constant Steady States

In this section, we will discuss the local and global stability of constant steady states. In the remaining part of this paper, we always assume that , which guarantees that system (2) has a unique positive constant steady state .

Recall that under the homogeneous Neumann boundary condition has eigenvalues and . Let be the eigenspace corresponding to with multiplicity . Let , , be the normalized eigenfunctions corresponding to . Let , where is an orthonormal basis of . Since , then and .

Theorem 3. For system (2),(a)the trivial constant steady state is unstable;(b)if and , the positive constant steady state is locally asymptotically stable.

Proof. The linearization of (2) at a constant solution can be expressed by , where , ,For each , is invariant under the operator , and is an eigenvalue of the operator on if and only if it is an eigenvalue of the matrixTherefore, the stability is reduced to consider the characteristic equation:where , and .(a)If , then . Obviously, if , for , one of the eigenvalues is . Therefore, is unstable.(b)If , then . It is easy to check that and for all , if and . Therefore, is locally asymptotically stable.
Now, we give the result of global stability of for system (2), which implies that plants and wrack will be spatially homogeneously distributed as time tends to infinity.

Theorem 4. Assume that and ; then is globally asymptotically stable for (2); that is, for any initial values , ,

Proof. If , from Theorem 3, we can conclude that is locally asymptotically stable. Since and , then . From Theorems 1 and 2, for any , and , , , and satisfyInequalities (20) show that and are a pair of coupled upper and lower solutions of system (2) as in the definition in [25, 26] because the nonlinearities in (2) are mixed quasimonotone. Therefore, it is clear that there exists such that, for any ,We define two iteration sequences and as follows: for ,where and . Then for ,and there exist and such thatTherefore, , , , , andSimplifying (25) yieldsSubtracting the first equation of (26) from the second equation, we haveIf we assume thatSubstituting (28) into (26),Hence, the equationhas two positive roots and . Equation (30) can be written as follows:Since and , (31) cannot have two positive roots. Hence, , and, consequently, . Then, from the results in [25, 26], the solution of system (2) satisfiesTherefore, from the above analysis, it can be concluded that the constant steady state is globally asymptotically stable for (2) if and . This completes the proof.

3. Nonconstant Steady States

In this section, we discuss the existence and nonexistence of nonconstant steady states of (2). The corresponding steady-state problem of (2) is the elliptic system:

3.1. A Priori Estimates

Now, we first cite two important lemmas. For notational convenience, we denote .

Lemma 5 (maximum principle [27]). Suppose that .(a)Assume that satisfiesIf , then .(b)Assume that satisfiesIf , then .

Lemma 6 (Harnack inequality [28]). Let be a positive solution to , where , satisfying the homogeneous Neumann boundary condition. Then there exists a positive constant , such that

Theorem 7 (upper bounds). For any positive solution of (33),

Proof. LetApplying Lemma 5, we can getThenTherefore,

Theorem 8 (lower bounds). For any positive constant , there exists a positive constant such that every positive solution of (33) satisfies and if .

Proof. Applying Lemma 5 again, we can getThenNote that and from (40). Equation (43) implies that for some positive constant .
Let . Then if . The Harnack inequality shows that there exists a positive constant such thatCombining (45) with , we can get for some positive constant . It follows from (44) that . This completes the proof.

3.2. Nonexistence of Nonconstant Steady States

In this subsection, we show the nonexistence of nonconstant positive solutions of (33) by the effect of large diffusivity. For convenience, we denote

Theorem 9. Let be a fixed positive constant. Then there exists a positive constant such that (33) has no nonconstant positive solutions if and .

Proof. Assume that is a positive solution of (33). Let thenMultiplying the first equation in (33) by and integrating over , Similarly, we multiply the second equation in (33) by to obtainAdding the above two equations, we can getFrom the PoincarƩ inequality, we can obtainSince , there exists a sufficiently small such that . Let ; then we can conclude that and if and . This completes the proof.

3.3. Existence of Nonconstant Steady States

In this subsection, we discuss the existence of nonconstant positive steady-state solutions of (33), which guarantees the existence of stationary patterns [29ā€“32], when the diffusion coefficients , vary and the other parameters are kept fixed. From Theorem 3, we can easily conclude that is locally asymptotically stable if . In this case, there might be no nonconstant positive steady-state solutions. Therefore, we will restrict our discussion to .

For simplicity, we denoteThen, and (33) can be written as follows:Therefore, is the solution of (33) if and only if it satisfieswhere is the inverse of with the homogeneous Neumann boundary condition. By simple computation, we can get

It is known that is an eigenvalue of on if and only if is an eigenvalue of the following matrix:Obviously,LetThen . Ifthen has two real roots defined byDenote , , and let be the multiplicity of .

To calculate the index of at , we introduce the following lemma.

Lemma 10 (see [33]). Suppose for all , then whereIn particular, if for all , then .
From Lemma 10, in order to calculate the index of at , we need to determine the range of for which .

Theorem 11. If , for some , and is odd, then there exists a positive constant such that (33) has at least one nonconstant positive solution for all .

Proof. Since , it follows that if is large enough, then (60) holds and . Furthermore, If , there exists such thatFrom Theorem 9, we know that there exists such that (33) with and has no nonconstant positive solution. Let be large enough such that , then there exists such thatIn the following, we will prove that, for any , (33) has at least one nonconstant positive solution. By way of contradiction, assume that the assertion is not true for some . By using the homotopy invariance of the topological degree, we can derive a contradiction.
Fixing , for , we defineand consider the following problem:Therefore, is a nonconstant positive solution of (33) if and only if it is a solution of (68) with . It is obvious that is the unique positive constant solution of (68). For any , is a nonconstant positive solution of (68) if and only if it is a solution of the following problem:Our earlier analysis has shown that (69) has no nonconstant positive solution for , and we have assumed that there is no such solution for at . It is obvious thatwhere . From (65) and (66), we have Because is odd, Lemma 10 yieldsFrom Theorems 7 and 8, there exist positive constants and such that the positive solutions of (69) satisfy on for all .
SetThen for all and . By the homotopy invariance of the Leray-Schauder degree [34],Since both equations and have the unique positive solution in , thenThis contradicts (74), and the proof is complete.

4. Bifurcation

4.1. Hopf Bifurcation

In this section, we study the existence of periodic solutions of system (2) by analyzing Hopf bifurcation from the positive constant steady state , and we will show the existence of spatially homogeneous and nonhomogeneous periodic orbits. For convenience, in this section, we assume that all eigenvalues are simple and that the space is one-dimensional spatial domain , where is a positive constant. In fact, this discussion also holds for the generic class of domains in higher dimensions. Let . In the following, we use as the main bifurcation parameter. To look for a possible Hopf bifurcation value, we rewrite the following necessary and sufficient condition from [35, 36].(A1)There exists such that And for the unique pair of complex eigenvalues near the imaginary axis , .

For convenience, (2) can be translated into the following system by the translation of and , while still letting and denote and respectively. This yieldsNow, we consider the linearization near of the above system:The characteristic equation of is where and

We will identify the Hopf bifurcation values which satisfy condition (A1).

(a) WhenthenObviously, if (82) and (83) hold, then and for any . And for any , if . This corresponds to the Hopf bifurcation of spatially homogeneous periodic solution. Furthermore, satisfying (82) is the unique case in which a Hopf bifurcation of spatially homogeneous periodic solution emerges for any . From the above discussion, we can easily get the following theorem.

Theorem 12. Assume that . If (82) and (83) hold, then system (2) undergoes a spatially homogeneous Hopf bifurcation.

In the following, we look for spatially nonhomogeneous Hopf bifurcation for . From Theorem 3, if and , then is locally asymptotically stable. In this case, there cannot be any Hopf bifurcation from . Hence, any potential spatially nonhomogeneous Hopf bifurcation point can make hold. In the following, we assume that .

(b) When , let . Then

Since , there exist and , such that , where is a spatially homogeneous Hopf bifurcation point.

Obviously,Ifthen , and there exist and , such that . Note that , and therefore attains its maximum at , and . DefineThen, for and , we define to be the roots of . These points satisfy

Clearly, and for . It is obvious that . Now it is only necessary to consider whether for all and . From (81), if , then for all . Summarizing the above analysis, we can obtain the main results in this section.

Theorem 13. Assume that , are defined as in (87). If (86) and hold, then for any , there exist points satisfyingsuch that system (2) undergoes a Hopf bifurcation at or , and the bifurcation periodic orbit near can be parameterized in the form of for and a small positive constant. Then where is the corresponding time frequency, is the corresponding spatial eigenfunction, and is the corresponding eigenvector:where is defined as in (78) and is the identity map. Moreover,(1)the bifurcation periodic orbits from are spatially homogeneous, which coincide with the periodic orbits of the corresponding ODE system;(2)the bifurcation periodic orbits from are spatially nonhomogeneous.

The stability and bifurcation direction of the spatially homogeneous periodic orbits bifurcating from can be determined by calculating the normal form according to [35]. Because the calculation process is lengthy and complex, it will not be discussed here, and the interested reader may refer to [35] for details. Because the positive constant steady state is unstable, has at least one pair of eigenvalues with positive real part. Therefore, the spatially nonhomogeneous periodic orbits bifurcating from are all unstable.

4.2. Turing Bifurcation

In this section, we will study Turing bifurcation of the positive constant steady state of system (2). Mathematically speaking, Turing instability occurs when a stable equilibrium point is driven unstable by the local dynamics and diffusion of species. In the following, we consider small space- and time-dependent perturbations for of system (2):where , are small enough and is the wave number. Substituting (92) into (2), we can obtain its characteristic matrix:Correspondingly, and are the roots of the following characteristic equation:where

The roots of (94) yield the dispersion relation:

It is obvious that the stability conditions for of the corresponding ODE system are given by and . It is clear that . Therefore, the stability of in system (2) changes with the sign of . It is easy to establish that for , where

If has positive values, then the range of instability for can be obtained, which is called Turing space. In this case, the positive constant steady state of system (2) is unstable. Consequently, Turing patterns emerge.

5. Numerical Simulations

In this section, we present some numerical simulations to verify the feasibility of the theoretical results obtained in previous sections. In these numerical simulations, we consider one-dimensional spatial domain . The reaction-diffusion system (2) contains six parameters: , , , , , and . If we choose them satisfying the conditions of Theorem 3 (b), then is locally asymptotically stable, as shown in Figure 1. Therefore, the model develops to a homogeneous vegetation state despite strong heterogeneity introduced by the initial conditions (Figure 1). Moreover, if we choose them satisfying the conditions of Theorem 12, then system (2) undergoes a spatially homogeneous Hopf bifurcation and generates spatially homogeneous periodic solutions, as shown in Figure 2. In this case, the plants and wrack are distributed uniformly over the whole space while changing periodically over time (Figure 2). However, if the values of parameters , , and are varied and the other parameters are fixed, from Theorem 13, system (2) undergoes a spatially nonhomogeneous Hopf bifurcation and generates spatially nonhomogeneous periodic solutions, as shown in Figure 3. Under the circumstances, the plants and wrack are heterogeneously distributed in space with high density on the edge and cyclical variations over time (Figure 3). In both cases, they change periodically with time. However, the former is spatially homogeneous regardless of location, while the latter is spatially nonhomogeneous and its distribution is strongly dependent on its spatial location. Furthermore, if the six parameters , , , , , and are chosen suitably, then system (2) can generate nonconstant positive steady-state solutions, as shown in Figure 4. In this case, the plants and wrack will eventually become distributed periodically in space but time-independent (Figure 4). According to the above Turing bifurcation analysis, if , , , , , and as in [8], we can obtain Turing patterns [37]. The numerical results about Turing patterns for model (2) can be found in [8].

In Figure 5, we use and as parameters to plot the bifurcation diagram of model (2). From Theorem 3, it is known that is always unstable. In region I, is globally asymptotically stable. All positive orbits converge to this constant steady state no matter where the initial point is. That is to say, the model will develop to a homogeneous state to any perturbations as time tends to infinity. In region II, is locally asymptotically stable. In this case, both the plants and wrack will be spatially homogeneously distributed to small spatiotemporal perturbations. Region III is Turing space described above, in which many Turing patterns exist. In this region, the plants and wrack have various spatial distributions. Regions IV and V are unstable regions. In region V, spatially nonhomogeneous Hopf bifurcation may occur. Under certain conditions, nonconstant positive steady-state solutions of system (2) may exist in regions III, IV, and V.

6. Conclusions and Discussion

In this paper, we have investigated a diffusive plant-wrack model with homogeneous Neumann boundary conditions theoretically and numerically. Mathematical analysis is used to study dissipation, persistence properties, and local and global stability of constant steady states. To prove global stability, a new method based on the upper and lower solution method is applied. Furthermore, the conditions of existence and nonexistence of nonconstant positive steady-state solutions of the reaction-diffusion system have been derived. Moreover, it has been found that system (2) has no nonconstant positive steady-state solutions if the diffusion coefficients , are large enough. According to Hopf bifurcation analysis, the existence of periodic solutions of system (2), including spatially homogeneous and nonhomogeneous periodic solutions, has been investigated. Moreover, the conditions for Turing instability have also been obtained by Turing bifurcation analysis.

The results presented here indicate that system (2) can exhibit rich dynamics, not only regular spatial pattern formation as described in [8]. From Figures 2 and 3, we can see clearly that the decay rate of wrack and the diffusion constants , have an important effect on the spatial distribution of plants and wrack. Different values of , , and can lead to different spatial distributions of vegetation. Fixing , , and as in Figures 2 and 3, if and , meet appropriate conditions, instability and spatial oscillations may occur, which causes a nonhomogeneous spatial vegetation distribution. Many field studies have shown that oscillations of population density occur in reality [38, 39], and Figures 2 and 3 indicate that the decay rate of wrack and diffusion can indeed result in oscillations with different amplitudes of biomass over time. From [8], it is known that a scale-dependent inhibitory effect of Carex stricta on its own growth, which is strongest at some distance from the tussock center, is key to explaining the formation of regular tussocks in space. In wetland ecosystems, the decay rate of wrack reflects variations in mortality. Moreover, differences in diffusion can lead to changes of biomass in space, which further enhance the scale-dependent inhibitory effect. Therefore, variations of , , and can create different spatial vegetation structures. Obviously, the amplitude of the plant biomass is larger than that of wrack in Figures 2 and 3 (note the difference in the scale of the -axis). Combining differences in amplitude with scale-dependent inhibition, it can be concluded that scale-dependent inhibition has more influence on plants than wrack.

By comparing Figures 1 and 4, it is apparent that different diffusion coefficients and can also lead to different distributions in space. In both cases, although the model ultimately develops to a homogeneous steady state, the spatial distributions of vegetation are not the same. The distribution is spatially homogeneous in Figure 1, but nonhomogeneous in Figure 4. That is to say, scale-dependent inhibition cannot affect the spatial distribution of vegetation in Figure 1 but has an important influence in Figure 4. From Figure 4, the same conclusion can be drawn as from Figures 2 and 3 that scale-dependent inhibition has a greater impact on plants than wrack.

Using the model proposed by van de Koppel and Crain [8], we have presented the first systematic theoretical analysis of the wetland ecosystem, which further deepen and expand the study of dynamics in the wetland ecosystem. Many studies [8, 40ā€“43] have suggested that tussock formation and the development of spatial self-organization structure within vegetation may have important consequences for community structure, especially for species interactions and biodiversity. Our results indeed coincide with these studies. All these results may provide a better understanding of the dynamics of wetland ecosystems.

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 Zhejiang Provincial Natural Science Foundation (nos. LY14C030006 and LY13A010010).