††thanks: This article has been accepted for publication in Physical Review E,
link:
https://journals.aps.org/pre/accepted/23078R21G221d52f82055f69e43409606ae97775c

# Effect of non-local grazing on dry-land vegetation dynamics

Mrinal Kanti Pal<EMAIL_ADDRESS>Swarup Poria  Department of
Applied Mathematics, University of Calcutta, 92 APC Road, Kolkata-700009,
India

###### Abstract

Dry-land ecosystem has turned into a matter of grave concern, due to growing
threat of land degradation and bioproductivity-loss. Self-organized vegetation
patterns are a remarkable characteristic of these ecosystems; apart from being
visually captivating, patterns modulate the system-response to increasing
environmental stress. Empirical studies hinted that herbivory is one the key
regulatory mechanism behind pattern formation and overall ecosystem
functioning. However most of the mathematical models have taken a mean-field
strategy to grazing; foraging has been considered to be independent of spatial
distribution of vegetation. To this end, an extended version of the celebrated
plant-water model due to Klausmeier, has been taken as the base here. To
encompass the effect of heterogeneous vegetation distribution on foraging
intensity and subsequent impact on entire ecosystem, grazing is considered
here to depend on spatially weighted average vegetation density, instead of
density at a particular point. Moreover, varying influence of vegetation at
any location over gazing elsewhere, is incorporated by choosing suitable
averaging function. A comprehensive analysis demonstrates that inclusion of
spatial non-locality, alters the understanding of system dynamics
significantly. The grazing ecosystem is found to be more resilient to
increasing aridity than it was anticipated to be in earlier studies on non-
local grazing. The system-response to rising environmental pressure is also
observed to vary depending on the grazer. Obtained results also suggest
possibility of multi-stability due the history-dependence of system-response.
Overall, this work indicates that the spatial heterogeneity in grazing
intensity has a decisive role to play in the functioning of water-limited
ecosystems.

##  I Introduction

Dry-lands, consisting of arid and semi-arid region and the dry subtropics,
make up roughly 41% world’s land mass and support approximately 38% human
population Reynolds  _et al._ (2007). Almost 10 to 20% of these areas face
acute land degradation and unfortunately this percentage is expected to grow
due to global warming Reynolds  _et al._ (2007); Gowda  _et al._ (2018). Dry-
land ecosystem has become a key priority for ecologists due to rising concerns
about desertification and biodiversity loss, and corresponding effect on
ecosystem functioning. Plants and vegetation constitute the underlying energy
base of all trophic levels and thus play a crucial role in keeping these water
limited ecosystems afloat.

Varying water stress often results in triggering self-organization of spatial
vegetation patterns von Hardenberg  _et al._ (2001). Since the pioneering work
Ref. Macfadyen (1950) in 1950, a substantial number of studies have documented
large scale spatial vegetation patterns in dry-lands via aerial photographs
and satellite images Deblauwe  _et al._ (2012); Tlidi  _et al._ (2018);
Fernandez-Oto  _et al._ (2019a). Apart from regular labyrinthine grass
patterns in arid or semi-arid landscapes, irregular patterns like groves
within grasslands or spots of uncovered ground within a grass matrix, have
been observed worldwide Deblauwe  _et al._ (2008); Martínez-García  _et al._
(2014); Fernandez-Oto  _et al._ (2014); Escaff  _et al._ (2015); Fernandez-Oto
_et al._ (2019b); Parra-Rivas and Fernandez-Oto (2020); Clerc  _et al._
(2021). These patterns are considered to be the key to understand the
processes responsible for desertification, and these understandings can be
utilized to combat the catastrophic effect of climate change Dawes and
Williams (2016); Mau  _et al._ (2015). The emerging characteristics of the
self-organization mechanisms govern the system engineering at ecosystem levels
like primary or secondary production, flexibility against increasing
environmental pressure, stability Bastiaansen  _et al._ (2018); Guttal and
Jayaprakash (2007). Due to the vast spatial and temporal scale of formation of
these patterns, mathematical models have become the primary tool in studying
vegetation dynamics (see Ref. Zelnik  _et al._ (2013) for review). Several
researchers attributed this self-organization of patterns to different
processes: positive feedback between vegetation biomass and water infiltration
HilleRisLambers  _et al._ (2001), competition among vegetation patches for
ground-water due to uptake by roots von Hardenberg  _et al._ (2001); Gowda
_et al._ (2014), non-local water uptake by laterally extended roots and
enhancement of root system with biomass growth Gilad  _et al._ (2004); Meron
_et al._ (2007), or plant-plant interaction only Lefever and Lejeune (1997);
Tlidi  _et al._ (2020). To incorporate the dispersal of plants and water
movement, several models based on partial differential equations have been
proposed HilleRisLambers  _et al._ (2001); Gilad  _et al._ (2007); Messaoudi
_et al._ (2020). One such model that have been a subject of several extension
in last two decade is Klausmeier model Klausmeier (1999) (which will be
detailed in next section) due to its lucidity and fundamental nature Eigentler
and Sherratt (2018, 2019); Sun  _et al._ (2021); Xue  _et al._ (2020); Kabir
and Gani (2022). This plant-water model is basically a reaction-diffusion-
advection system in which the water-uptake feedback loop and spatial
displacement of plant-water, have been taken care of in a minimalistic
approach.

Foraging by herbivores has long been identified as one of the principle
influencer on the physiognomy, structure, and functioning of vegetation,
ranging from landscape scale to a single plant systems Hobbs (1996); Noy-Meir
(1975). Grazing ecosystems supply large amounts of consumable protein, however
too much anthropogenic foraging activities makes ecosystem increasingly
vulnerable to stressful environments. Several studies Mueller  _et al._
(2014); Rietkerk and van de Koppel (1997) have emphasized that vegetation
patterns and foraging by herbivores are interlinked and are sensitive to
degradation under extreme environmental events, like droughts. However, most
of the prevailing spatial models on dry-land ecosystem have considered the
loss of vegetation due to herbivore grazing in a marginalized way; grazing
have been thought to be proportional to the vegetation density and modeled by
adding a constant term to the plant senescence. But, empirical studies have
revealed that foraging by herbivores depends on several factors like spatial
distribution of vegetation, quality of forage, behavior of grazer Focardi  _et
al._ (1996). A greater portion of grazers gets attracted by places with higher
vegetation concentration, thereby resulting in in-homogeneous grazing
pressure. Recently, in an elegant approach, Siero et al. Siero (2018); Siero
_et al._ (2019) have derived a grazing term that incorporates the underlying
effect of spatially heterogeneous vegetation distribution over foraging and
implemented it in the generic plant-water model by Klausmeier Klausmeier
(1999).

In Ref. Siero (2018); Siero  _et al._ (2019), Siero et al. have used mean-
density dependent response term to showcase the dependence of grazing pressure
at any position on vegetation elsewhere. This approach intrinsically presumes
herbivory at any particular location to depend equally on vegetation
everywhere. However, in real scenario the distance between grazer and
vegetation also plays a decisive role Siero (2018); Gordon and Prins (2008);
Worton (1989). Navigation to the forage, i.e. detecting and traveling to the
food item, depends on the characteristic of the grazer, like sight and
olfactory cues Gordon and Prins (2008); Roese  _et al._ (1991). Hence, the
grazing strategy of a herbivore will be more dependent on availability of
vegetation within a certain range, rather than the whole domain. To this end,
in this current work, the influence of vegetation over grazing at any
particular location is considered to vary with the intermediate distance. This
work primarily focuses on addressing two interconnected questions: (1) What
impact does this change bring to the self-organization of patterns? (2) Does
it have any effect on the system-response towards increasing aridity? Here the
mean density is replaced by a spatially weighted average density, using
convolution integral with a normalized weight function, that estimates the
utilization of space for grazing as a function of the distance. Moreover, a
generalized nature is maintained throughout the analysis, so that further
modification can be made by using data-driven weight functions.

In Sect. II.1, the working model is proposed after describing the Klausmeier
model in detail and then a thorough mathematical analysis is carried out in
Subsection II.2.2. In Sect. III, numerical simulations are carried out to
evaluate the resilience of the system toward varying environmental pressure.
In a parsimonious approach, two particular types of weight functions are
chosen here for having a comparative understanding and differentiating the
obtained results from earlier studies Siero (2018); Siero  _et al._ (2019).
Finally, all the findings are ecologically interpreted in Sect. IV and few
potential extensions to this model are suggested.

