Spatiotemporal patterns induced by cross-diffusion on vegetation model

1.   Introduction

Turing pattern is a kind of pattern formation caused by local instability of the system. In 1952, Turing ^[1] pointed out in his paper that diffusion will make steady-state unstable and eventually form a Turing pattern. Some people ^[2,3] have discovered many similarities between the dynamical instability of the uniform and the equilibrium phase transition when the system is far away from the thermodynamic equilibrium. In recent decades, the phenomenon of desertification has become more and more serious. Based on this, the exploration of the stability of the ecosystem has become a hot spot^[4,5,6].

Patterns of vegetation, well known for one of the most attractive and interesting characters of aerial images, are ubiquitous in semi-arid regions, in which water is regarded as the limiting resource for plant growth. In fact, the existence of vegetation is affected by many factors, such as overgrazing, deforestation and so on^[7]. Some models of vegetation consider the local effect between vegetation and water. It is clear that water is the source of life which is limited in arid and semi-arid areas and has a great impact on resulting in the formation of vegetation patterns. Hence, the study of relationship between the water and vegetation can play a certain guiding role in the stability of vegetation structure.

For the formation of vegetation patterns in semi-arid areas, many various models have been established. R. Lefever ^[8] established a univariate model (dimensionless model) that only included vegetation biomass:

where $u$ represents the number of vegetation or vegetation density, $\mu$ represents the mortality of vegetation, $A$ describes the sensitivity of mutual inhibition and promotion between vegetation, and $L$ represents the ratio of the spatial distance between vegetation mutual promotion and inhibition, model (1.1) reveals the relationship between the short-distance mutual promotion of individual vegetation and vegetation communities and the long-distance competition for resource mechanisms.

In 1999, Klausmeier ^[9] established a mathematical model containing two variables (vegetation and water) for the formation of regular vegetation patterns in semi-arid areas:

among them, $N$, $W$ respectively represent the number or density of vegetation and the density of water. The results showed that the regular vegetation pattern in the semi-arid area is caused by the traveling wave instability of the reaction-diffusion convection equation. The establishment of this model has provided much convenience for the study of vegetation patterns in semi-arid areas.

Later, some authors ^[10] introduced the soil water diffusion feedback existing between vegetation roots and soil water into the klausmeier model:

where the term $\frac{ \gamma w}{(1+\sigma w)}{n}$ describes plant growth at a rate that grows linearly with $w$ for dry soil. The spread of plants is modeled by the diffusion term $\nabla^{2}$. Model (1.3) contains a source tern $p$ representing precipitation and a loss term $-(1-\rho n)w$ representing evaporation. We ignore the evaporation reduced by vegetation evaporation ($\rho > 0$). Local uptake of water by plants (mostly transpiration) modeled by the tern $-w^2n$. $\beta$ describes positive feedback effects of water and biomass. $n$ represents plant population density and $w$ represents groundwater concentration. This nonlinear equation reflects the influence of mutual promotion and competition. The spatial term simulates the spread of vegetation and the transfer of water, and the cross-diffusion term explains the absorption of water by the roots of the vegetation.

In this paper, we ignore the evaporation reduced by vegetation and mainly study the properties of simplified Hardenberg's model with reaction diffusion when $\rho = 0$ as follows:

with Neumann boundary condition:

where $\Omega = \left[ {0, l} \right] \times \left[ {0, l} \right]$, $l$ is a positive bounded constant which gives out the size of the system in the directions of $x$ and $y$, $\varrho$ is the outward unit normal vector of the boundary $\partial \Omega$.

The rest of this paper is organized as follows: In Section 2, we analyze the stability of semi-arid ecosystems model and give the positive steady state existence conditions in order to explore the model's linear stability and Turing bifurcation conditions. Then in Section 3, we derive the amplitude equations of this model and consider the selection of Turing pattern close to the onset $\beta = \beta_{T}$ given later. Next in Section 4, based on the amplitude equations derived in Section 3, we analyze the conditions in different situations of Turing pattern steady-state. And in Section 5, the relaxation time describing the standard time of system to settle down is given to better the numerical simulations. At last, conclusions are shown from aspects of mathematics and biology, respectively, in Section 6.

