Nonlocal delay gives rise to vegetation patterns in a vegetation-sand model

1.   Introduction

Desertification is land degradation in dry semi-humid arid, semi-arid, and arid areas resulting from various factors, including climatic variation and human activities ^[1]. This process can bring about significant negative consequences for society and the economy, which include reduced agricultural output, destruction of property, and heightened health and safety risks ^[2,3]. Numerous ecosystems in arid and semi-arid regions are experiencing accelerated desertification ^[4,5]. In areas with frequent movement of sand caused by wind, desertification is primarily attributed to the displacement of sand, which leads to the depletion of vegetation cover ^[6,7].

However, regarding the research on vegetation modeling, more scholars focus on vegetation water modeling and less on vegetation-sand modeling. They believe that rainfall rate can induce transitions among bare soil state, vegetation pattern state, and homogeneous vegetation state ^[8,9,10,11]. In 1999, Klausmeier ^[12] first proposed a model with two variables, vegetation biomass and water density, for forming regular vegetation patterns in semi-arid areas. The results showed that the regularity of vegetation patterns in semi-arid regions is the traveling wave instability of the reaction-diffusion convection equation. The proposed model has facilitated the study of vegetation patterns in semi-arid areas. Meron et al. ^[13] presented a continuum model for vegetation patterns in water-limited systems. The model predicted transitions from bare soil at low precipitation to homogeneous vegetation at high rainfall through intermediate states of spot, stripe, and gap patterns.

In the past century, with significant growth and widespread occurrence of soil desertification and numerous sand transport and deposition events, the relevance of aeolian sand as a significant environmental concern has only emerged in the last two decades ^[14]. There is a growing recognition among people that vegetation patch formation using wind as a driver is important when dealing with environmental issues ^{[15,16,17,18,19]}. Zhang et al. ^[20] presented a model based on the two variables of vegetation cover and the height of aeolian sand accumulation

where $\frac{{\partial S}}{{\partial T}}$ and $\frac{{\partial V}}{{\partial T}}$ are the accumulation rate of sand and the growth rate of vegetation respectively, $K_0+MV(1-\frac{V}{V_0})-NS$ represents deposition by vegetation, $A_1\frac{\partial S}{\partial X}$ and $A_2\frac{\partial V}{\partial X}$ represent advection by prevailing wind and dispersal by prevailing wind respectively, $HV(1-\frac{V}{Vm})$ is vegetation growth, $PS\frac{V}{C+V}$ represents vegetation destroyed by sand, and $D_1(\frac{\partial^2 S}{\partial X^2}+\frac{\partial^2 S}{\partial Y^2})$ and $D_2(\frac{\partial^2 V}{\partial X^2}+\frac{\partial^2 V}{\partial Y^2})$ represent sand diffusion in all directions and dispersal in all directions respectively.

Zhang et al. ^[20] studied the effect of wind-sand on the spatial distribution of vegetation in the windy sand environment. The influence of the prevailing wind on the growth rate of the two variables is modeled as an advection term in the model, while the effect of other winds is modeled as a diffusion term. Many areas may not have significant prevailing winds. We can obtain the following model

Zhang et al. ^[20] considered that the aggregation effect of vegetation on sand decreases with the increase of vegetation biomass. Different from Zhang et al. ^[20], this paper considers the mechanism of aggregation effect of the linear action term, and the model (1.1) can be rewritten as follows

For ease of analysis, let

The dimensionless system (1.2) is obtained as follows

where $s$ and $v$ are the height of the sand accumulation layer and vegetation coverage respectively, $h$ describes the intrinsic growth rate of vegetation cover, $v_m$ represents the potential maximum value of vegetation coverage, $p$ is the coefficient of destructive effects by sand burial, $c$ is the half-saturation constant of sand capacity, $d$ denotes the ratio of the diffusivity of vegetation to that of sand, the term $ps\frac{v}{1+cv}$ models local destruction of vegetation by sand, and $\Delta$ is a standard Laplace operator. The vegetation-sand system reflects the influence of sand on vegetation.

Nonlocal interaction is one of the critical mechanisms of vegetation pattern formation. Many scholars have researched the nonlocal effects of vegetation models ^{[21,22,23,24,25]}. However, no one has yet studied the nonlocal interaction of sand and vegetation. Fine aeolian sand particles are deposited on the surface of leaves and stems, which will affect the life process of vegetation photosynthesis to a certain extent, and the deposited sand will reduce the water content of the soil and inhibit the growth of vegetation. In addition to the location of the vegetation itself, the aeolian sand in the vicinity and even the whole study area will destroy the development of the foliage. During the transportation of the aeolian sand, it will impact and abrase the stems and leaves of the vegetation, directly damaging the vegetation tissue and affecting its productivity. The accumulation of a large amount of aeolian sand may even bury the vegetation. Hence, it is interesting to study the nonlocal interaction between vegetation and sand. This paper introduces a nonlocal interaction in model (2.1) to better account for the effect of aeolian sand on vegetation growth throughout the study area.