##  II Model Description and Analyses

###  II.1 Model description

This section first reviews the original model introduced by Klausmeier, which
is then extended by incorporating non-local grazing to the equation for
vegetation. Then, Subsection II.2 presents the systematic mathematical
analysis of the modified system.

For modeling dry-land vegetation patterns on sloped terrain, Klausmeier
Klausmeier (1999) proposed a reaction-diffusion-advection system:

$\displaystyle\begin{cases}\begin{split}&\frac{\partial W}{\partial T}=A-LW-
RWN^{2}+V\frac{\partial W}{\partial X_{1}},\\\ \\\ &\frac{\partial N}{\partial
T}=RJWN^{2}-MN+D\nabla^{2}N,\end{split}\end{cases}$ (1)

where $\nabla^{2}=\frac{\partial^{2}}{\partial
X_{1}^{2}}+\frac{\partial^{2}}{\partial X_{2}^{2}}$ is the Laplacian operator.
$W(T;\vec{X})$ and $N(T;\vec{X})$ are the density of surface water and plant
biomass respectively at location $(T;\vec{X})\in{(T>0)\times\mathbb{R}^{2}}$.
The first equation of system (1) depicts the dynamics of water; here the
source term $A$ is the uniform rate at which water is added via precipitation.
The term $-LW$ accounts for the water-loss because of evaporation and the
advective term $V\frac{\partial W}{\partial X_{1}}$ symbolizes the downhill
movement of water along the sloping ground. The term $RWN^{2}$ corresponds to
the water-uptake by the plant-roots, here the appearance of non-linearity is a
consequence of the positive feedback between plant growth and water seepage.
The second equation presents vegetation dynamics, where $RJWN^{2}$ presents
the plant growth; $J$ being the yield of vegetation biomass per unit
consumption of water. The term $-MN$ specifies the loss of vegetation due to
natural death and grazing of plants by herbivores. Lastly the diffusion
$D\nabla^{2}N$ is for modeling the spatial spread of vegetation by means of
seed dispersal or clonal growth.

While considering banded vegetation on a slope, the Klausmeier model (1) fails
to produce stationary patterns in flat land. To this end several researchers
Zelnik  _et al._ (2013); Kealy and Wollkind (2012); Siteur  _et al._ (2014)
have omitted the advective term (i.e. $V=0$). Further they had added a soil
water diffusion term $E\nabla^{2}W$ to the first equation of system (1) to
encapsulate the movement of surface water due to the spatial heterogeneity in
infiltration rate Rietkerk  _et al._ (2002) :

$\displaystyle\begin{cases}\begin{split}&\frac{\partial W}{\partial T}=A-LW-
RWN^{2}+E\nabla^{2}W,\\\ \\\ &\frac{\partial N}{\partial
T}=RJWN^{2}-MN+D\nabla^{2}N.\end{split}\end{cases}$ (2)

In model system (1) and (2), senescence of vegetation biomass is assumed to be
independent of the spatial distribution of vegetation and grazing by
herbivores has been considered as included in this linear term $MN$. But
empirical studies Focardi  _et al._ (1996); Fryxell  _et al._ (2004) on semi
arid region tells that grazing depends on spatial distribution of vegetation;
herbivores get attracted more to superior forage. This leads to spatially
heterogeneous grazing pressure unlike system (2) where it is constant. To
incorporate this phenomena, Siero et al. Siero (2018); Siero  _et al._ (2019)
considered model (2) in one spatial dimension and modified it by adding a
density dependent response function $\tilde{G}$. Moreover, the grazing
pressure at any spatial point depends on vegetation density not only at that
point, but also elsewhere. Taking this non-locality in account they have taken
$\tilde{G}$ to be a function of mean vegetation density $\tilde{N}$, where

$\displaystyle\tilde{N}:=\frac{1}{|\Omega|}\int_{\Omega}N(z)dz,$ (3)

$|\Omega|$ denotes the length of the spatial domain $\Omega\in\mathbb{R}$
under consideration. Three different type grazing functions have been used,
the first type coincides with existing notion of linear mortality in
Klausmeier model, i.e. grazing pressure is considered to be constant:
$\tilde{G}=\tilde{M}_{loc}$. Type II is for sustained grazing (e.g. livestock
farming), foraging is optional in this case as food shortage will be
reimbursed by supplementary food; i.e. demographic response is kept at a
constant level all the time. Considering that herbivores maintain a saturating
functional response (Holling type II), the grazing function is given by
$\tilde{G}=\tilde{M}_{sus}/(\tilde{K}_{sus}+\tilde{N})$, where
$\tilde{K}_{sus}$ is the half persistence level. The third type is for natural
scenario where grazing is obligatory for survival, only a section of grazers
will be able to survive by acquiring sufficient amount of food and rest will
disappear. Approximating this demographic response by a sigmoid function
(Holling type III), the resulting grazing pressure becomes
$\tilde{G}=\tilde{M}_{nat}\tilde{N}/(\tilde{K}_{nat}^{2}+\tilde{N}^{2})$. The
resulting model is given by:

$\displaystyle\begin{cases}\begin{split}&\frac{\partial W}{\partial T}=A-LW-
RWN^{2}+E\nabla^{2}W,\\\ \\\ &\frac{\partial N}{\partial
T}=RJWN^{2}-(M+\tilde{G}(\tilde{N}))N+D\nabla^{2}N,\end{split}\end{cases}$ (4)

where $\nabla^{2}=\frac{\partial^{2}}{\partial X^{2}}$ is the one-dimensional
Laplacian operator. One noteworthy fact is that the non-local term appears
non-linearly in the system. For ease of mathematical analysis, system (4) is
non-dimensionalised with the following substitutions Siero (2018); Siero  _et
al._ (2019):

$\displaystyle w=WJ\sqrt{\frac{R}{L}},\qquad n=N\sqrt{\frac{R}{L}},\quad
x=X\sqrt{\frac{L}{D}},\qquad t=TL,$ $\displaystyle e=\frac{E}{D},\qquad
a=\frac{AJ}{L}\sqrt{\frac{R}{L}},\qquad m_{0}=\frac{M}{L},\qquad
m_{loc}=\frac{\tilde{M}_{loc}}{L},$ $\displaystyle
m_{sus/nat}=\frac{\tilde{M}_{sus/nat}}{L}\sqrt{\frac{R}{L}},\qquad
K_{sus/nat}=\sqrt{\frac{R}{L}}\tilde{K}_{sus/nat}.$

The dimensionless model is given by

$\displaystyle\begin{cases}\begin{split}&\frac{\partial w}{\partial t}=a-w-
wn^{2}+e\frac{\partial^{2}w}{\partial x^{2}},\\\\[5.0pt] &\frac{\partial
n}{\partial t}=wn^{2}-(m_{0}+G(\tilde{n}))n+\frac{\partial^{2}n}{\partial
x^{2}}.\end{split}\end{cases}$ (5)

where

$\displaystyle G(\tilde{n})=\begin{cases}\begin{aligned} &m_{loc}&&\text{for
local grazing (type I)},\\\\[5.0pt]
&\frac{m_{sus}}{K_{sus}+\tilde{n}}&&\text{for sustained grazing (type
II)},\\\\[5.0pt]
&\frac{m_{nat}\tilde{n}}{K_{nat}^{2}+\tilde{n}^{2}}&&\text{for natural grazing
(type III)}.\end{aligned}\end{cases}$ (6)

In this study, system (5-6) is taken as base model. One inherent assumption in
this model is that, grazing by herbivores at any particular spatial point
depends equally on vegetation densities at all other points. But, ecological
studies on browsing and grazing have showed that the rate of consumption of
vegetation biomass depends on the perceptual abilities (e.g. sight, olfaction)
of the animal Gordon and Prins (2008); Roese  _et al._ (1991). When the
observational area is vast (which is the case for most of the empirical
studies on dry-land vegetation dynamics), the grazing strategy of a herbivore
will be more dependent on availability of vegetation nearby, rather than the
whole domain. To capture this ecological fact more realistically, the
definition of mean vegetation density (3) is modified:

$\displaystyle\tilde{n}(x,t)=\int\rho(|x-x^{\prime}|)n(x^{\prime},t)dx^{\prime},$
(7)

