Analyzing vegetation pattern formation through a time-ordered fractional vegetation-sand model: A spatiotemporal dynamic approach

1.   Introduction

The vegetation-sand model is an essential theoretical framework that serves to elucidate the complex and often counterintuitive interactions between vegetation and sandy environments, particularly in the context of desertification and ecosystem degradation. As desertification continues to emerge as one of the most pressing environmental challenges of the 21st century, this model offers a critical lens through which scientists and policymakers can understand, predict, and mitigate the adverse effects of this global phenomenon ^[1,2,3]. To understand the dynamics of desertification and conditions under which it occurs, sophisticated models must be developed to simulate the interactions between vegetation, soil, and climatic factors ^[4,5,6,7]. Previously, mathematical models have been developed as tools to explore how vegetation growth, distribution, and decline are influenced by environmental variables such as soil moisture, water, wind, and precipitation ^{[6,7,8,9,10,11]}. These models also emphasize the critical role of vegetation in stabilizing sandy soils and mitigating desertification.

In arid and semi-arid regions, vegetation patterns are primarily influenced by water redistribution, which results from the varying infiltration rates of different vegetation types ^[6,7,8,9]. Studying the water plant model allows us to understand the mechanisms behind vegetation pattern formation and the impacts of various factors——including the vegetation growth rate, natural mortality rate, grazing, and precipitation——on vegetation distribution ^{[11,12,13,14,15]}. In sandy and windy areas, plant growth is significantly affected by the movement of sand particles, whereas seed dispersal is influenced by wind. Therefore, the effects of sand and wind in these environments must be considered ^[5]. The emergence of vegetation patterns in desert ecosystems has attracted significant interest from ecologists and researchers, as these patterns provide important insights into the dynamics of desert environments, plant adaptation, and ecological processes. Prior studies have focused on how wind influences the spatial distribution of vegetation, resulting in patterns such as stripes, spots, or gaps in ecosystems ^[1,2,3,4,5]. However, theoretical studies using partial differential equations had not been conducted until Zhang et al. ^[10] incorporated two variables——vegetation coverage rate and wind——into a model to explore the formation mechanisms of sand-vegetation patterns. Using this model, Zhang et al. identified several key processes on the interaction between vegetation and sand particles, including vegetation growth, sand particle deposition, vegetation diffusion, and horizontal sand transport. In addition, they presented vegetation-sand models for two specific species. The introduction of reaction-diffusion equations to study vegetation patterns in ecosystems, especially in desert environments, has led to significant advancements in theoretical ecology ^[16,17,18]. The research conducted by Zhang et al. has likely improved our understanding of how mathematical models characterize the dynamics of plant populations and their spatial distribution. In this model, the influence of the prevailing wind on the growth rates of the two variables was represented as an advection term, whereas the impacts of other winds were represented as diffusion terms. It is important to note that many regions lack significant dominant winds ^[18]. Unlike Zhang et al. ^[10,16,17], Li et al. ^[18] investigated the aggregation mechanism of linear action terms, which facilitated the reformulation of the vegetation-sand model as follows:

where $\frac{\partial S}{\partial T}$ and $\frac{\partial V}{\partial T}$ are the accumulation rate of sand and growth rate of vegetation, respectively, $K_0\; +\; MV\; -\; NS$ represents deposition by vegetation, $HV(1\; -\; \frac{V}{V_m})$ 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 vegetation dispersal in all directions, respectively. Following dimensionless treatment, the resulting model is expressed as

where

Research has also been conducted on the reaction-diffusion model of vegetation in sandy environments (1.2). Li et al.^[18] investigated the occurrence of diffusion instability and its impact on vegetation patterns due to Turing instability. Subsequently, the bifurcation theory was employed to extend the local partial bifurcation to the global bifurcation, deriving the conditions for determining the direction of bifurcation. A study of double eigenvalues was conducted using spatial decomposition techniques and the implicit function theorem. Zhang et al.^[17] proposed a vegetation-sand model that incorporates non-local interactions characterized by integral terms involving kernel functions. The instability of the Turing diagram was analyzed, leading to the derivation of its stability conditions. In addition, the amplitude equation at the critical value of the Turing bifurcation was derived using a multiscale method. In summary, although most prior studies on vegetation-sand models have focused on integer derivatives^{[10,16,17,18]}, there is growing interest in exploring more advanced mathematical frameworks, such as fractional derivatives, to better capture the complexities of ecological dynamics. The succession of arid ecosystems can span long periods, sometimes on the scale of hundreds of years, and succession processes in each region may vary with respect to climate, soil, and other regional factors. Owing to the locality of integer-order derivatives, there are some limitations in describing the succession. Fractional derivatives are more suitable for description than integer derivatives owing to their memorability and nonlocality^[19,20,21]. Significant progress has been made in research on time-fractional reaction-diffusion models, particularly in the fields of ecology, biology, and applied mathematics^{[22,23,24,25]}, as these models provide a more nuanced understanding of complex systems in which memory effects and nonlocal interactions play crucial roles. Similar to other ecological models, we aim to replace the integer-order derivative vegetation-sand model with a fractional-order derivative model to more accurately and comprehensively represent ecological dynamics. The fractional-order derivative model captures memory effects, simulates nonlocal interactions, and enhances predictive capabilities, thereby providing valuable insights into the complex relationship between vegetation and sand dynamics. This transformation is essential for advancing ecological research and enhancing management practices in response to environmental challenges. The pattern behavior of the fractional model is significantly more complex than that of the integer order ^[26,27,28]. There are several types of fractional time derivatives, including the Caputo, Riemann-Liouville, and Marchaud derivatives ^[20,29]. Among these, the Caputo derivative is widely used in fractional calculi. It is particularly useful in modeling processes that exhibit memory effects and nonlocal behavior, making it applicable in various scientific and engineering disciplines. In this paper, we briefly discuss the definition and basic properties of Caputo derivatives that were subsequently used in this study.