In reality, time delays are inevitable and substantially impact the dynamics ^{[26,27,28,29,30]}. It takes time for sand to be transported by wind to destroy vegetation, which shows that in the spatial scope of our study, the whole process of vegetation destroyed by sand will be affected by nonlocal effects and time delay. In the modelling of complex systems portraying vegetation-sand interactions, the spatially weighted average of nonlocal delay including the entire study area, first proposed by Britton, can be used to solve the problem that the study object is not at the current location of the study area at the current moment at the previous moments ^[31]. This reaction-diffusion equation with delay including a spatially weighted average over the entire study region is called the nonlocal reaction-diffusion equation. Nonlocal reaction-diffusion equations have been widely used in the fields of infectious diseases and predation by predation among animals, but in the field of eco-vegetation most of them are focused on the study of local and global existence and asymptotic approximation of the model solutions, as well as wavefront solutions and periodic wavefront solutions ^{[32,33,34,35,36,37]}. The use of nonlocal delay reaction-diffusion equations to study the vegetation patchwork in arid and semi-arid regions is a new research direction in ecological vegetation science, which can truly portray the actual situation of nonlocal action of vegetation ^{[38,39,40,41,42,43]}, and its theory about combating desertification still needs to be improved. Therefore, this paper establishes a vegetation-sand model with nonlocal effects in arid and semi arid regions, analyzes the model for Turing instability, reveals the influence of vegetation model parameters on vegetation patchwork, and provides theoretical basis for specific measures to combat land desertification.

Studying the effect of aeolian sand diffusion intensity on vegetation is the main objective of this paper, which compares the effect of two mechanisms on changes in vegetation pattern structure based on nonlocal delays. The rest of this paper is organized as follows. Section 2 derives a nonlocal delayed vegetation-sand model with soil-sand diffusion. Section 3 discusses the stability of the equilibria and the conditions for the emergence of the Turing pattern. We address a weakly nonlinear multiple scales analysis and obtain the amplitude equation for the Turing pattern in Section 4. In Section 5, numerical simulations are shown to verify the theoretical results.

2.   A mathematical model derivation

In semi-arid environments, vegetation and sand produce nonlocal interaction. In addition to the vegetation-covered areas, aeolian sand will have corresponding impacts on the vegetation throughout the study area. It will erode stems and leaves, damage plant tissues, and affect their productivity, among other things. To better describe the process, we create the following vegetation-sand system with nonlocal delay

For a detailed explanation of the kernel function ^[44,45]

The second and third of the above three kernels are referred to as weak and strong generic delay kernels respectively. The first kernel is the appropriate choice for giving a model with a discrete time delay, that is to say, the delay effect only involves the data exactly $\tau$ time units ago, $H(t)$ is the weak generic delay kernel for the Neumann problem and the expression is as follows

The kernel function $Q(x, y, t)H(t)$ represents the weight of eolian sand reaching the current position at any position before time $t.$ The expression of nonlocal delay is as follows

$Q(x, y, t)$ is the solution of the following system in a bounded domain

and

where $\sigma$ is the Dirac delta function and $D > 0$. According to (2.2) and the properties of $\sigma$ function, the above equation can be derived as follows

Based on the above derivation, (2.1) can be transformed into the following form

System (1.3) always has one bare sand state equilibrium point $K_0 = (1, 0).$ When $v\neq 0,$ we have

Let

The number of equilibrium points in system (1.3) depends on the relationship between above parameters.

(ⅰ) When the parameters satisfy the condition $B^2-4AC > 0, B-\sqrt{B^2-4AC} > 0,$ Eq (2.4) has two positive roots

Therefore, there are $s_{11} = 1+v_{11}, \; s_{12} = 1+v_{12}$. It can be obtained that system (1.3) has two uniformly vegetated equilibrium points $K_{11} = (s_{11}, v_{11}), \; K_{12} = (s_{12}, v_{12})$.

(ⅱ) When the parameters satisfy the condition $B^2-4AC = 0, B > 0,$ Eq (2.4) has a unique positive root, $v_2 = -\frac{B}{2A}$ and $s_2 = 1+v_2.$ Therefore, system (1.3) has an equilibrium point $K_2 = (s_2, v_2).$

(ⅲ) When the parameters satisfy the condition $B^2-4AC < 0,$ Eq (2.4) does not have any real roots, indicating that there is no equilibrium for system (1.3).