where $\rho(x)$ is a normalized kernel function (i.e. $\int\rho(|y|)dy=1$)
accounting for the weighted mean vegetation density. The use of absolute value
$|x-x^{\prime}|$ in definition (7) ensures spatial isotropy of the kernel
function. In a parsimonious representation, here it is assumed that the
influence of vegetation density at any position, over the grazing decreases
with the distance from the grazer; that is why the weight function $\rho$ is
considered to be a monotonically decreasing function of distance in this
study. In the next subsection II.2, at first the model (5-6-7) will be
analyzed mathematically in a general setting and then a Gaussian distribution
function will be used as a particular case for numerical simulation here.
Moreover, periodic boundary conditions will be taken to lessen the effect of
the boundary and to mimic an infinite domain.

###  II.2 Model analysis

####  II.2.1 Homogeneous steady states

This subsection summarizes the existence of homogeneous steady states
(hereafter, HSSs) of system (5-6). For a homogeneous steady state,
$\frac{\partial}{\partial t}w(t,x)=\frac{\partial}{\partial t}n(t,x)=0$ and
$\frac{\partial^{2}}{\partial x^{2}}w(t,x)=\frac{\partial^{2}}{\partial
x^{2}}n(t,x)=0$. Using these in system (5-6),

$\displaystyle a-w-wn^{2}=0,$ (8a) $\displaystyle n(wn-m_{0}-G(\tilde{n}))=0.$
(8b)

Hence there are two possibilities for HSS: a bare state ($\mathcal{B})$ with
$w=a,n=0$ everywhere; and homogeneously vegetated state ($\mathcal{V})$ with
$w=\frac{a}{1+n_{*}^{2}},n=n_{*}>0$, where $n_{*}$ satisfies

$\displaystyle wn_{*}-m_{0}-G(\tilde{n}_{*})=0.$ (9)

It is noteworthy that there can be more than one homogeneously vegetated state
($\mathcal{V})$ depending on parameter values and type of grazing function.
The main focus of this article is on the sustained and natural grazing type,
as the local grazing case have already been studied extensively Sherratt and
Lord (2007).

##### Local grazing

As $G(\tilde{n})=m_{loc}$ is constant in local grazing, it is basically same
as the pre-existing models that considers linear mortality. Apart from the
bare state ($\mathcal{B})$, there also exist two homogeneous steady states
$(w_{S},n_{S})$ and $(w_{N},n_{N})$, when $a>2(m_{0}+m_{loc})$ Sherratt
(2010),

$\displaystyle
w_{S}=\frac{2(m_{0}+m_{loc})^{2}}{a-\sqrt{a^{2}-4(m_{0}+m_{loc})^{2}}},$
$\displaystyle
n_{S}=\frac{a-\sqrt{a^{2}-4(m_{0}+m_{loc})^{2}}}{2(m_{0}+m_{loc})};$
$\displaystyle
w_{N}=\frac{2(m_{0}+m_{loc})^{2}}{a+\sqrt{a^{2}-4(m_{0}+m_{loc})^{2}}},$
$\displaystyle
n_{N}=\frac{a+\sqrt{a^{2}-4(m_{0}+m_{loc})^{2}}}{2(m_{0}+m_{loc})}.$

These two HSSs originate from a saddle-node bifurcation at
$a=2(m_{0}+m_{loc})$. Linear stability analysis against homogeneous
perturbation Sherratt (2005) reveals that the bare state $\mathcal{B}$ is
always stable (for $m_{0}+m_{loc}>0$, which is always true). Moreover,
$(w_{S},n_{S})$ is unstable (actually saddle) and $(w_{N},n_{N})$ is a node
Siteur  _et al._ (2014). But for ecologically meaningful parameter values for
semi-arid environment, $(w_{N},n_{N})$ is stable to homogeneous perturbation.
The location and stability of steady states against spatially uniform
perturbation, has been described in Figure 1. Furthermore, the linear
stability analysis against spatially heterogeneous perturbation Siteur  _et
al._ (2014) shows that the steady state $(w_{N},n_{N})$ can lose its stability
through Turing bifurcation; although the bare state $\mathcal{B}$ always
remains stable.

##### Sustained grazing

In this case, equation (9) takes the form

$\displaystyle wn-m_{0}-\frac{m_{sus}}{K_{sus}+\tilde{n}}=0.$ (10)

Now for uniformly vegetated state $\mathcal{V}$, using the normalization
condition of the kernel, we have

$\displaystyle\tilde{n}=\int\rho(|x-x^{\prime}|)n(x^{\prime},t)dx^{\prime}=n_{*}\int\rho(|x-x^{\prime}|)dx^{\prime}=n_{*}.$
(11)

Substituting $w=\frac{a}{1+n_{*}^{2}},n=n_{*}$ into equation (10) and using
(11), a cubic equation for $n_{*}$ is obtained:

$\displaystyle\begin{split}m_{0}n_{*}^{3}+(m_{0}K_{sus}-a+m_{sus})n_{*}^{2}+(m_{0}-aK_{sus})n_{*}+m_{0}K_{sus}&\\\
+m_{sus}=0&.\end{split}$ (12)

Only the positive solutions of (12) correspond to the physically admissible
$\mathcal{V}$ steady states. As all the parameters are non-negative, the
number of sign changes in the sequence of this cubic polynomial’s coefficients
is either zero or two. Moreover, there are no sign change in (12) when
$a<min(m_{0}K_{sus}+m_{sus},\frac{m_{0}}{K_{sus}})$, hence Descartes’ rule of
signs implies that equation (12) will have no positive solution in that case.
Therefore, for sufficiently small level of precipitation, the bare state
$\mathcal{B}$ will be the only HSS. Numerical simulation shows that for $a$
greater than some threshold value, there actually exists two branches of
homogeneously vegetated steady state $\mathcal{V}$ (See Figure 1).

(a)

(b)

(c)

Figure 1: Existence of HSSs and their stability against spatially homogeneous
perturbation: chained blue line( ) is the bare state $\mathcal{B}$ (always
stable); solid black line ( ) and dashed red line ( ) are the stable and
unstable branch of $\mathcal{V}$, respectively. The * denotes the critical
transition from one state to another, known as Tipping point. Parameters are
same as Ref. Siero  _et al._ (2019):
$e=500,m_{0}=0.225,m_{loc}=0.225,m_{sus}=m_{nat}=1.5,K_{sus}=K_{nat}=0.3$. (a)
Local Grazing, (b) Sustained grazing, (c) Natural grazing

##### Natural grazing

Proceeding in a similar fashion like sustained grazing, a quartic polynomial
for $n_{*}$ is obtained:

$\displaystyle\begin{split}m_{0}n_{*}^{4}+(m_{nat}-a)n_{*}^{3}+(K_{nat}^{2}+1)m_{0}n_{*}^{2}+(m_{nat}-aK_{nat}^{2})n&\\\
+m_{0}K_{nat}^{2}=0.&\end{split}$ (13)

Depending on model parameters, this quartic equation may or may not have
positive solution. If the model parameters are such that
$a<min(m_{nat},\frac{m_{nat}}{K_{nat}^{2}})$, then the sequence of
coefficients of the polynomial (13) will have no sign change; i.e. there will
be no uniformly vegetated state $\mathcal{V}$ for low precipitation level.
Numerical methods are used for finding roots of this quartic equation which
show that, for the parameter choice Ref. Siero  _et al._ (2019) (which was
based on earlier studies and empirical data), equation (13) has two positive
zeros when $a$ values are greater than some critical threshold. Consequently
for sufficiently large values of $a$, there will be two uniformly vegetated
states $\mathcal{V}$ apart from the bare state $\mathcal{B}$ (See Figure 1).

####  II.2.2 Stability analysis

Linear stability analysis is a broadly utilized tool to acquire the temporal
evolution of small perturbations to the homogeneous steady states of the
system. When given perturbation gets amplified over time, the steady state
will be unstable; but if the perturbation decays with time, taking the system
back to the homogeneous stationary state, then that steady state will be
stable. For notational convenience, let $(w_{eq},n_{eq})$ be a homogeneous
stationary state (it can be bare state $\mathcal{B}$ or vegetated state
$\mathcal{V}$). If a small perturbation is given to the HSS $(w_{eq},n_{eq})$,
the resulting density will be

$\displaystyle\begin{bmatrix}w(x,t)\\\
n(x,t)\end{bmatrix}=\begin{bmatrix}w_{eq}\\\
n_{eq}\end{bmatrix}+\begin{bmatrix}\epsilon_{1}\\\
\epsilon_{2}\end{bmatrix}\psi(x,t)$ (14)