The remainder of this paper is organized as follows. Section 2 presents a fractional vegetation-sand model, along with the definitions and tools applicable to fractional derivatives. In Section 3, we perform a comprehensive analysis of the temporal behavior of the solution in a diffusion-free scenario. This analysis identified a Hopf bifurcation, which signifies a substantial alteration in the system's dynamic behavior. Section 4 discusses the effects of spatial diffusion and fractional derivatives, as well as analyzing the existence of a Turing bifurcation. In Section 5, several numerical simulations are presented to validate the mathematical results. The significant outcomes are summarized in the concluding section.

2.   Fractional model and preliminaries

Fractional derivatives can account for the history of a system state, as the current state of a system depends on present as well as past variables ^[20,28]. Fractional derivatives offer advantages over traditional integer-order derivatives in certain contexts owing to the ability to capture memory effects and nonlocal behaviors in systems ^[30,31]. We briefly discuss the definition and basic properties of the Caputo derivatives ^[29] that were used throughout this study. First, we define the Caputo derivative as follows.

Definition 2.1. ^[20,29,32] The Caputo derivative of order $\delta \in (0, 1)$ for functional $g: [0, \infty) \rightarrow \mathbb{R}$ can be written in the following form:

where $\Gamma$ is Euler's gamma function. Let us now define the following nonlinear fractional-order system:

where $\delta \in (0, 1)$; $J \in \mathbb{R}^{n \times n}$; $g \in C^1(\mathbb{R}^n, \mathbb{R}^n)$; and $Dg(0) = 0$.

The Caputo fractional-order derivative is a generalization of traditional derivatives that allows for non-integer order differentiation, making it particularly useful for modeling complex systems that exhibit memory and hereditary properties ^{[32,33,34,35]}. Many biological processes exhibit memory effects, where the current state depends not only on the present conditions but also on past states. The Caputo fractional-order derivative captures these memory effects, allowing for more accurate modeling of phenomena such as population dynamics, drug absorption, and ecological interactions. The Caputo fractional-order derivative can be applied to models of population growth and decline, where the rate of change may depend on the history of population sizes or environmental factors. This can lead to more nuanced predictions about species interactions and ecosystem stability.

The stability analysis of fractional-order differential equations requires an examination of the roots of characteristic equations derived from the system. The complexity of fractional derivatives necessitates careful consideration of stability criteria specific to fractional systems. Recognizing these differences is crucial for the effective modeling and control of systems characterized by fractional dynamics. In fractional-order derivative systems, the concept of stability differs from that in differential systems, as outlined by the following theorem.

Theorem 2.1. ^[20,28] Assuming that $E^*$ is the equilibrium point of a fractional differential system, the linear system (2.2) is derived by linearizing the system around this equilibrium point. The system (2.2) is locally stable if all eigenvalues $\lambda$ of $M$ satisfy $|arg(\lambda)| > \frac{\delta\pi}{2}$$(0 < \delta < 1)$; conversely, the system (2.2) is unstable if $|arg(\lambda)| < \frac{\delta\pi}{2}$ for some values of the eigenvalues $\lambda$.

Let us now define the Hopf bifurcation of fractional-order systems using the parameter $\delta$. Hopf bifurcation is an important dynamic behavior that typically occurs in nonlinear systems when the stability of the equilibrium points changes along with system parameters, leading to the emergence of periodic solutions. The study of Hopf bifurcations in fractional-order systems is particularly significant because the dynamic characteristics of these systems differ from those in traditional integer-order systems. We explore how to analyze this phenomenon using appropriate mathematical tools and methods, and clarify the conditions and features of Hopf bifurcation in fractional-order systems.

Notation 2.1. ^[20,30,31] Assuming that $E^*$ is an equilibrium point of system (2.2), when the parameter $\delta$ takes the value $\delta = \delta_h$, system (2.2) undergoes Hopf bifurcation around $E^*$, provided that the following conditions hold:

(i) The Jacobian matrix of system (2.2) associated with the point $E$, has a pair of complex conjugate eigenvalues $\lambda_{1, 2} = a \pm ib$; these become purely imaginary at $\delta = \delta_h$.

