Coupled Effects of Turing and Neimark-Sacker Bifurcations on Vegetation Pattern Self-Organization in a Discrete Vegetation-Sand Model

Abstract

Wind-induced vegetation patterns were proposed a long time ago but only recently a dynamic vegetation-sand relationship has been established. In this research, we transformed the continuous vegetation-sand model into a discrete model. Fixed points and stability analyses were then studied. Bifurcation analyses are done around the fixed point, including Neimark-Sacker and Turing bifurcation. Then we simulated the parameter space for both bifurcations. Based on the bifurcation conditions, simulations are carried out around the bifurcation point. Simulation results showed that Neimark-Sacker bifurcation and Turing bifurcation can induce the self-organization of complex vegetation patterns, among which labyrinth and striped patterns are the key results that can be presented by the continuous model. Under the coupled effects of the two bifurcations, simulation results show that vegetation patterns can also be self-organized, but vegetation type changed. The type of the patterns can be Turing type, Neimark-Sacker type, or some other special type. The difference may depend on the relative intensity of each bifurcation. The calculation of entropy may help understand the variance of pattern types.

1. Introduction

Vegetation patterns have been widely observed in arid and semi-arid areas [1,2,3,4,5,6,7]. There have been many studies on the exploration of pattern self-organization mechanisms [8,9,10,11,12,13]. This research has contributed to the understanding of vegetation competition when resources are limited and the process by which robust land is converted into desert because of climate change or anthropogenic disturbances [14,15]. The relationship between self-organization and the production of entropy can find its origin in the pioneering works in [16,17]. Since then many studies have explored the relationship in a variety of fields [18,19]. The calculation of entropy can be used for pattern classification [20], therefore, the investigation of the mechanisms of self-organization will enhance the understanding of entropy. Besides, due to the complexity of the ecological scale, investigation of vegetation patterns has never stopped.

Vegetation bands lying perpendicular to the prevailing wind direction have been observed, showing dying trees in the windward edge and seedling regrowth on the leeward edge of each band [21,22,23]. While [24,25] suggested that in Jordan wind might be the initiating and driving factor of banded patterns and the accumulation of wind-blown material around isolated plants might act as a nucleus for the development of vegetation arcs. Dynamical models for these wind-induced vegetation patterns have been rare. Recently Zhang et al. [26] proposed a vegetation-sand model using partial differential equations to interpret the mechanism of vegetation self-organization in windy sandy environments [27,28,29,30,31]. In the model of Zhang et al. [26], wind is considered as prevailing and non-prevailing (defined as winds from all other directions except the prevailing wind direction) and each has different effects on vegetation dispersal and sand transportation. The effects of prevailing wind on sand and vegetation are modeled by advection terms, while the effects of non-prevailing wind are modeled as diffusion terms. Simulations showed that without the prevailing wind labyrinth vegetation patterns can be obtained, and with the prevailing wind banded vegetation patterns perpendicular to the prevailing wind direction can be obtained.

In Nature, wind events and sand movements are discrete. Wind collection and recording are also discrete. Therefore, we think the discretization of the continuous model is necessary. Discrete models have provided advantages in revealing complex nonlinear dynamics. A classic example is the discrete logistic model that exhibits period-doubling cascade and a route to chaos [32]. In contrast, the continuous logistic model provides a simple “S” form curve, and never demonstrates the above dynamic complexity. Through a variety of bifurcations, discrete models can generate periodic orbits, invariant circles, periodic windows, chaotic behavior and so on. More importantly, as births and deaths of bionts in predator-prey systems are discrete events, continuous models only make sense for very large populations. Other researchers propose that discrete models reveal the discontinuous properties (such as a patchy environment or a fragmented habitat) of predator-prey systems [33]. Discrete models may exhibit new dynamical behaviors. For example, Han et al. found that Turing instability and Turing patterns can occur in a simple discrete competitive Lotka-Volterra system rather than a continuous one [34].

In the discrete model, focusing on the time scale, Neimark-Sacker bifurcation can occur without spatial terms. Neimark-Sacker bifurcation is one of the most important bifurcations in discrete dynamics. Through Neimark-Sacker bifurcation, a system stable states bifurcate from a fixed point to an invariant circle. It should be noted that with spatial terms such as advection and diffusion, Turing bifurcation can also occur. Naturally, a question arises as to what the dynamics will be when Turing and Neimark-Sacker bifurcation occur at the same time, and how they influence the self-organization of vegetation patterns. To the authors’ knowledge, there have been no simulations of the coupling effects of Turing and Neimark-Sacker bifurcations on vegetation patterns.

In this research, we will transform the vegetation-sand model [26] to a spatially and temporally discrete model. Then, focusing on the stable fixed point, bifurcation analysis including Neimark-Sacker bifurcation and Turing bifurcation will be carried out. Bifurcation diagram and phase portrait will be shown. Based on bifurcation analysis, parameter values will be selected around each bifurcation. Given these parameters, simulations will be carried out to show how patterns transform from Turing bifurcation to Turing-Neimark-Sacker bifurcation (Turing bifurcation and Neimark-Sacker bifurcation occur at the same time). Different transformation processes will be simulated. Finally, we will discuss the possible reason that cause the different transformations.