where $|\epsilon_{1}|\ll 1$ and $|\epsilon_{2}|\ll 1$. Now plugging equation
(14) into equation (7) and using the normalization condition of the kernel
function,

$\displaystyle\tilde{n}(x,t)$
$\displaystyle=\int\rho(|x-x^{\prime}|)(n_{eq}+\epsilon_{2}\psi(x^{\prime},t))dx^{\prime}$
$\displaystyle=n_{eq}+\epsilon_{2}\int\rho(|x-x^{\prime}|)\psi(x^{\prime},t)dx^{\prime}.$
$\displaystyle\text{So,}\quad G(\tilde{n})$
$\displaystyle=G\bigl{(}n_{eq}+\epsilon_{2}\int\rho(|x-x^{\prime}|)\psi(x^{\prime},t)dx^{\prime}\bigr{)}.$

As the perturbation considered here is very small, the function $G$ can be
expanded using Taylor series:

$\displaystyle
G(\tilde{n})=G(n_{eq})+\epsilon_{2}G^{\prime}(n_{eq})\int\rho(|x-x^{\prime}|)\psi(x^{\prime},t)dx^{\prime}+O(\epsilon_{2}^{2})$
(15)

where
$G^{\prime}(n_{eq})=\bigl{(}\frac{dG}{d\tilde{n}}\bigr{)}_{\tilde{n}=n_{eq}}$.
Now substituting (14) into the first equation of system (5) and neglecting the
higher order terms in $\epsilon_{1}$ and $\epsilon_{2}$ we have,

$\displaystyle\begin{split}\epsilon_{1}\frac{\partial\psi}{\partial
t}=\epsilon_{1}e\frac{\partial^{2}\psi}{\partial
x^{2}}-\epsilon_{1}(1+n_{eq}^{2})\psi-\epsilon_{2}2w_{eq}n_{eq}\psi&\\\
+(a-w_{eq}-w_{eq}n_{eq}^{2})&.\end{split}$ (16)

As $(w_{eq},n_{eq})$ is homogeneous steady state of system (5), we have
$a-w_{eq}-w_{eq}n_{eq}^{2}=0$. Hence (16) becomes

$\displaystyle\epsilon_{1}\frac{\partial\psi}{\partial
t}=\epsilon_{1}e\frac{\partial^{2}\psi}{\partial
x^{2}}-\epsilon_{1}(1+n_{eq}^{2})\psi-\epsilon_{2}2w_{eq}n_{eq}\psi.$ (17)

Substituting (14) and (15) into the second equation of system (5) and
neglecting higher order terms in $\epsilon_{1}$ and $\epsilon_{2}$ we have,

$\displaystyle\epsilon_{2}\frac{\partial\psi}{\partial t}=$
$\displaystyle\epsilon_{1}n_{eq}^{2}\psi+\epsilon_{2}\biggl{\\{}\frac{\partial^{2}\psi}{\partial
x^{2}}+\bigl{(}-m_{0}+2w_{eq}n_{eq}-G(n_{eq})\bigr{)}\psi$ $\displaystyle-
n_{eq}G^{\prime}(n_{eq})\int\rho(|x-x^{\prime}|)\psi(x^{\prime},t)dx^{\prime}\biggr{\\}}$
(18)

Applying Fourier transform to this linearized integro-differential system
(17-II.2.2) for the evolution of $\psi$, we have

$\displaystyle\begin{bmatrix}\epsilon_{1}\\\
\epsilon_{2}\end{bmatrix}\frac{\partial\hat{\psi}(k,t)}{\partial
t}=L\begin{bmatrix}\epsilon_{1}\\\ \epsilon_{2}\end{bmatrix}\hat{\psi}(k,t),$
(19)

where

$\displaystyle
L=\begin{bmatrix}-ek^{2}-1-n_{eq}^{2}&-2w_{eq}n_{eq}\\\\[17.22217pt]
n_{eq}^{2}&-k^{2}-m_{0}+2w_{eq}n_{eq}-G(n_{eq})\\\
&-n_{eq}G^{\prime}(n_{eq})\hat{\rho}(k)\end{bmatrix}.$ (20)

$k$ is known as wave-number, $\hat{\psi}(k,t)=\int exp(ikx)\psi(x,t)dx$ is the
Fourier transform of the perturbation and similarly $\hat{\rho}(k)$ is the
Fourier transform of the kernel.

Now, considering $\hat{\psi}(k,t)\propto exp(\lambda(k)t)$, equation (19)
becomes

$\displaystyle\bigl{(}L-\lambda(k)I\bigr{)}\begin{bmatrix}\epsilon_{1}\\\
\epsilon_{2}\end{bmatrix}=\begin{bmatrix}0\\\ 0\end{bmatrix}.$ (21)

Non-zero solution of the linear system (21) exists if and only if

$\displaystyle Det(L-\lambda(k)I)=0,$

which yields a quadratic equation in $k^{2}$ for the linear growth rate
$\lambda$:

$\displaystyle\lambda^{2}-Trace(L)\lambda+Det(L)=0.$ (22)

Equation (22) is known as the dispersion relation. Linear stability of the
HSSs depends on the essential spectrum (which consists of the two solutions
$\lambda_{1},\lambda_{2}$ of the dispersion relation) of the linearized
system. $(w_{eq},n_{eq})$ is linearly stable if
$\operatorname{Re}(\lambda_{1})\leq 0$ and $\operatorname{Re}(\lambda_{2})\leq
0$, and unstable if either these two inequalities gets altered.

Here stability against both spatially uniform perturbation and heterogeneous
perturbation will be considered.

##### Spatially homogeneous perturbation

Note that, the case of spatially homogeneous perturbation corresponds to the
wave-number $k=0$. Then the linearized matrix (20) takes the form,

$\displaystyle\begin{bmatrix}-1-n_{eq}^{2}&-2w_{eq}n_{eq}\\\\[17.22217pt]
n_{eq}^{2}&-m_{0}+2w_{eq}n_{eq}-G(n_{eq})-n_{eq}G^{\prime}(n_{eq})\end{bmatrix}.$
(23)

If $(w_{eq},n_{eq})$ is taken to be the bare state $\mathcal{B}(a,0)$, then
(23) becomes,

$\displaystyle L_{\mathcal{B}}=\begin{bmatrix}-1&0\\\\[17.22217pt]
0&-m_{0}-G(0)\end{bmatrix}.$ (24)

From (24), it can be concluded that dispersion relation (22) for $\mathcal{B}$
will always have negative roots, i.e. bare state $\mathcal{B}$ is linearly
stable against spatially homogeneous perturbations for all three grazing
types.

Now, if we consider $(w_{eq},n_{eq})$ to be the uniform vegetated state
$\mathcal{V}$, then it will satisfy condition (9). So equation (23) takes the
form:

$\displaystyle
L_{\mathcal{V}}=\begin{bmatrix}-1-n_{eq}^{2}&-2w_{eq}n_{eq}\\\\[17.22217pt]
n_{eq}^{2}&w_{eq}n_{eq}-n_{eq}G^{\prime}(n_{eq})\end{bmatrix}.$ (25)

Unlike the case for bare state, the sign of the real parts of
$\lambda_{1},\lambda_{2}$ (which are actually the eigenvalues of matrix
$L_{\mathcal{V}}$) depend on the value of $w_{eq},n_{eq}$ and the grazing
function $G$. Using numerical schemes, the existence and linear stability of
these HSSs have been evaluated for all the three types of grazing (Figure 1).
Observe that up to a certain threshold value of precipitation $a$, bare state
$\mathcal{B}$ is the only HSS available. But at some critical value of $a$,
saddle-node bifurcation occurs and a stable and an unstable (saddle) branch of
uniformly vegetated state $\mathcal{V}$ appear.

##### Spatially heterogeneous perturbation

Putting $w_{eq}=a$ and $n_{eq}=0$ in matrix (20), it can be observed that the
linearized matrix for bare state $\mathcal{B}$ always have negative
eigenvalues, i.e. bare state $\mathcal{B}$ is stable against heterogeneous
perturbations also.