(ii) $\digamma_{1, 2}(\delta_h) = 0$, where $\digamma_{i}(\delta) = \frac{\delta\pi}{2}- \mathop{\min}\limits_{i = 1, 2}|\mathrm{arg} (\lambda_{i})|$.

(iii) $\frac{\partial\digamma_{1, 2}(\delta)}{\partial\delta}|_{\delta = \delta_h}\neq0.$

Notation 2.2. ^[24] Let $0 = \mu_1 < \mu_2 < \cdots < \mu_k < \rightarrow \infty$ be the eigenvalues of the Laplacian operator $-\Delta$ on $\Omega$ under no-flux boundary conditions, and $E(\mu_i)$ be the space of the eigenfunctions corresponding to $\mu_i$ in $C^1(\Omega)$. We define the following space decomposition:

(i) $X_{i j} = \left\{c \times \phi_{i j}: c \in \mathbb{R}^{2}\right\}$, where $\left\{\phi_{i j}\right\}$ is an orthonormal basis of $E\left(\mu_{i}\right)$ for $j = 1, \ldots, \; \mathit{\text{dim}} E(\mu_i)$.

(ii) $X = \left\{w \in\left[C^{1}(\bar{\Omega})\right]^{2}: \frac{\partial u}{\partial \eta} = \frac{\partial v}{\partial \eta} = 0\right. on \left.\partial \Omega\right\}$, hence $X = \bigoplus_{i = 1}^{\infty} \bigoplus_{j = 1}^{\operatorname{dim} E\left(\mu_{i}\right)} X_{i j}$, where $w = (u, v)^T$ and $\eta$ is the outward unit normal vector along the boundary of the domain.

The use of fractional-order time models in vegetation-sand systems enhances traditional ecological models by incorporating fractional derivatives related to time. These models are particularly effective at capturing memory effects and nonlocal interactions that can occur in ecological dynamics. The following sections outline the key components, formulas, and examples of the fractional-order time model specifically developed for vegetation-sand systems. In these models, time derivatives are substituted with fractional derivatives, typically expressed in Caputo form. The order of derivatives, commonly represented as $\delta$ under the condition $0 < \delta < 1$, indicates the level of memory within the system. To elucidate the relationship between vegetation and sand in arid ecosystems, the following fractional-order models were used:

where $\Delta$ denotes the Laplacian operator, $\Omega$ is a bounded domain in $\mathbb{R}^2$ with a smooth boundary $\partial\Omega$, and $\eta$ is the outward unit normal vector along $\partial\Omega$. The homogeneous Neumann boundary conditions indicate that there is no flux across the boundary.

System (2.3) always has one bare sand state equilibrium point $E_0 = (1, 0).$ When $\nu\neq0$, we have

The number of equilibrium points in the system (2.4) depends on the relationship between the above parameters ^[18].

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

Therefore, we have $s_{1} = 1+v_{1}$ and $s_{2} = 1+v_{2}$. It can be observed that the system (2.3) has two uniformly vegetated equilibrium points $E_{1} = (s_{1}, v_{1})$ and $E_{2} = (s_{2}, v_{2})$.

(ⅱ) When the parameters satisfy the condition $B^2- 4AC = 0, B > 0$, Eq (2.4) has a unique positive root, $v _{3} = -\frac B{2A}$ and $s_3 = 1+v_{3}$. Therefore, the system (2.3) has the equilibrium point $E_3 = (s_3, v_3).$

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

3.   Effect of time-fractional derivative

Before discussing the time-fractional spatiotemporal sand-vegetation model, it is essential to understand the corresponding time-fractional ordinary differential equation (FODE) model. The FODE serves as a foundational building block for the extension to more complex spatial models, such as reaction-diffusion equations. In this section, we investigate the behavior of the following temporal model using the time-fractional derivative

3.1.   Local stability analysis

We sequentially present the dynamic behavior of the model system around each steady state using the following theorems:

Theorem 3.1. The bare ground equilibrium $(1, 0)$ is locally asymptotically stable when $p > h$ but becomes unstable when $p\leq h$.

Proof. To analyze the stability of the equilibrium point $(1, 0)$, system (3.1) is linearized in its vicinity to derive the Jacobian matrix

The eigenvalues of the corresponding coefficient matrix are the roots of the characteristic equation expressed as

where the roots of this characteristic equation are

Two eigenvalues are $\lambda_1 = 2(h-p)$ and $\lambda_2 = -2$. If $p > h$ then $arg{\lambda_1} = arg{\lambda_2} = \pi$ which implies $arg{\lambda} > \frac{\delta\pi}{2}$. Hence the equilibrium point $(1, 0)$ is locally asymptotically stable when $p > h$. For $p\leq h$, we find $arg{\lambda_1} = 0$, $arg{\lambda_2} = \pi$, and consequently $arg{\lambda_2} < \frac{\delta\pi}{2}$ which implies $(1, 0)$ is unstable.

In order to discuss the stability of the internal equilibrium point, system (3.1) is linearized near the equilibrium point $E^*: = (s^*, v^*) = (s_j, v_j)(j = 1, 2, 3)$ to obtain the Jacobin matrix

