Patterns in a semiarid ecosystem are important because they directly and indirectly affect ecological processes, biodiversity, and ecosystem resilience. Understanding the causes and effects of these patterns is critical for long-term land surface management and conservation efforts in semiarid regions, which are especially sensitive to climate change and human-caused disturbances. It is known that there is a regular connection between the vegetation and the living species in a habitat since some animals evolved to live in a semiarid ecosystem and rely on plants for food. In this work, we have constructed a coupled mathematical model to connect the water resource, vegetation and living organisms and have investigated how the soil surface affects the resulting patterns for the long term. This study contributes to a better understanding of ecological patterns and processes in semiarid environments by shedding light on the complex interaction mechanisms that depend on the structure of semiarid ecosystems. The findings provide further critical insight into the influence of efforts for improving ecosystem resilience and adjusting to the challenges posed by climate change and human activities. Vegetation is an important component of an ecosystem because it helps energy and nutrient cycles, improves air quality, reduces flooding and soil erosion, and many other mechanisms [1,2]. Vegetation survives only on sunlight and groundwater minerals from the soil. In general, the availability of sunlight is not a problem for the vegetation, but sometimes the groundwater levels cause problems in semiarid ecosystems due to irregular rainfall [3]. This causes further food limitations for vegetation-dependent herbivores because they maintain a healthy balance by preventing vegetation's overgrowth [4]. It is known that diversity can be preserved through life-history trade-offs between growth rate and competitive ability [5]. This is the prediction of the intermediate disturbance hypothesis, which was first developed as a spatially homogeneous theory. Unlike theories that presume 'well-mixed' disturbance, patch dynamics predicts species diversity in relation to the intensity and magnitude of newly generated gaps [6]. In addition, many plants rely on herbivores, such as bees, to help them in pollination, and they enhance plant biodiversity at high productivity but have the reverse impact at low productivity [7]. As a result, low rainfall has a significant effect on regulating several food chains and food webs. Researchers have been using mathematical models as tools to capture and predict their behaviours in the long term. For instance, pattern formation in ecology helps in understanding the evolution of various living organisms in conjunction with non-living components. It is one of the key prerequisites that allows differently adapted species to coexist on a regional scale [8]. Many ecological theories have been developed to understand complex ecological systems and their interactions based on heterogeneous stationary and dynamic patterns [9]. Moreover, some important developments have been taking place in the context of complex systems that directly impact complex ecological systems studies [10,11]. Among others, research has also been progressing on vegetation patterns under climate change and water cycle pathways and their related factors, such as land surface coupled to other components of climate, through coupled climate models reflecting these systems' complexity [12][13][14][15]. Rich dynamics and associated patterns of different components of climate systems, including water and soil, have also been analyzed with a wide range of simplified mathematical models [16][17][18][19][20]. In general, the state of an ecosystem responds smoothly to continuous trends to a gradual change in external conditions, such as climate, groundwater reduction, habitat fragmentation, and many more [21][22][23]. But sometimes, a crucial situation arises when a drastic change happens in an ecosystem where the reversal is almost impossible. Such situations are often studied in the context of tipping phenomena and points of no return in ecosystems (e.g., [24] and references therein). There may be some indications or even no warning at all behind such drastic change [23,25]. This sudden dramatic shift of nature is called a catastrophic shift, and the critical threshold is known as 'catastrophic bifurcation' [26,27]. Some examples of this type of drastic change include the invasion of new species, a variety of human interventions, lightning-caused wildfires, etc. [28,29]. Heavy rain or extreme drought can also contribute to such climate events [15]. Therefore, preventing catastrophic change before its occurrence is essential for long-term ecosystem functioning. The major challenges to understanding such dramatic change are inefficient data and the infancy of mathematical theory. A number of indicators have been proposed prior to this shift [28], but to the best of our knowledge, all the underlying mathematical theory has been discussed only through saddle-node bifurcation [30]; however, this article is intended to explain a new underlying mechanism (transcritical bifurcation) for the regime shift in the context of vegetation pattern. Several mathematical models show long transient behaviours before converging to the final attractor [30], and they do not produce any clear trend in their properties. These long transients generally occur in a system close to different local and global bifurcation thresholds [27]. Sometimes, they do not converge to long-term stable solutions like constant or cyclic dynamics but persist for as many as tens of generations. In addition, a cyclic synchronous ecological system may lead to extinctions in the long run due to unexpected population fluctuations. For instance, a disturbance in two predators and one prey synchronously fluctuating population model may experience an overall predator populations boom, leading to a decrease in the prey populations, resulting in a sudden drop in the predator populations, and acting as a ripple effect throughout the ecosystem. On the other hand, long transients in ecological systems maintain biodiversity, making the ecosystem resilient and sustaining itself in the longer run [31]. Therefore, we intended to study long transients in the context of vegetation patterns. The soil surface of the earth is not uniform for all locations, and vegetation depends on it. Generally, the soil water travels through the ground in two ways: flow over the soil surface and underground diffusion [32], characterised by the soil surface.
In this work, we have added water diffusion into the existing model of water-plant dynamics in a semiarid ecosystem [33] and studied the role of the soil surface in a semiarid ecosystem. In addition, various animals and insects adapted to live in semiarid ecosystems; they depend on plants as their food source. Therefore, it is interesting to see how introducing herbivore species in the modified water-plant model affects the overall system's dynamics and enables us to understand the ecosystem better in the context of sustainability. As a baseline, we first consider the scenario where herbivore species are absent and study water-plant dynamics with and without water flow and its diffusion. Then, we discuss the effect of water flow and water diffusion in the presence of herbivores. As discussed, the temporal model without herbivores has at most three equilibrium points: water-only equilibrium (u 0 , v 0 ), which exists and is linearly stable for all the parameter conditions, and two coexisting equilibrium points can originate through a saddle-node bifurcation, and they are (u 1 , v 1 ) and (u 2 , v 2 ). These two coexisting equilibrium points exist for  > 2, and it is true for the parameter values given in Vegetation growth and spatial distribution are influenced by the climate, the physical environment, and human activity [15]. Many mathematical models predict such patterns in the presence of species mobility. For instance, a linearly stable coexisting equilibrium point for a temporal model may be unstable under a heterogeneous perturbation in the presence of diffusion terms [34]. In [33], the author has considered the water flow along the x-direction and the diffusion present only in the plant equation. In this case, the Turing instability analysis shows that the stable coexisting steady-state remains stable under heterogeneous perturbation when the water flow (a x ) is absent [34], and also for the flow up to a certain level (suppose a c x ), which depends on the values of  and . Therefore, the regular pattern formation is impossible on flat ground (a x = 0) and for the amount of water flows less than a c x . But, above the critical threshold of water flow a c x , the coexisting steady-state becomes unstable, and the system produces irregular transient solutions initially due to the heterogeneous perturbation of the initial conditions, and vegetation takes time to spread over the landscape with moving along the downhill. At a large time, the irregular transient solution changes into regular waves, and the vegetation patches keep the same distance between them. However, in the presence of water diffusion, the system produces Turing patterns for some diffusion coefficient irrespective of the water flow, as discussed in the earlier section. Furthermore, sufficient water flow induces stationary patterns to form strips propagating through the spatial domain. Diversity is an essential component of a system's persistence and adaptability in the face of environmental change [5], and here we discuss it by introducing the herbivores' interference into the existing water-plant model. After introducing the herbivores' dynamics, the resulting three-species model preserves all its equilibrium points with zero herbivores as the extra component. It has the water-only steady-state (u 0 , v 0 , w 0 ) = (, 0, 0), which is linearly stable. The steady-states (u 1 , v 1 , 0) and (u 2 , v 2 , 0) exist for  > 2, and (u 2 , v 2 , 0) is always unstable. The stability of (u 1 , v 1 , 0) depends on the feasible existence of the coexisting equilibrium point (u * , v * , w * ). This coexisting equilibrium point bifurcates through a transcritical bifurcation at (u 1 , v 1 , 0) [see the blue surface in Fig. 1(a)], and the coexisting equilibrium point gains stability and remains stable between blue and magenta surfaces. It loses stability through a sub-critical Hopf bifurcation (magenta surface in Fig. 1(a)), and the wateronly state remains the only linearly stable state. As the system is bistable, a trajectory can go to either of these equilibrium points depending on the initial conditions belonging to the appropriate basins of attraction. We first study the spatio-temporal dynamics of the model when water diffusion is present. We fix the temporal parameter values  = 2.8125,  = 0.8, and  = 0.05 and take  as the bifurcation parameter. The system produces a heterogeneous distribution for three components for  = 0.1 for different values of d u and d w . Figure 2 depicts the schematic diagram of patterns and the homogeneous solutions for the plants (v) for different combinations of d u and d w , which matches with theoretical investigations. Here, the silico trials show that the hot-spot pattern for the water corresponds to the cold-spot patterns for the plants and herbivores. In reality, the water negatively correlates with plants and herbivores because excess water damages the vegetation. The groundwater diffusion is generally low as it moves through the porous spaces between particles of unconsolidated soil or through networks of fractures and solution openings in consolidated rocks, which varies with locations [35]. Here, cold-spot patterns exist for the lower water diffusion on flat ground, and labyrinthine and hot-spot patterns form for higher diffusion rates [see Fig. 2]. In addition, negative correlations are also observed between the organisms. For zero groundwater diffusion, the non-trivial homogeneous steady state is linearly stable for the temporal parameter values used in Fig. 2, and no stationary pattern can be obtained [see Fig. 2]. The coexisting steady-state satisfies the Turing instability conditions when the water diffusion parameter crosses the critical Turing threshold [34], resulting in Turing patterns, e.g., see Fig. 3(a). Furthermore, this critical diffusion threshold can be obtained for the water by considering zero diffusion for the herbivores, as observed in The same type of moving strips for the water-plant model is observed for the herbivore-included model without groundwater diffusion.
Here, the groundwater flow makes the homogeneous steady-state unstable but requires a minimum speed. Though this critical speed depends on the other parameter values, the speed variation is slight compared to other parameter changes. In addition, these regular stripes can move diagonally for considering the same amount of water flow along x and y directions [e.g., see Fig. 4]. The ratio of water flow in the x and y directions determines the angle of moving regular strips. The non-homogeneous stationary patterns in Fig. 3(a) for the system in the presence of water diffusion no longer stay in one place in the presence of downhill water flow. They start moving from one side to another side [see Figs. 3(b)-(d)], even diagonally, when the flow is in both x and y directions [similar to Fig. 4], the Turing stability conditions are satisfied for each case. Sometimes, the water flow is insufficient to deform the Turing patch structure (e.g., the symmetric structure of the spots or labyrinthine); however, a good amount of water flow can force them to align in forming regular strips, and they move from one end to another. The moving speed and the distance between the two stripes depend on the water flow rate. Patterns may be stationary on a short-time scale, but they can be nonstationary over a longer time frame due to different pattern-forming mechanisms depending on the intensity of interaction and diffusivity. Sometimes, they oscillate regularly within a confined range or fluctuate irregularly, e.g., a regular moving stripe solution Fig. 3(d In contrast, the temporal model predicts the existence of the water-only steady state. The diffusivity for water and herbivores has to be small for chaotic solutions so that they do not spread out as quickly, leading to stronger local interactions. Furthermore, these solutions are observed in a narrow range of diffusion parameter values, which cannot be observed without water diffusion. Here, both diffusion parameters have limits; beyond that, the system produces either non-homogeneous stationary patterns or goes for extinction. Nevertheless, these chaotic solutions transform into regular moving stripes in the presence of sufficient groundwater flow [see Fig. 5(c)]. As discussed earlier, the temporal model predicts catastrophic shifts for the parameter values in the Hopf unstable domain. This devastating shift can be seen with small diffusion coefficients and the absence of groundwater flow. For instance, the spatial model settles to the water-only steady state for d u = d w = 1. In addition, it produces chaos for some parametric combinations of the diffusion parameters in the absence of groundwater flow [see Fig. 5(a)], and hence, slowly diffusing water can revive plants and herbivores from extinction. The spatial model produces the Turing and non-homogeneous stationary patterns for higher diffusion rates for the water and herbivores (e.g., d u = 500 and d w  100) on the right side of the blue surface in Fig. 1(a) in the absence of groundwater flow. These patterns are robust because they remain stationary by shifting the background parameter value  left-to-right or right-to-left of the magenta surface in Fig. 1(a). To verify these parametric shifts, we first simulate the spatial model till t = 2500 for  = 0.1 [see Fig. 6(a)]. Then, the final non-homogeneous stationary solution is considered as the initial condition for  = 0.11, and we rerun the simulation extended to t = 3000 [see Fig. 6(b)]. We call this shifting forward shifting. In this case, the parameter values  = 0.1 and  = 0.11 are on the magenta surface's left and right in Fig. 1(a). For this shifting of the parameter , the temporal model experiences catastrophic shifts; however, the spatial model leads to a coexistence scenario. On the other hand, the backward shifting of the parameter  from 0.11 to 0.1 also produces the same type of result. Furthermore, the system depicts the same behaviours for forward and backward shiftings for the parameters  and . The dynamics of the spatial model differ in this shifting for lower diffusion rates. The spatial model depicts the homogeneous stationary solution in the region left to the magenta surface in Fig. 1(a) and either a chaotic solution or water-only stationary solution in the region right to the magenta surface in Fig. 1(a) depending on the diffusion parameters. The forward shifting shows extinctions for the plants and herbivores. But the backward shift settles the system either to the coexisting steady state or the water-only steady state. In this case, the solution settles at the coexisting steady-state after the shifting when the system produces chaos before shifting. Otherwise, this shifting can not save the plants and herbivores from extinction. Furthermore, the stripe solutions can exist for such parameter values in the presence of water flow. The construction of sustainable development in ecological civilization is a major concern nowadays. Ecosystems often face disruptions due to natural changes or external forces (e.g., climate, human intervention), particularly semiarid ecosystems. Researchers have been trying to understand these interferences using mathematical models and predict their effect on the ecosystem's functioning, and here is no exception. We have developed a mathematical model comprising water, plant, and herbivore interactions. Studying different parametric ranges enables us to find the vivid possible outcomes (e.g., coexistence, catastrophic shift, stationary patterns, long transients, etc. ), for which it can capture the realistic scenarios of the organisms in semiarid ecosystems. The parameter values have been estimated based on the actual data. The silico trials of the temporal model reveal several exciting results, such as herbivores' over-consumption of plants can cause vegetation collapse and dependencies (e.g., herbivores), excessive rainfall reduces plant growth, etc. These extinctions can be understood by the model described here in the region right to the magenta surface in Fig. 1(a). The plant and herbivore can coexist in between the blue and magenta surfaces in Fig. 1(a) with respect to all of the bifurcation parameters: water supply, plant consumption by the herbivores, and the conversion factor.
In all three cases, restoring the coexisting state from extinction is impossible in the temporal setup, leading to a catastrophic shift. This catastrophic shift occurs due to the sub-critical Hopf bifurcation in the three-component model. No existing literature has identified the subcritical Hopf bifurcation behind such catastrophic change. Furthermore, the modified temporal model always has the water-only equilibrium point, which is linearly stable, and it does not have the coexisting equilibrium point for small parameter values of the bifurcation parameters (water-supply, consumption and conversion rates) [see Fig. 1]. Vegetation abundance plays a pivotal role in semiarid ecosystems by influencing local climate conditions. An increment in the water supply () helps coexist because agricultural production depends on water, which is observed in the considered model. Generally, the water assists positive plant growth, and more plant abundance supports the herbivores to survive. However, intensive rainfall can trigger critical transitions by destabilizing the system and altering the balance between vegetation growth, water availability and grazing by the herbivores, which leads to a significant yield loss. Furthermore, above a limit of water supply, depending on the other parameters, the ecosystem collapses. Semiarid ecosystems host diverse plant and animal communities that interact with each other and the environment. However, disruptions in these interactions can trigger critical transitions by destabilizing the ecological network [36]. This can be predicted through the considered model by increasing the consumption rate of the plants by the herbivores [see Fig. 1]. The model shows that the herbivore species require a minimum food supply for survival. On the other hand, a higher consumption rate causes plant extinction due to over-exploitation, which further affects their dependencies, e.g., herbivores, and in the end, both species die [see Fig. 1]. It also shows that a higher conversion rate () causes extinction, which proves that the herbivore cannot convert the maximum energy from the plant to power its life processes and to build more body tissues. This further validates that energy is lost whenever it is transferred from one trophic level to the next, called an energy pyramid, which decreases by going to higher trophic levels. In the last few decades, theoretical and empirical research has suggested that spatial heterogeneity strongly influences the distribution of organisms within their habitat and how populations colonize new areas and interact within a landscape [37]. The dynamics of a spatial model remain similar to the temporal model for lower diffusion rates but can be completely different for higher diffusion rates. We have shown different types of Turing stationary patterns in the absence of water flow in Figs. 2 and 3(a). Furthermore, Fig. 2 shows that water diffusion is important in the resulting non-homogeneous distributions. In an ecosystem, the ground surface is also an essential factor. Water stands on flat ground, which primarily restricts plants' growth by blocking the movement of oxygen into the soil and carbon dioxide out of the soil. The system produces the Turing pattern only in the presence of water diffusion [see Fig. 3(a)]. This shows groundwater diffusion is crucial for the plant population and overall system. Also, we have demonstrated that the Turing patterns can move from one spatial point to another in the presence of water flow [see Fig. 3(b)], depending on the flow speed and its direction. An increase in the water flow speed can alter the geometry of patches. Still, in the end, they form regular strips moving from one end to another [see Fig. 3(d)]. Therefore, these travelling strips are also Turing patterns but are not stationary. Thus, the Turing spatial structures further show the ecosystem's overall resilience. The temporal model predicted excessive rainfall could destroy the plant and herbivore community. This extinction scenario can occur for the spatial model for higher rainfall without groundwater flow and low or zero water diffusion parameter values. But, for a higher diffusion, the system shows nonhomogeneous stationary patterns, which are not Turing patterns [34]. These non-homogeneous stationary patterns can be sustained through forward and backward shifts as these parametric shifts can come at any time due to natural or human interventions. This further shows that the system is resilient in high water diffusion. Therefore, the spatial patterns can buffer the system against fluctuations, preventing rapid and catastrophic shifts. Furthermore, the water flow drives the system towards regular stripe patterns, which are more resilient, and it can be applied to other types of models where only diffusion is present. It is also true that diffusion can induce irregular distribution in a spatial model [38,39]. In the considered model, the temporal counterpart suggests the extinction of plant and herbivore species beyond the Hopf bifurcation [see Fig. 1]; however, the spatial model on the flat ground shows irregular time-varying distributions for some combinations of the diffusion parameters for the water and herbivores. This shows that self-organization allows the organisms' diversity and saves them from extinction. In addition, this chaos converted into regular strips moving from one end to another in the presence of water flow. Generally, low water diffusion occurs when the ground surface contains excessive rocks, and low diffusion occurs for herbivores when they are ill. Below these low diffusion parameters, plants and herbivores die; hence, the chaos introduced by the sub-critical Hopf bifurcation can be considered as the early warning signal of the catastrophic shift [40][41][42]. Different regular patterns exhibit self-organizing properties that make them resilient to some environmental changes, while irregular patterns are more susceptible to catastrophic shifts under abrupt environmental conditions. Interestingly, the three-component spatio-temporal model proposed here can predict more realistic scenarios in an ecosystem than the currently existing literature. The silico trials show that the presence of herbivores in an ecosystem can cause catastrophic shifts; however, sufficient water diffusion or water flow can sustain the ecosystem in the long run. This shows that diversity can accelerate the extinction rate of species if we do not take precautions.
Nevertheless, the modified model can be applied to the other types of herbivores, and in that case, the parameter values related to the herbivore species have to be calibrated. Additionally, our findings also set the stage for future studies exploring different natural and man-made factors in ecosystems, such as the spreading of fires, diseases, pollution, and introduction of invasive species [43]. Researchers have been working on a plant-water model in a semiarid environment [33], and it has been studied due to its simplicity and rich mathematical analysis. This spatio-temporal model is given by: where u(x, y, t) and v(x, y, t) denote the water and the plant densities, respectively. The uniform water supply rate is k 1 , which is lost due to evaporation at the rate of k 2 u. k 3 is the maximum rate of water consumption by plants, k 4 is the yield of plant biomass per unit of water consumed, and k 5 is the plant mortality rate. The parameter a 1 represents the water downhill speed, and d 2 is the plant dispersal diffusion rate. In this work, we introduce herbivore dynamics into the existing model (1), and these herbivores fed the vegetation. To maintain the model's simplicity, the Holling type II functional response is considered with a constant death rate of the herbivore. Along with this, we also include water diffusion in the model [44][45][46], which can capture the movement of surface water induced by spatial differences in infiltration rate [47]. In addition, the water downfall in the x and y directions is considered separately in the model to capture a wide range of realistic scenarios. Including all these factors in Klausmeier's model [33], we obtain: where u(x, y, t), v(x, y, t) and w(x, y, t) denote the densities of water, plants, and herbivores, respectively. Here, d 1 is the water diffusion, a 2 is the water flow in the y-direction, k 6 represents the rate of plant consumption by herbivores, k 7 is the half-saturation constant, k 8 is the conversion parameter where the plant biomass converts into herbivore biomass, and k 9 is the herbivore consumption rate at zero population growth. We non-dimensionalize this model by substituting where The habitat's area is large in an ecosystem; the periodic boundary conditions for all three variables are the most suitable in this case [33]. All the numerical simulations are done for the model on a 200 x 200 domain (in dimensional terms, 10,000 m 2 ) [33]. The initial conditions for each of the simulations are considered as a heterogeneous perturbation around the homogeneous steadystate: , and w 0 (x, y) = w * + o w ij with o = 10 4 , where  u ij ,  v ij , and  w ij are spatially uncorrelated Gaussian white noise terms. Parameter choice: The parameter values involved in the original model are given in [33]. We consider the species ostriches as herbivores in the model. Ostriches feed the plant biomass; in captivity, they are typically fed between 3 and 4 lbs (1.3 kg to 1.8 kg) of food daily ( In addition, for the spatial dimension, they need at least 1/3 of an acre (approximately 1349 m 2 ) per pair of birds. Here, we assume that each bird needs 1349 m 2 of spatial habitat to live, the same as a pair of birds. Therefore, k 6 = (1.3 to 1.8)x365/1349 = 0.352 to 0.487 kg m 2 individual 1 The mortality rate of baby ostriches is higher than that of juvenile ostriches. Combining both, the mortality rates of the ostriches range from 15% to 50% per year, i.e., k 9 = 0.15 to 0.5 year 1 . The equilibrium points for the water-plant model are the solutions of the nonlinear equations:   u  uv 2 = 0, and uv 2  v = 0. This system has at most three solutions, and they are (u 0 , v 0 ) = (, 0), (u 1 , v 1 ) = ((   2  4 2 )/2, 2/(   2  4 2 )), and (u 2 , v 2 ) = (( +  2  4 2 )/2, 2/( +  2  4 2 )). The feasibility condition for the last two equilibrium points is  > 2. Therefore, the water-plant model has three equilibrium points at most. The linear stability analysis shows (u 0 , v 0 ) and (u 1 , v 1 ) are stable, and (u 2 , v 2 ) is unstable. For studying Turing instability around the stable equilibrium point (u 1 , v 1 ), we substitute u = u 1 + oa 1 exp (t + i(lx + my)) and v = v . Similarly, we find the equilibrium points for the system (3), which are the solutions of the nonlinear equations: Here, the case  =  is not true because we fix h = 1. The Turing instability condition for this model can be studied following the same procedure as outlined for the water-plant model. The authors declare no competing interests. The authors declare no competing interests.