In order to discuss the stability of the above equilibrium point, system (1.3) is linearized near the equilibrium point $K^* = (s_*, v_*)$ to obtain the Jacobin matrix

Let $H = h-\dfrac{2hv_*}{v_m}-\dfrac{ps_*}{(1+cv_*)^2}, N = \dfrac{pv_*}{1+cv_*}.$ Its corresponding characteristic equation is

where

We can see that if $H < 1$ and $H < N$, then the equilibrium $(s_*, v_*)$ is locally asymptotically stable for (1.3). The following discussions are all assumed that $H < 1$ and $H < N$.

3.   Linear analysis and Turing patterns

3.1.   Stability and Turing bifurcation analysis of the diffusive system without nonlocal delay

In this subsection, we study (1.3) without considering the nonlocal delay. We start with the local system in one-dimensional space $\Omega = (0, \pi)$

Define the real Sobolev space

and let the complex extension space of $X$:$X_C = X\bigoplus iX = \{x_1+ix_2\mid x_1, x_2 \in X\}.$ The linearization operator for system (3.1) at $(s_*, v_*)$ is

Under the chi-square Neumann condition, the eigenvalue of operator $-{\nabla}^2$ is $\mu_k = k^2$ $(k = 0, 1, 2, \cdots)$, which satisfies

and ${\rm cos} (kx)\; (k\in N)$ corresponds to the characteristic function of $\mu_k.$ Taking the sequence of functions $\{{\rm cos}(kx)\}_{k = 0}^\infty$, it is a standard Orthonormal basis of space $L^2(0, \pi).$ Let

be the eigenfunction of $L$ corresponding to eigenvalue $\psi$, where $L(\phi, \varphi)^T = \psi(\phi, \varphi)^T.$ By direct calculation, we have

and

Let the characteristic equation of $L_k$ be

where

Obviously, when $H \le 0,$ for any $k\ge 0,$ we have $D_k\ge 0$ and $T_k < 0.$ At this point, the positive equilibrium point $(s_*, v_*)$ of system (3.1) is locally asymptotically stable.

In the following, assume that $0 < H < 1$ and $H < N$. If $d\ge H,$ then for any $k\ge 0,$ we obviously have $D_k\ge 0$ and $T_k < 0,$ which implies that $(s_*, v_*)$ is locally asymptotically stable. If $d < H,$ then let

Note that the discriminant of the quadratic function $f(z) = z^2+(2H-4N)z+H^2$ is

Thus, $f(z) = 0$ exists two positive real roots

If $z_1 < z < z_2,$ then $f(z) < 0$ and $\Delta < 0$, which implies $D_k > 0$ for any $k > 0$ when $z_1 < d < z_2$. Let $z_1-H = 2N-2H-2\sqrt{N(N-H)}.$ Since $H < N$, $N-H < \sqrt{N}\sqrt{N-H}.$ We can get $z_1 < H.$ Obviously, $Z_2 > H.$ Note that $z_1 < H < z_2$ and $(s_*, v_*)$ is locally asymptotically stable when $d\ge H$. Hence, $(s_*, v_*)$ is locally asymptotically stable when $d > z_1$.

Theorem 3.1. Assume that $H < 1$ and $H < N$ hold. Then, the equilibrium point $(s_*, v_*)$ is locally asymptotically stable for system (1.3) if $H\le0$. When $H > 0$, the equilibrium point $(s_*, v_*)$ is locally asymptotically stable for system (1.3) if $d > z_1$, where $z_1 = 2N-H-2\sqrt{N(N-H)} > 0$.

If $0 < d < z_1,$ then $\Delta > 0.$ The equation

has two positive real roots

Let $F(d) = H-d+\sqrt{(d-H)^2-4(N-H)d}$. Then,