where $H: = h-\frac{2hv^{*}}{v_{m}}-\frac{ps^{*}}{(1+cv^{*})^{2}}$, and $N = \frac{pv^{*}}{1+cv^{*}}$.

Its corresponding characteristic equation is

where

If $H < 1$ and $H < N$, the equilibrium point $E^* = (s^*, v^*)$ is locally asymptotically stable for system (3.1) when $\delta = 1$. We will now analyze the local stability of the equilibrium point $(s_*, v_*)$ for the case when $0 < \delta < 1$.

Theorem 3.2. The system (3.1) shows stable behavior around the interior equilibrium point $E^*$ if any one of the following conditions hold:

(i) $H\leq1$ and $H < N$ (see the cyan region in Fig. 1 and the second quadrant of the $H-N$ coordinates),

(ii) $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$(see the magenta region in Figure 1).

Proof. The characteristic equation of the Jacobian matrix at the interior equilibrium $E^*$ is expressed as

To check the local stability of the model system (3.1) around this interior steady state $E^*$, let us consider the following cases:

Case (Ⅰ): $H\leq1$ and $N > H$.

For this case, we have three sub-case scenarios as follows:

Sub-Case (Ⅰa): $H = 1$ and $N > 1$.

In this case, we can obtain both eigenvalues as complex conjugates (as $\lambda_{1, 2} = \pm\sqrt{N-1}i).$ Therefore,

which satisfies the properties of Theorem 2.1. Consequently, the model system demonstrates stable behavior around the interior steady state $E^*$.

Sub-Case (Ⅰb): $H < 1$ and $H+\frac{(H-1)^2}{4}\geq N > H$. After performing algebraic manipulation, it is observed that both eigenvalues $\lambda_1$ and $\lambda_2$ are negative real numbers

From Theorem 2.1, we know that the model system (3.1) exhibits stable behavior around this internal steady state $E^*$.

Sub-Case (Ⅰc): $H < 1$ and $H+\frac{(H-1)^2}{4} < N$. For this case, the Jacobian matrix $J$ will have a pair of conjugate complex eigenvalues with negative real parts. Therefore

So, by Theorem 2.1, we can tell that the model system (3.1) exhibits stable behavior around this interior steady state $E^*$.

Considering all the aforementioned sub-cases, we conclude that the model system (3.1) demonstrates stable behavior near the interior steady state $E^*$ when $H\leq1$ and $H < N$.

Now, let us consider the second case.

Case Ⅱ: $H > 1$, $N > \frac{(H-1)^2}{4}+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$.

For this case, we obtain a pair of complex conjugate eigenvalues satisfying

and

Therefore,

which satisfy the following restrictions (see the magenta region in Figure 1):

Therefore,

According to Theorem 2.1, the model system (3.1) demonstrates stable behavior in the vicinity of the interior steady state $E^*$ when $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$.

Remark 3.1. If $H < 1$ and $H < N$, the equilibrium point $E^* = (s^*, v^*)$ is locally asymptotically stable for both integer-order ordinary differential models (i.e., $\delta = 1$) and fractional ordinary differential models (i.e., $0 < \delta < 1$). However, if $H = 1$ or $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$, fractional differential models exhibit stability, while integer differential models show instability. This suggests that fractional derivatives with respect to time improve the model's stability. In summary, the incorporation of fractional derivatives with respect to time significantly enhances the stability of ecological models like the vegetation-sand model.

Theorem 3.3. The system (3.1) exhibits unstable behavior around the interior equilibrium point $E^*$ if any one of the following conditions hold:

(i) $N < H\leq1$ (see the white region $U_1$ in Figure 1),

(ii) $H > 1$ and $N\leq\frac{1}{4}(H-1)^2+H$ (see the white region $U_2$ in Figure 1),

(ii) $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta > \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$(see the white region $U_3$ in Figure 1).

Proof. We will examine each case individually.

Case (Ⅰ): $H = 1$ and $N\leq1$: For this case, we have three sub-case scenarios as follows:

Sub-Case (Ⅰa): When $H = 1$ and $N < 1$, the eigenvalues $\lambda_2$ become positive real numbers and hence

Sub-Case (Ⅰb): When $H = 1$ and $N = 1$, both eigenvalues become zero and hence

Based on Theorems 2.1 and 3.2, we can determine that the model system (3.1) exhibits unstable behavior near this internal steady state $E^*$.

Case (Ⅱ): $H > 1$ and $N\leq\frac{1}{4}(H-1)^2+H$: Both eigenvalues become positive real numbers and hence

For this situation, according to Theorem 2.1, we can determine that the model system (3.1) exhibits unstable behavior near this internal steady state $E^*$.

Case (Ⅲ): $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta > \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$.

In this particular scenario, we derive a set of eigenvalues that are complex conjugates of each other, indicating that

and

Therefore,

Based on Theorem 2.1, it can be inferred that the model system described by Eq (3.1) demonstrates unstable dynamics in the vicinity of the internal steady state $E^*$.