So far, all the mathematical analysis have been carried out without specifying
the explicit form of kernel function. But, from the linearized matrix (20), it
is evident that the linear stability of the state $\mathcal{V}$ against
heterogeneous perturbation, depends on the Fourier transform of the kernel.
So, it directly depends on the structure of the kernel function. Without loss
of generality, it is assumed that, vegetation at any point have lesser
influence over the grazer the further they are; and after certain cut-off
distance, the influence is too insignificant to consider. Keeping this
assumption in mind, a cut-off Gaussian function Fuentes  _et al._ (2004) is
used here:

$\displaystyle\rho_{1}(x)=\frac{1}{\sigma\sqrt{\pi}\operatorname{erf}(w/\sigma)}exp\biggl{[}-\biggl{(}\frac{x^{2}}{\sigma^{2}}\biggr{)}\biggr{]}\bigl{\\{}\theta[w-x]\theta[w+x]\bigr{\\}}.$
(26)

Here, $\theta$ is the Heaviside function and $w$ denotes the effective radius
of action. The shape of the kernel function is characterized by its width
($\sigma$) and the cut-off length ($w$). For further analysis, the spatial
domain is considered to be finite: $[-L,L]$. The cut-off distance $w$ is
always within the bound of domain length (i.e. $w\leq L$). The term
$\frac{1}{\sigma\sqrt{\pi}\operatorname{erf}(w/\sigma)}$ is the normalization
factor, which ensures that the integral of $\rho_{1}(x)$ over the whole domain
equals to unity and it serves the purpose of being an averaging function. For
this choice of $\rho$, the Fourier transform will be:

$\displaystyle\hat{\rho}_{1}(k)=\frac{exp[-(\sigma
k/2)^{2}]}{2\operatorname{erf}(w/\sigma)}\biggl{[}\operatorname{erf}\biggl{(}\frac{w}{\sigma}+\frac{ik\sigma}{2}\biggr{)}+\operatorname{erf}\biggl{(}\frac{w}{\sigma}-\frac{ik\sigma}{2}\biggr{)}\biggr{]}.$
(27)

To have a comparative understanding of the role of the kernel function, a
second choice of $\rho$ is considered here:

$\displaystyle\rho_{2}(x)=\frac{1}{2w}\bigl{\\{}\theta[w-x]\theta[w+x]\bigr{\\}},$
(28)

and

$\displaystyle\hat{\rho}_{2}(k)=\frac{\sin(kw)}{kw}.$ (29)

Without loss of generality, both the kernels will be considered with $w=L$ for
numerical simulations. In this scenario, $\sigma$ is the parameter which
regulates the variation in dependence of grazing intensity over vegetation
elsewhere. One noteworthy fact is, for $w=L$, kernel choice $\rho_{2}$
actually corresponds to the uniform weight function, which is the case in Ref.
Siero (2018); Siero  _et al._ (2019). This enables us to compare results of
this study with the findings of Ref. Siero (2018); Siero  _et al._ (2019).

Unlike the case for homogeneous perturbation (i.e. $k=0$), here the growth
term $\lambda$ will be dependent not only on the model parameters, but also on
the wave number $k$ of the perturbation. As the saddle branch of steady state
(dashed red part in Figure 1) is already unstable against the homogeneous
perturbation, only the stable node (solid black part) will be considered here.
These stable nodes will get destabilized if the real part of one or both of
the solutions $\lambda(k)$ of dispersion relation (22) become positive. With
the variation in the bifurcation parameter, instability can happen in mainly
two ways: Hopf instability (when $Trace(L)$ goes to positive from being
negative); and Turing instability (when $Det(L)$ becomes negative from being
positive). Numerical simulation reveals that, for realistic choice of
ecological parameters Klausmeier (1999); Siero  _et al._ (2019), these stable
nodes loose stability only by Turing bifurcation, no Hopf instability is
observed in this model. When a uniformly vegetated state is driven out of
Turing stability region, spatio-temporal vegetation pattern appears.

(a)

(b)

Figure 2: Stability scenario against heterogeneous perturbation in $(a,k)$
space. The colored solid lines represent the boundaries of Turing prediction
region for different choice of kernel: cut-off Gaussian $\rho_{1}$ having
$\sigma=1$ ( ), $\sigma=8$ ( ), $\sigma=15$ ( ); and uniform kernel $\rho_{2}$
( ). The chained line ( ) denotes the lower bound of precipitation, below
which the steady state $\mathcal{V}$ ceases to exist . The dashed lines
(colored in aforesaid order) represents the mode of perturbations with largest
growth rate. Parameters are same as Figure 1. (a) Sustained grazing, (b)
Natural grazing

In Figure 2, values of $a$ at which the homogeneously vegetated state
($\mathcal{V}$) becomes Turing unstable, have been derived. To have a
comparative understanding, different structures of kernel function are
considered. The black chained line ( ) marks the lower bound of precipitation
for which uniform vegetated state ($\mathcal{V}$) persists. Inside the region
in $a-k$ space, bounded by the colored solid curves and to the right of the
black chained line ( ), the amplitude of given perturbation grows with time,
i.e. for every $(a,k)$ pair from this region, the dispersion equation (22)
will have at least one solution $\lambda$ having positive real part. This
region is termed as Turing prediction region, Siteur  _et al._ (2014); Siero
_et al._ (2015) as self-organized spatio-temporal patterns can be expected in
this region due to Turing type instability. Moreover, the most unstable mode
of the perturbation (i.e. the maximum of $Re(\lambda(a,k))$) for each of $a$
in Turing instability band is derived, and it is denoted by dashed lines
(colored in their respective order). When precipitation is getting reduced
over time, there will be no pattern formation initially; but at the maximum
precipitation value in the Turing prediction region (i.e. the intersection
points of dashed and solid curves of same color) the HSS looses stability via
Turing bifurcation and spatio-temporal patterns appear. These intersection
points will be denoted by $a_{T}$ hereafter. In other words, the stable branch
(black solid curve) in Figure 1 will end at $a_{T}$, long before reaching the
tipping point, which is basically the level of precipitation at which the non-
spatial version of the system shifts from vegetated state $\mathcal{V}$ to the
less desirable stable state $\mathcal{B}$. It can also be seen from Figure 2
that the pattern formation starts at less arid situation in case of natural
grazing, compared to the sustained grazing. Model runs have confirmed that
patterns generated at $a_{T}$ have a wave-number very close to the $k$ value
of these intersection points. Usually in-homogeneous patterns of vegetation,
generated when random perturbations are applied to the uniformly vegetated
state, are expected to have the mode close to the wave-number of the
perturbation that have largest growth rate. However, which of these wave-
numbers gets chosen at a particular precipitation level, is largely
unpredictable and in section III it will be seen that the history of
environmental variations has significant role to play in this selection of
wave-number.

##  III Response To Changing Precipitation Level

It is well-known that the stability characteristics of a system, acquired by
linear stability analysis, are generally local in nature, i.e. these
predictions holds to be true in a close vicinity to the steady state.
Therefore in this section, to have a complete understanding of the dynamical
scenario, system (5-6) will be studied with gradually varying precipitation
parameter. Earlier studies Siteur  _et al._ (2014); Sherratt (2013) on dry
land vegetation revealed that rising environmental pressure (e.g. degrading
precipitation level) drives the system to a coarsening cascade of transitions
to patterns with increasing wavelengths and for adequately low precipitation
level the system goes to bare state. Here, a spatial domain
$\Omega=[-250,250]$ (resembling 1 km Siero (2018)) with periodic boundary
conditions will be considered for numerical simulation. All other parameter
values (mentioned in Figure 1) are in accordance with previous studies
Klausmeier (1999). Initially for every simulation, the precipitation level is
kept at $a=3(\approx 800$ mm year${}^{-1})$, sufficient to support a uniformly
vegetated state. Then the precipitation level is set to decrease at a constant
rate $\frac{da}{dt}=-10^{-4}(\approx 0.1$ mm year${}^{-1})$ until $a=0$ is
reached Siero  _et al._ (2019). Moreover, to reduce numerical artifacts and
incorporate intrinsic noise sources that deterministic equations fail to
capture, spatially and temporally uncorrelated multiplicative uniformly
distributed noise of small amplitude $0.05\%$ is added to both the components
at every $\frac{1}{4}$ year Siteur  _et al._ (2014). Simulations are carried
out in MATLAB using the Backward Time Centered Space (BTCS) difference scheme.
To avoid discrepancy that may arise near the boundaries of $\Omega$ during
numerical evaluation of the mean density (7), both the kernel functions
$\rho_{1}$, $\rho_{2}$ are modified accordingly Pal and Poria (2021); Pal  _et
al._ (2020). Following earlier studies Siero  _et al._ (2019); Siero (2018),
in the numerical simulations vegetation density $n$ less than $10^{-6}$ has
been considered as zero.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Simulation results for system (5-6) with kernel $\rho_{1}$ (26)
having different $\sigma$ values: (a-b) $\sigma=1$, (c-d) $\sigma=8$, (e-f)
$\sigma=15$. The left side column is for Sustained Grazing (a-c-e) and the
second column (b-d-f) is for Natural Grazing. In all of these model runs,
starting from a homogeneously vegetated state at $a=3$, the precipitation
level is set to decrease at a rate $\frac{da}{dt}=-10^{-4}$