Recall that $H > z_1$ and $d < z_1$. It follows from $d < H < N$ that $d < 2N-H$. Hence, we can get $F^{'}(d) < 0$ and $\mu_+(d)$ is monotonically decreasing on $d.$ Define

Let $d\to 0^+$. Then, we have

Clearly, we see that $\Phi_1\cap \Phi_2\not = \varnothing$, which implies that the positive equilibrium $(s_*, v_*)$ of system (3.1) is unstable. So, we obtain the Turing instability of $(s_*, v_*)$ for $d$ is small.

Theorem 3.2. Assume that $0 < H < 1$ and $H < N$ hold. Then, there exists sufficiently small $\tilde{d}$ such that for $0 < d < \tilde{d}$, the positive equilibrium $(s_*, v_*)$ is Turing unstable for system (3.1).

3.2.   Stability analysis for the local system with nonlocal delay

In this section, we will study the conditions under which the equilibrium point $(s_*, v_*, n_*)$ is stable when the system (2.3) has no diffusion.

Begin with the local system

where $t > 0$. The linearized system of (3.3) at equilibrium $(s_*, v_*, n_*)$ is

where

We analyzed the stability of the equilibrium point $(s_*, v_*, n_*)$ in Table 3-1 ^[46].

The coefficients in Table 3-1 are shown below

Condition 1. When $H < 1$ and $H < N$, $P_1(0) > 0$ and $P_3(0) > 0.$

Condition 2. When $H < 1$ and $H < N$, $P_1(0)P_2(0)-P_3(0) > 0.$

Let $(s_*, v_*, n_*) = (1+\beta, \beta, 1+\beta)$ and $X_1 = \dfrac{2h\beta}{v_m}-h$. Then, we have

When $X_1\in(-\infty, X_{1, -})\cup(X_{1, +}, +\infty)$, we can calculate $P_1(0)P_2(0)-P_3(0) > 0.$

Hence, the two conditions in Table 3-1 are satisfied simultaneously, so $(s_*, v_*, n_*)$ is asymptotically stable to (3.3).

3.3.   Stability and Turing bifurcation analysis of the diffusive system with nonlocal delay

In this subsection, we study (2.3) as spatial diffusion increases. We analyze the characteristic equation and obtain the conditions of Turing bifurcation generation using the Routh-Hurwitz criterion.

It can be seen from the above analysis that the equilibrium $(s_*, v_*, n_*)$ is locally asymptotically stable when the vegetation-sand system does not consider spatial diffusion. When the vegetation-sand system is coupled with spatial diffusion, it can be concluded from the Routh-Hurwitz criterion that the equilibrium $(s_*, v_*, n_*)$ is stable under the following conditions

Different nonlocal effects correspond with different real part curves of eigenvalues. It is easy to see that in the appropriate parameter range, as the intensity of nonlocal effect gradually increases, the real part of eigenvalue gradually increases and Turing patterns also appear. Obviously, near the equilibrium point $(s_*, v_*, n_*)$, the real part of the eigenvalue of the linearized system (2.3) can be regarded as a function of wave number $k$. We show the dispersion relation under different nonlocal intensities in Figure 1. The necessary condition for Turing pattern in system (2.3) is that the non-diffusion system is stable, while the diffusion system is unstable. Then, we will analyze the conditions of Turing pattern in equilibrium $(s_*, v_*, n_*)$ from the following three conditions.

Let

where $k$ is the wavenumber, $\lambda$ is the perturbation growth rate in time $t$, and $i^2 = -1.$ The exponential solution (3.4) is substituted into the (2.3) and the following characteristic equation is obtained through calculation

The solution of the characteristic (3.5) is found and then the dispersion relation is as follows

where

We show the dispersion relation for different nonlocal effects in Figure 1. A necessary condition for the Turing pattern in (2.3) is that the system is stable in the absence of diffusion and unstable in the presence of diffusion. Therefore, we will analyze the Turing pattern's equilibrium conditions $(s_*, v_*, n_*)$ in the following three conditions.

Condition 1. $P_1(k) > 0.$

Since $k^2+dk^2+Dk^2 > 0$ and the equilibrium $(s_*, v_*, n_*)$ is stable without diffusion, Condition 1 always holds.

Condition 2. $P_3(k) > 0.$

Theorem 3.3. In Condition 2, the equilibrium point of the system is stable in the absence of diffusion and unstable in the presence of diffusion, which is a necessary and sufficient Condition for the Turing pattern. If the inequalities $f_2^2-3f_1f_3 > 0, \; u_{\rm min} = u_1 > 0,$ and $F_{\rm min} = F(u_1) < 0$ hold, then the Turing bifurcation could occur.

Proof. Let $P_3(k) = F(k^2)$ and $u = k^2.$ Then, $J(k)$ is equivalent to the following expression

where

Analyze the properties of the polynomial $F(u)$ as follows.

(ⅰ) $\lim \limits_{u \to +\infty} F(u) = +\infty.$

(ⅱ) We can obtain two extremum points by calculating the first partial derivative of (3.2). The forms of the two extreme points are as follows

(ⅲ) We can obtain $F(u)$ second-order deflection by calculating

Using the properties of a unary cubic polynomial, we have

Thanks to $f_0 > 0,$ we obtain that if

is satisfied, Turing bifurcation could occur. Because $u_{\rm min} = u_1 = \dfrac{-f_2+\sqrt{f^2_2-3f_1f_3}}{3f_3}$ represents the number of the wave, $u_{\rm min} = u_1 > 0,$ to ensure that $u_1, u_2$ is meaningful, thus $f_2^2-3f_1f_3 > 0.$

In Condition 2, if the inequalities $f_2^2-3f_1f_3 > 0, u_{\rm min} = u_1 > 0,$ and $F_{\rm min} = F(u_1) < 0$ hold, then the Turing bifurcation could occur. The proof is completed. Λ

Condition 3. $P_1(k)P_2(k)-P_3(k) > 0.$

Theorem 3.4. If the inequalities $g_2^2-3g_1g_3 > 0, u_{\rm min} = u_3 > 0,$ and $G_{\rm min} = G(u_3) < 0$ hold, then the Turing bifurcation could occur.

Proof. Let $G(k^2) = P_1(k)P_2(k)-P_3(k)$ and $u = k^2.$ Then, $G(u) = g_3u^3+g_2u^2+g_1u+g_0,$ where $g_3 = 2dD+d+D+d^2D+d^2+dD^2+D^2 > 0,$

Next, analyze the properties of the polynomial $G(u)$

(ⅰ) $\lim \limits_{u \to +\infty} G(u) = +\infty.$

(ⅱ) We can obtain two extremum points by calculating the first partial derivative of $G(u).$ The forms of the two extreme points are as follows

(ⅲ) We can obtain $G(u)$ second-order deflection by calculating

Using the properties of a unary cubic polynomial, we get

Thanks to $g_0 > 0,$ we obtain that if the following conditions are satisfied, Turing bifurcation could occur

Because $u_{\rm min} = u_3 = \frac{-g_2+\sqrt{g^2_2-3g_1g_3}}{3g_3}$ represents the number of the wave, $u_{\rm min} = u_3 > 0,$ to ensure that $u_3, u_4$ is meaningful, thus $g_2^2-3g_1g_3 > 0.$

Similar to the Condition 2, if the inequalities $g_2^2-3g_1g_3 > 0, u_{\rm min} = u_3 > 0,$ and $G_{\rm min} = G(u_3) < 0$ hold, then the Turing bifurcation could occur. The proof is completed.

Λ

4.   Amplitude equation analysis of Turing patterns

In this section, $D$ is the control parameter, and $D_T$ is the bifurcation threshold for the Turing pattern. $D$ is close to the initial $D_T,$ and the eigenvalues are about zero, corresponding to the slow-changing critical mode. In contrast, the off-critical methods relax quickly, so only the perturbation of $k$ around $k_T$ should be considered. We first calculate the critical wave number $k_T$ and bring the necessary $k_T$ into $P_3(k) = 0$ to obtain the bifurcation threshold $D_T$ of the Turing pattern.

According to Theorem 3.3, we get

Therefore, $D_T$ satisfies the following equality

where

To obtain the amplitude equation, the perturbation solution $(s-s_*, v-v_*, n-n_*)$ of (2.3) can be represented by $X = (s, v, n)^T.$ Then the linearized form of (2.3) near the uniform steady state $(s_*, v_*, n_*)$ is written as follows

where

According to ^[47], for $D$ sufficiently close to $D_T,$ the solutions of (2.3) can be expanded in a hexagonal planform as

in which $A_j$ and $\overline{A}_j$ are the amplitudes corresponding to the modes of $\mathit{\boldsymbol{{\rm k}}}_j$ and $-\mathit{\boldsymbol{{\rm k}}}_j.$ Here $\mathit{\boldsymbol{{\rm k}}}_j$ is the wave vector with $|\mathit{\boldsymbol{{\rm k}}}_j| = k_T.$ Let $W = (W_1, W_2, W_3)^T.$ System (4.1) can be expressed as follows

where

Since we calculate a sufficiently small neighborhood of the critical value of $D,$ we extend the bifurcation parameter $D$ at $D_T$

By expanding $X$ and the non-linear terms $W$ by the small parameter $\epsilon,$ we can obtain

and

with

The linear operator $L$ can be decomposed as

where

Separating the time scale into $T_1$ and $T_2$ by the small $\epsilon$ through the chain rule of differentiation, we have

with

Since the amplitude $A_j$ changes slowly, derivative $T_0$ does not affect on amplitude $A_j.$ We consider each time scale $T_i$ as an independent variable, and then the reciprocal of time can be obtained as follows

Substituting (4.3)–(4.11) into (4.2), we have

$\epsilon:$

${\epsilon}^2:$

${\epsilon}^3:$

As the linear operator of a vegetation-sand system at a critical point, $(s_1, v_1, n_1)^T$ corresponds to the linear combination of the eigenvectors in which the eigenvalue is $0.$ By calculating (4.12), we can write the corresponding modes of the three wave vectors separately. We have

where

and $c.c.$ represents the complex conjugate. $V_j(j = 1, 2, 3)$ is the amplitude of ${\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_j\cdot\mathit{\boldsymbol{{\rm r}}})$ which is under first-order perturbation. The perturbation higher-order term determines it.

For (4.13), we get that

According to the Fredholm solvability condition, the vector function of the right-hand side of (4.17) must be orthogonal with the zero eigenvectors of the operator $L^+_T$ to ensure the existence of the nontrivial solution of the (4.17), where $L^+_T$ is the adjoint operator of $L_T.$ The zero eigenvector of operator $L^+_T$ is determined by

where

Customarily, $G^j_s, G^j_v, G^j_n(j = 1, 2, 3)$ are the coefficients corresponding to ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_j \cdot\mathit{\boldsymbol{{\rm r}}})$ in $G_s, G_v$ and $G_n$ respectively. Therefore, substituting (4.15) into (4.17) and sorting out the coefficients of ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_1 \cdot\mathit{\boldsymbol{{\rm r}}}),$ ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_2 \cdot\mathit{\boldsymbol{{\rm r}}})$, and ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_3 \cdot\mathit{\boldsymbol{{\rm r}}}).$ Using the Fredholm solvability condition, we obtain the following system