2.   Stability analysis of semi-arid ecosystems model

In this section, we analyze the stability of the corresponding ODE model at first,

where

with Neumann boundary condition:

where $\Omega = \left[ {0, l} \right] \times \left[ {0, l} \right]$, $l$ is a positive bounded constant which gives out the size of the system in the directions of $x$ and $y$, $\varrho$ is the outward unit normal vector of the boundary $\partial \Omega$.

2.1.   Positive steady state existence conditions

We first consider about the existence of the solutions. For model (2.1), when it doesn't change with time, the second of equations

can be rewritten to the following:

Theorem 2.1. (1) If $(-\gamma+\sigma v) < 0$ holds, then the solution $w_{*}$ of Eq (2.4) must be positive.

(2) If $(-\gamma+\sigma v) < 0$, $v-\sigma < 0$ and $\sigma p-1 < 0$, then Eq (2.4) has the unique positive real solution. Furthermore, if $p(-\gamma+\sigma v)+v < 0$, then the system (2.3) has the unique positive equilibrium $(n_{*}, w_{*})$.

Proof. Denote

since the parameters are positive, so

thus it is obvious that $F(+\infty) < 0$. So there exists $w_{*}$ making $F\left(w_{*}\right) = 0$.

Denote

so

when $w > 0$, since $F(0) = p > 0$ and $F(+\infty) < 0$, so there is only one positive solution $w_{*}$ making $F\left(w_{*}\right) = 0$.

Since $F(p) = p^{3}(-\gamma+\sigma v)+p^{2}(-\sigma+v)+\sigma p^{2}$ and $p(-\gamma+\sigma v)+v < 0$, then $F(p) < 0$, and since $F^{\prime}(w) < 0$ is proved in (2) when $w > 0$, so $p > w_{*}$, notice that $p-w_{*}-w_{*}^{2} n_{*} = 0$ is equivalent to $n_{*} = \frac{p-w_{*}}{w_{*}^{2}}$, that is to say $n_{*} > 0$.

2.2.   Linear stability analysis and Turing bifurcation

Until now, we have already gotten the positive steady state $\left(n_{*}, w_{*}\right)$ of (2.3). And now we discuss about its stability, then analyze the conditions for model (2.2) in 2-D space to admit the Turing bifurcation ^{[11,12,13,14]}. The Jacobian matrix corresponding to this equilibrium is

where

Denote

then, we get the linear perturbations equation:

Therefore, the eigenvalues of linearized operator can be derived by discussing roots of following series of equations. Denote

where $k^2 = l^2+m^2,$ and $l = 1, 2, 3, \cdot \cdot \cdot, \ \ m = 1, 2, 3, \cdot \cdot \cdot \;$.

If there is no diffusion, $k$ is equal to $0$, and we get

So $\left(n_{*}, w_{*}\right)$ is stable.

Next we focus on the diffusion terms:

where

When the diffusion term is considered, the case is different. The sign of the $Det_{k}$ is dominated by the cross diffusion parameter $\beta$. So we should discuss the cases following to make the Turing instability occur.

Denote

We assume that

(N0) $\beta > \beta_{00}$,

(N1) $\beta > \beta_{01}$,

(N2) $\beta > \beta_{02}$,

and the cases are as follows:

(1) If $\beta = 0$, then ${Det}_{k} > 0$, the positive equilibrium $\left(n_{*}, w_{*}\right)$ is stable, the Turing instability will not take place;

(2) If $\beta > 0$, and then we will choose the conditions of $\beta$, furthermore:

(i) If (N0) holds, then the symmetry axis of equation $Det_{k} = 0$ is positive, that is to say, we can have $k$ to get $(Det_k)_{\min }$.

(ii) If (N1) holds, then $(Det_k)_{\min } < 0$.

(iii) If (N2) holds, then $k_{\min}^{2} > \frac{1}{2}$.

Condition (N1) guarantees that $(Det_k)_{\min } < 0$, where $k_{\min }^{2} = k_{T}^{2}.$ Moreover, condition (N2) guarantees that the minimal point $k_{min}^{2} > 0.5.$