(a)

(b)

Figure 4: Simulation results for system (5-6) with kernel $\rho_{2}$ (28). The
left side column is for Sustained Grazing (a) and the second column (b) is for
Natural Grazing. Starting from a homogeneously vegetated state at $a=3$, the
precipitation level is set to decrease at a rate $\frac{da}{dt}=-10^{-4}$

To showcase the role of kernel function in the system-response (to decreasing
precipitation level), model runs have been performed with both cut-off
Gaussian $\rho_{1}$ and uniform kernel $\rho_{2}$. Different shapes of kernel
$\rho_{1}$ (i.e. $\rho_{1}$ having different widths $\sigma$) have been used
in these simulations, results of three of which are shown in Figure 3. It can
be seen from these figures that for both type of grazing, as the precipitation
reduces, patterns of vegetation appears just after the initial uniformly
vegetated state loses its stability through Turing bifurcation. This
transition of stability happens at the critical precipitation level $a_{T}$,
mentioned in the previous section. Moreover, wavelengths of these initial
patterns matches very closely with the estimation made by linear stability
analysis. For example, in case of $\sigma=8$ with sustained grazing, the
bifurcation point in Figure 2 is $(a_{T},k_{T})=(2.954,0.395)$, so the
critical wavelength would be $2\pi/k_{T}\approx 15.91$, hence the resulting
pattern is expected to have approximately 31 ridges in the computational
domain [-250,250], which is exactly the case in Figure 3. Furthermore, Figure
2 reveals that the $a_{T}$ value for $\sigma=1$ is much greater than $a=3$
(which is the starting point of the simulation), so in this case pattern
emerges right from the beginning (Figure 3-Figure 3). But as the width of the
weight function $\rho_{1}$ increases, this critical precipitation level
$a_{T}$ decreases. Figure 3-Figure 3 reflect the same phenomena, a
homogeneously vegetated state persists until precipitation level reaches
$a_{T}$. However, the rate of lowering of $a_{T}$ falls off with increment in
$\sigma$ and after a certain value (which is $\sigma\approx 15$ for aforesaid
parameter choice) no more significant change in the value of $a_{T}$ occurs.
With further decline in precipitation the wavelengths of these patterns
remains constant initially, and then undergo number of sudden transitions to
patterns with larger and larger wavelength. And after a sufficiently low value
of precipitation all of these vegetation patches go extinct and a bare desert
state is attained. This desertification threshold varies with the width
($\sigma$) of kernel $\rho_{1}$, as $\sigma$ rises the threshold value of
precipitation for desertification also rises, for both type of grazing.

In Figure 4 simulations have been performed for $\rho_{2}$ with the same
setting as Figure 3. Like the case for $\rho_{1}$, as the precipitation
reduces, the initial homogeneously vegetated state undergoes a Turing
bifurcation resulting in formation vegetation patterns, and then after a
number of pattern transitions a regime shift to desert state occurs. In case
of sustained grazing this regime shift is to a state with no vegetation, but a
little amount of vegetation continues to exist (the thin strips from $a\approx
2$ to $a\approx 0.54$ in Figure 4) for natural grazing, before going to a
fully degraded state finally. As expected, these observations matches exactly
with the findings of Ref. Siero (2018); Siero  _et al._ (2019).

A qualitative comparison between these kernel choices shows that, unlike the
kernel $\rho_{2}$, in case of Gaussian kernel $\rho_{1}$, patterned vegetated
state goes on till a precipitation level that is very low, in case of
sustained grazing also. Moreover, shift toward a completely degraded state
happens for much higher value of precipitation in case of $\rho_{2}$ compared
to $\rho_{1}$. One more interesting fact is, as the width of kernel $\rho_{1}$
grows, the characteristics of the system become more and more similar to the
system with uniform kernel $\rho_{2}$.

(a)

(b)

Figure 5: Spatially averaged vegetation density vs precipitation corresponding
to the previous simulations, with cut-off Gaussian $\rho_{1}$ having
$\sigma=1$ ( ), $\sigma=8$ ( ), $\sigma=15$ ( ), $\sigma=25$ ( ); and uniform
kernel $\rho_{2}$ ( ). (a) Sustained Grazing and (b) Natural Grazing

In Figure 5, the spatial averages (i.e. biomass per unit area) of the
vegetation biomass are plotted for the aforementioned simulation runs. For
both kernel choices, declining precipitation level results in reduction of
vegetation biomass. This reduction is not gradual, rather we can see sudden
changes in the density curve (Figure 5) whenever the patterns undergoes
wavelength adaptations in Figure 3-Figure 4; and for sufficiently low
precipitation these curves abruptly shift to a state with no vegetation. This
abrupt transition comes much earlier for kernel $\rho_{1}$ that has higher
width. Furthermore, the curves for simulations with cut-off Gaussian with
higher width $\sigma$, coalesces with the curve for $\rho_{2}$.

We have also carried out numerical simulations similar to those in Figure
3-Figure 4 but with a higher rate of precipitation decay
$(\frac{da}{dt}=-10^{-2},\frac{da}{dt}=-10^{-3})$ to better understand how the
system-response is dependent on the rate of variation of precipitation level.
Like in the earlier cases, here too as the precipitation decays, the system
shifts from an initial vegetated state to patterned states and then to a
complete desert state after a number of pattern rearrangements and wavelength
adaptations. However, the system is found to be undergoing wavelength
adaptations with increasing step-size for growing rates of variation in $a$;
i.e., when the rate of change is higher, a much higher proportion of patches
become extinct while undergoing wavelength adaptation compared to the cases in
Figure 3-Figure 4. Due to this, desertification occurs at precipitation levels
where stable patterned states previously existed at earlier instances. This
infers the possibility of the rate-dependent transition, but a detailed study
of this phenomenon is out of the scope of this work.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Simulation results for system (5-6) with kernel $\rho_{1}$ (26)
having different $\sigma$ values: (a-b) $\sigma=1$, (c-d) $\sigma=8$, (e-f)
$\sigma=15$. The left side column is for Sustained Grazing (a-c-e) and the
second column (b-d-f) is for Natural Grazing. Initial vegetation distribution
is taken from Figure 3 at precipitation where very few vegetation patches are
left and then the precipitation level is set to increase at a rate
$\frac{da}{dt}=10^{-4}$

Earlier studies Siero  _et al._ (2019); Sherratt (2013) have demonstrated that
extended Klausmeier model (2) manifests hysteresis phenomena: numerical run
for decreasing precipitation yielded patterns with shorter wavelength than
patterns for increasing precipitation. To know the dependence of system-
response on its history, simulations have been carried out with the same
setting of Figure 3-Figure 4, but this time with precipitation increasing at a
rate $\frac{da}{dt}=10^{-4}$. Following previous studies Siero  _et al._
(2019); Siteur  _et al._ (2014), simulation is performed up to a precipitation
level $a=3.5$. Vegetation distributions of Figure 3-Figure 4, where a very
little amount of vegetation is left (i.e. at the precipitation level just
before vegetation density $n$ becomes lesser than $10^{-6}$), are taken as the
initial conditions. The simulation results (Figure 6) for kernel $\rho_{1}$
having different width, show that with increasing precipitation, the initial
vegetation distribution continues for some time and number of vegetated
patches remains the same. Further increment in precipitation results in a
sequence of transitions to patterns with shorter and shorter wavelength, and
ultimately a regime shift to a state with homogeneous vegetation takes place.
In case of $\sigma=1$, the system is unable to recover fully and no
homogeneously vegetated state appears till $a=3.5$ for both type of grazing.
But as the width $\sigma$ increases, the final transition from patterned state
to uniformly vegetated state happens (Figure 6-Figure 6). This regime shift
occurs for slightly higher level of precipitation for natural grazing, which
reciprocates the earlier observations for decreasing precipitation. Moreover,
this final regime shift occurs at precipitation level which is greater than
the Turing bifurcation threshold $a_{T}$ in Figure 3. Furthermore, as expected
the simulation results for uniform kernel $\rho_{2}$ are in line with Ref.
Siero  _et al._ (2019); not shown here for the sake of brevity. It is observed
that, unlike Ref. Siero  _et al._ (2019), the restoration of vegetation under
increasing precipitation can be seen not only for natural grazing, but in case
of sustained grazing also. Another noteworthy fact is, the wavelength of these
patterns at any given precipitation level is larger with respect to their
counterpart in Figure 3. This can be understood simply by counting the number
of ridges in Figure 6 at any particular value of $a$ and comparing them with
their corresponding result in Figure 3. So at any specific precipitation
level, there are multiple possible stable states at which the system may
reside, depending on its history.