3.2.   Hopf-bifurcation analysis

In this section, we investigate the Hopf-bifurcation conditions analytically in a fractional-order system. Hopf bifurcation is a critical phenomenon in dynamical systems that occurs when a system's stability changes, leading to the emergence of periodic solutions (oscillations) as parameters vary. The coexistence equilibrium point for ordinary differential equations, that is the model (3.1) with $\delta = 1$, loses its stability through Hopf bifurcation when $H = 1$ and $N > H$. Considering $\delta$ as the bifurcation parameter and $\delta_h$ as the Hopf bifurcation threshold, the Hopf bifurcation condition $\delta = \delta_h$ applies to the coexistence equilibrium point $E^*$ of Eq (3.1), where $0 < \delta < 1$. Now we can state the Hopf-bifurcation condition for the coexistence equilibrium point $E^*$ for the model (3.1) following the approach reported in ^[26,28].

Let us assume $E^*$ is an unstable equilibrium point of the fractional-order system (3.1) such that the characteristic equation has two complex conjugate eigenvalues with positive real parts and hence the arguments of the eigenvalues are lying within $(0, \frac{\delta\pi}{2})$. Mathematically we can say that $E^*$ is unstable as

We can identify a threshold $\delta$, denoted as $\delta_h$, such that $\digamma(\delta) = 0$. We demonstrate that $E^*$ exhibits a Hopf bifurcation as $\delta$ crosses the threshold $\delta_h$.

Theorem 3.4. The system (3.1) undergoes Hopf bifurcation around the interior steady state $E^*$ at

where $H > 1$ and $N > \frac{1}{4}(H-1)^2+H$.

Proof. For $H > 1\; \; and \; \; N > \frac{1}{4}(H-1)^2+H$, the two roots of the characteristic equation (3.5) are complex conjugates with positive real parts. Hence

for some $\delta$. Now we can choose $\delta$ such that $\digamma(\delta) = 0$ and the root of this equation is denoted by $\delta_h$ where

and

Therefore, all the conditions for the occurrence of Hopf bifurcation are satisfied.

4.   Spatiotemporal behavior of the reaction-diffusion system

The spatiotemporal behavior of a time-fractional sand-vegetation system is an intriguing area of study that combines concepts from fractional calculus, reaction-diffusion dynamics, and ecological modeling. This system typically describes how vegetation interacts with the sand substrate over time and space, considering memory effects inherent in the dynamics. Depending on the parameters, this can result in various spatial patterns, including uniform distributions, clumping, or even oscillatory patterns as a result of feedback mechanisms between vegetation and sand. Analyzing the stability of equilibrium points and identifying bifurcations (such as Turing bifurcations and Hopf bifurcations) can reveal conditions under which periodic or spatially structured solutions emerge. For instance, as parameters change, the system may transition from stable equilibrium to oscillatory behavior or pattern formation.

We start with the spatiotemporal dynamics model in one-dimensional space $\Omega = (0, l\pi)$:

The linearization operator for system (4.1) at $(s_*, v_*)$ is

According to Notation 2.1, the eigenvalue of operator $-\frac{\partial^2 }{\partial x^2}$ is $\mu_k = \frac{k^2}{l^2}\left(k = 0, 1, 2, \cdots\right)$, which satisfies

and $\cos(\frac{kx}{l})\left(k\in N\right)$ corresponds to the characteristic function of $\mu_k$.

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

The roots of this characteristic equations are

It is well known that the positive homogeneous steady state $E^*$ is locally stable if the roots of the characteristic Eq (4.2) verify $|arg(\lambda_{1, 2})| > \frac{\delta\pi}{2}$.

Theorem 4.1. Assume that $H\leq1$ and $H < N$ hold. If $H\leq0$, then the equilibrium point $E^*$ is locally stable for system (4.1). When $H > 0$, the equilibrium point $E^*$ is locally stable for system (4.1) if $d > d_{*}$, where

Proof. Clearly, if $H\leq0$, then we have $Det_k\geq0$ and $Tr_k < 0$ for all non-negative positive integers. Therefore, both eigenvalues are negative real numbers, indicating that $|arg(\lambda_{1, 2})| > \frac{\delta\pi}{2}$.

Assume that $0 < H < 1$ and $H < N.$ If $d\geq H$, then for any $k\geq0$, we obviously have $Det_k > 0$ and $Tr_k < 0$, which implies that $E^*$ is locally asymptotically stable.

If $d < H$, then let

for all positive integers $k$. From the above inequality, we get $d^2+(2H-4N)d+H^2 < 0$.

Note that the discriminant of the quadratic function

Thus, $\theta(d) = 0$ and there exists two positive real roots

If $d_* < d < d^*$, then $\min Det_k > 0$ for any $k > 0$. Let $d_*-H = 2N-2H-2\sqrt{N(N-H)}$. Since $H < N, N-H < \sqrt{N}\sqrt{N-H}$. We can get $d_* < H$. Obviously, $d^* > H$. Note that $d_* < H < d^*$ and $E^*$ is locally asymptotically stable when $d\geq H$. Hence, $E^*$ is locally asymptotically stable when $d > d_*.$