We know that $\beta = \beta_{02}(d)$ increases monotonically in $d$ and intersects with $\beta = \beta_{01}$ at the point $d = d_0$, where $d_{0} = 2 \sqrt{\delta\left[n_{*}\left(1+2 w_{*} n_{*}\right)+\frac{w_{*}^{2}}{\left(1+\sigma w_{*}\right)^{2}} \gamma n_{*}\right]}$.

We take

then we get following conclusion:

Lemma 2.1. Suppose (N0) holds, then (N1) and (N2) hold if and only if $\beta > \beta_{B}$.

Denote

then $Det_{k} = 0$ when $\beta = \beta_{*}(k, d)$.

Lemma 2.2. Assume that (N0) holds, function $\beta = \beta_{*}(k, d)$ has following properties:

(a) It will reach the minimum $\beta = \beta_{01}$ at $d = d_ m(k)$, and $d = d_ m(k)$ decreases monotonically as $k$ increases, where

(b) It increases monotonically in $d$ and $k$ when $d > d_ m(k)$.

In Figure 1, a graph is present of functions $\beta = \beta_{01}$, $\beta = \beta_{02}(d)$ and $\beta = \beta_{*}(k, d)$, $d > 0$, $k = \sqrt{2}, \sqrt{5}, \sqrt{8}, ...$, and in Figure 2, we present a graph of Turing bifurcation line with different $k$-modes, which illustrates the results of Lemmas 2.1 and 2.2.

Theorem 2.2. Assume that (N0) holds:

(1) For any given $k_{i}$, where $k_{i}^2 = l^2+m^2, l = 1, 2, 3, ...$, $m = 1, 2, 3, ...,$ and $k_{i}$ is ranked orderly according to the increasing of $l^2+m^2$, like $k_{0} = \sqrt{2}$, $k_{1} = \sqrt{5}$, $k_{2} = \sqrt{8}$, $k_{3} = \sqrt{10}$, $k_{4} = \sqrt{13}$, $k_{5} = \sqrt{17}$, $...$, we have follows:

(i) When $\beta = \beta_{*} \triangleq \beta_{*}\left(k_{1}, d\right)$, $d \in(d_{k_{1}, k_{1}+1}, d_{k_{1}-1, k_{1}})$, if $k = k_1$, characteristic equation $D_{k} = 0$ has a simple real eigenvalue $\lambda = \lambda(k_{1}, \beta)$ with $\lambda(k_{1}, \beta_{*}) = 0$, $\left.\frac{d D_{k_{1}}(\lambda, \beta)}{d \beta}\right|_{\lambda = 0, \beta = \beta_{*}} < 0$, and all other roots of $D_{k}\left(\lambda, \beta_* \right)$ are possession of negative real parts.

(ii) At $\beta = \beta_{*} \triangleq \beta_{*}(k_{1}, d)$, system (2.2) undergoes $k_{1}$-mode Turing bifurcation at $(n_{*}, w_{*})$.

(2) The relationship between $\beta$ and $\beta_{*}$ indicates the occurrence of Turing instability:

(i) $\beta < \beta_{*}(d)$, $d > 0$ is the critical curve to produce Turing instability.

(ii) When $\beta = \beta_{*}(d)$, $d > 0$, there is no Turing instability in system (2.2) with the asymptotically stable equilibrium $(n_{*}, w_{*})$.

(iii) When $\beta > \beta_{*}(d)$, $d > 0$, the diffusion we add makes the equilibrium $(n_{*}, w_{*})$ unstable.

(3) On the critical curve $\beta = \beta_{*}(d)$, $d > 0$, $k_{i}$-mode Turing bifurcation occurs when $d \in(d_{k_{i}, k_{i}+1}, d_{k_{i}-1, k_{i}})$, and $(k_{i}, k_{i+1})$-mode Turing-Turing bifurcation occurs when $d = d_{k_{i}, k_{i}+1}$, $i \in N^*$, which means that one Turing bifurcation curves with wave numbers $k_{i}$ and the other one with $k_{i+1}$ are intersected with each other.

