Cost-benefit analysis for the installation of cogeneration CSP technology in Cyprus

1.   Introduction

The coexistence of trees and grasses on the tropical grassland has always been the topic of theoretical and empirical research. One of the landscape characteristics of savanna is the coverage of grass layer and tree layer, especially the coverage of community, and one of the mechanisms causing this landscape is random repeated disturbance. In other words, the fire disturbance suppresses the specific stage of tree life cycle, limits the density of upper trees, and generates an open and unsaturated tree canopy, which promotes the growth and development of grasses ^[1,2,3,4]. In Savanna, fire is a frequent disturbance that has a "inhibitor effect" on the growth of trees and grasses ^[5,6]. Moreover, these two species can adapt to the local fire system, and they will be able to survive and reproduce in repeated fires ^[7]. However, too frequent or too infrequent fires may lead to species extinction ^[8,9]. To analyze the fire effect on the coexistence balance of trees and grasses, ecologists have established many types of tree-grass models with fire feedback, and some of which assume that fire was periodicity, unpredictable, and determined by external factors ^[10,11]. Others believe that fire come from the internal feedback of vegetation pattern, and the probability of fire is a predictable disturbance affected by many factors such as vegetation and its rich characteristics ^[4,12,13,14]. For areas covered by grasses and trees (including grassland and woodland), and the fire issue is fundamental. This is because fire plays a key role in controlling the generation mechanism of tree-grass spatial patterns, which is a global problem ^[15,16]. However, most scholars analyze the fire effect in vegetation pattern dynamics of tree-grass by collecting observations and statistical analysis, but there are few in-depth theoretical analysis.

In grassland ecosystem, the formation area of tree-grass pattern is broad, and it needs a long evolution process (decades or even longer). For this long-term and large-scale spatial evolution mechanism of tree-grass, considering the propagation and diffusion characteristics of vegetation is one of the key factors to explore the generation of tree-grass spatial pattern in grassland ecosystem. Therefore, the reaction-diffusion theory is suitable for the dynamic analysis of vegetation pattern ^[17]. Meanwhile, many researches agree that the reaction-diffusion system is one of the basic equations describing the movement of nature, and their application range covers almost many subjects ^[18,19,20], such as predator-prey model ^[21,22,23], infectious disease model ^[24,25], population competition model ^[26,27] and neuron oscillator model ^[28,29], and when control parameter of the model cross the threshold value (bifurcation occur), a qualitative change in the behavior of a system will be induced. In order to understand the underlying ecological mechanisms in the process of tree-grass coexistence, explore the effect of fire disturbance on grassland desertification, learn from the analysis results of bifurcation in these studies, and expand the dynamic method of vegetation pattern research, which has important guiding significance for delaying savanna desertification and promoting the restoration of tree-grass vegetation. There is little theoretical analysis on the dynamic characteristics of tree-grass model with the fire frequency as control parameter, and there are also few numerical simulations to verify and predict the distribution of tree-grass density pattern. In the following analysis, we will use the bifurcation dynamics theory and method to analyze the effect of fire frequency on tree-grass coexistent equilibrium from the perspective of nonlinear dynamic system theory.

The other organizational structures of this paper are as follows: In section 2, a universal tree-grass competition model under fire feedback is selected, and the tree-grass model with time-delay and diffusion term is established based on it. In section 3, we mainly analyze the stability of the equilibrium point of the tree-grass model, the Hopf bifurcation induced by time delay, and the sufficient conditions for the Turing bifurcation caused by diffusion based on the theorem of differential equations and bifurcation theory. In section 4, the amplitude equation of tree-grass model with time delay and diffusion is solved by multi-scale analysis method ^[17]. Taking the delay as the control parameter, and the stability of the Turing pattern was analyzed. In section 5, the sensitivity analysis of fire frequency to each variable of tree-grass model is discussed. Based on the analysis, a series of numerical simulations on the tree-grass pattern are carried out, and compared with the spatial structure of tree-grass density distribution on Hulunbuir grassland, which is obtained through field observation. It is found that the theoretical analysis results are consistent with the actual grassland ecosystem.

2.   The model

We represent the dynamics of the savanna system with a tree-grass model (2.1) as the theoretical model of this paper ^[6].

By introducing the dimensionless transformation as follows:

then the system (2.1) becomes:

where $ B_1 = \frac{r_w}{r_g} $, $ \delta(f) = c\frac{g}{r_{g}} $ ($ c $ is a dimensionless parameter that describes the strength of the vegetation feedback on fire, which in this case is driven by grass biomass $ G $ of a given pyrogenicity). $ W $, $ G $ represent the biomass of trees and grasses respectively. $ r_w $, $ r_g $ represent the change rate of the biomass of trees and grasses respectively. $ k_w $, $ k_g $ represent the maximum biomass(carrying capacity) of trees and grasses that can be supported in a given area respectively. $ M_w $, $ M_g $ are constant, which represent the loss rate of the biomass of trees and grasses in fire respectively. $ \delta(f) $ represents the fire frequency. $ d w/d t $, $ d g/d t $ represent the rate of change of state variables $ w $ and $ g $ in the range $ [0, \; 1] $ respectively. The growth term $ {B_1}w(1 - w) $ in the first equation of model (2.2) indicates that trees are planted where there are no trees, and grasses does not inhibit their growth. The growth term $ g(1 - g - w) $ in the second equation of model (2.2) indicates that grasses is planted in a place where neither trees nor grasses, that is, the increase of grass biomass is not only inhibited by its own accumulation, but also decreases with the increase of tree biomass, where $ \delta(f){M_w}w $, $ \delta(f){M_g}g $ represent the disappearance terms of tree and grass biomass caused by fire frequency respectively.

For the convenience of calculation, let $ B_1 = B $, $ w = u $, $ g = v $, $ M_{w} = a $, $ M_{g} = b $, $ \delta(f) = \delta $, then the system (2.2) becomes:

we have $ a > 0 $, $ b > 0 $, $ \delta > 0 $ and $ B > 0 $ for the biological implications. The establishment of model (2.3) reflects the two characteristics of coexistence balance between upper and lower species in vegetation spatial^[6]. One is asymmetric competition, where trees shade the surface vegetation(grasses), but the surface vegetation(grasses) has relatively little influence on the tree layer, resulting in the gradual replacement of the grass layer by trees without interference (see the direction of change indicated by the black solid line in Figure 1). It shows that in grassland ecosystem, when the interference of fire frequency is not considered, trees are in a dominant position and at the expense of grass layer, and the competitive effect of grasses is not enough to prevent the rapid supplement of tree biomass ^[3,30,31]. The other is the different sensitivity to fire feedback, compared with grass layer, when the fire intensity is higher, the damage degree of tree layer is more serious(see width of brown solid line in Figure 1). This is because the tree and grass species in the core area of fire are burned on a large scale. Due to the slow growth and the long regeneration cycle of trees, it can not recover within limited interval of fires. When considering the interference of fire frequency, grasses as the fuel of fire determines the fire intensity (see the positive feedback of the red dotted line in Figure 1). Since the canopy of trees is higher than the heart of the fire, the water content of trees is higher than that of grasses, and trees are considered to be the factor to inhibit the occurrence of fire in the early stage of fire(see the negative feedback of the blue dotted line in Figure 1). Trees can also switch from inhibiting the fires to promoting them if they occur frequently. Note: The model structure is shown in Figure 1, which is improved according to the reference^[12].

For the model (2.3) under the coupling of the fire, tree species and grass species may diffuse in the direction of regions rich in resources, or may diffuse in areas with few competitors because of the asymmetric competition between species, or may diffuse with the influence of external environmental factors such as water flow. Based on the above complex ecological process, It is reasonable to study the influence of diffusion on the evolution trend of vegetation pattern by using the bifurcation theory of reaction-diffusion equation (use Laplace operator $ \Delta $ to describe the diffusion characteristics of trees and grasses).

In addition, there may be months or even years between two consecutive fires, during which tree species and grass species will regrow and affect each other ^[6]. Therefore, a series of effects of fire frequency cycle (tempering time) on tree-grass density distribution need to be considered. It is worth noting that the delay effect caused by fire frequency is one of the key parameters affecting the dynamic change of vegetation landscape pattern in Savanna ^[6]. The vegetation density distribution before the occurring time of fire directly affects the next state of vegetation growth, that is, a series of consequences caused by the time delay "$ \tau $" of fire frequency can not be ignored. Based on the above considerations, this paper introduces the effect of time delay $ \tau $ of fire frequency to the reduction item $ a\delta u(t) $, $ b\delta v(t) $ respectively, then the system (2.3) in the following reaction-diffusion form:

under Neumann boundary condition

with the initial condition:

where $ u(x, y, t) $, $ v(x, y, t) $ represent the state variables of trees and grasses with respect to time $ t $ and spatial position $ (x, y) $, respectively. Let $ \Omega = (0, L)\times (0, L) $, which represents a bounded domain with the smooth boundary $ \partial \Omega $ in two-dimensional space $ {\bf{R}}^{2} $, where $ L $ represents a positive boundary constant. $ \Delta = \frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2} $ is a Laplacian operator in space $ \Omega $. $ \delta a u(x, y, t-\tau) $ and $ \delta b v(x, y, t-\tau) $ with delay term $ \tau $ in model (2.4) represent the loss of biomass caused by the duration of fire, respectively. $ d_{11}, \; d_{22} $ represent the passive diffusion coefficients of $ u(x, y, t), \; v(x, y, t) $ respectively, and we have $ d_{11} > 0, \; d_{22} > 0 $ for biological sense. $ d_{12}, \; d_{21} $ represent the cross-diffusion coefficients of $ u(x, y, t), \; v(x, y, t) $ respectively. From the definition of cross-diffusion in reference ^[32], then the coefficient $ d_{12} $ indicate the influence of $ v(x, y, t) $ density to $ u(x, y, t) $ density. That is to say, $ u(x, y, t) $ is repelled from $ v(x, y, t) $ for $ d_{12} > 0 $, and $ u(x, y, t) $ is attracted to $ v(x, y, t) $ for $ d_{12} < 0 $. When the role of $ u(x, y, t) $ and $ v(x, y, t) $ switched, then $ d_{21} $ has the same meaning. In addition, we assume $ d_{11}d_{22}-d_{12}d_{21} > 0 $, which indicates that the flux of the respective densities in the spatial domain depends more strongly on their own density than on the other ^[33,34].

For the tree-grass model (2.4), the two main explanations for the coexistence balance of trees and grasses are either competition for resources or different sensitivity to external disturbances, such as sensitivity to fire frequency^[35]. From the perspective of resource competition, trees and grasses are in a competitive relationship. In order to avoid competition with the grasses (trees), the trees(grasses) will choose to migrate from the high density area to the low density area of the grasses (trees), it was corresponded to $ d_{12} > 0\; (d_{21} > 0) $.

From the perspective of sensitivity to fire frequency, grasses is the fuel of fire, while trees is more damaged by fire ^[12]. Under the disturbance of fire frequency, the proportion between trees and grasses, and the distribution of vegetation density will change ^[3,4,35]. From the sensitivity to fire disturbance, if the fire frequency is within the threshold range, because the adult trees have certain fire resistance, and grasses will sacrifice its own biomass as the fuel to support the continuous combustion of fire. Then the tree layer is in a dominant position, can quickly regrow from the infrequent occurrence of fire, and expand the range of shade trees at the expense of the grass layer, thus slowly reducing grass biomass. Therefore, in the process of vegetation restoration under infrequent fire interference, the grass population will be gradually replaced by trees, and the Savanna gradually expands to the forest, that is the tree population will spread from the low-density area to the high-density area of grasses, which was corresponded to $ d_{12} < 0 $.

With the increase of fire frequency, tree and grass biomass will continue to decrease, and more frequent fires will lead to the extinction of tree species. However, grass species can recover rapidly from limited tempering interval of fire due to the fast growth cycle, and tree species are no longer in a dominant competitive position due to the slow growth and metabolism (compared with seedlings, grasses is in a dominant competitive position) ^[12]. Finally, the grass population grow vigorously due to short recover circle and lack of inhibition from dominant competitors (trees), and delay the regeneration of tree population. Therefore, in the process of vegetation restoration under frequent fire interference, the tree population will be gradually replaced by grasses, and the Savanna gradually expands to the grassland, that is the grass population will spread from the low-density area to the high-density area of trees, which was corresponded to $ d_{21} < 0 $. Therefore, regardless of $ d_{12} < 0 $ or $ d_{21} < 0 $, it has practical ecological significance.

3.   Stability analysis and Turing instability conditions

Set $ f(u(t), v(t)) = B u(t)(1-u(t)-a \delta) $, $ g(u(t), v(t)) = v(t)(1-u(t)-v(t)-b \delta) $. Next, solve the equations $ f(u(t), v(t)) = 0 $ and $ g(u(t), v(t)) = 0 $, then the corresponding ODE system (2.3) always has three boundary equilibrium points $ E_0 = (0, 0) $, $ E_1 = (u_{1}, 0) $, $ E_2 = (0, v_{2}) $ and one coexistence equilibrium point $ E_* = (u_{*}, v_{*}) $. In order to satisfy the biological significance, the positive equilibrium condition $ (H_{0}) $ holds, where $ u_1 = 1-a \delta/B $, $ v_2 = 1-b \delta $, $ u_* = u_1 $ and $ v_* = a \delta/B-b \delta $.

3.1.   Stability analysis of equilibrium points

Let $ \tilde{u}(x, y, t) = u_{*}-u(x, y, t) $, $ \tilde{v}(x, y, t) = v_{*}-v(x, y, t) $, which are substituted into Eq (2.4) and removing the sign "$ \sim $", the Taylor expansion at $ E_* = (u_{*}, v_{*}) $ is as follows:

then the linearization of the system (3.1) is:

the matrices are:

where $ a_{11} = B(1-2u_{*}) $, $ a_{21} = -v_{*} $, $ a_{22} = 1-u_{*}-2v_{*} $, $ a_{13} = -a\delta $, $ a_{24} = -b\delta $, and the nonlinear term of the system (3.1) is

The characteristic equation corresponding to Eq (3.2) is:

where $ p_1^k = -a_{11}-a_{22}+d_{11} k^2+d_{22} k^2 $, $ p_2^k = \left(d_{11} k^2-a_{11}\right) \left(d_{22} k^2-a_{22}\right)-d_{12}k^2 \left(d_{21} k^2-a_{21}\right) $, $ p_3^k = -a_{13}-a_{24} $, $ p_4^k = a_{24}\left(a_{11}-d_{11} k^2\right)+a_{13}\left(a_{22}-d_{22} k^2\right) $ and $ p_5 = a_{13} a_{24} $.

Under the condition $ (H_0) $, assume $ a_{11} < 0, \; a_{22} < 0 $ hold, then we obtain the condition

where $ \delta _{0} = \frac{B}{a} $. Then under the condition $ (H_1) $, the characteristic equation (3.3) with $ k = 0, \; \tau = 0 $ becomes:

then we have

The corresponding ODE system (2.3) has an unique asymptotically stable positive equilibrium point $ E_{*} $. Next, using the same method to judge the properties of the other three boundary equilibrium points $ E_{0}, \; E_{1}, \; E_{2}. $ The characteristic eigenvalues of the system (2.3) at $ E_{0} $ are $ {\lambda _1}({E_0}) = B - af $ and $ {\lambda _2}({E_0}) = 1 - bf $, then we have $ {\lambda _1}({E_0}) > 0, \; \; {\lambda _2}({E_0}) > 0 $. The characteristic eigenvalues of the system (2.3) at $ E_{1} $ are $ {\lambda _1}({E_1}) = - {\lambda _1}({E_0}) $ and $ {\lambda _2}({E_1}) = v_{*} $, then we have $ {\lambda _1}({E_1}) < 0, \; \; {\lambda _2}({E_1}) > 0 $. The characteristic eigenvalues of the system (2.3) at $ E_{2} $ are $ {\lambda _1}({E_2}) = {\lambda _1}({E_0}) $ and $ {\lambda _2}({E_2}) = - {\lambda _2}({E_0}) $, then we have $ {\lambda _1}({E_2}) > 0, \; \; {\lambda _2}({E_2}) < 0 $.

In summary, the properties of the four balance points $ E_{0}, \; E_{1}, \; E_{2} $, $ E_{*} $, and the property of eigenvalues of the system in the neighborhood of these equilibrium points are summarized in Table 1. The main results and conclusion are as below.

Lemma 3.1. For the system (2.3), if the condition $ ({H}_1) $ holds, then

$ Case \; i $. The tree-grass coexistence equilibrium $ E_{*} = (u_{*}, \; v_{*}) $ is asymptotically stable.

$ Case \; ii $. The wooded equilibrium $ E_{1} = (u_{1}, \; 0) $ is a unstable saddle.

$ Case \; iii $. The grassland equilibrium $ E_{2} = (0, \; v_{2}) $ is a unstable saddle.