##  IV Discussion and Conclusion

This study aims to evaluate the influence of the intrinsic spatial non-
locality in herbivore grazing over the response of dry-land ecosystem to
changing environmental conditions. While modeling vegetation dynamics in arid
ecosystem, Ref. Siero (2018); Siero  _et al._ (2019) have elegantly
incorporated a mean-density dependent grazing response; here it is further
modified by using a Gaussian weight function for deriving the mean vegetation
density. This modification is driven by the simple assumption that a grazer at
any location will be more influenced by the forage nearby than that of further
away. It has been observed that such a simple ecological aspect of herbivory
effects the dynamics of ecosystem significantly.

Results of this study clearly reconfirm that apart from the water and nutrient
availability, herbivory is also a major player in the functioning of dry-land
ecosystem. Simulations in section III show that the system responds to the
change in precipitation by regulating vegetation biomass through self-
organization of patterns. With increasing aridity, system starting from a
homogeneously vegetated state ($\mathcal{V}$) first shifts to a patterned
state due to Turing bifurcation and then undergoes a sequence of pattern
adaptations to patterns with larger and larger wavelength, followed by a
regime shift to completely barren desert state. In non-spatial models, with
declining precipitation the system was thought to stay in stable uniformly
vegetated state until a tipping point (marked by $*$ in Figure 1) is reached,
where a critical transition to the alternate stable state ($\mathcal{B}$)
happens Rietkerk  _et al._ (2021). However, findings of this study suggest
that the steady state ($\mathcal{V}$) with uniform vegetation cover, loses
stability against spatial effects through Turing bifurcation, at much higher
precipitation level than that of the tipping point and vegetation persists
through rearrangement of patterns for precipitation level beyond the tipping
point, which reiterates the significance of self-organized patterns for
maintaining productivity in dry-land ecosystem Siteur  _et al._ (2014);
Bastiaansen  _et al._ (2020). One of the new phenomena observed in this study
is that the system response to increasing environmental pressure heavily
depends on the characteristic of grazer, for example perceptual abilities
(e.g. sight, olfaction). From Figure 2 it can be noticed that when the width
$\sigma$ of the weight function ($\rho$) is small (i.e. scenario where grazer
relies mainly on the vegetation present in a close vicinity), the shift from
homogeneous vegetation cover to patterned state happens at relatively higher
precipitation level. Moreover, for both sustained and natural grazing,
complete desertification comes at more arid condition for the case with
smaller width. From comparison of Figure 3 and Figure 4 it is evident that the
ultimate critical transition to barren desert state takes place at more dry
condition in case of Gaussian kernel, i.e. the grazing ecosystem is more
resilient to increasing aridity than it was considered to be in previous
studies Siero (2018); Siero  _et al._ (2019). Another interesting fact
observed in this study is that the decline in average vegetation density in
response to the increasing aridity is not gradual, rather sudden oscillations
can be observed in Figure 5 whenever the patterns rearrange themselves by
adapting wavelength. Observations in Ref. Siero  _et al._ (2019) suggest that
in case of human controlled grazing, the final regime shift occurs at a high
precipitation level (Figure 4); however a very small amount of vegetation
continues to exist for far more arid situation in case of natural grazing
(Figure 4). In contrast, this study exhibits resilience to decreasing
precipitation in both type of grazing (Figure 3). Another noteworthy fact is
that as width of the Gaussian influence function increases, the dynamical
behavior of the system become increasingly alike to the system with uniform
kernel $\rho_{2}$ (which resembles the model in Ref. Siero (2018); Siero  _et
al._ (2019)). This can be easily explained by comparing the shape of weight
functions $\rho_{1}$ and $\rho_{2}$. With growing $\sigma$ the bell shape of
Gaussian function gets more and more flattened, hence the weights all over the
range becomes almost equal, which resembles the case for uniform kernel.

Our findings also imply that the response of patterned ecosystems to
environmental variation depends not only on the magnitude of the variation but
also on the rate at which conditions change. Similar type of results had been
observed in Ref. Siteur  _et al._ (2014). When the rate of change in
precipitation level is rapid, the adaptation process is less gradual. Such
circumstances prevent vegetation patches from rearranging and the system
shifts abruptly to a completely degraded state for less arid conditions
compared to the scenarios with a slow rate of change. So to tackle the
increasing environmental pressure efficiently, it is also necessary to
identify the critical rates of change in environmental conditions.

This work demonstrates possibility of restoration in grazing ecosystem, which
is in line with previous studies Siteur  _et al._ (2014); Sherratt (2013).
Model runs with increasing precipitation, exhibit complete recovery to
uniformly vegetated state at higher precipitation level, preceded by a number
of pattern transitions to patterns with shorter and shorter wavelengths.
Unlike Ref. Siero  _et al._ (2019), improving environmental condition results
in restoration of vegetation not only for natural grazing, but also for
sustained grazing case. Furthermore, qualitative comparison between Figure 3
and Figure 6 reveals that although the rate of change in precipitation are
same, the system response to increasing and decreasing precipitation are
contrasting. This infers to the history-dependence of system response under
stressful environment (better known as Hysteresis) and possibility of multi-
stability in the grazing ecosystem, which reciprocates the observations of
earlier studies Siero  _et al._ (2019); Siteur  _et al._ (2014).

Summarizing, this study reveals how the incorporation of herbivore-grazing as
function of distance from the grazer, can substantially alter the ecosystem
response to changing environmental condition. Apart from providing novel
insights into grazing ecosystem, this study also reconfirms the findings of
recent studies in arid ecosystem. Moreover, outcomes of this work necessitates
inclusion of spatial non-locality while modeling herbivory in vegetation
systems. However, it must be mentioned that this study is based on a
phenomenological approach, intended to integrate the scale-dependent feedback
and spatially non-local herbivore grazing. It can further provide theoretical
framework for future data acquisition on herbivore foraging using
sophisticated techniques and satellite images. Assimilating the biological
behavior-driven movement of herbivores with realistic data support while
choosing weight function would be a research problem worth pursuing.
Furthermore, in order to understand the transitions among qualitatively
different patterns in the presence of herbivore grazing, the two dimensional
version of this model needs to be explored.

###### Acknowledgements.

This work is supported by UGC, Govt. of India under NET-JRF Scheme. We are
grateful to Jasashwi Mandal and Swarnendu Banerjee for providing technical
support. We would also like to thank anonymous reviewers for their insightful
suggestions.