where

Applying orthogonal condition

we can obtain that

where

The system of (4.20) is an amplitude equation with first-order perturbation. We can find that the second-order coefficient of the equation is more significant than zero. Then, we can get amplitude $V_j (j = 1, 2, 3)$ diverges. So, introducing higher-order perturbation terms into the solution of (4.13) is necessary, such as ${\rm exp}(2i\mathit{\boldsymbol{{\rm k}}}_j \cdot\mathit{\boldsymbol{{\rm r}}})$, ${\rm exp}(i(\mathit{\boldsymbol{{\rm k}}}_1-\mathit{\boldsymbol{{\rm k}}}_2)\cdot\mathit{\boldsymbol{{\rm r}}})$, and so on.

with

where

For (4.14), we get that

Substituting the solutions (4.15) and (4.22) of the upper two order perturbation equation into (4.14) and sorting out the coefficients of ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_1 \cdot\mathit{\boldsymbol{{\rm r}}}),$ and ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_2 \cdot\mathit{\boldsymbol{{\rm r}}})$, and ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_3 \cdot\mathit{\boldsymbol{{\rm r}}})$ by using the Fredholm solvability condition, we obtain the following system

where

and

Applying orthogonal condition

we can obtain that

The amplitude $A_j$ is the coefficient of ${\rm{ }}{\rm exp}(i\mathit{\boldsymbol{{\rm k}}}_j \cdot\mathit{\boldsymbol{{\rm r}}}) (j = 1, 2, 3)$ at each level. We have