Proof. Firstly, we know directly that

When $d \in(d_{k_{i}, k_{i+1}}, d_{k_{i-1}, k_{i}}),$ $\lambda = 0$ is a root of $D_{k_{i}}(\lambda, \beta_{*}).$ And $Det_{k_{i}} = 0$, $Det_{k} > 0$, where $k \neq k_{i}$, $k = \sqrt{2}, \sqrt{5}, \sqrt{8}, ...\;$. From a direct calculation, we obtain that

thus $\lambda = 0$ is a simple root.

Next, we show that the transversal condition is valid. Let $\lambda = \lambda(k, \beta)$ be the root of $D_{k}(\lambda, \beta)$ satisfying $\lambda(k_{i}, \beta_{*}) = 0,$ then

With the formulas above, we can obtain this:

Remark 2.1. Now we have three conclusions as follows:

(1) $\beta < \beta_{01}$. The constant critical value $\beta_{01}$ does not change with diffusion.

(2) $\beta < \beta_{B}(d)$, $d > 0$. The critical value $\beta_{B}(d)$ is determined by diffusion rather with mode number $k$.

(3) $\beta < \beta_{*}(d)$, $d > 0$. Both diffusion and mode number $k$ have influence on the critical value $\beta_{*}(d)$. And equilibrium $(n_{*}, w_{*})$ is asymptotically stable.

We can see that both terms (1) and (2) in Remark 2.1 are sufficient conditions, which Turing instability does not take place under, and term (3) in Remark 2.1 is not only sufficient but also necessary.

Remark 2.2. The critical curve of Turing instability $\beta = \beta_{*}(d)$, $d > 0$ is called the first Turing bifurcation curve, one with the corresponding characteristic equation having no root with positive real part. It is a piecewise smooth curve, with piecewise point called Turing-Turing bifurcation point $T_{k_{i}.k_{i+1}}$, see in Figure 2. Notice that wave number $k$ and diffusion coefficient $d$ affect the expression of the first Turing curve, so stable spatial pattern can be found with wave number k, where $k = \sqrt{2}, \sqrt{5}, \sqrt{8}...\;$.

Remark 2.3. Focusing on Lemma 2.1 and Figure 2, It is clear that if diffusion ratio $\beta$ is relatively constant, wave number $k$ of spatial pattern is affected by the diffusion coefficient $d$. Smaller the diffusion coefficient $d$ is, larger the wave number $k$ is.

3.   Derivation of amplitude equation

We haven't determined the selection of Turing pattern based on the above discussion. In this section, we will take the selection of Turing pattern into consideration with free condition in 2-D space. $\beta$ is close to the onset $\beta_{T}$, the eigenvalues are around zero corresponding to the critical modes which vary slowly, meanwhile the off-critical modes relax quickly, so only perturbations with $k$ around $k_{T}$ should be considered. Since amplitude equations dominate the dynamics of the system, we analyze the stability of different patterns near the onset using the amplitude equations. There are two methods to drive the coefficients of amplitude equations: One is symmetrical analysis and the other is standard multiple-scale analysis ^[15,16]. The critical value is

The Taylor series expansion is

and then

is considered, and then substituted into system (2.2). The terms $O(\rho^4)$ are omitted for the sake of convenience, then we get a system which approaches system (2.2) as follows:

and by substituting the specific expressions into Eq (3.2), we get

The solutions of model (3.2) can be expanded to

where $c.c$ stands for the complex conjugate.

Denote

and

where $N_{1}, N_{2}$ are given in Appendix. Then, model (3.2) can be changed to as the follows:

where

For this system (3.5), the parameter $\beta$ behavior near the onset $\beta_{T}$ is analyzed. In this way, $\beta$ is expanded in the following term:

where $\varepsilon$ is considered as a small parameter. The variable $u$ and the nonlinear term $N$ are expanded according to $\varepsilon$ so that we obtain the follows:

where $N_{3}$, $N_{4}$ are given in Appendix.

And $L$ can be expanded as follows:

where