$ Case \; iv $. The desert equilibrium $ E_{0} = (0, \; 0) $ is a unstable source. It is to be noted that the existence of $ E_{*} $, $ E_{1} $ or $ E_{2} $ destabilizes the desert equilibrium $ E_{0} $.

Obviously, $ {\lambda _1}({E_*}) < 0 $, $ {\lambda _2}({E_*}) < 0 $ when the condition $ ({H}_0) $, $ ({H}_1) $ hold. This implies that the coexistence equilibrium $ E_* $ of (2.3) is locally asymptotically stable. From the above discussion and Lemma 3.1, we have the following conclusion.

Theorem 3.1. If $ ({H}_0) $, $ ({H}_1) $ hold, the coexistence equilibrium $ E_* $ of (2.3) is locally asymptotically stable.

The above result indicates that $ E_* $ attracts all feasible solutions under the condition $ ({H}_0), ({H}_1) $. Biologically, it means that tree species and grass species coexist steadily in the same area for a long time.

3.2.   Hopf bifurcation analysis

In this section, we want to study the Hopf bifurcation of Eq (2.4) caused by time delay $ \tau $. Now multiply both sides of the system (3.3) by $ e^{\lambda\tau} $, then:

suppose Eq (3.5) has a pair of pure imaginary roots, that is $ \lambda = \pm \mathrm{i} w $, $ (w > 0) $, $ \mathrm{i}^{2} = -1 $. Substituting $ \lambda = \mathrm{i} w $ into (3.5), and separating the real and imaginary parts, we have:

the direct calculation results are as follows:

Noting that $ \text{cos}^{2}(\omega \tau)+\text{sin}^{2}(\omega \tau) = 1 $, then:

and $ w $ is the positive real root of (3.8), in which:

From (3.7), we have:

According to the existence theorem of Hopf bifurcation ^[36,37], let $ \tau $ be the Hopf bifurcation parameter, $ \tau _l $ is the critical value of Hopf bifurcation. If $ \text{Re}\lambda = 0, \; \text{Im}\lambda \neq 0 $, then Hopf bifurcation occurs when the system (3.5) satisfies the transversal condition $ \text{Re}\left(\frac{d \lambda }{d \tau }\right)_{\tau = \tau _{\min }}\neq 0 $. Set $ \tau _{\min} = min \left\{\tau _l\right\} $, then we have following theorem:

Theorem 3.2. Let $ k = 0, $ if the conditions $ (H_1) $ and $ \mathit{\text{Re}}\left(\frac{d \lambda }{d \tau }\right)_{\tau = \tau _{\min }}\neq 0 $ hold, then

$ Case \; i $. When $ \tau \in \left[0, \tau _{\min }\right) $, the system (2.4) is locally asymptotically stable near the equilibrium point $ E_{*} $.

$ Case \; ii $. When $ \tau = \tau _{\min } $, the system (2.4) undergoes Hopf bifurcation the equilibrium point $ E_{*} $.

Next, we do a simple numerical simulation to verify the existence of Hopf bifurcation. Under the condition $ (H_{1}) $, fixed parameters $ a = 1.5 $, $ b = 0.5 $, $ B = 1.2 $ and $ \delta = 0.3 $. By solving Eqs (3.8) and (3.9), we get a unique positive real root $ w = 0.34 $, the critical value $ \tau_{min} = 6.8 $ for Hopf bifurcation, the equilibrium point $ E_* = \left(u_*, v_*\right) = (0.6625, 0.225) $, and $ {{d\lambda } / {d\tau }} = 0.00768 - 0.0387{\text{i}} $, then $ {\mathop{\rm Re}\nolimits} {({{d\lambda } / {d\tau }})_{\tau = {\tau _{\min }}}} > 0 $. With the change of time $ t $, we obtain the phase diagrams of the system (2.4) without diffusion according to the different values of delay $ \tau $, see Figure 2. In Figure 2 $ (a) $ and $ (b) $, selected the time delay value $ \tau = 2.8, \; \tau = 4.87 $, belong to the region $ \left[0, \tau _{\min }\right) $, and the system is asymptotically stable nearby $ E_{*} $. In Figure 2 $ (c) $, select the time delay value $ \tau = 6.87 $, belong to the region $ \tau > \tau_{min} $, then the unstable periodic solution appears, and the Hopf bifurcation occurs nearby $ \tau = \tau_{min} $.

3.3.   Turing instability analysis

In this section, we mainly analyze the sufficient conditions of Turing instability. In order to analyze the effect of $ \tau $ on the Turing instability of the system (2.4), we develop $ u(x, y, t-\tau) $ and $ v(x, y, t-\tau) $ into power series^[31], and only the first order term left as follows:

Substituting (3.10) into the system (2.4), then we get:

After calculation, we have:

where $ A_{11}(\tau) = \frac{a_{11}+a_{13}}{a_{13}\tau +1} $, $ A_{21}(\tau) = \frac{a_{21}}{a_{24} \tau +1} $, $ A_{22}(\tau) = \frac{a_{22}+a_{24}}{a_{24} \tau +1} $, $ D_{11}(\tau) = \frac{d_{11}}{a_{13} \tau +1} $, $ D_{12}(\tau) = \frac{d_{12}}{a_{13} \tau +1} $, $ D_{21}(\tau) = \frac{d_{21}}{a_{24} \tau +1} $ and $ D_{22}(\tau) = \frac{d_{22}}{a_{24} \tau +1} $. All the coefficients of the system (3.11) with the cross-diffusion are related to the time delay $ \tau $, so we only need to analyze this new system (3.11). Define:

then the characteristic equation corresponding to the linearization of the system (3.11) is as follows:

where $ \Gamma = D_{11}(\tau)D_{22}(\tau)-D_{12}(\tau)D_{21}(\tau) $, $ L_1 = A_{22}(\tau)D_{11}(\tau)-A_{21}(\tau)D_{12}(\tau)+A_{11}(\tau)D_{22}(\tau), $ and $ L_{2} = A_{11}(\tau)A_{22}(\tau) $.

In order to analyze the conditions of Turing instability, we assume that the condition $ (H_2^1):\; a > b, \; 0 < B < \frac{a}{a+b}, \; \tau < \frac{1}{a \delta } $ holds when the loss rate of tree biomass $ a $ is greater than that of grass biomass $ b $ (for $ a > b $), and the condition $ (H_2^2):\; 0 < a < b, \; 0 < B < \frac{a}{2 b}, \; \tau < \frac{1}{b \delta } $ holds when the loss rate of tree biomass $ a $ is less than that of grass biomass $ b $ (for $ a < b $). Combine the condition $ (H_2^1) $ with $ (H_2^2) $ means the following condition holds:

since $ DET_{k} $ is a curve function of the wave number $ k $. The necessary condition for Turing bifurcation also needs to satisfy the conditions:

When the system (3.11) satisfies the condition $ (H_{1})-(H_{4}) $, it means that there is always a wave number $ k_{c} $ that leads to Turing instability. That is to say, when $ k = k_c\neq 0 $ the Turing instability occurs if $ \text{Re}(\lambda) = 0, \; \text{Im}(\lambda)\neq0 $. Next, we want to find the critical control parameter $ \tau_{c} $ and the critical wave number $ k_{c} $ for Turing instability. From $ DET_{k} = 0 $ and the roots $ k^{2} $ are:

according $ \Lambda = L_1^2-4\Gamma L_2 = 0 $, the critical wave number $ k_c^2 = \frac{L_1}{2 \Gamma} $.

Next, we will discuss the stability of the system (3.11) with passive diffusion but without cross-diffusion (for $ d_{12} = d_{21} = 0 $). Then the system (3.11) becomes:

the general linear analysis shows that the sufficient condition for the generation of Turing pattern is to satisfy conditions $ (H_{2})-(H_{4}) $. However, in the characteristic equation corresponding to the system (3.13) without cross-diffusion, there have

no matter how the wave number $ k^{2} $ is chosen, the curve $ DET_k|_{\left(d_{12} = d_{21} = 0\right)}\geq L_{2} > 0 $, which is contradicted with the condition $ (H_{4}) $. Therefore, when $ d_{12} = d_{21} = 0 $, the Turing pattern of the system (2.4) will not appear, it means that the model (2.4) is always asymptotically stable without cross-diffusion at the positive equilibrium point $ E_{*} $. Therefore, cross-diffusion is essential and indispensable to the Turing instability analysis of the system (2.4). On the other hand, since the coefficients of the system (3.11) are all related to the time delay $ \tau $, the Turing instability of the system (2.4) is also related to the time delay $ \tau $ and the fire frequency $ \delta $. If $ \delta $ is chosen as Turing bifurcation parameter, from $ DET_k = 0 $ we get the value $ \delta_{c} $ as:

where

Theorem 3.3. If conditions $ (H_1)-(H_{4}) $ hold, and the wave number $ k^{2} = m^{2}+n^{2}, \; (m, n)\in \mathbb{N}_{+}^{0} $. Then the system (2.4) undergoes $ k $-mode Turing instability if any of the following conditions is true.

$ Case \; i $. The wave number $ k^{2} $ and Turing bifurcation parameter $ \delta $ are satisfy $ k_c^2 < k^2 < k_+^2 $ and $ \delta _c < \delta < \delta _0 $ respectively.

$ Case \; ii $. The fire frequency $ (\tau, \delta)\in \mathit{\text{Region}}\left\{\eta_1 < 0, \; \eta_2 > 0, \; \eta_3 > 0, \; \eta_4 > 0\right\} $.

Remark 3.1. For the case $ i $ in the Theorem 3.3, we make the following assumptions for the Turing instability analysis of the system (2.4) for convenient. That is, the Turing instability conditions are $ (H_{1})-(H_{4}) $, which would mean that both the condition $ \delta > \delta_{c} $ and the condition $ 0 < \delta < \delta _0 $ hold. Therefore, the range of Turing instability induced by fire frequency is summarized as $ \delta _c < \delta < \delta _0 $.

Remark 3.2. For the case $ ii $ in the Theorem 3.3, we make the following assumptions for the Turing instability analysis of the system (2.4) for convenient. That is, if the time delay $ \tau $ is chosen as Turing bifurcation parameter, set curves $ \eta_1 = TR_0 $, $ \eta_2 = DET_0 $, $ \eta_3 = L_1 $, $ \eta_4 = \Lambda $, and these curves are functions of $ \tau $ and $ \delta $. In addition, these curves will have to satisfy the Turing instability conditions $ (H_{1})-(H_{4}) $, which can be simplified to the region $ \left\{\eta_1 < 0, \; \eta_2 > 0, \; \eta_3 > 0, \; \eta_4 > 0\right\} $.

In order to further determine the conditions of Turing instability caused by the combined effects of diffusion and time delay, the case $ i $ and case $ ii $ in Theorem 3.3 are discussed through numerical simulation. If the conditions $ (H_{2})-(H_{4}) $ hold, the parameters are taken as $ a = 0.95 $, $ b = 1.1 $, $ B = 0.5 $, $ d_{11} = 0.0004 $, $ d_{12} = 0.005 $, $ d_{21} = -0.0001 $ and $ d_{22} = 0.00495 $.

Discussing the relationship between $ k $ and $ \delta $ in Case $ i $ of theorem 3. Set $ \tau = 0.1 $, then we get $ k_{c}^{2} = 135 $, $ \delta_{c} = 0.37697 $ and $ \delta_{0} = 0.52631 $. According to Theorem 3.3, Turing instability occurs in the range $ \delta _c < \delta < \delta _0 $. We discretize $ \delta $ and observe the relationship between $ DET_{k} $ and $ k^{2} $, as shown in Figure 3. The system (2.4) is asymptotically stable and Turing instability does not occur for $ \delta < \delta_{c} $, see the curve $ \delta_{1}, \delta_{2} $. The system (2.4) occurs Turing instability for $ \delta_{c} < \delta < \delta_{0} $, see the curve $ \delta_{4}, \delta_{5} $, and the shaded areas is the Turing instability region. The numerical simulation is consistent with the case $ i $ of Theorem 3.3.

Discussing the relationship between $ \tau $ and $ \delta $ in case $ ii $ of theorem 3, the curves of $ (\tau, \; \delta) $ are expressed as follows:

the Turing instability occurs in the region of conditions $ \eta_1 < 0, \; \eta_2 > 0 $, $ \eta_3 > 0 $ and $ \eta_4 > 0 $. By solving the above inequality equations, we obtain two simplified Turing instability regions with respect to $ (\tau, \; \delta) $, which are $ \text{Region I:}\left\{0.37351 < \delta < 0.52632, \; 0 < \tau \leq 1.7273\right\} $ and $ \text{Region II:}\left\{ 1.7273 < \tau < 2.4339, \; 0.37351 < \delta < \frac{0.9091}{\tau }\right\} $. By merging the $ \text{Region I} $ and the $ \text{Region II} $, the range of Turing instability is obtained as shown in Figure 4, where the gray region represents $ \text{Region I} $ and the blue region represents $ \text{Region II} $.

Next, we simulate the trends of curves $ \eta _1 $, $ \eta _2 $, $ \eta _3 $ and $ \eta _4 $ in these two regions respectively. In Figure 5, selected $ \delta = 0.49 $ and we are easy to find the Turing instability occurs for $ 0 < \tau < 1.34669 $, where Figure 5$ (a) $ represents the trend of curves $ \eta_1 $, $ \eta_2 $, $ \eta_3 $, $ \eta_4 $ with $ \tau > 0 $. Figure 5$ (b) $ represents the truncated part of 5$ (a) $ with $ 0 < \tau < 1.34669 $, where the shaded area is the Turing instability region. In Figure 6, selected $ \delta = 0.4 $ and we are easy to find the Turing instability occurs for $ 0 < \tau < 1.7408 $, where Figure 6$ (a) $ represents the solution of curves $ \eta_1 $, $ \eta_2 $, $ \eta_3 $ and $ \eta_4 $ with $ \tau > 0 $. Figure 6$ (b) $ represents the truncated part of Figure 6$ (a) $ with $ 0 < \tau < 1.7408 $ and the shadow areas is the Turing instability region. The numerical simulation is consistent with the case $ ii $ of Theorem 3.3.

From the above discussions, the Turing instability of system (2.4) and the conditions for the generation of Turing pattern are obtained, but the generation and variation of the Turing pattern are still a problem. Therefore, in the following discussion, we will derive the amplitude equation for the Turing pattern starting nearby $ \tau = \tau_{c} $, which explains the stability of different forms of Turing pattern and the structural transformation between them ^[17].

4.   The amplitude equation of Turing patterns

In this section, we will use multi-scale analysis to derive the amplitude equations of the Turing pattern of the system (2.4) near the case $ k^{2} = k^{2}_{c} $ and $ \tau = \tau_{c} $, and refer to Ouyang Qi^[17] for the detailed information of the derivation process and the relevant conclusions of the pattern. We develop the variables $ u(x, y, t-\tau) $, $ v(x, y, t-\tau) $ into power series of $ \tau $ and only remain the first order terms. The effect of diffusion is always considered, and let $ u(x, y, t-\tau) = u(t-\tau) $, $ v(x, y, t-\tau) = v(t-\tau) $. Then the abstract form of system (2.4) at the equilibrium $ E_{*} $ is:

The Turing pattern of the system (4.1) can be represented by three wave vectors $ {\bf{k}}_{1} $, $ {\bf{k}}_{2} $, $ {\bf{k}}_{3} $ with an included angle of $ 120^{o} $, which are satisfied $ {\bf{k}}_{1}+{\bf{k}}_{2}+{\bf{k}}_{3} = {\bf{0}} $ and $ |{\bf{k}}_{i}|^{2} = {\bf{k}}_c^{2} $, but with different directions. We still use the symbol $ {\bf{U}} = (u, v, \hat{u}, \hat{v})^{T} $ stands for the solution of the system (4.1), and rewrite it at the critical value of Turing bifurcation:

where $ {\bf{r}} = (x, y) $ is the spatial row vector in two-dimensional space $ \Omega $, $ x, y $ is spatial variable, $ c.c. $ denotes the complex conjugate roots of the former term, then the system (4.1) can be rewired as:

where $ \mathfrak{L} $ is linear operator, $ \mathcal{G} $ is nonlinear operator and we have $ \mathfrak{L} = \mathfrak{L}_{c}+\left(\tau_{c}-\tau\right)\mathfrak{M} $ with

where $ a_{31} = \frac{1}{\tau_{c}}, \; \mathfrak{m}_{11} = -\frac{\partial a_{31}}{\partial \tau}|_{\tau = \tau_{c}} = \frac{1}{\tau \tau_{c}} $. The system (4.2) occurs the Turing bifurcation close to the critical values $ {k^2}{\text{ = }}k_c^2 $ and the time delay $ \tau{\text{ = }}{\tau _c} $. We introduce an additional parameter $ \kappa $, expand the bifurcation parameter $ \tau $ near the critical value $ \tau_{c} $ as follows:

where

According to the chain rule of differentiation, we use the method of multi-scale analysis and introduce multiple time scales $ T_{0} = t, \; T_{1} = \kappa t, \; T_{2} = \kappa^{2}t $, and the derivative of time should be transformed into:

Substituting (4.3), (4.4) into the system (4.2), which can be expanded in orders of small perturbation $ \kappa $, then we get the following results:

where $ {\bf{0}} $ denotes the $ 2\times1 $ null matrix. Because $ \mathfrak{L}_{c} $ is the symbol of linear operator of the system (4.2) at the critical point, $ \left(u_1, v_1, \hat{u}_{1}, \hat{v}_{1}\right)^T $ is the linear combination of the eigenvectors correspond to eigenvalue zero. For the first order term of $ \kappa $, we solve Eq (4.6) yields:

where $ z_{1} = z_{2}z_{3} $, $ z_2 = k_c^2 \tau _c+1 $, and $ z_3 = \frac{d_{12} k_c^4 \tau _c+d_{12} k_c^2}{-a \delta +B k_c^2 \tau _c-2 B u_* k_c^2 \tau _c-2 B u_*+B-d_{11} k_c^4 \tau _c-d_{11} k_c^2} $. $ {\bf{\Upsilon}}_{j} $ is the amplitude of mode $ \text{Exp}\left(\mathrm{i}{\bf{k}}_{j}{\bf{r}}\right) $ under the first order perturbation, and its form is determined by the higher order perturbation term. For the second order term of $ \kappa $, we solve Eq (4.7) yields:

From Fredholm solubility condition, the system (4.10) must be orthogonal to the zero eigenvector of the adjoint operator $ \mathfrak{L}_{c}^{*} $ corresponding to $ \mathfrak{L}_{c} $, then the system has a non-trivial solution. The zero eigenvector of $ \mathfrak{L}_{c}^{*} $ is:

Let $ \Theta_{u}^{j}, \; \Theta_{v}^{j}, \; \Theta_{\hat{u}}^{j}, \; \Theta_{\hat{v}}^{j} $ represent the coefficients of $ \Theta_{u}, \; \Theta_{v}, \; \Theta_{\hat{u}}, \; \Theta_{\hat{v}} $ corresponding to $ \text{Exp}(\mathrm{i}{\bf{k}}_{j}\cdot{\bf{r}}) $ in the system (4.9). According to the orthogonality condition:

where $ c_1 = \frac{-a \delta +B k_c^2 \tau _c-2 B u_* k_c^2 \tau _c-2 B u_*+B-d_{11} k_c^4 \tau _c-d_{11} k_c^2}{\left(k_c^2 \tau _c+1\right) \left(d_{21} k_c^2+v_*\right)} $, $ c_2 = -\frac{a \delta \tau _c}{k_c^2 \tau _c+1} $ and $ c_3 = \frac{b c_1 c_2}{a} $, then we get:

The system (4.6) is the amplitude equation under the first-order perturbation. While the coefficient of the second-order term is greater than zero, and the amplitude $ {\bf{\Upsilon}}_{j} $ is divergent function, so we need to introduce the higher-order perturbation to saturate with it. According to the orthogonal condition and the formula in the reference^[19,20], the system (4.10) is calculated as follows:

and we get

where

For the third order of $ \kappa $, we set:

Similarity, we let $ \mathcal{X}_{u}^{j}, \; \mathcal{X}_{v}^{j}, \; \mathcal{X}_{\hat{u}}^{j}, \; \mathcal{X}_{\hat{v}}^{j} $ represent the coefficients of $ \mathcal{X}_{u}, \; \mathcal{X}_{v}, \; \mathcal{X}_{\hat{u}}, \; \mathcal{X}_{\hat{v}} $ corresponding to $ \text{Exp}(\mathrm{i}{\bf{k}}_{j}\cdot{\bf{r}}) $ in the system (4.13). According to the orthogonality condition, we have:

where

Let $ {\bf{E}}_j = {\bf{E}}_j^u = z_3 {\bf{E}}_j^v = z_2 {\bf{E}}_j^{\hat{u}} = z_1 {\bf{E}}_j^{\hat{v}} $ be the coefficients of $ \text{Exp}(\mathrm{i}{\bf{k}}_{j}\cdot{\bf{r}}) $, then:

Multiplying (4.11) and (4.14) by $ \kappa $ and $ \kappa^{2} $, respectively, and considering together with the variables in (4.5) and (4.15), the amplitude equation corresponding to $ {\bf{E}}_{j} $ is obtained:

where $ \mu = \frac{\tau-\tau_{c}}{\tau_{c}} $, $ \varsigma_{0} = \frac{\tau \left(c_1 z_2+c_2 z_3+c_3+z_1\right)}{\left(z_2-1\right) \left(c_2 z_3+c_3\right)} $, $ h = \frac{\tau \left(-2 B z_1^2-2 c_1z_2 \left(z_1+z_2\right)\right)}{z_1 \left(z_2-1\right)\left(c_2 z_3+c_3\right)} $, $ g_1 = \frac{G_2^1 \tau }{z_1^2\left(z_2-1\right) \left(c_2 z_3+c_3\right)} $, $ g_2 = \frac{G_2^2 \tau }{z_1^2\left(z_2-1\right)\left(c_2 z_3+c_3\right)} $.

Next, we analyze the stability of the amplitude equation. Each amplitude in system (4.16) can be expressed as $ {\bf{E}}_j = \theta _j\text{Exp}\left(\mathrm{i}\rho _j \right), \; j = 1, 2, 3 $, and $ \theta_{j} = |{\bf{E}}_j| $, $ \rho_{j} $ represents the corresponding phase angle. Substituting $ {\bf{E}}_j $ into Eq (4.16) and separating the real and imaginary parts yields four differential equations of the real variables as follows:

where $ \rho = \rho _1 +\rho _2 +\rho _3 $. Next, the linear stability of the amplitude equation (4.17) and the Turing pattern are analyzed. According to the reaction-diffusion pattern dynamics theory^[17], the properties are summarized as follows:

Proposition 4.1. The amplitude equation (4.17) satisfies

$ Case \; i $. If the system is in steady state $ (\rho > 0) $, the sum of pattern phase can only take steady states $ \rho = 0 $ and $ \rho = \pi $. In which, the solution of $ \rho = 0 $ corresponding to pattern $ H_{0} $ is stable for $ h > 0 $, and the solution of $ \rho = \pi $ corresponding to pattern $ H_{\pi} $ is stable.

$ Case \; ii $. The spot pattern are $ H_{0} $ and $ H_{\pi} $, which is introduced by subcritical bifurcation of the equation. The stripe pattern is $ S_{p} $, which is introduced by the supercritical bifurcation of the equation.

$ Case \; iii $. If the mode equation has a steady-state solution, the coefficients of the second order term should satisfy $ g_{1} > 0, \; g_{2} > 0, $ and $ \frac{g_{2}}{g_{1}} > 1 $.

$ Case \; iv $. The bifurcation parameter of the amplitude equation is $ \mu $, the formula of the control parameters are:

Proposition 4.2. If $ h\neq0 $ holds, and the relationship between the critical values is $ {\mu _1} < {\mu _2} < {\mu _3} < {\mu _4} $, then the system (4.17) has five kinds of stationary solutions, as shown below.

$ Case \; i $. The stationary state: $ \theta _1 = \theta _2 = \theta _3 = 0 $. The stationary state is always stable for $ \mu < \mu _2 = 0 $.

$ Case \; ii $. Stripe pattern: $ {\theta _1} = \sqrt{\frac{\mu}{g_{1}}}, \; {\theta _2} = {\theta _3} = 0 $. The stripe pattern is always stable for $ \mu > \mu _{3} = \frac{h^{2}g_{1}}{(g_{1}-g_{2})^{2}} $.

$ Case \; iii $. When $ \theta^{-} = \frac{|h|\pm \sqrt{h^{2}+4\left(g_{1}+2g_{2}\right)\mu}}{2\left(g_{1}+2g_{2}\right)} $, we have two kinds of spot pattern $ H_{0} $ and $ H_{\pi} $ exist for $ \mu > \mu _1 $, where The spot pattern $ H_{\pi} $ is always unstable for $ h > 0 $. In which, for the stationary state $ \theta^{-} = \frac{|h|-\sqrt{h^{2}+4\left(g_{1}+2g_{2}\right)\mu}}{2\left(g_{1}+2g_{2}\right)} $, then the corresponding pattern will never stable, and for the stationary state $ \theta^{+} = \frac{|h|+\sqrt{h^{2}+4\left(g_{1}+2g_{2}\right)\mu}}{2\left(g_{1}+2g_{2}\right)} $, then the corresponding pattern is stable for $ {\mu _1} < \mu < {\mu _4} $