Together with (4.9), (4.10), (4.20), and (4.25), the amplitude equations about $A_j$ are as follows

with

where $l_1, l_2, l^+_2, l^+_3, G_1, F_1, H, G_2, F_2, G_3, {\rm and}\; F_3$ are defined as (4.16), (4.18), (4.19), (4.21), (4.23), and (4.24).

Next, we study the stability of the amplitude equation of the vegetation-sand model. We can construct different types of pattern structures in Turing space, and a stable Turing pattern structure is a steady-state solution for the corresponding (4.26). Each amplitude of (4.26) can be decomposed into a corresponding mode ${\gamma}_i = |A_i|$ and phase angle ${\theta}_i.$ By substituting $A_i = {\gamma}_i exp(i{\theta}_i)$ into (4.26) and separating the real part from the imaginary part, we can get

where $\theta = {\theta}_1+{\theta}_2+{\theta}_3.$ The following conclusions can be drawn for (4.27).

Theorem 4.1. For homogeneous reaction-diffusion systems, assume that $Z_1 > 0$ and $Z_2 > 0.$

(i) A uniform steady state solution $(\boldsymbol{O})$: ${\gamma}_1 = {\gamma}_2 = {\gamma}_3 = 0.$ If $\phi < {\phi}_2 = 0,$ then $(\boldsymbol{O})$ is stable. If $\phi > {\phi}_2 = 0,$ then $(\boldsymbol{O})$ is unstable.

(ii) A stripe pattern diagram $(\boldsymbol{S})$: ${\gamma}_1 = \sqrt{\frac{\phi}{Z_1}} \neq 0,$ ${\gamma}_2 = {\gamma}_3 = 0$ with $\phi > {\phi}_2 = 0.$ If $\phi > {\phi}_3 = \frac{{\eta}^2Z_1}{{(Z_2-Z_1)}^2},$ then $(\boldsymbol{S})$ is stable. If $\phi < {\phi}_3,$ then $(\boldsymbol{S})$ is unstable.