Separating the dynamical behavior of the system (3.5) according to different time scale is regarded as the core of the standard multiple-scale analysis scale. We just need to separate the time scale (i.e. $T_{0} = t$, $T_{1} = \varepsilon t$, $T_{2} = \varepsilon^2 t$). Each $T_{i}$ corresponds to the dynamical behaviors of the variables, whose scales are $\varepsilon^{-1}$, an independent variable. So the derivative with respect to time converts to the following term:

For model (3.5), we consider the following result:

Substituting Eqs (3.6)–(3.9) into Eq (3.5) and expanding it according to different orders of $\varepsilon$, we can obtain three equations as follows:

where $F_{1}^{n}$, $F_{1}^{w}$ are given in Appendix.

where $F_{2}^{n}$, $F_{2}^{w}$ are given in Appendix.

For the first order of $\varepsilon$, as $L_{T}$ is the linear operator of the system near the onset, $(n_{1}, w_{1})^T$ is linear combination of the eigenvectors that corresponds to the eigenvalue zero. Solving the solution of the first order of $\varepsilon$, we can obtain

where

$\left|k_{j}\right| = \left|k_{T}\right|$, and $\left|w_{j}\right|$ is the amplitude of the mode $exp(ik_{j}r)$ if the system is under the first-order perturbation. The perturbational term of the higher order determines its form. For the second order of $\varepsilon$, we obtain Eq (3.13) above.

First we directly solve the equation with the solutions of Eq (3.15) which we have already calculated and by the comparing with the corresponding coefficients. We will get the following results:

Secondly, according to the Fredholm condition, the vector function of the right-hand side of Eq (3.13) must be orthogonal with the zero eigenvectors of operator $L_{T}^+$, so that the existence of the nontrivial solution of this equation can be ensured. $L_{T}^+$ is the adjoint operator of $L_{T}$, it is as follow:

In this system, the zero eigenvectors of operator $L_{T}^+$ are

The orthogonality condition is

that is

where $R_{1}$, $R_{2}$ are given in Appendix.

For the third order of the $\varepsilon$ in Eq (3.14), as with its second step of the second order, we have

that is

where $B_1$–$B_{11}$ are given in Appendix.

Expanding the above equation, we have

where $C_1$–$C_6$ are given in Appendix.

Denote $A_i = A_i^v$, and then, we expand amplitude $A_{1}$ in the following form:

Substituting it into Eq (3.18), the amplitude equation corresponding to $A_{1}$ can be obtained as follows:

where $C_1$–$C_6$ are given in Appendix.

4.   Analysis of Turing pattern steady-state

Eq (3.20) can be decomposed to mode $\rho_{i} = |A_{i}|$ with a corresponding phase angle $\varphi_{i}$. Then, four differential equations of the real variables are obtained as follows^[15,16,17]:

where $\varphi = \varphi_{1}+\varphi_{2}+\varphi_{3}$.

The dynamical system (4.1) has four kinds of solutions:

(i) The stationary state $(0, 0, 0)$ is stable;

(ii) Stripe pattern, given by

and exist only when $\mu$ is of the same sign as $g_{1}$;

(iii) Hexagonal pattern, given by

Denote

the existence condition of $\rho_{1}, \; \rho_{2}, \; \rho_{3}$ is $\mu > \mu_{1}$;

(iiii) Mixed structure pattern, given by

where $g_{2} > g_{1}$.

Without discussing solution (i) which is trivial all time, we firstly take the stability of stripe pattern into consideration. Setting $\rho_{1} = \rho_{0}+\delta \rho_{1}$, $\rho_{2} = \delta \rho_{2}$, $\rho_{3} = \delta \rho_{3}$, substituting them into Eq (4.1), and then linearizing, we get

The characteristic equation is

so the eigenvalues are

denote

then we get stable stripe pattern when $\mu > \mu_{2}$ and $\mu > \mu_{3}$.

Secondly, we focus on the stability of hexagonal pattern. Substituting $\rho_{i} = \rho_{0}+\sigma_{i}$, where $i = 1, 2, 3$, into Eq (4.1), then linearizing, we get

where

The character equation is

so the eigenvalues are

denote

only when