$ Case \; iv $. There is spot-stripe mixture pattern when $ {g_2} > {g_1} $, one is $ \theta_1 = \frac{|h|}{g_{2}-g_{1}} $, the others are $ {\theta _2} = {\theta _3} = \sqrt {\frac{\mu-g_1\theta_1^2}{g_{2}+g_{1}}} $, and it is always unstable for $ \mu > {g_1}\theta _1^2 $.

5.   System identification

In this subsection, we make sensitivity analysis on the system (2.4) with $ k = 0 $, and discuss the influence of fire frequency $ \delta $ to the density distribution of trees and grasses. Then we studied the effect of fire frequency $ \delta $ to the Turing pattern generation of tree-grass coexistence system (2.4) on the basis of numerical simulation by Matlab, which were based on the above sensitivity analysis results. For the tree-grass model (2.4), we are very interested in the self-organizing variation of its Turing pattern, so we chose the Neumann boundary condition(zero flux), which means that external input is not considered.

5.1.   Sensitivity analysis

Sensitivity analysis is a method to judge the sensitivity of change in a state or output variable of nonlinear dynamic model to system environmental parameters. Therefore, sensitivity analysis of tree-grass model parameters is a very important step in the generation of vegetation pattern. The model (2.2) involves many parameters, among which the fire frequency $ \delta $ has a great effect on tree-grass coexistence. and the delay effect $ \tau $ of fire frequency also cannot be ignored.

In order to verify whether the system (2.4) can truly reflect the influence of $ \delta $ on the density distribution of trees and grasses, we will make sensitivity analysis on the model (2.4) without considering the effect of spatial diffusion. Refer to the parameters value shown in Table 2 and the other parameters are $ \tau = 1.23 $, $ B_{1}(B) = r_{w}/r_{g} = 0.533 $, $ \delta = 1.7 $, Time step $ \Delta t = 2500 $. Reducing mathematical expression by simplifying variables, the sensitivity equation of model (2.4) with it$ k = 0 $ is:

where the sensitivity of trees and grasses to fire frequency are set as:

To eliminate influence of dimension the relative sensitivity of trees and grasses to fire frequency are $ \delta S_{1}(t)/ u(t) $, $ \delta S_{2}(t)/v(t) $. We will use the $ RK-4 $ method to carry out sensitivity analysis for the system (5.1), and the variation trend of the sensitivity and relative sensitivity of trees $ u(t) $ and grasses $ v(t) $ to fire frequency was obtained, as shown in Figures 7 and 8.

In Figures 7(a) and 8(a), both the sensitivity $ S_{1} $ and the relative sensitivity $ \delta S_{1}(t)/ u(t) $ to $ \delta $ showed a trend of less than $ 0 $, which means that the biomass of trees decreases as the increases of fire frequency. According to discrete data of the program, when $ t = 20 $, we have $ S_{1} = -0.185 < 0 $ and $ \delta S_{1}(t)/ u(t) = -0.4617 $. This indicates that when the fire frequency $ \delta $ increases by $ 10\% $, which will reduce the biomass of trees by $ 4.617\% $. From the perspective of biology, model (2.4) regard adult trees and young trees as the tree species, and believe that the extinction rate of tree species is not only inhibited by its growth, but also directly related to fire disturbance, so the occurrence of fire will lead to the death of some trees, especially young trees. Some scholars have found that adult trees are hardly affected by the incipient fire, because the water content of trees is higher than that of grasses, and the bark protects the crown and stem higher than the flame area ^[38]. When the frequency of fire increases, adult trees will gradually become the fuel for fire and support its continuous burning. However, the tree growth cycle is too long, resulting in the tree species can not recover in time, and the tree biomass will decrease steadily, interspersed with relatively rare supplementary events ^[31]. If the fire frequency exceeds the threshold, it will eventually lead to the extinction of tree species, and the extreme phenomenon of complete grassland or desertification will occur in the forest steppe ecotone.

In Figures 7(b) and 8(b), both the sensitivity $ S_{2} $ and the relative sensitivity $ \delta S_{2}(t)/ v(t) $ to $ \delta $ showed a trend of being less than $ 0 $, then greater than $ 0 $, and finally stabilized in this region. According to discrete data of the program, we divided the above change trend into two stages: when $ t \in [0, \; \; 1.7] $, we have $ {S_2}(t) \le 0 $, $ \delta S_{2}(t)/ v(t)\le 0 $, and when $ t > 1.7 $ we have $ {S_2}(t) > 0 $, $ \delta S_{2}(t)/ v(t) > 0 $. In which, when $ t = 0.8 $, the sensitivity and relative sensitivity of grass species to fire frequency $ \delta $ reach the minimum value $ {S_2}(t) = - 0.004254 $ and $ \delta S_{2}(t)/ v(t) = -0.02987 $. In other words, the fire frequency $ \delta $ increases by $ 10\% $, which will reduce the biomass of grass species by $ 0.2987\% $. When $ t{\text{ = 25}} $, the sensitivity and relative sensitivity of grass species to fire frequency are $ {S_2}(t) = 0.1362 $, $ \delta S_{2}(t)/ v(t) = 0.9762 $, and finally reach a stable equilibrium. In other words, the fire frequency $ \delta $ increases by $ 10\% $, the biomass of grasses will increase by $ 9.762\% $. From the perspective of biology, another state variable in the model (2.4) is grass species, whose extinction rate is not only inhibited by its own growth, but also directly related to the growth of tree species and fire feedback. The grass species is the main fuel for fire, and it needs enough grass biomass accumulation to burn continuously, so grasses will be significantly reduced in the early stage of fire. However, the short growth cycle of grasses makes it able to the recovery partial growth within the interval of fire frequency rapidly. In addition, with the increase of fire frequency, the biomass of trees always decreases at a stable rate. When the dominant competitor(trees) decreases, the grass biomass begins to increase slowly due to the lack of inhibition by tree species, and finally reaches a stable balance. The above ecological phenomena exactly correspond to the change trend that grass biomass first decreases and then increases with the increase of fire frequency.

In conclusion, the theoretical analysis of the sensitivity and relative sensitivity of tree-grass biomass to fire frequency $ \delta $ accords with the practical ecological significance. The analysis shows that the regeneration capacity of trees, grasses and the fire frequency $ \delta $ are the key to predict the coexistence and balance of species in the Savanna, and the model (2.4) without diffusion can well reflect the different effects of fire frequency on trees and grasses. Next, we take the time delay $ \tau $ as the Turing bifurcation parameter, and carry out numerical simulation on the spatial Turing pattern introduced by diffusion in the system (2.4) to verify the theoretical analysis.

5.2.   Numerical simulation and pattern structure

In this section, the Matlab program is used to verify the influence of fire frequency $ \delta $ and time delay $ \tau $ on the Turing pattern of the model (2.4). For the tree-grass model, we are very interested in the self-organizing changes of its Turing pattern, so we select Neumann boundary conditions in specific numerical simulation, which means there is no external input.

According to the significance of cross-diffusion $ d_{21} $ of the tree-grass model (2.4) in the second part $ 2 $, we performed numerical simulation with the two aspects of cross-diffusion coefficient $ d_{21} > 0 $ and $ d_{21} < 0 $, respectively. The numerical algorithm of reaction-diffusion system is often complicated and requires tedious Matlab programming. For a more comprehensive theoretical analysis of the effects of $ \delta $ and $ \tau $, we use the different programs to approximate the solutions of reaction-diffusion model (2.4) and (4.1). Among them, the semi-implicit difference method (with respect to time) is used for $ d_{21} > 0 $, and the finite difference method is used for $ d_{21} < 0 $.

The case of $ d_{21} > 0 $

We mainly analyze the influence of the time delay $ \tau $ on the Turing pattern of the system (2.4) with the semi-implicit difference method, which is an improvement on the method of Garvie ^[39]. The method involves approximations at the current time level $ t_{n} $ and the previous time level $ t_{n-1} $, introduce delay $ \tau $ in the iteration process, and finally get the sparse, banded and linear system of algebraic equations. Take $ \Omega \subset {\bf{R}}^{2} $, we carry out the simulation algorithm on a discrete square region with sides of $ [400\times400] $, the distance between adjacent lattices is $ \triangle h = 0.25 $, and the time step is $ \triangle t = 10^{-6} $. The diffusion coefficients are $ d_{11} = 0.4 $, $ d_{12} = 0.5 $, $ d_{21} = 0.396\; (d_{21} > 0) $ and $ d_{22 } = 0.495 $. Let $ a_{1} = a \delta $, $ b_{1} = b \delta $ and the other related parameters reference to Table 2, then we have $ {E_*} = (0.675, 0.309) $, $ \tau _c = 1.23 $ and $ k_c^2 = 884 $.