(iii) Two hexagonal layouts $(\boldsymbol{H}_\boldsymbol{0}, \boldsymbol{H}_\mathbf{\pi})$: ${\gamma}_1 = {\gamma}_2 = {\gamma}_3 = \frac{|\eta| \pm \sqrt{{\eta}^2+4(Z_1+2Z_2)\phi}}{2(Z_1+2Z_2)}$ with $\phi > {\phi}_1 = -\frac{{\eta}^2}{{4(Z_1+2Z_2)}^2}.$ If $\phi < {\phi}_4 = \frac{(2Z_1+Z_2){\eta}^2}{(Z_2-Z_1)},$ then the solution ${\gamma}^+ = \frac{|\eta| + \sqrt{{\eta}^2+4(Z_1+2Z_2)\phi}}{2(Z_1+2Z_2)}$ is stable and the other solution ${\gamma}^- = \frac{|\eta| - \sqrt{{\eta}^2+4(Z_1+2Z_2)\phi}}{2(Z_1+2Z_2)}$ is unstable.

(iv) A mixed structure solution $(\boldsymbol{MS})$: ${\gamma}_1 = \frac{|\eta|}{Z_2-Z_1}, {\gamma}_2 = {\gamma}_3 = \sqrt{\frac{\phi-Z_1{\gamma}^2_1}{Z_1+Z_2}}$, which exists under the condition $\phi > {\phi}_3$ and $Z_2 > Z_1,$ and the solution is always unstable.

It should be pointed out that when ${\tau}_0 > 0,$ this means that if $\eta > 0,$ then the stable hexagon pattern is $H_0(\theta = 0)$ and $H_{\pi}(\theta = \pi)$ is unstable; if $\eta < 0,$ then the stable hexagon pattern is $H_{\pi}(\theta = \pi)$ and $H_0(\theta = 0)$ is unstable.

Theorem 4.2. In the interval where Turing pattern is generated, other parameters are fixed and $D$ and $\tau$ are set as variable parameters.

(i) When $\phi_2 < \phi < \phi_3$, the system (2.3) appears spot pattern.

(ii) When $\phi_3 < \phi < \phi_4$, according to different initial conditions, the system (2.3) appears spot pattern or stripe pattern.

(iii) When $\phi > \phi_4$, the spot pattern in the system (2.3) is transformed to the stripe pattern.

5.   Numerical simulations

In this section, we give some numerical simulations to support the theoretical results obtained in the above section. Choose a region of size $30 \times 30$ whose boundary satisfies Neumann boundary conditions. We set the time region as [0, 96] and the time step as $\Delta t = 0.002$. The initial value is the random perturbation at the equilibrium point $(s_*, v_*, n_*)$. The simulation program runs until the main features of the pattern do not seem to change.

(1) We use numerical simulation to verify the above theoretical analysis. Selecting the values of different parameters $h, p, v_m, c, d$, and $D$, we can calculate the values of $\eta$, $Z_1$, $Z_2$, $\phi_1$, $\phi_2$, $\phi_3$, $\phi_4$, and $\phi$ according to the expression of the amplitude equation coefficients in Section 4. In order to observe the pattern structure of vegetation, we fix the other parameters and select three sets of parameters that differ only in the value of $\tau$. We continue to keep $d = 0.01, h = 2.8, vm = 4, p = 3.9, c = 3$, and $D = 3.$ When $\tau = 0.35,$ the spot pattern appears in Figure 2(a); As $\tau$ decreases, the vegetation pattern structure shows mixed pattern of stripe and spot pattern (see Figure 2(b)); Figure 2(c) transforms into stripe pattern when $\tau = 0.31.$ Choose a region of size $90 \times 90$ whose boundary satisfies Neumann boundary conditions. The corresponding results are shown in Fig 2. Specifically, when the first set of parameter values is taken, $\phi$ is between $\phi_2$ and $\phi_3$, and system (1.3) presents a spot pattern (see Figure 2(a)); when the second set of parameter values is taken, $\phi$ is between $\phi_3$ and $\phi_4$, the spot pattern loses its stability and the stripe pattern begins to appear, showing a mixed pattern (see Figure 2(b)); when the third set of parameter values is taken, $\phi$ is greater than $\phi_4$, and the spot pattern disappears to a striped pattern (see Figure 2(c)). All our simulations are carried out on the basis of Table 1 at $\tau$.

(2) Let $h = 24.56, p = 17, v_m = 10, c = 0.2, \tau = 10$. The equilibrium $(s_*, v_*, n_*)$ is locally asymptotically stable for (3.3). See Figure 3. Let $h = 26, p = 17, v_m = 10, c = 0.2, \tau = 10$. The equilibrium $(s_*, v_*, n_*)$ is unstable for (3.3). See Figure 4. Choose the initial conditions $(s_0, v_0, n_0)$ = (2.1, 2.6, 2.1).