Theorem 4.2. Assume that $0 < H\leq1$ and $H < N$ hold. Then, there exists sufficiently small $\tilde{d}$ such that for

the positive equilibrium $E^*$ is Turing unstable for system (4.1).

Proof. Suppose that $0 < H\leq1$ and $H < N$ hold. Then

If system (4.1) produces Turing instability, then we must have

for some positive integer $k$. If $Det_k < 0$, then the characteristic equation ($4.2$) has two real eigenvalues and one of them is greater than zero. So, then $|arg(\lambda_{1, 2})| < \frac{\delta\pi}{2}$. At this point, the positive equilibrium point $E^*$ of system (4.1) is unstable.

The following derivation applies to the interval of $d$ for which $Det_k < 0$. If $Det_k < 0$, then we have

which are equivalent to

By using the above inequalities (4.7) and Theorem 4.1, a necessary condition for the instability of $E^*$ is

It can be easily seen that when $0 < d < d_*$, the equation

has two roots

Thus, in order to get the instability of $E^*$, we must have ${(k^-)}^2 < {k}^2 < {(k^+)}^2$ for some positive integer $k$. If there exists some $k$ between $k^-$ and $k^+$, then it must satisfy $k^+-k^-\geq1$, which results in

According to the inequality above, we obtain $(\frac{1}{l^2}+1)d+2\sqrt{N-H}\sqrt{d}-H\leq0$. Here, we define a function

It is easy to check that $\Upsilon(\sqrt d) = 0$ has a unique positive root

and $\Upsilon(\sqrt d)\leq0$ when $0 < d \leq \tilde{d}$.

Remark 4.1. In Theorems 4.1 and 4.2, we examined the influence of the diffusion coefficients $d(: = \frac{D_2}{D_1})$, and identified the region of diffusion-induced instability. Specifically, for pattern formation to occur, the sand must diffuse significantly faster than the vegetation, i.e., $D_2 < < D_1$.

Theorem 4.3. Assume that $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$ hold. Then, the system (4.1) shows stable behaviour around the interior equilibrium point $E^*$ if any one of the following conditions hold:

(i) $d > \max\{d_*, \hat{d}\}$, where $\hat{d}: = l^2(H-1)-1$,

(ii) $d_* < d < \hat{d}$ and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$ for all positive integer $k$.

Proof. Assume that $H > 1$ holds. The following derivation applies to the interval of $d$ for which $Tr_k < 0$ for all $k$. If $Tr_k < 0$, then $Tr_k = 0$ has two opposable roots $k_1 = -\sqrt{\frac{l^2(H-1)}{1+d}}$ and $k_2 = \sqrt{\frac{l^2(H-1)}{1+d}}$. Because $k$ is a positive integer, there is no $k$ that satisfies $Tr_k > 0$ when $k_2 < 1$, whereas some $k$ satisfies $Tr_k > 0$ when $k_2\geq1$. According to the above analysis, we can get

Assume that $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$ hold. The following derivation applies to the interval of $d$ for which $Det_k < 0$ for all $k$. According to Theorems 4.1 and 4.2, we get

According to Eqs (4.12) and (4.13), it is easy to get that the equilibrium point $E^*$ is stable when $d > \max\{d_*, \hat{d}\}$.

If $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$, and $d_* < d < \hat{d}$ hold, then we obtain a pair of complex conjugate eigenvalues satisfying

and $Re (\lambda _{1}) = Re (\lambda _{2}) = Tr_k > 0$, which satisfy the following restrictions (see Figure 1a):

Therefore,

Based on the above discussion and Theorem 3.1(ii), the following theorem provides the obtained results.

Theorem 4.4. Assume that $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4(N-H)-(H-1)^{2}}}{H-1})$ hold. Then, the system (4.1) shows the following types of instability behavior around the interior equilibrium point $E^*$:

(i) If $\tilde{d} > d > \hat{d}$, then the system (4.1) shows Turing instability.

(ii) If $d < min\{\hat{d}, \tilde{d}\}$, then the system (4.1) shows Turing-Hopf instability.

(iii) If $d^* < d\leq\hat{d}$ and $\delta > \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$ for some positive integer $k$, then the system (4.1) shows Hopf instability.

Theorem 4.5. Assume that $H > 1$, $N > \frac{1}{4}(H-1)^2+H$, and $d^* < d\leq\hat{d}$ hold. Then the system (4.1) undergoes Hopf bifurcation at $\delta_h = \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$ for some positive integer $k\geq0$.

5.   Numerical simulation

In this section, we will conduct a numerical investigation into the impact of time-fractional derivatives on both the temporal and spatiotemporal behavior of the solution.

5.1.   Numerical simulation of a time-fractional ordinary differential model

Selecting the values of different parameters $h$, $p$, $v_m$, and $c$, we can calculate the values of $H$, $N$, $v_*$, and $s_*$ according to the expression of the equilibrium solution in Section 2.