Let $ \tau = \epsilon_{\tau} \triangle t $ and $ t_{\epsilon} = \epsilon \triangle t $. First, we denote $ \hat{f}_{i} $ and $ \hat{g}_{i}, \; (i = 1, 2) $ as the discrete functions corresponding to $ f_{i} $ and $ g_{i} $, respectively. Discretizing system (2.4) as follows:

where $ i, j = 0, 1, \cdots, \frac{L}{\triangle h} $. Denote

and we get

where the constant matrices $ {\bf{B}}_{1} $ and $ {\bf{B}}_{2} $ can be referred to ^[39].

According to the formula (4.18), we have $ \mu _1 = -0.292578 < \mu_{ 2} = 0 < \mu_{ 3 } = 5.66835 < \mu_{ 4 } = 22.7874 $, $ h = 138.065 > 0 $ and $ g_{2}/g_{1} \approx 2.02011 $. From Proposition 4.1, the pattern $ {H_0}(\rho = 0) $ is always stable, and from Proposition 4.2, the expression of the stripe patter $ {S_p} $, the spot pattern $ {H_\pi} $ and $ {H_0}^{\pm} $ are as follows:

Then the bifurcation diagram of the amplitude equation (4.17) is as follows: Then the bifurcation diagram of the amplitude equation (4.17) is as follows:

Figure 9 shows the bifurcation of the model (4.17) for $ h > 0 $, in which the amplitude is a function of the control parameter $ \mu $, which is used to describe the bifurcation of stripe pattern $ S_p $ and spot pattern $ H_{0}^{\pm} $, respectively. when $ \mu > {\mu _1} $, the pattern $ {H_0} $ and the uniform state are always stable, then there exists a first bistable region $ {\mu _1} < \mu < {\mu _2} $ with respect to $ \mu $. When the parameter gradually increase from $ {\mu _1} < 0 $ to $ {\mu _2} = 0 $, the steady state begins to lose stability, and part of the spot pattern in the region $ {\mu _2} < \mu < {\mu _3} $ deforms into the stripe, that is, the stripe pattern $ S_p $ appears but it is unstable. With the change of parameters, the stripe pattern $ S_p $ tends to be stable when $ \mu > {\mu _3} $. At this time, both pattern $ {H_0} $ and stripe pattern $ S_p $ are stable, and there is a second bistable region $ {\mu _3} < \mu < {\mu _4} $ with respect to $ \mu $. When $ \mu > {\mu _4} $, the spot pattern $ {H_0} $ begins to lose its stability. According to Proposition 4.1, 4.2, then we have the following proposition:

Proposition 5.1. Under the Turing conditions, the time delay $ \tau $ is chosen as the critical value of Turing bifurcation with $ \tau = \tau_{c}(1+\mu) $ and $ k^2 = m^2+n^2 $, $ (m, n) \in \mathbb{N}_{0}^{+} $. Then the Turing instability of system (2.4) occurs near $ \tau = \tau_{c} $, and we have the following conclusions.

$ Case \; i $. When $ \tau \in\left(\left(1+\mu _1\right) \tau _c, \; \tau _c\right) $, which is correspond to the interval $ \mu \in(\mu _1, \; \mu _2) $, and there exist stable spot pattern $ H_{0} $ and stable stationary state.

$ Case \; ii $. When $ \tau \in\left(\tau _c, \; \left(1+\mu _3\right) \tau _c\right) $, the spot pattern $ H_{\pi} $ appears and is always unstable for $ \tau > \tau _c $, which is correspond to the interval $ \mu > \mu_{2} $.

$ Case \; iii $. When $ \tau \in\left(\tau _c, \; \left(1+\mu _3\right) \tau _c\right) $, the stripe pattern $ S_p $ appears, which is unstable in interval $ \tau \in \left(\tau _c, \; \left(1+\mu _3\right) \tau _c\right) $ and stable in interval $ \tau \in \left(\left(1+\mu _3\right) \tau _c, \; \left(1+\mu _4\right) \tau _c\right) $.

$ Case \; iv $. When $ \tau \in \left(\left(1+\mu _3\right) \tau _c, \; \left(1+\mu _4\right) \tau _c\right) $, the spot pattern $ H_{0}^{+} $ and the stripe pattern $ S_p $ are all stable, which are correspond to the interval $ \mu \in(\mu _3, \; \mu _4) $.

$ Case \; v $. When $ \tau \in \left(\left(1+\mu _4\right) \tau _c, \infty\right) $, the spot pattern $ H_{0}^{+} $ is no longer stable, while the stripe pattern $ S_p $ is still stable, which are correspond to the interval $ \mu > \mu _4 $.

In the numerical simulation, the evolution process of different types of trees $ u(x, y, t) $ and grasses $ v(x, y, t) $ patterns can be observed. Here, we only give the tree species density distribution pattern $ u(x, y, t) $, see Figures 10 and 11.

In Figure 10, take $ \tau = 1.2 $, from Proposition 5.1 we have $ \tau \in ({\tau _c}(1{\text{ + }}{\mu _1}), \; {\tau _c}) = (0.87, 1.23) $, which is equivalent to $ \mu = -0.02439\in\left(\mu_{1}, \; \mu_{2}\right) $. It shows the evolution process of the spatial patterns of tree species $ u(x, y, t) $ after $ 0 $ iteration, $ 80,000 $ iterations and $ 300,000 $ iterations, in which the spot pattern $ H _0 $ occupied the whole domain eventually, and we find that this pattern needs a long iteration to stabilize. Starting from Figure 10$ (a) $ with a homogeneous state $ E_{*} $, and the random interference will cause more and more noticeable spots appears eventually, see Figure 10$ (b)(c) $. In fact, we should pay attention to the situation that $ \mu = -0.02439 $ is very close to the critical value $ \mu_{2} $. According to the analysis above, there should be only $ H _0 $ spot patterns under this situation. In other words, the numerical simulations illustrate well the theoretical analysis. From the viewpoint of biology, in the case of hot spot patterns, the tree species $ u(x, y, t) $ lies in the isolated cycle region with high-density, and the remainder region is composed of low-density.

In Figure 11, take $ \tau = 2.02 $, from Proposition 5.1 we have $ \tau \in ({\tau _c}(1{\text{ + }}{\mu _3}), \; {\tau _c}(1{\text{ + }}{\mu _4})) = (1.82, 2.93) $, which is equivalent to $ \mu = 8.9178\in\left(\mu_{3}, \; \mu_{4}\right) $. It shows the evolution process of the spatial patterns of the tree species $ u(x, y, t) $ species after $ 0 $ iteration, $ 60,000 $ iterations and $ 300,000 $ iterations, in which the spot pattern $ H _0 $ and the stripe pattern $ S_{p} $ occupied the whole domain eventually. Through random disturbance and long time iteration, the spatial pattern is finally stable in the coexistence of spot and stripe, which is well illustrated by the numerical simulation, see Figure 11$ (a) $. When the control parameter $ \tau $ is far away from the Turing critical value $ \tau_{c} $, some spot patterns gradually evolve into stripe patterns, see Figure 11$ (b) $ and $ (c) $. Theoretically, there is a bistable state between spot pattern and stripe pattern. With the increase of $ \tau $, the pure stripe pattern does not exist, because the amplitude equation can only explain the spatial pattern near the critical value $ \tau_{c} $.

The case of $ d_{21} < 0 $.

We mainly analyze the influence of the time delay $ \tau $ on the Turing pattern of the system (4.1) with the finite difference method. Suppose the boundary is a discrete square region $ \Omega = [2\times2] $, and the distance between adjacent lattices is $ \triangle h = 0.02 $. The diffusion coefficients are $ d_{11} = 0.0004 $, $ d_{12} = 0.005 $, $ d_{21} = -0.0001\; (d_{21} < 0) $, $ d_{22} = 0.00495 $, and the other parameters are $ B = 0.5 $, $ a_{1} = 0.95 $, $ b_{1} = 1.1 $ and $ (u_{*}, \; v_{*}) = (0.63, \; 0.392) $, then $ \tau_{c} = 1.8553 $, $ k_{c}^{2} = 330 $. Refer to the formula 4.18, we have $ \mu _1 = - 0.0004 < \mu_{ 2} = 0 < \mu_{ 3 } = 0.2415 < \mu_{ 4 } = 0.6918 $, $ h = 138.065 > 0 $ and $ g_{2}/g_{1} = 1.156 > 1 $.