(3) Let $h = 3.3333, p = 2.8245, v_m = 10, c = 3.1111, d = 100$, and $(s_0, v_0) = (8.0248+1.337{\rm cos}3x, 7.0248+1.337{\rm cos}3x).$ The conditions $0 < H < 1, H < N$, and $d > z_1$ hold true. By Theory 3.1, the equilibrium $(s_*, v_*)$ is stable for system (3.1). See Figure 5. Take $h = 3.3333, p = 2.1234, v_m = 10, c = 1.1111, d = 0.00061$, and $(s_0, v_0) = (5.1533+1.337{\rm cos}1.7x, 4.1533+1.337{\rm cos}1.7x).$ Then, $d \to 0^+$ hold. Hence, by Theory 3.2, the equilibrium $(s_*, v_*)$ is Turing unstable for system (3.1). See Figure 6. Let $h = 3.3333, p = 2.8245, v_m = 10, c = 1.1111, d = 1$, and $(s_0, v_0) = (3.1213+1.337{\rm cos}1.7x, 2.1213+1.337{\rm cos}1.7x).$ System (3.1) can induce spatially inhomogeneous Hopf bifuracation. See Figure 7.

(4) Let $d = 0.01, h = 2.8, v_m = 4, p = 3.9, c = 3, D = 1$. Figure 8 shows the succession of vegetation pattern when $\tau = 0.01$. It can be seen from the figure, as time goes on, the vegetation patterns evolves from an initially uniform distribution to gradually uneven clustering, forming gapped patterns, and eventually reaching stability where it no longer changes with time.

(5) Let $d = 0.01, h = 2.8, v_m = 4, p = 3.9, c = 3, D = 1$. Figure 9 shows the vegetation pattern in the final stable state with different delays. As $\tau$ gradually increases, we observed that the vegetation patterns formed gapped, clustered, and ring-like patterns, with vegetation transitioning from dense to gradually sparse. This reflects the nonlocal delayed impact of aeolian sand on vegetation. The greater the delay, the higher the degree of vegetation destruction, resulting in a sparser distribution.

(6) In Figure 10, we consider the effect of the destruction by sand burial $p$ on vegetation pattern structure. Let $d = 0.01, h = 2.8, v_m = 4, c = 3, D = 1, \tau = 0.01.$ When the parameter $p$ is small, the gap pattern appears (Figure 10(a)); as the parameter $p$ continues to increase, the spot pattern appears, showing the spot stripe mixed pattern (Figure 10(c)); when the parameter $p$ further increases, the vegetation patches form a ring-like pattern.

6.   Conclusions

Aeolian sand vegetation patterns have been observed and described for a long time, but there are few explorations on their formation mechanisms. To explore the ecological mechanism of its formation, to enrich the scenarios of vegetation competition for resources in different ecosystems, and to deepen the understanding of desertification and vegetation restoration process, this paper mainly investigates the spatial and temporal complexity of vegetation in arid and semi-arid regions from the aspects of vegetation system modelling and theoretical analysis. First, a vegetation model with nonlocal effects is constructed, the equilibrium state of the model is calculated, linear analysis is carried out at the equilibrium state, and the spatio-temporal dynamics of the model are investigated. Second, the theoretical basis of vegetation for combating desertification is proposed. Finally, the effect of the nonlocal delay on the vegetation sand model is shown through numerical simulation.

The effect of nonlocal delay on vegetation patterns is as follows: a certain range as $\tau$ increases, the density of vegetation decreases, which reflects that the ability of sand to destroy vegetation increases as $\tau$ increases, so the nonlocal delay can change the structure of vegetation pattern. In terms of theoretical analysis, Turing space was obtained using linearized analysis for the model, but the pattern structure and its stability could not be determined. Amplitude equations that can describe the dynamical behaviour of the system near the destabilisation point of the vegetation model were derived using nonlinear analysis, both to describe the different parameter spaces corresponding to different patterns and to determine whether the corresponding patterns are stable or not, which makes up for the shortcomings of the linear analysis method.

The improved vegetation-sand dynamics model in this paper is theoretical. The theoretical model is widely accepted in the field of vegetation pattern, which explains the main mechanism of the vegetation pattern formation process and promotes the research process of nonlocal delay of the vegetation-sand model, but it is still difficult to explore the more complex vegetation-sand model, such as the addition of the advection term and so on, to a certain extent. For studying the landing of sand particles and the discontinuity of vegetation, a discrete model may be more suitable. Discrete models can use discrete event systems to describe the state changes of sand particles and vegetation, including the movement of sand, landing, and the frequency of sand impact on vegetation. We will also focus on studying discrete models of vegetative-sand on networks ^[48,49].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The work is supported by the National Natural Science Foundation of China (61872227, 12061081

Conflict of interest

The authors declare there is no conflict of interest.