all the eigenvalues for $\rho_{+}$ are negative, which means that we have the stable hexagon pattern; while the eigenvalues $s_{1}$, $s_{2}$ for $\rho_{-}$ are positive which means that we have the unstable hexagon pattern. And the solutions (iiii) are always unstable, so that we do not make a discussion about it.

In next section, we will simulate the following solutions:

(1)

(2)

(3)

where if the $\varphi = 0$ in system (4.1), the $H_{0}$ mode hexagon is stable when $\mu\in(\mu_{1}, \mu_{4})$ and unstable when $\mu > \mu_{4}$ while the $H_{\pi}$ mode hexagon is unstable when $\mu > \mu_{2}$, and if the $\varphi = \pi$ in system (4.1), the $H_{0}$ mode hexagon is unstable when $\mu\in(\mu_{2}, \mu_{4})$ and stable when $\mu > \mu_{4}$ or $\mu\in(\mu_{1}, \mu_{2})$ while the $H_{\pi}$ mode hexagon is stable when $\mu > \mu_{2}$.

5.   Numerical simulations

In this section, extensive numerical simulations of the spatially extended model (1.4) in two-dimensional space are performed, with the qualitative results shown here. The space and time of the problem need to be discretized in order to solve differential equations with the help of computer. A discrete domain with $M*N$ lattice sites whose space between the lattice points is $\Delta h$ is built up to solve the continuous problem corresponding to the reaction-diffusion system in two-dimensional space. The time evolution is in a discrete process, with steps of $\Delta t$, and can be solved by the Euler method. Then, we calculate the Laplacian describing diffusion by using finite differences in the discrete system. The Turing pattern obtained from numerical simulations of Eq (2.2) satisfies initial and boundary conditions. We choose the time step $\Delta {t} = 0.0005$ and a system size of $50*50$ with the space step $\Delta {x} = 0.50$, $\Delta {y} = 0.50$. We keep $\gamma = 1.6$, $\delta = 5.2$, $p = 0.2958$, $\sigma = 1.6$, $v = 0.2$. By calculating, we analyze the Turing pattern of the system that we have already discussed about in following steps:

(1) First of all, we calculate the control parameters $\beta$ and the important control parameters $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, $\mu_{4}$ by the parameters given above and give a picture to show that result directly;

(2) Secondly, we calculate the relaxation time $\tau$ whose expression is given later, which means the time of the system to nearly settle down, and then we give the evolution pictures of $\mu$ in different intervals, such as $[\mu_{1}, \mu_{2}]$, $[\mu_{2}, \mu_{3}]$ and $[\mu_{3}, \mu_{4}]$ as well.

5.1.   Parameters calculate

The values of control parameter $\beta$ is corresponding to parameters $\mu_{1} = -0.00467$, $\mu_{2} = 0$, $\mu_{3} = 0.07179$, $\mu_{4} = 0.30022$, and $\beta_{T} = 7.3738$ is the critical value of Turing bifurcation. Then we present the figures below to show the different solutions' stability regarded to $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, $\mu_{4}$.

In Figure 3, we take the amplitude as a function of the control parameter $\mu$ which is used to describe the bifurcation of the stripe pattern $S_{p}$ and the spot pattern $H_{0}$ or $H_{\pi}$ respectively.

5.2.   Numerical analysis of evolutions

Firstly, by the analyze in Section 3, we obtain the relaxation time $\tau$ as follow:

then we get 4.15492 after substituting the parameters into the formula of $\tau$. It is equivalent to 8310 iterations when we choose the time step as $\Delta {t} = 0.0005.$

Secondly, we give the evolution pictures with $\mu$ in different intervals, such as $[\mu_{1}, \mu_{2}]$, $[\mu_{2}, \mu_{3}]$ and $[\mu_{3}, \mu_{4}]$. And different types of patterns are observed in the process of numerical simulations. The distributions of water species $w$ and vegetation species $n$ are always in the same type. Here, we only presert the distribution of the population density pattern of vegetation species $n$.

5.2.1.   Uniform state and hexagon pattern