From Proposition 4.1, 4.2, the spot pattern $ H_{0} $ is always stable, the expressions of $ S_p $, $ H_{\pi} $ and $ H_{0}^{\pm} $ refer to (5.3), the theoretical analysis of Turing bifurcation for the control parameters $ \mu $ of the system (4.17) is consistent with subsection $ 5.2.1 $, and the bifurcation diagram refer to Figure 9. For the convenience of recording simulation results, we set $ u(x, y, t) = {u_1}(x, y, t), \; v(x, y, t) = {v_1}(x, y, t) $, $ \hat{u}(x, y, t) = {u_2}(x, y, t) $ and $ \hat{v}(x, y, t) = {v_2}(x, y, t) $. Since variables $ u_2(x, y, t) $ and $ v_2(x, y, t) $ are obtained by expanding the power series with delay $ \tau $ in variables $ u_1(x, y, t) $ and $ v_1(x, y, t) $ respectively. Therefore, in the simulation results, the pattern $ u_2(x, y, t) $ and $ v_2(x, y, t) $ can reflect the pattern of variables $ u_1(x, y, t) $ and $ v_1(x, y, t) $ from another Angle.

Figure 12 shows the density pattern of variables $ u_{1}(x, y, t), v_{1}(x, y, t), u_{2}(x, y, t) $ and $ v_{2}(x, y, t) $ with the time steps $ \Delta t = 5\times10^{-10} $, $ m = 11 $, $ n = 15 $ at the time $ T = 1000 $. Selected $ \tau = 1.8555 $ belong to the region $ \left(\tau_c\left(1+\mu _1\right), \; \tau _c\right) = (1.8553, \; 1.8556) $, that is, the control parameter $ \mu \in \left(\mu _1, \; \mu _2\right) $. At this time, there exist two kinds of pattern: one is a stable spot pattern, the other is a stationary state, and the numerical simulation results agree well with theory analysis.

Figure 13 shows the density pattern of variables $ u_{1}(x, y, t), v_{1}(x, y, t), u_{2}(x, y, t) $ and $ v_{2}(x, y, t) $ with the time steps $ \Delta t = 5\times10^{-10} $, $ m = 11 $, $ n = 15 $ at the time $ T = 1000 $. Selected $ \tau = 2.2556 $ belong to the region $ \left(\left(1+\mu _3\right)\tau _c, \; \left(1+\mu _4\right)\tau _c\right) = (2.12962, \; 2.63433) $, that is, the control parameter $ \mu \in \left(\mu _3, \; \mu _4\right) $. At this time, there exist a stable mixture pattern with coexistence of spots and stripes, and the numerical simulation results agree well with theory analysis.

In order to more clearly observe the influence of time delay $ \tau $ on the spatial self-organization structure change of trees and grasses, take $ \tau = 2.3556 $ in $ \mu \in \left(\mu _3, \mu _4\right) $ and simulate the snapshot graph of population density of variables $ {u_1}(x, y, t) $, $ {u_2}(x, y, t) $ at time $ T = 50 $, $ T = 200 $, $ T = 500 $ and $ T = 1000 $, respectively, See Figures 14 and 15. In this region, the pattern structure of $ {u_1}(x, y, t) $, $ {u_2}(x, y, t) $ are similar to $ {v_1}(x, y, t) $, $ {v_2}(x, y, t) $ respectively, and all of them always present a stable spot-stripe mixture pattern, but only the slightly change in value, so the simulated pattern of $ {v_1}(x, y, t) $ and $ {v_2}(x, y, t) $ are ignored. From the biological point of view, with the increase of the frequency of fire cycle $ \tau $, which means that fires are getting less frequent. To some extent, the fire will promote the coexistence balance of trees and grasses, and the density distribution of tree-grass will present the spot-stripe mixture pattern, and finally reach the coexistence balance. The numerical simulation results are consistent with the theoretical analysis.

5.3.   Tree-grass spatial pattern in Hulunbuir

From the analysis of tree-grass simulation pattern in part 5.2, when the generation conditions of Turing pattern are true, the range of environmental parameters is determined by theoretical analysis. With the change of control parameters $ \mu $ (related to fire frequency $ \delta $) near the critical point, the tree-grass density distribution shows a pattern of hexagonal(spot), stripe and spot-stripe mixture. Hulunbuir is located in eastern Inner Mongolia in China, in a semi-arid and semi humid area, and it has a rich spatial patterns of the forest-steppe ecotone. Go deep into the forest-steppe ecotone in Hulunbuir and observe the spatial self-organization structure of tree-grass vegetation distribution by UAV, as shown in Figures 16 and 17.

Figure 16 shows tree-grass vegetation pattern is located at $ 120^{o}38^{'}\; E $ and $ 49^{o}05^{'}\; N $ on July $ 10 $, $ 2021 $, which is taken by UAV. By observing the vegetation pattern in the forest-steppe ecotone, it is verified that the spatial self-organizing structure of tree-grass presents spot pattern.

Figure 17 shows tree-grass vegetation pattern is located at $ 120^{o}31^{'}\; E $ and $ 49^{o}05^{'}\; N $ on June $ 30 $, $ 2021 $, which is taken by UAV. By observing the vegetation pattern in the forest-steppe ecotone, it is verified that the spatial self-organizing structure of tree-grass presents spot-stripe mixture pattern.

From the comparison, it is found that the simulation results of Turing pattern in theoretical analysis are consistent with the spatial patterns of tree-grass density distribution in Hulunbuir ecotone, such as the simulation of Turing spot pattern in Figure 10, Figure 12 correspond to the actual space-construction of tree-grass density distribution in Figure 16, the simulation of Turing spot-stripe mixture pattern 11, 13, 14 and 15 correspond to the actual space-construction of tree-grass density distribution in Figure 17. Although the vegetation patterns obtained by numerical simulation are likely to be overly regular and change over time, so that it is difficult to find the tree-grass vegetation spatial structure in the actual background which is highly similar to it, however, these types of tree-grass pattern (spot, stripe and spot-stripe mixture) do exist in the forest-steppe ecotone. Therefore, the dynamic properties of vegetation pattern induced by Turing bifurcation can predict the biological characterization of spatial structure of tree-grass density distribution.

6.   Conclusions

In this paper, we study a tree-grass competition model (2.4) with diffusion term and delay effect $ \tau $ caused by fire frequency. We mainly discuss the Hopf bifurcation induced by delay, the Turing bifurcation induced by the diffusion, and the dynamic behaviors of the system near the bifurcation critical value mathematically. By analyzing the characteristic equation corresponding to the linearization, the sufficient conditions of Turing instability are obtained. Based on multi-scale analysis, the amplitude equation of the system (2.4) is derived, which helps us to determine the selection and competition of Turing pattern. The influence of trees and grasses to fire frequency $ \delta $ was obtained through sensitivity analysis, which showed that changes of $ \delta $ affected pattern generation of tree-grass, and the rationality of the tree-grass model (2.4) is confirmed. Finally, a series of numerical simulation experiments are carried out to try to understand the characteristics of tree-grass vegetation pattern under the biological background. With the effect of fire frequency, the density distribution of trees and grasses will eventually present different types of spatial steady-state patterns (include spot, stripe and spot-stripe mixture), which is consistent with the spatial structure of tree-grass density distribution in Hulunbuir ecotone, China (see section 5.3 for details).

In the spatial structure of vegetation, the connection and fusion of pattern of tree-grass density distribution are caused by the mathematical mechanism of bifurcation in the model (2.4). At present, we mainly focus on the relationship between the Turing pattern and bifurcation parameters, and find that the Turing pattern is highly dependent on these parameters (include fire frequency and cross-diffusion). For our study, the fire frequency was taken as the control parameter, so the Turing pattern has a great relationship with the fire frequency, moreover the direct factor affecting fire frequency is biomass diffusion (cross-diffusion). Generally speaking, if we can establish the corresponding relationship between state variables and environmental parameters of tree-grass model, then it can provide certain guidance for predicting the regularity of fire frequency feedback intensity on the coexistence and balance of trees and grasses by comparing the pattern structures, so as to reduce unnecessary loss of human resources.

Acknowledgments

This research was supported by Doctoral Scientific Research Foundation of Inner Mongolia Minzu University (No.BS648).

Conflict of interest

The author declares that they have no competing interests.