2. A Discrete Vegetation-Sand Model

The continuous vegetation-sand model was proposed by [26] and is shown as:

$$\begin{gathered} \frac{\partial S}{\partial t}=k_0+mV\left(1-\frac{V}{V_0}\right)-nS-a_1\frac{\partial S}{\partial x}+D_1\left(\frac{{\partial }^{2}S}{\partial x^2}+\frac{{\partial }^{2}S}{\partial y^2}\right), \\ \frac{\partial V}{\partial t}=hV\left(1-\frac{V}{V_m}\right)-pS\frac{V}{C+V}-a_2\frac{\partial V}{\partial x}+D_2\left(\frac{{\partial }^{2}V}{\partial x^2}+\frac{{\partial }^{2}V}{\partial y^2}\right) \end{gathered}$$

(1)

The model is considered in the dimensional form rather than the non-dimensional form. This is because there are few parameters, and each parameter has its own ecological significance. Here S (cm) represents the height of sand deposition. Parameter k₀ (cm·d⁻¹) is the accumulation rate of air sand without vegetation. m (cm·d⁻¹) is the coefficient representing the trapping effect of vegetation. Parameters k₀ and m are affected by sand flow intensity. The trapping effects of vegetation decrease with the increase of vegetation cover until a certain value (near 100%) of vegetation cover (saturation constant of vegetation trapping effects), and V₀ (%) is two times of this certain value. n (d⁻¹) is the coefficient representing the decrease of accumulation rate with the increase of sand height. a₁ (m·d⁻¹) is the coefficient representing the translation of sand dunes by wind. D₁ (m²·d⁻¹) is the diffusion coefficient of sand without strong unidirectional wind. V (%) represents vegetation coverage, and h (d⁻¹) represents the intrinsic growth rate of vegetation. V_m (%) is the potential maximum of vegetation coverage. p (cm⁻¹·d⁻¹) is the coefficient of destruction effect by sand. C (%) is a constant representing how sand tolerance increases with the increase of vegetation cover. a₂ (m·d⁻¹) is the advection coefficient representing the dispersal of vegetation by wind. D₂ (m²·d⁻¹) is the diffusion coefficient representing the dispersal of vegetation without strong unidirectional wind. x and y are space axes, and t is time. Wind direction is along the positive x direction. When there is no prevailing wind, the effects of wind are modeled as diffusion terms. When there is a prevailing wind, its effects can be modeled by adding an advection term to the diffusion term as shown in Equation (1).

In this research, the above continuous model will be transformed to a discrete model (Equation (2)). We consider the model on a N × N lattice, and the two variables can be expressed as S_(i,j,t) and V_(i,j,t) (i,j ∈ {1, 2, 3, …, N} and t ∈ Z⁺), that represent the height of sand deposition and the vegetation coverage in lattice (i,j) at time t, respectively. According to the former research works of [33,35,36,37], there are two stages, reaction stage and diffusion stage, when we discretize the continuous model (Equation (1)). The spatial dispersal stage, that are advection and diffusion, is considered firstly as:

$$\{\begin{array}{c}S{\prime }_{(i,j,t)}=S_{(i,j,t)}+\frac{\tau }{d}{a}_1{\nabla }_d{S}_{(i,j,t)}+\frac{\tau }{d^2}{D}_1{\nabla }_d^2{S}_{(i,j,t)}\\ V{\prime }_{(i,j,t)}=V_{(i,j,t)}+\frac{\tau }{d}{a}_2{\nabla }_d{V}_{(i,j,t)}+\frac{\tau }{d^2}{D}_2{\nabla }_d^2{V}_{(i,j,t)}^{}\end{array}.$$

(2)

where S_(i,j,t) and V_(i,j,t) are the transitional variables of S and V after one step of dispersal stage. τ (days) and d (meters) are the time step and space step respectively. ${\nabla }_d$ denotes the discrete form of advections. ${\nabla }_d^2$ denotes the discrete form of the Laplacian operator:

$$\begin{gathered} {\nabla }_d{S}_{(i,j,t)}=S_{(i+1,j,t)}-S_{(i,j,t)} \\ {\nabla }_d{V}_{(i,j,t)}=V_{(i+1,j,t)}-V_{(i,j,t)} \\ {\nabla }_d^2{S}_{(i,j,t)}=S_{(i+1,j,t)}+S_{(i-1,j,t)}+S_{(i,j+1,t)}+S_{(i,j-1,t)}-4S_{(i,j,t)} \\ {\nabla }_d^2{V}_{(i,j,t)}=V_{(i+1,j,t)}+V_{(i-1,j,t)}+V_{(i,j+1,t)}+V_{(i,j-1,t)}-4V_{(i,j,t)}. \end{gathered}$$

(3)

Boundary conditions are set as periodic conditions. Then we consider the reaction stage:

$$\{\begin{array}{c}{S}_{(i,j,t+1)}=f_1(S{\prime }_{(i,j,t)},V{\prime }_{(i,j,t)})=S{\prime }_{(i,j,t)}+\tau f(S_{(i,j,t)},V_{(i,j,t)})\\ V_{(i,j,t+1)}=g_1(S{\prime }_{(i,j,t)},V{\prime }_{(i,j,t)})=V{\prime }_{(i,j,t)}+\tau g(S_{(i,j,t)},V_{(i,j,t)})\end{array},$$

(4)

in which:

$$\{\begin{array}{c}f(S,V)=k_0+mV\left(1-\frac{V}{V_0}\right)-nS\\ g(S,V)=hV\left(1-\frac{V}{V_m}\right)-pS\frac{V}{C+V}\end{array},$$

(5)

Equations (2) and (4) including both diffusion and reaction stages are defined as our discrete model.

As the previous work in [26] has been done on the condition of equilibrium in the continuous model, we focus on one of the most important situations, which has three equilibria under the same parameter condition. The three equilibria are: one boundary equilibrium (S_b,V_b), one unstable saddle point (S₂,V₂), and one locally asymptotical stable node (S₁,V₁). These equilibria are also fixed points in the discrete model, and can be expressed as:

$$\left(S_b,V_b\right)=\left(\frac{k_0}{n},0\right),$$

(6)

$$\left(S_1,V_1\right)=\left(\frac{k_0}{n}+\frac{m{b}_4}{n}\left(1-\frac{b_4}{V_0}\right),b_4\right),$$

(7)

$$\left(S_2,V_2\right)=\left(\frac{k_0}{n}+\frac{m{b}_5}{n}\left(1-\frac{b_5}{V_0}\right),b_5\right),$$

(8)

in which:

$$\begin{gathered} b_1=\left(\frac{h}{p{V}_m}-\frac{m}{n{V}_0}\right),b_2=\left(\frac{m}{n}-\frac{h}{p}+\frac{hC}{p{V}_m}\right),b_3=\left(\frac{k_0}{n}-\frac{hC}{p}\right), \\ {b}_4=\frac{-b_2+\sqrt{b_2^2-4b_1{b}_3}}{2b_1},b_5=\frac{-b_2-\sqrt{b_2^2-4b_1{b}_3}}{2b_1}, \end{gathered}$$

(9)

In [26], stability analysis shows that (S₁,V₁) is the only locally asymptotically stable interior equilibrium. Therefore focusing on (S₁,V₁), we have the Jacobian matrix:

$$J(S,V)=\left(\begin{array}{cc}1-n\tau & \frac{\tau m(-2V+V_0)}{V_0}\\ -\frac{\tau pV}{C+V}& 1+\tau (h-\frac{2hV}{V_m}-\frac{pSC}{{(C+V)}^2})\end{array}\right),$$

(10)

and eigenvalues:

$${\lambda }_{1,2}=\frac{t{r}_0\pm \sqrt{t{r}_0{}^2-4{\Delta }_0}}{2}$$

(11)

In which (S₁,V₁) is the interior fixed point, and:

$$\begin{array}{l}t{r}_0=2-\tau (n-h+\frac{2h{V}_1}{V_m}+\frac{p{S}_{1}C}{{(C+V_1)}^2})\\ {\Delta }_0=1+\tau (h-n-\frac{2h{V}_1}{V_m}-\frac{p{S}_{1}C}{{(C+V_1)}^2})+{\tau }^2(-nh+\frac{2nh{V}_1}{V_m}+\frac{np{S}_{1}C}{{(C+V_1)}^2}+\frac{mp{V}_1(-2V_1+V_0)}{(C+V_1)V_0})\end{array}.$$

(12)

Therefore (S₁,V₁) is locally asymptotically stable only when|λ₁| < 1, and |λ₂| < 1, which yields:

$$\{\begin{array}{c}{\Delta }_0<1\\ -(1+{\Delta }_0)<t{r}_0<1+{\Delta }_0\end{array},$$

(13)

We set most parameter values unchanged in this research: h = 0.2 d⁻¹; V_m = 100%; p = 0.045 cm⁻¹·d⁻¹; m = 0.2 cm·d⁻¹; n = 0.07 d⁻¹; V₀ = 200%; C = 10%; D₂ = 0.01 m²·d⁻¹. And other parameters such as k₀, D₁, a₂, a₁, τ and h₁ will vary.

3. Bifurcation Analysis

Turing bifurcation is most commonly used in the study of vegetation self-organized patterns, and there are many studies focusing on Turing type patterns [38,39,40]. Spatial symmetry breaking can be induced by Turing bifurcation, leading to the formation of patterns that are stationary in time and oscillatory in space. While there is another type of symmetry breaking, temporal symmetry breaking. Temporal symmetry breaking can be induced by Neimark-Sacker bifurcation, giving rise to states that are homogeneous in space and oscillatory in time. When temporal and spatial symmetry breakings take place simultaneously, the discrete model will generate the patterns that are oscillatory in both space and time. In this section, the parameter conditions of Neimark-Sacker bifurcation, Turing bifurcation and both bifurcations will be obtained. Bifurcation analysis will be focused on the fixed point (S₁,V₁). The condition of Neimark-Sacker bifurcation will be calculated and bifurcation diagram will be shown. Then Turing bifurcation will be re-investigated as the model has been transformed to a discrete model. And based on the two bifurcation conditions, we choose two key parameters τ and a₁ and calculate the parameter space that satisfy both bifurcations. In this section, we will only show the key results of Neimark-Sacker and Turing bifurcation analysis. The complete calculation process of each bifurcation condition can be seen in the Appendix A.

3.1. Neimark-Sacker Bifurcation Analysis

According to [41], assume that U₀(μ) is the asymptotically stable fixed points of map:

$$U_0\to F\left(U_0,\mu \right),U_0\in R^2,\mu \in R^1$$

(14)

and its eigenvalues are conjugate $\lambda \left(\mu \right)$, $\overline{\lambda }\left(\mu \right)$. If:

$$\left|\lambda \left({\mu }_0\right)\right|=1,{\lambda }^j\left({\mu }_0\right)\ne 1,j=1,2,3,4$$

(15)

$$\frac{d}{d\mu }\left(\left|\lambda \left({\mu }_0\right)\right|\right)\ne 0,$$

(16)

$$a=-\mathrm{Re}\left(\frac{\left(1-2\overline{\lambda }\right){\overline{\lambda }}^2}{1-\lambda }{\xi }_{11}{\xi }_{20}\right)-\frac{1}{2}{\left|{\xi }_{11}\right|}^2-{\left|{\xi }_{02}\right|}^2+\mathrm{Re}\left(\overline{\lambda }{\xi }_{21}\right)\ne 0,$$

(17)

in which:

$${\xi }_{20}=\frac{1}{8}\left(\left(F_{\tilde{x}\tilde{x}}-F_{\tilde{y}\tilde{y}}+2G_{\tilde{x}\tilde{y}}\right)+i\left(G_{\tilde{x}\tilde{x}}-G_{\tilde{y}\tilde{y}}-2F_{\tilde{x}\tilde{y}}\right)\right),$$

(18)

$${\xi }_{11}=\frac{1}{4}\left(\left(F_{\tilde{x}\tilde{x}}+F_{\tilde{y}\tilde{y}}\right)+i\left(G_{\tilde{x}\tilde{x}}+G_{\tilde{y}\tilde{y}}\right)\right),$$

(19)

$${\xi }_{02}=\frac{1}{8}\left(\left(F_{\tilde{x}\tilde{x}}-F_{\tilde{y}\tilde{y}}-2G_{\tilde{x}\tilde{y}}\right)+i\left(G_{\tilde{x}\tilde{x}}-G_{\tilde{y}\tilde{y}}+2F_{\tilde{x}\tilde{y}}\right)\right),$$

(20)

$${\xi }_{21}=\frac{1}{16}\left(\left(F_{\tilde{x}\tilde{x}\tilde{x}}+F_{\tilde{x}\tilde{y}\tilde{y}}+G_{\tilde{x}\tilde{x}\tilde{y}}+G_{\tilde{y}\tilde{y}\tilde{y}}\right)+i\left(G_{\tilde{x}\tilde{x}\tilde{x}}+G_{\tilde{x}\tilde{y}\tilde{y}}-F_{\tilde{x}\tilde{x}\tilde{y}}-F_{\tilde{y}\tilde{y}\tilde{y}}\right)\right).$$

(21)

then the Map (14) undergoes Neimark-Sacker bifurcation when μ = μ₀,

Our discrete model can be transformed into a map:

$$\left(\begin{array}{c}S\\ V\end{array}\right)\to \left(\begin{array}{c}S+\tau (k_0+mV(1-\frac{V}{V_0})-nS)\\ V+\tau (hV(1-\frac{V}{V_m})-pS\frac{V}{C+V})\end{array}\right)$$

(22)

Considering Map (22), (S₁,V₁) is also the asymptotically stable fixed point. Then the map undergoes Neimark-Sacker bifurcation at the fixed point (S₁,V₁), when τ = τ₀, in which:

$${\tau }_0=\frac{n-h+\frac{2h{V}_1}{V_m}+\frac{p{S}_{1}C}{{(C+V_1)}^2}}{-nh+\frac{2nh{V}_1}{V_m}+\frac{np{S}_{1}C}{{(C+V_1)}^2}+\frac{mp{V}_1(-2V_1+V_0)}{(C+V_1)V_0}}$$

(23)

and τ₀ satisfies the following conditions:

$${\tau }_0(n-h+\frac{2h{V}_1}{V_m}+\frac{p{S}_{1}C}{{(C+V_1)}^2})\ne 2,3.$$

(24)

$$\begin{array}{l}a=-\frac{AN-MB}{32a_{12}^2{\beta }^2({(1-\alpha )}^2+{\beta }^2)}-\frac{1}{32a_{12}^2{\beta }^2}({(G_{1\tilde{w}\tilde{w}}+G_{1\tilde{z}\tilde{z}})}^2+{(G_{2\tilde{w}\tilde{w}}+G_{2\tilde{z}\tilde{z}})}^2)\\ -\frac{1}{64a_{12}^2{\beta }^2}({(G_{1\tilde{w}\tilde{w}}-G_{1\tilde{z}\tilde{z}}-2G_{2\tilde{w}\tilde{z}})}^2+{(G_{2\tilde{w}\tilde{w}}-G_{2\tilde{z}\tilde{z}}+2G_{1\tilde{w}\tilde{z}})}^2)\\ +\frac{1}{16a_{12}\beta }(\alpha (G_{2\tilde{w}\tilde{w}\tilde{z}}+G_{2\tilde{z}\tilde{z}\tilde{z}})+\beta (G_{2\tilde{w}\tilde{w}\tilde{w}}+G_{2\tilde{w}\tilde{z}\tilde{z}}))\ne 0\end{array}$$

(25)

in which:

$$\begin{gathered} M=\beta \left(3-4\alpha \right)\left({\alpha }^2-{\beta }^2\right)-2\alpha \beta \left(\left(1-\alpha \right)\left(1-2\alpha \right)-2{\beta }^2\right), \\ N=\left(1-3\alpha +2{\alpha }^2-2{\beta }^2\right)\left({\alpha }^2-{\beta }^2\right)+\left(6\alpha -8{\alpha }^2\right){\beta }^2, \\ A=(G_{1\tilde{w}\tilde{w}}+G_{1\tilde{z}\tilde{z}})(G_{1\tilde{w}\tilde{w}}-G_{1\tilde{z}\tilde{z}}+2G_{2\tilde{w}\tilde{z}})-(G_{2\tilde{w}\tilde{w}}+G_{2\tilde{z}\tilde{z}})(G_{2\tilde{w}\tilde{w}}-G_{2\tilde{z}\tilde{z}}-2G_{1\tilde{w}\tilde{z}}), \\ B=(G_{2\tilde{w}\tilde{w}}+G_{2\tilde{z}\tilde{z}})(G_{1\tilde{w}\tilde{w}}-G_{1\tilde{z}\tilde{z}}+2G_{2\tilde{w}\tilde{z}})+(G_{1\tilde{w}\tilde{w}}+G_{1\tilde{z}\tilde{z}})(G_{2\tilde{w}\tilde{w}}-G_{2\tilde{z}\tilde{z}}-2G_{1\tilde{w}\tilde{z}}). \end{gathered}$$

(26)

Figure 1 shows the variations of V versus parameter τ when the parameter values satisfy the Neimark-Sacker bifurcation condition. Note that this bifurcation diagram has no visible periodic windows. When τ < 16.6997 days, the fixed point is asymptotically stable. When τ = 16.6997 days, the system starts to bifurcate around a fixed point. As the value of τ increases, the stable states of the system in the phase plane (V,S) experience several stages, such as invariant circle (τ = 16.7013 days as shown in Figure 1b), period-7 (τ = 16.755 days as shown in Figure 1c), and then invariant circle again (τ = 16.8667 days as shown in Figure 1d). Then the stable states of the system may go through several multi-period and invariant circle stages.

3.2. Turing Bifurcation Analysis

According to [42], assume that U₀(μ) is the fixed points of one reaction-diffusion model:

$$U_{t+1}=v_{1}f(U_t,\mu )+v_2{\nabla }_d^2{U}_t+v_3{\nabla }_d{U}_t,$$

(27)

In which, ${\nabla }_d$ and ${\nabla }_d^2$ denotes the discrete form of advection and Laplacian operator respectively. f is smooth. v₁, v₂, v₃ are coefficients. λ(μ) is the eigenvalues of the model at U = U₀(μ). When:

• In the absence of diffusion and advection, U₀(μ) is asymptotically stable;
• With diffusion and advection, max(λ(μ)) > 0;

then the model undergoes Turing bifurcation.

In the discrete models, eigenvalues of ${\nabla }_d$ and ${\nabla }_d^2$ are difficult to obtain. But according to the method in [43], the eigenvalues of ${\nabla }_d$ can be obtained as:

$${\lambda }_{kl}^1=2\sin {\varphi }_k\exp \left(({\varphi }_k-\frac{\pi }{2})i\right).$$

(28)

where ${\varphi }_k=\left(k-1\right)\pi /N$, $k\in \left\{1,2,3,\cdots ,N\right\}$ and i represents the complex unit, that is $i=\sqrt{-1}$. While the eigenvalues of ${\nabla }_d^2$ can be obtained as:

$${\lambda }_{kl}=4\left({\sin }^2{\varphi }_k+{\sin }^2{\varphi }_l\right),$$

(29)

in which ${\varphi }_k=\left(k-1\right)\pi /N$, ${\varphi }_l=\left(l-1\right)\pi /N$, and $k,l\in \left\{1,2,3,\cdots ,N\right\}$.

In this research, we have known that (S₁,V₁) is the asymptotically stable fixed point of our discrete vegetation-sand model (Equations (2) and (4)). The eigenvalues of the model are:

$${\lambda }_{\pm }\left(k,l,\tau \right)=-\frac{1}{2}tr\left(k,l,\tau \right)\pm \frac{1}{2}\sqrt{tr{\left(k,l,\tau \right)}^2-4\Delta \left(k,l,\tau \right)},$$

(30)

in which:

$$tr\left(k,l,\tau \right)=t{r}_0(\tau )+\frac{\tau }{d}\left(a_{11}(\tau )a_1+a_{22}(\tau )a_2\right)a_1{\lambda }_{kl}^1+\frac{\tau }{d^2}\left(a_{11}(\tau )D_1+a_{22}(\tau )D_2\right){\lambda }_{kl},$$

(31)

$$\Delta \left(k,l,\tau \right)={\Delta }_0\left(\tau \right)\left(1+\frac{\tau }{d}{a}_1{\lambda }_{kl}^1-D_1\frac{\tau }{d^2}{\lambda }_{kl}\right)\left(1+\frac{\tau }{d}{a}_2{\lambda }_{kl}^1-\frac{\tau }{d^2}{D}_2\delta {\lambda }_{kl}\right).$$

(32)

Here tr₀(τ), ∆₀(τ), a₁₁(τ), and a₂₂(τ) are denoted for reminding that they are dependent on τ. Based on the two eigenvalues, define:

$${\lambda }_m\left(k,l,\tau \right)=\max \left(\left|{\lambda }_{+}\left(k,l\right)\right|,\left|{\lambda }_{-}\left(k,l\right)\right|\right),$$

(33)

$${\lambda }_{mm}\left(\tau \right)=\underset{k=1,l=1}{\overset{N}{\max }}{\lambda }_m\left(k,l,\tau \right).$$

(34)

where λ_mm(τ) represents the maximal value of absolute modulus of both eigenvalues in (33). Therefore, the threshold condition for the occurrence of Turing bifurcation requires λ_mm(τ) = 1. When λ_mm(τ) > 1, Turing instability occurs; when λ_mm(τ) < 1, the discrete system stabilizes at the homogeneous states.

Turing bifurcation can be shown through the variations of eigenvalues λ(k,l) as shown in Figure 2. In Figure 2a, we can see that the effects of the perturbation numbers k and l are symmetric. Thus we let k = l, and we can obtain the variations of eigenvalues versus l as shown in Figure 2b. When there is no perturbation, the system remains at the fixed point. When the advection coefficient a₁ = 0.005 m·d⁻¹, the eigenvalues of the system remains at less than 1 with the increase of perturbation number l, but when the diffusion coefficient a₁ increases more than a_1c, the eigenvalues will exceed 1 with the increase of l and Turing bifurcation occurs.

3.3. Parameter Space for Neimark-Sacker and Turing Bifurcation

Based on the bifurcation analysis above, parameters satisfying both Turing and Neimark-Sacker bifurcations will be selected to investigate the coupled effects of both bifurcations. Before that, we select two parameters which are parameter τ and a₁, to show the parameter space for both bifurcations. τ is defined as the time step, it can reflect the frequency of wind events. However, as wind frequency is varying naturally (here we consider the real wind frequency, rather than the measured or recorded data), therefore, the variation of τ is meaningful and the four spatial parameters, D₁, D₂, a₁ and a₂, can reflect the values of wind speed, non-prevailing wind and prevailing wind, respectively. The previous study [26] gave that a₁/a₂ = D₁/D₂, and we mainly investigate prevailing wind. We can see that the variation of parameter a₁ (parameter a₂ will vary accordingly) represents the variation of prevailing wind speed. Note that wind, especially prevailing wind is the most important driving force in the process of vegetation pattern self-organization. Wind speed and frequency are the two key indexes of wind. Therefore, simulations on vegetation pattern self-organization will also be carried out through the variation of parameter τ and a₁.

Figure 3 shows the parameter space that satisfies Turing and Neimark-Sacker bifurcations around bifurcation point τ₀. Blue area and green area mean parameter values satisfy only Neimark-Sacker bifurcation and only Turing bifurcation respectively. In Figure 3a, τ₀ = 16.6997 days. And in Figure 3b, τ₀ = 16.0428 days. Given τ ∈ [0.98τ₀, 1.04τ₀] in both Figure 3a,b, we can see that Turing bifurcation area in Figure 3b is bigger than that in Figure 3a. And this will cause differences in the self-organization of patterns (will be shown in the next section). Turing-Neimark-Sacker bifurcation is defined as when parameter values satisfy both Turing and Neimark-Sacker bifurcations. Note that, we define the dark green area in Figure 3 as Turing-Neimark-Sacker bifurcation area. The condition for this area is:

$$\tau >{\tau }_0,{\lambda }_m\left(1,1,\tau \right)>1,{\lambda }_{mm}\left(\tau \right)>{\lambda }_m\left(1,1,\tau \right),$$

(35)

4. Coupled Effects of Turing and Neimark-Sacker Bifurcations on Vegetation Pattern Self-Organization

On the basis of the bifurcation analysis and the calculation of the parameter space for both Turing and Neimark-Sacker bifurcations, we will investigate the coupled effects of Turing and Neimark-Sacker bifurcation on vegetation pattern self-organization. Numerical simulations will be carried out around the bifurcation point τ₀. With the increase of time scale τ, system will go through only Turing bifurcation to Turing-Neimark-Sacker bifurcation. Note that the color scale in all the following figures are automatically selected by the software MATLAB 7.12.0 (2011a) command “pcolor” according to the minimum and the maximum values of vegetation cover at each time. Precisely, the blue area represents low vegetation coverage, while the red area represents high vegetation coverage. Sometime, the minimum value of vegetation coverage reaches zero, which means vegetation dies out in the vegetation-sand system and only bare sand can be seen. Besides, after the simulation of vegetation pattern self-organization in each figure, the values of vegetation coverage will be mapped to the integers in the domain [1, 256] (which represents the 256 colors in the color spectrum), and then the Shannon entropy of each figure will be calculated.

Firstly, vegetation patterns induced by only Neimark-Sacker bifurcation are shown in Figure 4. Given parameters k₀ = 3.095 cm·d⁻¹, h = 0.2 d⁻¹, V_m = 100%, p = 0.045 cm⁻¹·d⁻¹, m = 0.2 cm·d⁻¹, n = 0.07 d⁻¹, V₀ = 200%, C = 10%, D₁ = 0.01 m²·d⁻¹, D₂ = 0.01 m²·d⁻¹, a₁ = 0.01 m·d⁻¹, a₂ = 0.01 m·d⁻¹, d = 2 m, Neimark-Sacker bifurcation occurs when τ = τ₀ = 16.6997 days. When τ > 16.6997 days, spatially heterogenous vegetation patterns can be induced by the Neimark-Sacker bifurcation. With the increase of τ, wavelength of vegetation patterns becomes smaller, and the whole pattern becomes more complex as shown in Figure 4b. The type of the spatially heterogeneous vegetation patterns is complex and irregular, and difficult to define. In this research, we define it as the Neimark-Sacker type patterns. The values of entropy in Figure 4a,b are calculated as 7.7066 and 7.7708, respectively.

Figure 5 shows the transformation of vegetation patterns from under Turing bifurcation to under Turing-Neimark-Sacker bifurcation. Simulations are carried out around the Neimark-Sacker bifurcation point τ₀ = 16.6997 days. When τ < τ₀, Turing bifurcation occurs and Figure 5a,b show that regular striped vegetation patterns can be formed. Note that the amplitude of vegetation patterns in Figure 5a,b is not very big. When τ > τ₀, Neimark-Sacker bifurcation and Turing bifurcation both occur, and vegetation patterns are changed. As shown in Figure 5c–e, we can see that pattern wavelength becomes smaller, stripes become narrower, curvy and irregular. Compared with the patterns in Figure 4, the vegetation patterns in Figure 5c–e are like stripes coupled with the Neimark-Sacker type patterns, but they are more similar to Neimark-Sacker type patterns. The values of entropy in Figure 5a–e are calculated as 5.8924, 7.8776, 7.5626, 7.5309 and 7.4856, respectively.

Figure 6 shows another series of vegetation patterns with the increase of time scale τ. When τ < τ₀ = 16.0428 days, Turing bifurcation occurs and Figure 6a,b show that regular striped vegetation patterns with large amplitude can be formed. When τ > τ₀, Neimark-Sacker bifurcation occurs. However the striped vegetation patterns seem unchanged as shown in Figure 6c–e. Obviously, these patterns are Turing type. The values of entropy in Figure 6a–e are calculated as 7.3101, 7.5809, 7.0303, 7.5018 and 7.3576, respectively. From Figure 5 and Figure 6, we can see that the coupled effects of Turing-Neimark-Sacker bifurcations can be different.

Figure 7 and Figure 8 will show two special situations of Turing-Neimark-Sacker bifurcations. From previous study, we know that without prevailing wind, labyrinth patterns can be obtained, while with prevailing wind, striped patterns can be formed. These two patterns will be shown in Figure 7 and Figure 8, respectively, and the effects of Turing-Neimark-Sacker bifurcations will be simulated.

Figure 7a,b show that labyrinth patterns can be induced by Turing bifurcation, but Neimark-Sacker bifurcations destroy the self-organizations and induces spots as shown in Figure 7c–e. With the increase of parameter τ, vegetation spots become less and more isolated. The values of entropy in Figure 7a–e are calculated as 7.1693, 6.9560, 3.5169, 2.1530 and 2.0576, respectively.

Figure 8a,b show that regular striped vegetation patterns can be induced by Turing bifurcation and Turing-Neimark-Sacker bifurcations induce the breakdown of the self-organizations as shown in Figure 8c–e. Vegetation stripes can also be formed, but the number and the width of vegetation stripes drastically decreases. With the increase of parameter τ, vegetation stripes become less and shorter. The values of entropy in Figure 8a–e are calculated as 7.0886, 6.3024, 2.3306, 4.1024 and 2.9138, respectively.

The different simulation results in Figure 5, Figure 6, Figure 7 and Figure 8 need to be discussed. From the comparison between Figure 5a and Figure 6a, we can see that the intensity of Turing bifurcation is different. Although we cannot assess the density of Turing bifurcation quantitatively, we can tell the relative intensity from the amplitude of the patterns. Clearly, the amplitude of vegetation patterns in Figure 6a is bigger than that in Figure 5a. Therefore, when Neimark-Sacker bifurcation also occurs, the coupled effects of Turing and Neimark-Sacker bifurcations will make Figure 5a more like Neimark-Sacker type patterns, and make Figure 6a more like Turing type patterns. Simulations in Figure 7 and Figure 8 show quite different results. From Figure 7 and Figure 8, we can see that vegetation patterns induced by Turing-Neimark-Sacker bifurcations are neither Turing type (especially Figure 7) nor Neimark-Sacker Type patterns. It seems that the coupled effects of Turing and Neimark-Sacker bifurcations make the degradation process of vegetation fiercer. Therefore, labyrinth patterns are broken down into spots, vegetation stripes become less. It is necessary to point out that simulations in Figure 7 and Figure 8 may depend on the special coupling of Turing and Neimark-Sacker bifurcations.

In the real world, vegetation bands are not perfectly straight. There have been several explanations on this issue. Some researchers think that this may be due to the stochastic micro-topography. There is another interpretation in [26], that the self-organization process takes a very long time (about 30,000 days). Before the regular bands are formed, the values of some parameters may change. For example, advection coefficients a₁ and a₂, as the speed and the direction of prevailing wind are changing all the time. Therefore, the patterns in the real world are curvy bands. In this research, we propose another possible explanation. We choose the parameter τ to be the varying parameter because the variation of time scale reflects the variation of wind frequency well. As wind frequency is changing all the time, the variation of τ becomes meaningful. Then the variation of τ causes the occurring of Neimark-Sacker bifurcation. Neimark-Sacker bifurcation and Turing bifurcation may occur at the same time. The coupled effects of Neimark-Sacker and Turing bifurcations make vegetation patterns more irregular and similar to real vegetation patterns in nature. Therefore, the irregular patterns of vegetation may be induced by the coupled effects of Neimark-Sacker and Turing bifurcation.

From the calculated values of entropy in all figures, we can see that the values of entropy don’t vary much in Figure 4, Figure 5 and Figure 6, but the values of entropy in Figure 7 and Figure 8 decrease drastically when Turing and Neimark-Sacker bifurcations both occur. The decrease of entropy reflects the decrease of information in the figures. This is consistent with what we can observe directly from the figures, as most vegetation dies out due to the coupled effects of Turing and Neimark-Sacker bifurcations, leaving only little vegetation.

5. Conclusions

The continuous model of wind-induced vegetation patterns is quite new and attractive. In this research, we transform it to a discrete model. The discretization of the continuous vegetation-sand model is due to the consideration that wind events and sand movements are both discrete. Based on the discrete model, this research considered Neimark-Sacker bifurcation, which was not considered in the previous study [26]. Special and complex vegetation patterns can be self-organized (Figure 4) when parameter values satisfy Neimark-Sacker bifurcation. Simulation results also show that labyrinth (as shown in Figure 7a,b) and striped vegetation patterns (as shown in Figure 6 and Figure 8a,b) can be obtained, which are the key results in [26]. Besides, with the increase of time scale τ, we simulated how vegetation patterns transform from Turing bifurcation to Turing-Neimark-Sacker bifurcation (Figure 5, Figure 6, Figure 7 and Figure 8). Especially the transformations from Figure 5a–e, the transformations from Figure 7a–e and the transformations from Figure 8a–e reveal the coupling effects of Turing and Neimark-Sacker bifurcations.

From all the bifurcation analysis, simulations and discussion, we can conclude that:

(a)

The method of discretization provides a new scenario in the study on wind-induced vegetation patterns. It preserves the patterns that can be obtained in [26].

(b)

After discretization, the variation of time scale becomes possible. This variation helps to understand pattern self-organization under different ecological scales.

(c)

With the discrete vegetation-sand model, we investigated the coupled effects of Turing and Neimark-Sacker bifurcations. Under both bifurcation conditions, the type of simulated patterns depends on the intensity of each bifurcation. Sometime one bifurcation effect dominates the self-organization of patterns, while specially they couple with each other and lead to pattern variation, which can also be supported by the variance of entropy in Figure 7 and Figure 8. However, the coupled effects of Turing and Neimark-Sacker bifurcations are so complex that new methods may be needed to assess the intensity of each bifurcation.