We take $\mu = -1.08725 * 10^{-6}\in(\mu_{1}, \mu_{2})$. The pattern of uniform stationary solution occupies the whole domain eventually at the end. Figure 4 shows the evolution process of the spatial pattern of species $n$ after 0 iterations, 500 iterations, 8310 iterations, 20000 iterations, 300000 iterations, 400000 iterations and 500000 iterations.

5.2.2.   Hexagon pattern

We take $\mu = 0.0717635\in(\mu_{2}, \mu_{3})$. At this time, the whole domain is occupied with $H_{0}$ hexagon pattern eventually. Figure 5 shows the evolution process of the spatial pattern of species n after the same numbers before. And the relation time is also 8310 iterations. From the Figure 5, we can find when the pattern settle down, the stripes decrease gradually into nonexistent while the spots are still in our sight.

5.2.3.   Hexagon and strip pattern

We take $\mu = 0.261217\in(\mu_{3}, \mu_{4})$. At this time, $H_{0}$ hexagon pattern and strip are coexistence in the whole domain in the end. Figure 6 shows the evolution process of the spatial pattern of species $n$ after the same iterations as above and with the same relaxation time 8310 iterations too, and we can find the coexistence of the $H_{0}$ hexagon pattern and the strips clearly as the theoretical analysis prejudged.

5.2.4.   Strip pattern

We take $\mu = 0.300274 > \mu_{4}$. In this way, $H_{0}$ hexagon pattern is dismissing while the strips are substitute for the $H_{0}$ hexagon in the whole domain eventually. In Figure 7 forms squares consist of horizontal and vertical strips finally, and they are both shown with the same iteration numbers of evolution and the same relaxation time 8310 iterations as above. Comparing with the analysis, there should be only stripe pattern under this situation. Evidently the both series of the pictures are consistent with the theoretical analysis very well too.

6.   Conclusions

For a semi-arid ecosystem model, a great amount of works discuss the stability and bifurcation based on Lyapunov method, but some of the studiers don't prefer to consider the diffusion effects, especially the cross diffusion. In this paper, we take it into consideration that a semi-arid ecosystems system is with cross diffusion effect in 2-D spatial domain. And some conditions of Turing instability to the model (1.4) is given with the help of Turing's bifurcation theory. We study both the amplitude equations and the stability of different patterns in detail. We also introduce the conception "relaxation time", with which we could find the standard to show the time when the system nearly settles down.

Furthermore, we also focus on the biology significance of the model (1.4). In the natural world, two species are in relationships of both coexistence and competition. The specie "water" may be recognized as a restrained survival source by the other specie "vegetation", meanwhile, the diffusion of the vegetation make a influence of the water since the water is absorbed by the vegetation root. This phenomenon are caused by cross-diffusion. For the sake of this, the effect of cross-diffusion on the pattern structure is considered and the results indicate that cross-diffusion is a critical term to the Turing spatial pattern formulation. Without cross diffusion, the instability will not take place in the system refer to Section 2 discussion. Besides, the numerical results are well consistent with the theory. When we take the stability of different patterns into consideration, such as spots, stripes, we find it is changing with the $\mu$ refer to Section 4 discussion. Biologically, it means that if we randomly plant some trees, it will occurs in four situations: (1) If $\mu < \mu_{1}$, the number of trees will probably decrease. So the region will change into a desert; (2) If $\mu_{2} < \mu < \mu_{3}$, the pattern of trees is more like some spots called small bushes; (3) If $\mu_{3} < \mu < \mu_{4}$, the pattern of the trees will change into the state of coexistence of both stripes and spots. So the gaps between two groups of the vegetation is either lines or labyrinths; (4) If $\mu > \mu_{4}$, the pattern of tree groups will be like some stripes. Based on the descriptions of the four situations, we would optimise the $\mu$, dominated by parameters physically, larger than $\mu_{2}$, so that the desert will probably not increase, especially when the pattern is strip, the $\mu$ should be lager than $\mu_{3}$ in case it will degenerate into nothing. By these useful decisions, not only a constructive way to ensure the stationary pattern formation of the environment is presented, but also an interesting usage of Turing pattern is provided.

Conflict of interest

The authors declare that they have no conflicts of interest.

Appendix

where