(ⅰ) The parameters are set at $h = 2.6$, $p = 1.2$, $v_m = 26$, $c = 0.2$, and $\delta = 0.8$. By substituting the parameter values, we can easily obtain $v^* = 1.7191$, $s^* = 2.7191$, $H = 0.4493$, and $N = 1.5351$. According to the first case of Theorem 2.2, if $H\leq1$ and $H < N$, then the equilibrium $(s^*, v^*) = \; (2.7191, 1.7191)$ of system (3.1) is locally asymptotically stable (see Figure 2a). The observation that the numerical solution converges to the equilibrium point over time, despite a small perturbation in the initial conditions, aligns with the conclusion of Theorem 2.2.

(ⅱ) The parameters are set at $h = 2.6$, $p = 1.2$, $v_m = 30$, and $c = 0.45$. By substituting the parameter values, we can easily obtain $v^* = 4.6796$, $s^* = 5.6796$, $H = 1.0823$, and $N = 1.8081$. When $\delta = 1$, system (3.1) transforms into an integer-order model, and the equilibrium point is stable. However, when $\delta = 1$, based on the second case of Theorem 2.2, if $H > 1$ and $N > 1.0840$, then the equilibrium $(s^*, v^*) = (5.6796, 4.6796)$ of system (3.1) is locally asymptotically stable (see Figure 2b). Similar to the previous analysis, the observation that the numerical solution converges to the equilibrium point over time aligns with the conclusion of Theorem 2.2.

(ⅲ) In Figure 3a, the values $h = 26$, $p = 17$, $v_m = 10$, $c = 0.55$, $\delta = 0.8$, and the initial conditions $v_1 = 0.6115$ and $s_1 = 1.6115$ have been considered. In this case, the conditions $H > 1$, and $N\leq\frac{1}{4}(H-1)^2+H$ are satisfied, thus the interior equilibrium is unstable. The observation that a minor disturbance in the initial conditions leads to a shift in the equilibrium point of the numerical solution, as time progresses, is an important aspect of dynamic systems and aligns with the conclusions drawn from Theorem 2.2. In Figure 3b, the values $h = 26$, $p = 17$, $v_m = 10$, $c = 0.22$, $\delta = 0.75215$, and the initial conditions $v_1 = 0.6320$ and $s_1 = 1.6320$ have been considered. In this case, the conditions $H > 1$, and $N < \frac{1}{4}(H-1)^2+H$ are satisfied, thus the interior equilibrium is unstable, i.e., system (3.1) can induce Hopf bifuracation. Figure 3b illustrates the bifurcation branch of Hopf at the equilibrium point, demonstrating the emergence of periodic solutions in the model through numerical simulations, consistent with the theoretical results of Theorem 3.4.

5.2.   Numerical simulation of a time-fractional reaction-diffusion model

We begin by delineating the regions of stability and instability that align with Theorems 4.1 and 4.2. We present the stable critical diffusion coefficient as a function of the $N-d$ axis (see Figure 4) when $H = 0.4493$ and $N > 0.4493$. To be more precise, pattern formation occurs in region $A_3$ $($below the red curve $d = \tilde{d}(N)$$)$ and the system is stable in region $A_1$ $($above the black curve $d = d_*(N)$$)$. In the region $A_2 = \{(N, d): d_*(N) > d > \tilde{d}(N)\}$, it is unable to determine the stability of $E^*$.

To illustrate the results given by Theorems 4.1 and 4.2, we choose parameters $h = 2.6$, $p = 1.2$, $v_m = 26$, $c = 0.2$, $L = 10$, and $\delta = 0.8$. By incorporating the parameter values, we can easily obtain $v^* = 1.7191$, $s^* = 2.7191$, $H = 0.4493$, and $N = 1.5351$. By Theorems 4.1 and 4.2, we have $\tilde{d} = 0.0388$ and $d_* = 0.0905$. In fact, to obtain pattern occurrence, the sand must diffuse much faster than the vegetation (i.e., $D_2 \ll D_1$ in system (1.2)). According to Theorem 4.2(ⅰ), we can obtain bounded region $(0, 0.0388](\ni d)$ (see Figure 4) and Turing instability occurs in this region (see Figure 6a, b). By Theorem 4.2, there also exists an unbounded region $(0.0905, \infty)(\ni d)$ in which the equilibrium $(s^*, b^*)$ is stable (see Figure 5a, b).

To illustrate the results given by Theorems 4.3–4.5, we choose parameters $H = 2$ and $N > \frac{9}{4}$:

(ⅰ) We present the stable critical diffusion coefficient as a function of the $N-d$ axis when $L = 1.03$ (see Figure 7a). To be more precise, Turing instability occurs in region $D_3$ $($$\tilde{d} < d < \tilde{d}(N)$$)$, stability occurs in region $D_1$ $($above the black curve $d = d_*(N)$$)$, and Turing-Hopf instability occurs in region $D_4$ $($below the magenta curve $d = \hat{d}(N)$$)$. In the region $D_2 = \{(N, d): d_*(N) > d > \tilde{d}(N)\}$, it is unable to determine the stability of $E^*$.