## References

  * Reynolds  _et al._ (2007) J. F. Reynolds, D. M. S. Smith, E. F. Lambin, B. L. Turner, M. Mortimore, S. P. J. Batterbury, T. E. Downing, H. Dowlatabadi, R. J. Fernández, J. E. Herrick, _et al._ , science 316, 847 (2007).
  * Gowda  _et al._ (2018) K. Gowda, S. Iams, and M. Silber, Scientific reports 8, 1 (2018).
  * von Hardenberg  _et al._ (2001) J. von Hardenberg, E. Meron, M. Shachak, and Y. Zarmi, Physical Review Letters 87, 198101 (2001).
  * Macfadyen (1950) W. A. Macfadyen, The Geographical Journal 116, 199 (1950).
  * Deblauwe  _et al._ (2012) V. Deblauwe, P. Couteron, J. Bogaert, and N. Barbier, Ecological monographs 82, 3 (2012).
  * Tlidi  _et al._ (2018) M. Tlidi, M. Clerc, D. Escaff, P. Couteron, M. Messaoudi, M. Khaffou, and A. Makhoute, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20180026 (2018).
  * Fernandez-Oto  _et al._ (2019a) C. Fernandez-Oto, D. Escaff, and J. Cisternas, Ecological Complexity 37, 38 (2019a).
  * Deblauwe  _et al._ (2008) V. Deblauwe, N. Barbier, P. Couteron, O. Lejeune, and J. Bogaert, Global Ecology and Biogeography 17, 715 (2008).
  * Martínez-García  _et al._ (2014) R. Martínez-García, J. M. Calabrese, E. Hernández-García, and C. López, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20140068 (2014).
  * Fernandez-Oto  _et al._ (2014) C. Fernandez-Oto, M. Tlidi, D. Escaff, and M. Clerc, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20140009 (2014).
  * Escaff  _et al._ (2015) D. Escaff, C. Fernandez-Oto, M. Clerc, and M. Tlidi, Physical Review E 91, 022924 (2015).
  * Fernandez-Oto  _et al._ (2019b) C. Fernandez-Oto, O. Tzuk, and E. Meron, Physical review letters 122, 048101 (2019b).
  * Parra-Rivas and Fernandez-Oto (2020) P. Parra-Rivas and C. Fernandez-Oto, Physical Review E 101, 052214 (2020).
  * Clerc  _et al._ (2021) M. G. Clerc, S. Echeverría-Alar, and M. Tlidi, Scientific reports 11, 1 (2021).
  * Dawes and Williams (2016) J. Dawes and J. Williams, Journal of mathematical biology 73, 63 (2016).
  * Mau  _et al._ (2015) Y. Mau, L. Haim, and E. Meron, Physical Review E 91, 012903 (2015).
  * Bastiaansen  _et al._ (2018) R. Bastiaansen, O. Jaïbi, V. Deblauwe, M. B. Eppinga, K. Siteur, E. Siero, S. Mermoz, A. Bouvet, A. Doelman, and M. Rietkerk, Proceedings of the National Academy of Sciences 115, 11256 (2018).
  * Guttal and Jayaprakash (2007) V. Guttal and C. Jayaprakash, Journal of theoretical biology 248, 490 (2007).
  * Zelnik  _et al._ (2013) Y. R. Zelnik, S. Kinast, H. Yizhaq, G. Bel, and E. Meron, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371, 20120358 (2013).
  * HilleRisLambers  _et al._ (2001) R. HilleRisLambers, M. Rietkerk, F. van den Bosch, H. H. Prins, and H. de Kroon, Ecology 82, 50 (2001).
  * Gowda  _et al._ (2014) K. Gowda, H. Riecke, and M. Silber, Physical Review E 89, 022701 (2014).
  * Gilad  _et al._ (2004) E. Gilad, J. von Hardenberg, A. Provenzale, M. Shachak, and E. Meron, Physical Review Letters 93, 098105 (2004).
  * Meron  _et al._ (2007) E. Meron, H. Yizhaq, and E. Gilad, Chaos: An Interdisciplinary Journal of Nonlinear Science 17, 037109 (2007).
  * Lefever and Lejeune (1997) R. Lefever and O. Lejeune, Bulletin of Mathematical biology 59, 263 (1997).
  * Tlidi  _et al._ (2020) M. Tlidi, E. Berríos-Caro, D. Pinto-Ramo, A. Vladimirov, and M. G. Clerc, Physica D: Nonlinear Phenomena 414, 132708 (2020).
  * Gilad  _et al._ (2007) E. Gilad, J. von Hardenberg, A. Provenzale, M. Shachak, and E. Meron, Journal of Theoretical Biology 244, 680 (2007).
  * Messaoudi  _et al._ (2020) M. Messaoudi, M. G. Clerc, E. Berríos-Caro, D. Pinto-Ramos, M. Khaffou, A. Makhoute, and M. Tlidi, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 093136 (2020).
  * Klausmeier (1999) C. A. Klausmeier, Science 284, 1826 (1999).
  * Eigentler and Sherratt (2018) L. Eigentler and J. A. Sherratt, Journal of mathematical biology 77, 739 (2018).
  * Eigentler and Sherratt (2019) L. Eigentler and J. A. Sherratt, Bulletin of Mathematical Biology 81, 2290 (2019).
  * Sun  _et al._ (2021) G. Q. Sun, H. T. Zhang, J. S. Wang, J. Li, Y. Wang, L. Li, Y. P. Wu, G. L. Feng, and Z. Jin, Nonlinear Dynamics 104, 1677 (2021).
  * Xue  _et al._ (2020) Q. Xue, G.-Q. Sun, C. Liu, Z.-G. Guo, Z. Jin, Y.-P. Wu, and G.-L. Feng, Nonlinear Dynamics 99, 3407 (2020).
  * Kabir and Gani (2022) M. H. Kabir and M. O. Gani, Journal of Theoretical Biology 536, 110997 (2022).
  * Hobbs (1996) N. T. Hobbs, The Journal of Wildlife Management , 695 (1996).
  * Noy-Meir (1975) I. Noy-Meir, The Journal of Ecology , 459 (1975).
  * Mueller  _et al._ (2014) E. N. Mueller, J. Wainwright, A. J. Parsons, and L. Turnbull,  _Patterns of land degradation in drylands_ (Springer, 2014).
  * Rietkerk and van de Koppel (1997) M. Rietkerk and J. van de Koppel, Oikos , 69 (1997).
  * Focardi  _et al._ (1996) S. Focardi, P. Marcellini, and P. Montanaro, Journal of Animal Ecology , 606 (1996).
  * Siero (2018) E. Siero, Journal of theoretical biology 436, 64 (2018).
  * Siero  _et al._ (2019) E. Siero, K. Siteur, A. Doelman, J. van de Koppel, M. Rietkerk, and M. B. Eppinga, The American Naturalist 193, 472 (2019).
  * Gordon and Prins (2008) I. J. Gordon and H. H. T. Prins,  _The ecology of browsing and grazing_ , 195 (Springer, 2008).
  * Worton (1989) B. J. Worton, Ecology 70, 164 (1989).
  * Roese  _et al._ (1991) J. H. Roese, K. L. Risenhoover, and L. J. Folse, Ecological modelling 57, 133 (1991).
  * Kealy and Wollkind (2012) B. J. Kealy and D. J. Wollkind, Bulletin of mathematical biology 74, 803 (2012).
  * Siteur  _et al._ (2014) K. Siteur, E. Siero, M. B. Eppinga, J. D. Rademacher, A. Doelman, and M. Rietkerk, Ecological Complexity 20, 81 (2014).
  * Rietkerk  _et al._ (2002) M. Rietkerk, M. C. Boerlijst, F. van Langevelde, R. HilleRisLambers, J. van de Koppel, L. Kumar, H. H. Prins, and A. M. de Roos, The American Naturalist 160, 524 (2002).
  * Fryxell  _et al._ (2004) J. M. Fryxell, J. F. Wilmshurst, and A. R. Sinclair, Ecology 85, 2429 (2004).
  * Sherratt and Lord (2007) J. A. Sherratt and G. J. Lord, Theoretical population biology 71, 1 (2007).
  * Sherratt (2010) J. A. Sherratt, Nonlinearity 23, 2657 (2010).
  * Sherratt (2005) J. A. Sherratt, Journal of mathematical biology 51, 183 (2005).
  * Fuentes  _et al._ (2004) M. A. Fuentes, M. N. Kuperman, and V. M. Kenkre, The Journal of Physical Chemistry B 108, 10505 (2004).
  * Siero  _et al._ (2015) E. Siero, A. Doelman, M. B. Eppinga, J. D. Rademacher, M. Rietkerk, and K. Siteur, Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 036411 (2015).
  * Sherratt (2013) J. A. Sherratt, Ecological Complexity 14, 8 (2013).
  * Pal and Poria (2021) M. K. Pal and S. Poria, Chaos: An Interdisciplinary Journal of Nonlinear Science 31, 053108 (2021).
  * Pal  _et al._ (2020) S. Pal, M. Banerjee, and S. Ghorai, Applied Mathematical Modelling 79, 809 (2020).
  * Rietkerk  _et al._ (2021) M. Rietkerk, R. Bastiaansen, S. Banerjee, J. van de Koppel, M. Baudena, and A. Doelman, Science 374, eabj0359 (2021).
  * Bastiaansen  _et al._ (2020) R. Bastiaansen, A. Doelman, M. B. Eppinga, and M. Rietkerk, Ecology letters 23, 414 (2020).