(ⅱ) We present the stable critical diffusion coefficient as a function of the $N-d$ axis when $L = 1.2$ (see Figure 7b). Turing instability occurs in region $U_3$ $($$\hat{d}(N) < d < \tilde{d}(N)$$)$, Turing-Hopf instability occurs in region $U_6$$(d < min\{\hat{d}(N), \tilde{d}(N)\})$, the stable regions are $U_1$ $(d > \max\{d_*(N), \hat{d}(N)\})$ and $U_5$ $(d_*(N) < d < \hat{d}(N))$ when $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$, Turing-Hopf instability occurs in region $U_6$, Hopf instability occurs in region $U_5$ $(d < min\{\hat{d}(N), \tilde{d}(N)\})$ when $\delta > \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$, and in the regions $U_2$ and $U_4$, it is unable to determine the stability of $E^*$.

(ⅲ) We present the stable critical diffusion coefficient as a function of the $N-d$ axis when $L = 1.5$ (see Figure 7c). In this case, because $\hat{d}(N) > \tilde{d}(N)$, there is no Turing instability. The stable regions are $S_1$ and $S_2$ when $\delta < \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$, Turing-Hopf instability occurs in region $S_4$, Hopf instability occurs in region $S_2$ when $\delta > \frac{2}{\pi}\arctan(\frac{\sqrt{4Det_k-Tr_k^2}}{Tr_k})$, and in the region $S_3$, it is unable to determine the stability of $E^*$.

6.   Conclusions

In this paper, we introduce a fractional derivative to vegetation desert modeling with the objective of comprehensively describing and understanding the complex interactions between vegetation and sand, as well as their dynamic characteristics. Our study shows that the fractional characteristics of time series significantly affect the stability and evolution processes of the system. First, the analysis results indicate a nonlinear relationship between vegetation coverage and sand morphology. This relationship is not only influenced by environmental factors, but is also closely related to historical states. The incorporation of fractional derivatives with respect to time enhances the stability of vegetation-sand models by accounting for memory effects, and improving resilience to disturbances. These characteristics make fractional derivative models powerful tools for understanding complex ecological dynamics and for developing effective management strategies in the face of environmental challenges. Through numerical simulations, we observed differences in system behavior under different initial conditions, further validating the effectiveness of the fractional-order model in capturing complex dynamics. Furthermore, the research findings demonstrate that the fractional-order model can predict the feedback mechanisms between vegetation and sand more accurately. The figure illustrating the convergence of the numerical solution to the equilibrium point, even with small perturbations in initial conditions, underscores the stability of the model. This behavior aligns well with theoretical research, validating the underlying mathematical framework and enhancing the model's applicability in ecological management and conservation efforts. Such findings contribute to a deeper understanding of ecosystem dynamics and inform strategies for maintaining ecological balance in the face of change. This study enriches the theoretical foundation of spatiotemporal dynamic analyses and offers new ideas for future ecological models. The Turing instability analysis revealing that the diffusion rate of sand must be significantly faster than that of vegetation for vegetation patterns to emerge is a crucial insight. It underscores the importance of spatial dynamics in ecological interactions and provides valuable information for managing and restoring ecosystems in arid and semi-arid regions. By understanding these relationships, we can better address challenges related to desertification and promote sustainable land management practices. We hope that study will inspire scholars in related fields and promote the further exploration of ecosystem complexity. In summary, the fractional-order vegetation-sand model exhibits significant theoretical value and application potential in describing the spatiotemporal dynamics of ecosystems.

Future research should explore the impact of additional environmental factors on system dynamics to improve the model's applicability and predictive capabilities. In conclusion, further investigation of the vegetation-desert model incorporating fractional diffusion offers a valuable opportunity to deepen our understanding of vegetation dynamics and desertification processes. By incorporating memory effects and superdiffusion, along with conducting comprehensive stability analyses, researchers can obtain valuable insights to guide ecosystem management and conservation efforts. This research direction has substantial potential to address urgent environmental challenges associated with desertification and land degradation.

Author Contributions

Conceptualization, Y. Maimaiti; writing—original draft preparation, Y. Maimaiti, Z. Y. Lv, W. Zhang, and A. Muhammadhaji; writing—review and editing, Y. Maimaiti, W. Zhang, and A. Muhammadhaji; supervision, Y. Maimaiti; methodology, Y. Maimaiti. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

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

Acknowledgments

The authors are grateful to the anonymous referees, who carefully read the manuscript and made valuable comments and suggestions. The work is supported by the National Natural Science Foundation of China (Grant no: 12301639), Natural Science Foundation of Xinjiang Uygur Autonomous Region (Grant no: 2023D01C166), Talent Project of Tianchi Doctoral Program in Xinjiang Uygur Autonomous Region (Grant no: 5105240152k), and Open Project of the Key Laboratory of Applied Mathematics of Xinjiang Uygur Autonomous Region (Grant no: 2023D04045).

Conflict of interest

The authors declare there is no conflict of interest.