category: literaturenote citekey: grifofarequilibriumtravellingpulses2024 title: Far-from-equilibrium travelling pulses in sloped semi-arid environments driven by autotoxicity effects authors: "Grifò, Gabriele; Iuorio, Annalisa; Veerman, Frits" year: 2024 date: 2024-00-00 2024 doi: 10.1137/24M1669499 url: "https://arxiv.org/abs/2405.15602v1" zotero_key: I8ZJFIWD zotero_storage: 85WE3XA4 collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Grifò, Gabriele"

status: converted

Far-from-equilibrium travelling pulses in sloped semi-arid environments driven by autotoxicity effects

Gabriele Grif`o¹,² , Annalisa Iuorio³,4,^∗ , Frits Veerman⁵ ¹Department of Mathematical, Computer, Physical and Earth Sciences, University of Messina, V.le F. Stagno D'Alcontres 31, I-98166 Messina, Italy. ²Department of Engineering, University of Messina, C.da di Dio, I-98166 Messina, Italy ³Department of Engineering, Parthenope University of Naples,

Centro Direzionale - Isola C4, 80143 Naples, Italy ⁴Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria ⁵Mathematical Institute, University of Leiden, Einsteinweg 55, 2333 CC Leiden, The Netherlands

(∗) Corresponding author: annalisa.iuorio@uniparthenope.it

Abstract

In this work, an extension of the 1D Klausmeier model that accounts for the toxicity compounds is considered and the occurrence of travelling stripes is investigated. Numerical simulations are firstly conducted to capture the qualitative behaviours of the pulse–type solutions and, then, geometric singular perturbation theory is used to prove the existence of such travelling pulses by constructing the corresponding homoclinic orbits in the associated 4-dimensional system. A scaling analysis on the investigated model is performed to identify the asymptotic scaling regime in which travelling pulses can be constructed. Biological observations are extracted from the analytical results and the role of autotoxicity in travelling patterns is emphasized. Finally, the analytically constructed solutions are compared with the numerical ones, leading to a good agreement that confirms the validity of the conducted analysis. Numerical investigations are also carried out in order to gain additional information on vegetation dynamics.

Keywords: Travelling pulses, vegetation patterns, geometric singular perturbation theory, autotoxicity, reaction–diffusion–advection models

1 Introduction

One of the major challenges in the context of climate change is desertification, which poses a danger for the durability of dryland areas [1]. To predict and mitigate desertification catastrophes, considerable effort has been spent to study ecosystem response to rapid change, and to identify indicators that provide reliable information on ecosystem health and robustness [2, 3]. In this context, the occurrence and dynamics of vegetation patterns have been increasingly recognised as a key predictive indicator [4–7]. Recent studies [8, 9] show how vegetation patterns can be interpreted as a signature of ecosystem resilience, providing a gradual (in contrast to catastrophic) sequence of patterned states towards the desert one.

In arid and semi-arid environments, the self-organisation of biomass can mainly be attributed to local positive feedback between water and biomass [10]. However, vegetation patterns are also observed in non-arid (non water-limited) ecosystems, suggesting the presence of alternative mechanisms that may contribute to the emergence of spatial vegetation patterns. One of these mechanisms, called autotoxicity, is an example of plant-soil negative feedback that arises from the presence of soil-borne pathogens, changing the composition of soil microbial communities, and/or the accumulation of autotoxic compounds from decomposing plant litter [11, 12].

Although the emergence of stationary and oscillatory vegetation patterns is believed to be ubiquitous, these patterning phenomena tend to occur over long timescales and far away from populated areas, making on site experiments and detailed observations challenging to perform. Therefore, mathematical modelling plays a crucial role in analysing these pattern formation processes, providing analytical tools to assess an ecosystem's response to stressors. There is an extensive body of mathematical literature on reproducing spatially periodic vegetation patterns and predicting their spatial evolution, using various models [10, 13–24]. Among these models, the Klausmeier model [13] is one of the most simple systems still able to reproduce the occurrence of vegetation patterns as migrating bands. The Klausmeier model is a two-compartment reaction–advection–diffusion model describing water and biomass evolution in arid and semi-arid enviroments. However, this model does not take other processes into account that may influence vegetation dynamics. Therefore, several extensions to the Klausmeier model have been proposed, that include: (i) the emergence of stationary patterns on flat terrains [25–30], (ii) the influence of inertial effects [28–33], (iii) the presence of secondary seed dispersal [32–35], (iv) the occurrence of toxicity compounds [36–39], and (v) a finite soil carrying capacity [40–43].

To the best of our knowledge, the existing literature on water-biomass-toxicity systems focuses on flat semi-arid environments [36 –39]. We aim to eludicate the influence of autotoxicity on migrating vegetation bands in sloped semi-arid terrains, by means of a generalisation of the Klausmeier model. We study far–from–equilibrium solutions using geometric singular perturbation theory, in particular solutions that are bi-asymptotic to the desert state. We also present the results of numerical simulations that reflect our analytical results, and that allow us to extract additional information on the influence of autotoxicity on travelling vegetation pulses.

This paper is organised as follows. In Section 2, we introduce an extension of the original Klausmeier model, adding a toxicity component. In the same section, we present the results of numerical simulations that show the emergence of migrating vegetation bands. In Section 3, a scaling analysis is performed to identify the asymptotic scaling regime in which travelling vegetation bands can be found. We construct a pulse–type solution in the singular limit, and prove its persistence. In Section 4, we use numerical continuation methods to study the behaviour of the constructed solution for a variety of parameter ranges, to illustrate and extend the asymptotic analysis in Section 3. Future directions and concluding remarks are presented in Section 5.

2 The model

Following existing literature in which the influence of toxic compounds on vegetation dynamics is studied [18, 36 –38], we formulate a generalisation of the two-compartment Klausmeier model [13], incorporating the effect of autotoxicity. In particular, we add a ordinary differential equation (ODE) that describes the temporal evolution of a toxic compound. Additional kinetic terms describe the interaction between the biomass and the toxic compound. Moreover, assuming a ramp-like topography, we reduce the spatial dimension of the model to focus on the principal direction of the slope. The resulting model reads

$$\begin{split} W_t &= p - rB^2W - lW + \nu W_x, \ B_t &= cB^2W - (d + sT)B + DBB{xx}, \ T_t &= q(d + sT)B - (k + wp)T, \end{split} \tag{1}$$

in which W(x, t), B(x, t), and T(x, t) represent the surface water, biomass, and autotoxicity densities at location x ∈ Ω ⊂ R (the positive x-direction being uphill) and time t ∈ R +, respectively.

The motion of surface water W is dominated by anisotropic transport, modelled by an advection term with speed ν > 0, mimicking downhill water flow. The biomass B is assumed to diffuse isotropically with diffusion coefficient DB. In addition, the biomass density increases due to water availability with growth rate c and decreases both due to intrinsic mortality with rate d, and due to the presence of toxic compounds with sensitivity s. The surface water density is fed by precipitation p (assumed to be constant over the intrinsic time scale of our model), and decreases both due to evaporation or drainage with rate l and due to water uptake by biomass with rate r. The additional toxicity component T decreases due to natural decay with rate k, and is washed out by precipitation proportional to w. The toxicity increases due to the decomposition of biomass, proportional to q. For a detailed description of all the parameters of system (1), we refer to Table 1.

We nondimensionalise model (1) by applying the rescaling

$$\tilde{x} = \sqrt{\frac{l}{D_B}}x, \quad \tilde{t} = lt, \quad W = \frac{\sqrt{lr}}{c}U, \quad B = \sqrt{\frac{l}{r}}V, \quad T = \frac{ql\sqrt{l}}{(k+pw)\sqrt{r}}S.$$ (2)

The associated nondimensional parameters are defined as

Parameter	Description	Units
c	Growth rate of B due to water uptake	m4d −1kg−2
d	Death rate of biomass B	−1 d
k	Decay rate of toxicity T	−1 d
l	Water loss due to evaporation or drainage	d−1
p	Precipitation rate	kg d−1 m−2
q	Proportion of toxins in dead biomass	-
r	Rate of water uptake	m4d −1kg−2
s	Sensitivity of plants to toxicity T	m2d −1kg−1
w	Washing out of toxins by precipitation	m2kg−1
DB	Diffusion coefficient of biomass B	m2d −1
ν	Water advection speed	m d−1

Table 1: An overview of the parameters of model (1).

$$\mathcal{A} = \frac{cp}{l\sqrt{lr}}, \qquad \mathcal{B} = \frac{d}{l}, \qquad \mathcal{D} = \frac{l}{k+pw},$$

$$\mathcal{H} = \frac{sq\sqrt{l}}{(k+pw)\sqrt{r}}, \qquad \varepsilon = \frac{\sqrt{D_B l}}{\nu}.$$ (3)

Here, A is proportional to the precipitation rate p and hence can be interpreted (by fixing all other dimensional parameters) as its nondimensional counterpart. The parameter B (resp. D) represents the ratio between the (linear) biomass and water (resp. water and toxicity) loss. Analogous to A, the parameter H is directly proportional to the biomass sensitivity to autotoxicity s; therefore, it can be interpreted (considering all other parameters fixed) as the nondimensional version of this parameter. The parameter ε quantifies the ratio between the (normalised) diffusion speed of the biomass and the downhill advection speed of the surface water. As the downhill flow of water is usually significantly faster than the spread of biomass due to growth, we assume that ε is asymptotically small, i.e. 0 < ε ≪ 1. This parameter will play a pivotal role in the analysis in Section 3.

Using the rescaling (2) and the nondimensional parameters (3), system (1) becomes (dropping the tilde notation for the dependent variables ˜x and t˜)

$$U_t = \mathcal{A} - U - V^2 U + \varepsilon^{-1} U_x,$$

$$Vt = V^2 U - \mathcal{B}V - \mathcal{H}SV + V{xx},$$

$$\mathcal{D}S_t = \mathcal{B}V + \mathcal{H}SV - S.$$ (4)

Note that the original Klausmeier model [13] is a special case of model (4): taking H = 0 fully decouples the autotoxicity equation.

2.1 Preliminary observations and main result

The influence of autotoxicity can already be observed in the spatially homogeneous steady states of system (4). In addition to the trivial desert state

$$(U*, V, S_) = (\mathcal{A}, 0, 0),$$ (5)

system (4) admits two nontrivial vegetated states

$$(U{\pm}, V{\pm}, S{\pm}) = \left(\frac{\mathcal{A}}{1 + V{\pm}^2}, \frac{\mathcal{A} \pm \sqrt{\mathcal{A}^2 - 4\mathcal{B}(\mathcal{B} + \mathcal{A}\mathcal{H})}}{2(\mathcal{B} + \mathcal{A}\mathcal{H})}, \frac{\mathcal{B}V{\pm}}{1 - \mathcal{H}V{\pm}}\right)$$ (6)

which explicitly depend on the value of the toxicity coupling parameter H. In particular, these nontrivial spatially homogeneous steady states exist if and only if A > 2B H + √ 1 + H² , which provides a nontrivial condition on the original 'Klausmeier' parameters A and B. This condition reduces to the non-toxic Klausmeier condition A > 2B in the limit H = 0.

Existing theoretical and experimental works [44–47] show that vegetation patterns on sloped semi-arid terrain are predominantly characterised by pulse or wave train solutions that originate from the desert state. From a mathematical viewpoint, these patterns are determined by far–from–equilibrium dynamics that are beyond the reach of classical tools such as linear or weakly nonlinear stability analysis [48–50].

Direct numerical simulation of model (4) confirms this observation. Using Matlab® [51], we integrate (4) over a spatial domain of length L = 1000 with periodic boundary conditions, with constant initial data for both water U and autotoxicity S, i.e. U(x, 0) = 0.5 and S(x, 0) = 0 ∀x ∈ [0, L]. For the vegetation variable V , we take as initial data a Gaussian pulse centered in x = 300 with standard deviation σ = 2 5 and amplitude Vmax = 10. The model parameters are fixed (in accordance with [13, 38, 46]) as A = 1.2, B = 0.45, ε = 0.005, D = 4.5 and H = 1. In Figure 1(d), the space–time evolution of the biomass density clearly indicates a narrow peak travelling at constant speed through the domain. Figure 1 a)–c) depicts the spatial profile of this travelling pulse, highlighting the nonlinear and spatially multi-scale nature of this solution. In particular, the component profiles suggest a hierarchy of three asymptotically distinct spatial scales, which we denote as 'superslow', 'slow' and 'fast'. Note that this terminology is slightly misleading, as these scales are not connected to the temporal evolution of the pulse profile; rather, they indicate how the pulse profile changes as a function of x (for fixed t). Using these terms, V is observed to change on the fast spatial scale; the profile of S displays both slow and fast changes, and the profile of U encorporates both superslow and fast behaviour. This spatial scale separation will play a crucial role in the analysis in Section 3.

Inspired by the observations depicted in Figure 1, we investigate in this paper the existence of travelling pulses with a fixed wave speed C > 0. The existence of far–from–equilibrium pulse–type solutions has been studied in previous works, both in the context of specific models [40, 46, 52] and general classes of reaction-diffusion systems [53, 54], using geometric singular perturbation theory (GSPT). We follow an analogous approach for our model (4) that includes autotoxicity. As such, our work can be compared directly to [46], where GSPT was used to establish the existence of travelling pulses in (4) in the absence of autotoxicity, i.e. for H = 0. In comparison with [46], we show that autotoxicity not only introduces novel behaviour –which is, from a general viewpoint, to be expected when considering a model of higher complexity– but also simplifies the asymptotic analysis by avoiding the need to use geometric blow-up techniques, as was necessary in [46].

Our main analytical result is presented in the following Theorem:

Theorem 1. For 0 < ε ≪ 1 sufficiently small, there exists a unique θ⁰ > 0 such that there exists a travelling wave solution (U, V, S) (x, t) = (U, V, S) (x − Ct) of system (4) with wave speed C = A2θ 0 ε 1/³ + O(1).

The constant θ⁰ coincides with the one provided in [46, Theorem 1.1], i.e. θ⁰ ≈ 0.8615.

Figure 1: Spatial profiles for surface water U (a), biomass V (b) and toxicity S (c) of a travelling pulse, obtained by numerical integration of system (4) over a domain with length L = 1000 with periodic boundary conditions. The parameter values are A = 1.2, B = 0.45, ε = 0.005, D = 4.5, and H = 1.

3 Existence of travelling pulses

In this section, we apply GSPT to prove the existence of a constant speed travelling pulse solution to (4). Introducing a comoving frame variable ξ = x−Ct, system (4) can be rephrased as

$$\left(\mathcal{C} + \frac{1}{\varepsilon}\right) U_{\xi} + \mathcal{A} - U - UV^{2} = 0,$$

$$V{\xi\xi} + \mathcal{C}V{\xi} - \mathcal{B}V + UV^{2} - \mathcal{H}SV = 0,$$

$$\mathcal{C}\mathcal{D}S_{\xi} + \mathcal{B}V + \mathcal{H}SV - S = 0.$$ (7)

System (7) is equivalent to the first-order system

$$U_{\xi} = \frac{\varepsilon}{1 + \varepsilon C} \left( U - A + UV^{2} \right),$$

$$V_{\xi} = Q,$$

$$Q_{\xi} = BV - UV^{2} + \mathcal{H}SV - CQ,$$

$$\mathcal{CDS}_{\xi} = S - \mathcal{B}V - \mathcal{H}SV.$$ (8)

Inspired by [46] and by the different spatial scales visible in the pulse profile shown in Figure 1, we introduce the following asymptotic rescaling:

$$V = c^{-1}\varepsilon^{-2/3}v, Q = \varepsilon^{-1}q, \xi = c^{2}\varepsilon^{1/3}\tau,$$

$$C = \varepsilon^{-1/3}c, S = c^{-1}\varepsilon^{-2/3}s, \delta = \varepsilon^{2/3}c.$$ (9)

In these rescaled variables, system (8) becomes

$$\dot{u} = \frac{1}{1+\delta} \left( uv^2 + \delta^2 \left( u - \mathcal{A} \right) \right), \dot{v} = c^3 q, \dot{q} = \delta \mathcal{B} v - uv^2 + \mathcal{H} sv - c^3 q, \mathcal{D} \dot{s} = \delta s - \delta \mathcal{B} v - \mathcal{H} sv,$$ (10)

where ˙ = ^d dτ . Note that δ has taken over the role of ε as the asymptotically small parameter. Single stripe patterns travelling uphill are given by orbits of (10) that are homoclinic to the equilibrium (u, v, q, s) = (A, 0, 0, 0), which corresponds to the spatially homogeneous desert state (U∗, V∗, S∗) in (5).

A priori, it is not immediately clear that system (10) admits three (rather than two) asymptotic scales. Introducing the new variable

$$w = (1+\delta)u + q + v + \mathcal{D}s \tag{11}$$

will aid the separation of scales in the upcoming analysis. Using w in favour of u, system (10) is rephrased as

$$\dot{w} = \delta s + \frac{\delta^2}{1+\delta} (w - q - v - \mathcal{D}s - a),$$

$$\dot{v} = c^3 q,$$

$$\dot{q} = \delta \mathcal{B}v - \frac{1}{1+\delta} (w - q - v - \mathcal{D}s) v^2 + \mathcal{H}sv - c^3 q,$$

$$\mathcal{D}\dot{s} = \delta s - \delta \mathcal{B}v - \mathcal{H}sv,$$ (12)

where a = (1 + δ) A. In the following, we will refer to system (12) as the fast system. Our goal is to construct a homoclinic orbit to the equilibrium (w, v, q, s) = (a, 0, 0, 0).

The fast-slow structure of system (12) induced by the presence of the small parameter 0 < δ ≪ 1 allows us to construct the solution described in Theorem 1 as perturbations of singular orbits obtained by matching portions of critical manifolds (i.e. manifolds of equilibria of system (12) for δ = 0) with fast jumps along heteroclinic orbits arising in the corresponding layer problem. The reformulation of system (12) with respect to the slow and the superslow scale is described in Section 3.1, whereas the construction of singular homoclinic orbits is illustrated in Section 3.3. The persistence of singular orbits for 0 < δ ≪ 1 stated in Theorem 1 is proved in Section 3.4. Our theoretical findings are illustrated and extended by numerical results obtained via the continuation software AUTO [55] in Section 4.

3.1 Critical manifolds and reduced dynamics

We take δ = 0 in the fast system (12) to obtain the layer problem

$$\dot{w} = 0, \tag{13a}$$

$$\dot{v} = c^3 q,\tag{13b}$$

$$\dot{q} = -\left(w - q - v - \mathcal{D}s\right)v^2 + \mathcal{H}sv - c^3q,\tag{13c}$$

$$\mathcal{D}\dot{s} = -\mathcal{H}sv. \tag{13d}$$

The equilibria of the layer problem (13) determine the critical manifolds

$$\mathcal{M}^{(1)} = { v = 0, q = 0, s = 0 }, \tag{14a}$$

$$\mathcal{M}^{(2)} = { v = w, q = 0, s = 0 }, \tag{14b}$$

$$\mathcal{M}^{(3)} = { v = 0, q = 0 }. \tag{14c}$$

Geometrically, the one-dimensional critical manifolds M(1) and M(2) are lines, and M(1) is contained in the two-dimensional hyperplane M(3). For future reference, we denote the equilibria of the layer problem (13) as

$$p_1(w) = (w, 0, 0, 0), \quad p_2(w) = (w, w, 0, 0), \quad p_3(w, s) = (w, 0, 0, s),$$ (15)

allowing us to parametrise the critical manifolds as M(1) = S w∈R p1(w), M(2) = S w∈R p2(w) and M(3) = S (w,s)∈R² p3(w, s). We note that M(1) transversally intersects M(2) (and therefore M(3)) at the origin. The critical manifold M(3) is normally hyperbolic, since the equilibrium p3(w, s) is a saddle in the normal (v, q)-directions. The critical manifold M(2) is normally hyperbolic away from the origin, with eigenvalues −c 3 , − Hw D and w 2 in the normal directions. Since M(1) is fully contained in M(3), its persistence is related to the dynamics on M(3), which we will turn to next.

We observe that the hyperplane M(3) (14) is invariant under the flow of (12). In addition, M(1) ⊂ M(3) is itself also invariant under the flow of (12). In contrast to the other two critical manifolds, the critical manifold M(2) is not invariant under the flow of the full system (12). The flow on M(2) is to leading order given by

$$\dot{q} = \delta \mathcal{B} w,$$ (16a)

$$\mathcal{D}\dot{s} = -\delta \mathcal{B}w. \tag{16b}$$

Since s = 0 on M(2), the reduced dynamics here would be biologically unfeasible (s decreases and hence becomes negative). The critical manifold M(2) is therefore not going to play a role in the construction of the travelling pulse solutions (see Figure 2). Hence, we focus our

attention on the dynamics on M(3), which contain the dynamics on M(1) .

To study the reduced flow on the critical manifold M(3), we introduce σ = δτ . In this slow variable, system (12) reads

$$w_{\sigma} = s + \frac{\delta}{1+\delta} \left( w - q - v - \mathcal{D}s - a \right), \tag{17a}$$

$$\delta v_{\sigma} = c^3 q,\tag{17b}$$

$$\delta q_{\sigma} = \delta \mathcal{B}v - \frac{1}{1+\delta} \left( w - q - v - \mathcal{D}s \right) v^2 + \mathcal{H}sv - c^3 q, \tag{17c}$$

$$\delta \mathcal{D}s_{\sigma} = \delta s - \delta \mathcal{B}v - \mathcal{H}sv. \tag{17d}$$

The flow on M(3) is therefore given by

$$w_{\sigma} = s + \frac{\delta}{1+\delta} \left( w - \mathcal{D}s - a \right), \tag{18a}$$

$$\mathcal{D}s_{\sigma} = s. \tag{18b}$$

System (18) is itself a (linear) slow-fast system, which can be brought in Fenichel normal form by observing that

$$\frac{\mathrm{d}}{\mathrm{d}\sigma}(w - \mathcal{D}s) = \frac{\delta}{1+\delta}((w - \mathcal{D}s) - a). \tag{19}$$

First, we observe that the set { s = 0 }, representing M(1), is invariant under the flow of (18). At the same time, { s = 0 } is the critical manifold of the slow-fast system (18); the flow on M(1) is given by

$$w_{\zeta} = \frac{1}{1+\delta} \left( w - a \right) \tag{20}$$

in terms of the superslow variable ζ := δσ = δ 2 τ . The flow (20) has a single unstable equilibrium at w = a, which corresponds to the global desert state.

Second, the linear coordinate transform (19) inspires us to foliate the critical manifold M(3) by one-dimensional manifolds

$$\mathcal{M}^{(3),w_0} := \left{ p_3 \left( w, \frac{w - w_0}{\mathcal{D}} \right) : w_0 \in \mathbb{R}_0^+ \right}, \tag{21}$$

cf. (15). Note that M(3),w⁰ are precisely the (unstable) Fenichel fibres for the flow (18); the flow on M(3),w⁰ is independent of w0, and given by Ds^σ = s. The critical manifolds and the reduced flow on these critical manifolds is shown in Figure 2.

Figure 2: Slow (blue) and superslow (purple) dynamics on the critical manifolds M(i) , i = 1, 2, 3 (14) (v, s, w)-space for q = 0.

3.2 Layer problem

We recall that our sought-after travelling pulses correspond to solutions homoclinic to the equilibrium (w, v, q, s) = (a, 0, 0, 0). From the analysis carried out in Section 3.1, we see that the equilibrium (a, 0, 0, 0) lies on the invariant manifold M(3); however, the dynamics on M(3) are linear (18). Therefore, M(3) cannot contain a homoclinic orbit, and we conclude that any orbit homoclinic to (a, 0, 0, 0) must necessarily contain fast segments, indicating an excursion away from M(3)

In the singular limit, the flow on M(3) (18) is decomposed into a first superslow segment on M(1) from p1(a) = (a, 0, 0, 0) to a point (w ∗ , 0, 0, 0), followed by a slow transition on M(3),w^∗ from (w ∗ , 0, 0, 0) to p3(a, s^∗ ) = (a, 0, 0, s^∗ ) for particular choices of w ^∗ and s ∗ , see Figure 3. The fast transition closing the loop would be a heteroclinic orbit in the layer problem (13) from the equilibrium p3(a, s^∗ ) to p1(a).

Our goal therefore is to determine the value of the constants w ^∗ and s ^∗ by investigating the layer problem, in particular by ensuring that there exists a unique fast connection from p3(a, s^∗ ) to p1(a). A priori, obtaining this heteroclinic orbit in the three-dimensional layer problem (13) would be prohibitive; however, we will exploit the structure of the layer problem to show that the target heteroclinic lies close to the union of two heteroclinic connections to and from p2(a), see Figure 4. The existence of the pair of heteroclinic connections to and from p2(a) is determined in Lemma 2 and Lemma 3.

Lemma 2. There exists a unique c = c ^∗ > 0 for which system (13) admits a heteroclinic orbit from p2(a) = (a, a, 0, 0) to p1(a) = (a, 0, 0, 0). This heteroclinic orbit lies in the hyperplane (w, v, q, s) : w = a, s = 0 . In particular,

$$c^* = (a\,\theta_0)^{2/3} \tag{22}$$

with θ⁰ ≈ 0.8615.

Proof. The hyperplane { s = 0} is invariant under the flow of system (13). On this hyperplane, the layer problem reduces to

$$\dot{v} = c^3 q, \dot{q} = -(w - q - v) v^2 - c^3 q,$$ (23)

which coincides with the layer problem studied in [46, Proposition 2.2]. The existence of a unique value of c = c ^∗ = (a θ0) ²/³ with θ⁰ ≈ 0.8615 such that there is a heteroclinic orbit connecting the two equilibria p1(a) and p2(a) hence directly follows.

Lemma 3. There exists a unique s = s ^∗ = a D for which system (13) admits a heteroclinic orbit from p3(a, s^∗ ) = (a, 0, 0, s^∗ ) to p2(a) = (a, a, 0, 0). This heteroclinic orbit lies in the hyperplane (w, v, q, s) : w = a, w − v − q − Ds = 0 .

Proof. We fix w = a, and observe that the line S s p3(a, s) = M(3) ∩ {w = a} is invariant under the flow of the layer problem (13). For fixed s = ¯s, p3(a, s¯) is an equilibrium of the layer problem (13) with eigenvalues

$$\lambda_3^c(a,\bar{s}) = 0 \quad \text{and} \quad \lambda_3^{s,u}(a,\bar{s}) = -\frac{1}{2}c^{3/2}\left(-c^{3/2} \pm \sqrt{c^3 + 4\mathcal{H}\bar{s}}\right).$$ (24)

The corresponding centre, stable, and unstable eigenvectors are given by

$$\boldsymbol{\eta}^c = (0, 0, 1), \qquad \boldsymbol{\eta}^{s,u} = \left(-\lambda{s,u}, -\mathcal{H}\bar{s} + \lambda{s,u}, \frac{\mathcal{H}\bar{s}}{\mathcal{D}}\right).$$ (25)

Figure 3: Schematic representation of the slow (blue) and superslow (purple) dynamics of the travelling pulse. The dynamics take place on M(3) (18): the superslow segment is contained in M(1) (20), and the slow segment is contained in M(3),w⁰ (21). See also Figure 2.

Figure 4: Schematic representation of the fast dynamics associated to a travelling pulse solution to (12), in (v, q, s)-space for fixed w = a. The existence of the heteroclinic connection from p3(a, s^∗ ) to p2(a) follows from Lemma 3; the existence of the heteroclinic connection from p2(a) to p1(a) follows from Lemma 2. The invariant planes Π (27) and { s = 0 } are indicated in green.

Performing a coordinate change (v, q, s) = (0, 0, s¯) + wcη ^c + wsη ^s + wuη u , we obtain for the centre dynamics

$$\dot{w}_c = \frac{(\mathcal{D}(\bar{s} + w_c) - a)(\lambda_3^s(a, \bar{s})w_s + \lambda_3^u(a, \bar{s})w_u)^2}{\mathcal{D}},$$ (26)

from which we see that the plane {w^c = 0} is globally invariant if and only if ¯s = a ^D =: s ∗ . Reverting to the original coordinates, this plane transforms to

$$\Pi = { a - v - q - \mathcal{D}s = 0 }; \tag{27}$$

note that this is O(δ) close to the hyperplane {u = 0, w = a}, cf. (11).

We observe that both p2(a), p3(a, s^∗ ) ∈ Π. Since the eigenvalues λ 2 (a) of the equilibrium p2(a) are given by

$$\lambda_2^{(1)}(a) = -c^3, \quad \lambda_2^{(2)}(a) = -\frac{\mathcal{H}a}{\mathcal{D}}, \quad \lambda_2^{(3)}(a) = a^2,$$ (28)

the stable manifold of p2(a) is two-dimensional. The stable eigenspace of p2(a) is spanned by (−1, 1, 0) and −c ³D² , aDH, c3D − aH . Inspection of these eigenvectors reveals that the stable eigenspace of p2(a) is contained in the plane Π; hence, the two planes coincide. The stable manifold W^s p2(a) is therefore tangent to Π at p2(a). However, since both W^s p2(a) and Π are two-dimensional manifolds that are invariant under the flow of (13) and p2(a) is hyperbolic, it follows that W^s p2(a) and Π must coincide, as the stable manifold is (locally) unique.

The dynamics of (13) on Π = W^s p2(a) are given by

$$\dot{v} = c^3 q, \dot{q} = \frac{\mathcal{H}}{\mathcal{D}} (a - q - v) v - c^3 q.$$ (29)

The equilibrium (v, q) = (0, 0) of (29), representing p3(a, s^∗ ), is of saddle type. Its onedimensional unstable manifold lies per construction inside the stable manifold of p2(a), represented by the equilibrium (a, 0) of (29), which is a sink. It follows that a heteroclinic orbit from (0, 0) to (a, 0) in (29) is guaranteed to exist.

3.3 Singular orbit

We can combine the results obtained for the layer problem (13) from Lemma 2 and Lemma 3 with the analysis of the reduced problems presented in Section 3.1 to obtain a singular structure for the homoclinic orbit to the desert state (a, 0, 0, 0). In the following, for any w, w such that w < w we refer to subsets of M(1) and M(3),w⁰ using the notation

$$\mathcal{M}^{(i)}(\underline{w}, \overline{w}) = \left{ p_i(w) : w \in [\underline{w}, \overline{w}] \right}, \quad i = 1, 2,$$

$$\mathcal{M}^{(3), w_0}(\underline{w}, \overline{w}) = \left{ p_3(w, s(w)) : w \in [\underline{w}, \overline{w}] \right}.$$ (30)

Proposition 4. For any fixed value of a, B, D, H > 0 and given c = c ^∗ = (a θ0) 2/3 , there exists a unique singular skeleton orbit Φ⁰ homoclinic to p1(a) = (a, 0, 0, 0) corresponding to a singular solution to system (12), composed of segments of orbits of the layer problem (13) and the reduced problems (18) and (20). In particular, the singular solution Φ⁰ is given by the concatenation

$$\Phi0 = \mathcal{M}^{(1)}(0, a) \cup \mathcal{M}^{(3), 0}(0, a) \cup \phi{f, 1} \cup \phi_{f, 2}, \tag{31}$$

where ϕf,¹ and ϕf,² are the heteroclinic orbits determined in Lemma 2 and Lemma 3, respectively.

Proof. The concatenation of the two fast heteroclinic orbits ϕf,¹ and ϕf,² determines a fast heteroclinic connection from p3(a, s^∗ ) to p1(a), with s ^∗ = a D as given in Lemma 3; Lemma 2 determines the value of c ∗ (22).

From the definition of M(3),w⁰ (21), it follows that p3(a, s^∗ ) ∈ M(3),⁰ . From the fibration of M(3) along M(1), combined with the slow-fast structure of the flow (18) on M(3), it follows that the slow singular segment is composed of a superslow segment on M(1) from p1(a, 0, 0, 0) to the origin, followed by a slow segment on M(3),⁰ from the origin to p3(a, s^∗ ).

For a sketch of the singular skeleton Φ0, see Figure 5.

Figure 5: The singular skeleton orbit Φ⁰ (31) given in Proposition 4, obtained by matching superslow (purple), slow (blue) and fast (green) orbits of the reduced problems (20), (18) and the layer problem (13).

3.4 Persistence

The singular skeleton orbit Φ⁰ (31) provides the backbone for the main result of this paper: the existence of a homoclinic solution in system (12). The persistence of Φ⁰ for 0 < δ ≪ 1 is given in the following Theorem.

Theorem 5. Let 0 < δ ≪ 1 be sufficiently small, en let a, B, D, H be fixed O(1) and positive. There exists a unique value c = c ∗ δ for which a unique orbit Φ^δ to (12) exists that is homoclinic to p1(a) = (a, 0, 0, 0). Moreover, Φ^δ is O(δ) close to Φ⁰ (31), and c ∗ δ is O(δ) close to c ∗ (22).

Proof. Our goal is to show that the singular skeleton orbit Φ⁰ (31) persists for 0 < δ ≪ 1 as a solution to (12) that is O(δ) close to Φ0.

First, we focus on the dynamics on the manifold M(3), which is invariant under the flow of (12). The dynamics on M(3) are given by (18). System (18) is linear, and can thus be solved explicitly, yielding

$$w(\sigma) = a + c_1 \mathcal{D}e^{\frac{\sigma}{\mathcal{D}}} + (c_2 - c_1 \mathcal{D})e^{\frac{\delta\sigma}{1+\delta}}, \qquad s(\sigma) = c_1 e^{\frac{\sigma}{\mathcal{D}}}.$$ (32)

Equivalently, solving

$$\frac{dw}{ds} = 1 + \frac{\delta}{1+\delta} \frac{(w - \mathcal{D}s - a)}{s},$$

we obtain

$$w(s) = a + \mathcal{D}s + k_1(s(1+\delta))^{\frac{\delta \mathcal{D}}{1+\delta}} =: w_3(s; k_1).$$ (33)

Therefore, taking a neighbourhood of the point p3(a, s^∗ ) on M(3) as follows

$$\Sigma := {(w, v, q, s) : v = q = 0, w = a, |s - s^*| < \mu}$$ (34)

for µ sufficiently small, we have that the α-limit set of Σ solely consists of the point p1(a). The section Σ can thus be interpreted as the exit section of the intersection of the unstable manifold of p1(a) with M(3) given by

$$W_{\delta,\mu}^{u}(p_1(a)) := \left{ (w, v, q, s) : v = q = 0, w_3(s, k_1^+(\mu)) < w < w_3(s, k_1^-(\mu)), s > 0 \right}, \quad (35)$$

where w3(s, k^± 1 (µ)) = a + Ds 1 − a±µD Ds 1⁻ ^δ^D 1+δ ! , see also Figure 6. We remark that by definition W^u δ,µ(p1(a)) ∩ {w = a} = Σ.

Figure 6: Sketch of the dynamics on M(3) as described in the proof of Theorem 5. The teal set corresponds to W^u δ,µ(p1(a)) as defined in (35), the red section represents Σ as in (34), and the blue curve shows the set Γ defined in (36). The teal and blue points correspond to p1(a) and p3(a, s^∗ ), respectively.

The concatenation of the pair of heteroclinics determined in Lemma 2 and Lemma 3 provides a heteroclinic connection between p1(a) and p3(a, s^∗ ) via p2(a). We first observe that the argument used in the proof of Lemma 2 remains valid for other values of w, with associated speed value c = c ∗ ^w = (w θ0) 2/3 . Furthermore, the construction used in the proof of Lemma 3 also carries over for other values of w, as an invariant hyperplane containing a heteroclinic from p3(w, s^∗ ^w) to p2(w) can be found when s = s ∗ ^w = w D

Consequently, let us consider the neighbourhood of p1(a) given by M(1)(a − µD, a + µD) and track its singular fast flow backwards until it reaches again M(3). Our aim is to show that on this hyperplane W^s (M(1)(a − µD, a + µD)) intersects Σ transversely in the point p3(a, s^∗ ), which in turn will guarantee the persistence of a perturbed solution O(δ)-close to the singular one.

We start by observing that M(1)(a − µD, a + µD) admits a strongly stable foliation, which by construction we can follow backwards until M(2) - a manifold of normally hyperbolic equilibria of the layer problem (13). This set thus corresponds to the 2-dimensional manifold W^u (M(2)(a−µD, a+µD)). This manifold transversely intersects the 2-dimensional manifold W^s (M(2)(a−µD, a+µD)), which by construction (we remark here that w acts as a parameter

Figure 7: Sketch of the dynamics described in the proof of Theorem 5 in (v, s, w)-space, for fixed q = 0. The teal set corresponds to W^u δ,µ(p1(a)) as defined in (35), the red section represents Σ as in (34), and the blue solid (dashed) curve shows the set Γ as in Eq. (36) (M(2) in Eq. (14)). The teal, green, and blue points correspond to p1(a), p2(a), and p3(a, s^∗ ), respectively, and the green section reproduces the unstable eigenspace of a neighbourhood of M(3),⁰ (a − κ, a + κ)) as defined in (30).

for the layer problem) intersects M(3) in the curve given by

$$\Gamma = {(w, v, q, s) : v = q = 0, w \in [a - \mu \mathcal{D}, a + \mu \mathcal{D}], s \in [s^* - \mu, s^* + \mu] }.$$ (36)

The curve Γ intersects Σ in the unique point p3(a, s^∗ ). Thanks to Fenichel theory [56] and the Exchange Lemma [57, 58], analogously tracking M(1)(a−µD, a+µD) for 0 < δ ≪ 1 until M(3), will lead to an intersection Σ in a point p ∗ 3 (a, s^∗ δ ) which is O(δ)-close to p3(a, s^∗ ); the corresponding final piece of the orbit "returning" to w = a on M(3) is then given by

$$w3^*(s) = a + \mathcal{D}s \left( 1 - \left( \frac{s\delta^*}{s} \right)^{1 - \frac{\delta \mathcal{D}}{1 + \delta}} \right). \tag{37}$$

We note that for every fixed 0 < δ ≪ 1 we can choose µ such that |s ∗ ^δ−s ∗ | ≤ µ and w is nonneg-

ative in $$\mathcal{W}^{u}_{\delta,\mu}(p_1(a))$$ as long as $\mu$ is $\mathcal{O}(\delta \log(\delta))$ and satisfies $\mu \leq \frac{1}{\mathcal{D}} \left( \left( \frac{1+\delta}{\delta \mathcal{D}} \frac{\delta \mathcal{D}}{1+\delta} \right)^{\frac{1+\delta}{1+\delta(1-\mathcal{D})}} - a \right)$ .

Therefore, the orbit homoclinic to p1(a) perturbs to a solution Φ^δ for 0 < δ ≪ 1 which is O(δ) close to Φ0, as claimed.

Proof of Theorem 1. The statement of Theorem 1 directly follows from Theorem 5, observing that the travelling pulse solution (U, V, S) (x, t) to system (4) corresponds to the homoclinic orbit Φ^δ in system (12). Reverting the rescalings given in (9) with c = c ∗ (22), we obtain the leading order expression for the pulse wave speed

$$\mathcal{C} = \left(\frac{\mathcal{A}^2 \theta_0^2}{\varepsilon}\right)^{1/3} + \mathcal{O}(1).$$

4 Numerical continuation

In order to illustrate the validity of our analysis in Section 3, culminating in Theorem 1, and to provide insight into the role of the autotoxicity parameters H and D on the shape and dynamics of vegetation pulse patterns, we carry out a numerical analysis of system (8) by numerical continuation using the software AUTO [55], to determine the bifurcation structure of pulse solutions to system (8). The results of this numerical analysis are shown in Figures 8, 9 and 10. All solutions were numerically determined on a domain of size L = 1000 with periodic boundary conditions.

The first row of Figure 8 shows the relation between the wave speed C and the parameters A (panel (a)), H (panel (b)), and D (panel (c)). In particular, panel (a) shows a saddle-node bifurcation of travelling pulses; an equivalent bifurcation occurs in the toxicity-free version of (8) [46, Figure 8 (a)]. For the upper branch, an increase of the precipitation rate A leads to an increase in the migration speed C, which is in line with the leading order expression for C as a function of A as stated in Theorem 1, indicated as a red dashed curve in Figure 8 (a). A solution on the upper branch for A = 1.2 and C = 5.8, indicated by a black dot, is used as a starting point for the other continuations.

Panels (b) and (c) show that, for the chosen solution, the wave speed does not measurably depend on the value of H and D; this is again in line with the statement of Theorem 1. Note that, in all panels, ε = 0.005 is fixed. Apart from panel (a), this implies that δ = εC = 0.029 is fixed as well. The strong agreement with the statement of Theorem 1 suggests that δ = 0.029 is indeed sufficiently small, except near the saddle-node bifurcation point A = 0.45, C = 1.90. The pulse solutions on the lower branch of panel (a) are characterised by a U-profile asymptotically close to the equilibrium value U^∗ = A. Since these profiles have little ecological relevance and numerically turn out to be unstable, we have therefore excluded these solutions from our analysis; see also [46, Remark 1.3].

In the second and third row of Figure 8, we show the influence on the toxicity parameters H and D on the maximum value of the field variables U, V and S. An increase in H leads to a monotonic increase of the maximum value of the surface water density Umax and of the maximum value of the toxicity Smax, whereas the maximum value of the biomass Vmax decreases monotonically. This behaviour is consistent with the fact that an increment of H corresponds to a higher plant sensibility to toxicity, leading to a lower biomass amplitude. The effect of increasing D is reciprocal to that of H; this is consistent with the fact that an increase in D corresponds to a slower growth of the toxicity S (cf. System (4)), reducing the stress of vegetation induced by the toxicity.

In Figure 9, we show the profiles of the individual state variables U, V and S for three fixed values of the parameters H and D, corresponding to the coloured circles shown in Figure 8. In panels (a-c), we observe that an increase of the plant sensitivity to autotoxicity, measured

Figure 8: Bifurcation diagrams obtained by performing numerical continuation of system (4) with AUTO varying (a) A, (b, d-f) H, and (c, g-i) D. Panels (d-f, g-i) show the dependence of the field variables (U, V, S) on the above parameters, whereas panels (a-c) depict how the migration speed C is affected by A, H, and D, respectively. The red dashed line in panel (a) represents the migration speed derived analytically in Theorem 1.

The starting point for all the figures is the same as the one fixed in Figure 1 (A = 1.2, B = 0.45, ε = 0.005, D = 4.5 and H = 1); subsequently, the corresponding control parameter is varied. All the results in panels (b),(c),(d)-(i) are related to the same point in panel (a) (top branch with A = 1.2 and C = 5.8), indicated by a black dot, since A is always fixed in all the other panels and C is not affected by H and D.

by H, leads to a narrowing of the biomass pulse and a decrease of its amplitude, whereas the amplitude of the toxicity pulse increases. In contrast, slowing down the toxicity dynamics by increasing D (panels (d-f)) increases the with of the biomass pulse and lowers the amplitude of the toxicity component.

In order to compare the numerical solutions depicted in Figure 9 with the analytical results obtained in the asymptotic case 0 < δ ≪ 1 (see Figures 4-5), we plot in Figure 10 the profiles as

Figure 9: Profiles of the field variables (first column: U, second column: V , third column: S) in terms of the length-rescaled spatial variable ξ/L, for different values of H (panels (a-c)) and D (panels (d-f)) obtained by performing numerical continuation of system (4) with AUTO. Different colours correspond to different parameter configurations as shown in Figure 8.

homoclinic orbits projected onto (v, s, w)-space (first column) and onto (v, q, s)-space (second column), with colours corresponding to the solutions shown in Figure 8 and Figure 9. We observe that the value of H has little to no influence on the shape of the homoclinic, which corresponds to the observation that the singular skeleton described in Proposition 4 does not depend on H. In fact, the value of H only occurs in the eigenvalues of p² (28) and p³ (24), determining the dynamics in the neighbourhood of the stable/unstable manifolds of the equilibria, which disappear in the singular limit δ ↓ 0. However, the case H = 0 is special, see the Discussion.

In contrast to the lack of influence of H, the value of D does significantly change the shape of the homoclinic. Decreasing D (panels (c,d)) increases the value of s ^∗ = a D (cf. Lemma 3), shifting the location of the equilibrium p3(a, s^∗ ) further away from the q-axis. Note that the position of p2(a) (15) does not depend on D.

5 Discussion

In this manuscript, we have studied the influence of autotoxicity on far-from-threshold vegetation pulses in sloped semi-arid environments, extending the Klausmeier model with a third

Figure 10: Homoclinic orbits obtained by performing numerical continuation of system (4) with AUTO in (v, s, w)-space (panels (a,c)) and in (v, s, q)-space (panels (b,d)). Different colours correspond to different parameter configurations as shown in Figure 8.

component that models the autotoxicity density. After employing a suitable rescaling, we have used geometric singular perturbation theory to establish the existence of a homoclinic orbit to the desert state, whose spatial profile corresponds to numerically observed travelling pulses.

The influence of autotoxicity is quantified by two parameters: H, measuring biomass sensitivity to autotoxicity, and D, measuring the autotoxicity decay rate. We observe that, although the shape of the pulse profile (strongly) depends on H and D (cf. Figure 9), the pulse propagation speed does not depend on H or D, at least to leading order in δ (cf. Theorem 1 and Figure 8 (b)-(c)). That is, in the model studied in this paper, the presence of autotoxicity does not influence the observed dynamics. This is in direct contrast with the results of [38], where equivalent analytical techniques were used in the context of a similar model to conclude that autotoxicity has a direct influence on far-from-threshold pulse dynamics. In the same context of [38], equivalent conclusions were drawn in [36, 37]; in addition, the influence of autotoxicity on the dynamics of vegetation patterns has also been observed in other contexts [22, 23]. A crucial difference between the latter models and the model studied in this paper, is that in the original Klausmeier model as studied e.g. in [46] and in the current paper, water movement is driven by advection, whereas in the models studied in [22, 23, 36 –38], water movement is driven by diffusion. Generalising this meta-observation, one could hypothesise that autotoxicity does influence vegetation dynamics in environments where water movement is diffusion driven, but does not influence vegetation dynamics in environments where water

movement is advection driven (such as on sloped topographies). From a more general mathematical perspective, the effect of coupling an ODE to an existing reaction-advection-diffusion system on the dynamics of localised structures in that system clearly merits further research, ideally in the context of model classes, rather than using specific models.

Our numerical investigation with AUTO allows us to highlight similarities and differences with the analysis carried out in [46]. Analogously to [46], the constructed orbits correspond to the upper branch of the A − C bifurcation diagram shown in Figure 8(a), whereas the lower branch is associated to solutions whose U-component remains asymptotically close to U^∗ = A. When continuing in H up to the limit H = 0, AUTO correctly retrieves the single pulse solution constructed in [46]: in this case, the toxicity (S) dynamics are still nontrivial, but they decouple from the vegetation (V ) dynamics (see Figure 9). It is worth noting that, while the numerics straightforwardly connect pulse solutions of (8) (for nontrivial H) to the pulse solutions considered in [46] (for H = 0), the geometric construction behind these two solution families is substantially different. One major difference is that the variable w as introduced in (11) contains an s-component which does not appear in [46], cf. Figure 10. Another crucial observation is that the equilibria p² and p³ lose their hyperbolicity for H = 0, see (28) and (24). In particular, this means that the arguments employed in the proof of Theorem 5 cannot be used anymore, as these rely on the (normal) hyperbolicity of several geometric objects such as M(3) .

We can conclude that, from a purely mathematical perspective, the inclusion of autotoxicity desingularlises the pulse construction problem. Where in [46], the use of singular blowup techniques was needed to study the dynamics near the origin, the presence of the autotoxicity component s suffices to drive the dynamics away from the origin, see Figure 6. In fact, the minimum distance of the slowest orbit segment to the origin is O(δ log δ), cf. (33).

We did not address the stability of the travelling pulse as a solution to the full PDE (4), which provides an interesting topic of future research. Furthermore, the extension of (4) to include inertial effects [33] would provide a –to our knowledge– hitherto unexplored analytical case study for the dynamics of far-from-equilibrium patterns in hyperbolic reaction-transport models.

Acknowledgements. AI and GG are members of Gruppo Nazionale per la Fisica Matematica (GNFM), Istituto Nazionale di Alta Matematica (INdAM), which has partially supported this work. AI acknowledges support from a scholarship (borsa per l'estero) granted by INdAM to carry out a research stay at the University of Leiden and an FWF Hertha Firnberg Research Fellowship (T 1199-N). GG acknowledges support from MUR (Italian Ministry of University and Research) through PRIN2022-PNRR Project No. P2022WC2ZZ "A multidisciplinary approach to evaluate ecosystems resilience under climate change".

References

[1] United Nations Convention to Combat Desertfication. The Global Land Outlook. UNCCD, Bonn, second edition, 2022. ISBN 978-92-95118-53-9. URL https://www. unccd.int/resources/global-land-outlook/global-land-outlook-2nd-edition.

[2] A.G. de Soyza, W.G. Whitford, J.E. Herrick, J.W. Van Zee, and K.M. Havstad. Early warning indicators of desertification: examples of tests in the Chihuahuan Desert. Journal of Arid Environments, 39(2):101–112, 1998. doi: 10.1006/jare.1998.0391.
[3] G. Tirabassi, J. Viebahn, V. Dakos, H.A. Dijkstra, C. Masoller, M. Rietkerk, and S.C. Dekker. Interaction network based early-warning indicators of vegetation transitions. Ecological Complexity, 19:148–157, 2014. doi: 10.1016/j.ecocom.2014.06.004.
[4] K. Gowda, S. Iams, and M. Silber. Signatures of human impact on self-organized vegetation in the Horn of Africa. Scientific Reports, 8(1):3622, 2018. doi: 10.1038/s41598- 018-22075-5.
[5] R.M. May. Thresholds and breakpoints in ecosystems with a multiplicity of stable states. Nature, 269(5628):471–477, 1977. doi: 10.1038/269471a0.
[6] M. Rietkerk, F. van den Bosch, and J. van de Koppel. Site-specific properties and irreversible vegetation changes in semi-arid grazing systems. Oikos, 80:241–252, 1997. doi: 10.2307/3546592.
[7] P.M. Saco, M. Moreno-de las Heras, S. Keesstra, J. Baartman, O. Yetemen, and J.F. Rodr´ıguez. Vegetation and soil degradation in drylands: nonlinear feedbacks and early warning signals. Current Opinion in Environmental Science & Health, 5:67–72, 2018. doi: 10.1016/j.coesh.2018.06.001.
[8] R. Bastiaansen, A. Doelman, M.B. Eppinga, and M. Rietkerk. The effect of climate change on the resilience of ecosystems with adaptive spatial pattern formation. Ecology Letters, 43:414–429, 2020. doi: 10.1111/ele.13449.
[9] M. Rietkerk, R. Bastiaansen, S. Banerjee, J. Van De Koppel, M. Baudena, and A. Doelman. Evasion of tipping in complex systems through spatial pattern formation. Science, 374, 2021. doi: 10.1126/science.abj0359.
[10] E. Meron. Nonlinear Physics of Ecosystems. CRC Press, Boca Raton, first edition, 2015. doi: 10.1201/b18360.
[11] S. Mazzoleni, G. Bonanomi, F. Giannino, M. Rietkerk, S. Dekker, and F. Zucconi. Is plant biodiversity driven by decomposition processes? an emerging new theory on plant diversity. Community Ecology, 8(1):103–109, 2007. doi: 10.1556/ComEc.8.2007.1.12.
[12] S. Mazzoleni, G. Bonanomi, G. Incerti, M.L. Chiusano, P. Termolino, A. Mingo, M. Senatore, F. Giannino, F. Carten`ı, M. Rietkerk, and V. Lanzotti. Inhibitory and toxic effects of extracellular self-dna in litter: a mechanism for negative plant–soil feedbacks? New Phytologist, 205(3):1195–1210, 2015. doi: 10.1111/nph.13121.
[13] C.A. Klausmeier. Regular and irregular patterns in semiarid vegetation. Science, 284: 1826–1828, 1999. doi: 10.1126/science.284.5421.1826.
[14] R. HilleRisLambers, M. Rietkerk, F. van de Bosch, H.H.T. Prins, and H. de Kroon. Vegetation pattern formation in semi-arid grazing systems. Ecology, 82(1):50, 2001. doi: 10.1890/0012-9658(2001)082[0050:VPFISA]2.0.CO;2.
[15] M. Rietkerk, M.C. Boerlijst, F. van Langevelde, R. HilleRisLambers, J. van de Koppel, H.H.T. Prins, and A. de Roos. Self-organisation of vegetation in arid ecosystems. The American Naturalist, 160(4):534, 2002. doi: 10.1086/342078.
[16] J.A. Sherratt. An analysis of vegetation stripe formation in semi-arid landscapes. Journal of Mathematical Biology, 51:183–197, 2005. doi: 10.1007/s00285-005-0319-5.
[17] K. Siteur, E. Siero, M.B. Eppinga, J.D.M. Rademacher, A. Doelman, and M. Rietkerk. Beyond turing: The response of patterned ecosystems to environmental change. Ecological Complexity, 20:81–96, 2014. doi: 10.1016/j.ecocom.2014.09.002.
[18] F. Carten`ı, A. Marasco, G. Bonanomi, S. Mazzoleni, M. Rietkerk, and F. Giannino. Negative plant soil feedback explaining ring formation in clonal plants. Journal of Theoretical Biology, 313:153–161, 2012. doi: 10.1016/j.jtbi.2012.08.008.
[19] S. van der Stelt, A. Doelman, G. Hek, and J.D.M. Rademacher. Rise and fall of periodic patterns for a generalized Klausmeier-Gray-Scott model. Journal of Nonlinear Science, 23(7):39–95, 2013. doi: 10.1007/s00332-012-9139-0.
[20] A. Marasco, F. Giannino, and A. Iuorio. Modelling competitive interactions and plant–soil feedback in vegetation dynamics. Ricerche di Matematica, 69:553—-577, 2020. doi: 10.1007/s11587-020-00497-6.
[21] L. Eigentler and J.A. Sherratt. An integrodifference model for vegetation patterns in semi-arid environments with seasonality. Journal of Mathematical Biology, 81:875–904, 2020. doi: 10.1007/s00285-020-01530-w.
[22] A. Iuorio, M.B. Eppinga, M. Baudena, F. Veerman, M. Rietkerk, and F. Giannino. How does negative plant-soil feedback across lifestages affect the spatial patterning of trees? Scientific Reports, 2023. doi: 10.1038/s41598-023-44867-0.
[23] A. Iuorio, M. Baudena, M.B. Eppinga, F. Giannino, M. Rietkerk, and F. Veerman. Travelling waves due to negative plant-soil feedbacks in a model including tree life-stages. Mathematical Biosciences, 2023. doi: 10.1016/j.mbs.2023.109128.
[24] K. Spiliotis, L. Russo, F. Giannino, and C. Siettos. Numerical bifurcation analysis of turing and symmetry broken patterns of a vegetation pde model. Preprint, 2023. doi: 10.48550/arXiv.2303.13248.
[25] B.J. Kealy and D.J. Wollkind. A nonlinear stability analysis of vegetative turing pattern formation for an interaction-diffusion plant-surface water model system in an arid flat environment. Bulletin of Mathematical Biology, 74(4):803–833, 2012. doi: 10.1007/s11538- 011-9688-7.
[26] Y. Zelnik, S. Kinast, H. Yizhaq, G. Bel, and E. Meron. Regime shifts in models of dryland vegetation. Philosophical Transactions of the Royal Society A, 371(2004):20120358, 2013. doi: 10.1098/rsta.2012.0358.
[27] G.Q. Sun, L. Li, and Z.K. Zhang. Spatial dynamics of a vegetation model in an arid flat environment. Nonlinear Dynamics, 73(5):2207–2219, 2013. doi: 10.1007/s11071-013- 0935-3.
[28] G. Consolo and G. Grif`o. Eckhaus instability of stationary patterns in hyperbolic reaction-diffusion models on large finite domains. Partial Differential Equations and Applications, 3:57, 2022. doi: 10.1007/s42985-022-00193-0.
[29] G. Grif´o, G. Consolo, C. Curr´o, and G. Valenti. Rhombic and hexagonal pattern formation in 2D hyperbolic reaction-transport systems in the context of dryland ecology. Physica D, 449:133745, 2023. doi: 10.1016/j.physd.2023.133745.
[30] C. Curro, G. Grifo, and G. Valenti. Turing patterns in hyperbolic reaction-transport vegetation models with cross-diffusion. Chaos, Solitons & Fractals, 176:114152, 2023. doi: 10.1016/j.chaos.2023.114152.
[31] G. Consolo, C. Curro, G. Grifo, and G. Valenti. Oscillatory periodic pattern dynamics in hyperbolic reaction-advection-diffusion models. Physical Review E, 105:034206, 2022. doi: 10.1103/PhysRevE.105.034206.
[32] G. Consolo, G. Grif`o, and G. Valenti. Dryland vegetation pattern dynamics driven by inertial effects and secondary seed dispersal. Ecological Modelling, 474:110171, 2022. doi: 10.1016/j.ecolmodel.2022.110171.
[33] G. Grif`o. Vegetation patterns in the hyperbolic Klausmeier model with secondary seed dispersal. Mathematics, 11:1084, 2023. doi: 10.3390/math11051084.
[34] S. Thompson and G. Katul. Secondary seed dispersal and its role in landscape organization. Geophysical Research Letters, 36(2):L02402, 2009. doi: 10.1029/2008GL036044.
[35] S. Thompson, S. Assouline, L. Chen, A. Trahktenbrot, T. Svoray, and G. Katul. Secondary dispersal driven by overland flow in drylands: Review and mechanistic model development. Movement Ecology, 2:7, 2014. doi: 10.1186/2051-3933-2-7.
[36] A. Marasco, A. Iuorio, F. Carten`ı, G. Bonanomi, F. Giannino, and S Mazzoleni. Water limitation and negative plant-soil feedback explain vegetation patterns along rainfall gradient. Procedia Environmental Sciences, 19:139–147, 2013. doi: 10.1016/j.proenv.2013.06.016.
[37] A. Marasco, A. Iuorio, F. Carten`ı, G. Bonanomi, D.M. Tartakovsky, S. Mazzoleni, and F. Giannino. Vegetation pattern formation due to interactions between water availability and toxicity in plant-soil feedback. Bulletin of Mathematical Biology, 76:2866–2883, 2014. doi: 10.1007/s11538-014-0036-6.
[38] A. Iuorio and F. Veerman. The influence of autotoxicity on the dynamics of vegetation spots. Physica D, 427:133015, 2021. doi: 10.1016/j.physd.2021.133015.
[39] G. Consolo, G. Grif`o, and G. Valenti. Modeling vegetation patterning on sloped terrains: The role of toxic compounds. Physica D, 459:134020, 2024. doi: 10.1016/j.physd.2023.134020.
[40] R. Bastiaansen, P. Carter, and A. Doelman. Stable planar vegetation stripe patterns on sloped terrain in dryland ecosystems. Nonlinearity, 32(8):2759, 2019. doi: 10.1088/1361- 6544/ab1767.
[41] L. Eigentler. Species coexistence in resource-limited patterned ecosystems is facilitated by the interplay of spatial self-organisation and intraspecific competition. Oikos, 130: 609–623, 2021. doi: 10.1111/oik.07880.
[42] G. Consolo and G. Grif`o. Turing vegetation patterns in flat arid environments with finite soil carrying capacity. Ricerche di Matematica, 2023. doi: 10.1007/s11587-023-00783-z.
[43] E. Byrnes, P. Carter, A. Doelman, and L. Liu. Large amplitude radially symmetric spots and gaps in a dryland ecosystem model. Preprint, 2023. doi: 10.48550/arXiv.2208.13167.
[44] D.J. Tongway. Banded Vegetation Patterning in Arid and Semiarid Environments. Springer, New York, first edition, 2001. doi: 10.1007/978-1-4613-0207-0.
[45] D. Dunkerley. Banded vegetation in some australian semi-arid landscapes: 20 years of field observations to support the development and evaluation of numerical models of vegetation pattern evolution. Desert, 23(2):165–187, 2018. URL https://jdesert.ut. ac.ir/article\_69115.html.
[46] P. Carter and A. Doelman. Traveling stripes in the Klausmeier model of vegetation pattern formation. SIAM Journal of Applied Mathematics, 78:3213–3237, 2018. doi: 10.1137/18M1196996.
[47] L. Sewalt and A. Doelman. Spatially periodic multipulse patterns in a generalized Klausmeier–Gray–Scott model. SIAM Journal on Applied Dynamical Systems, 16(2):1113– 1163, 2017. doi: 10.1137/16M1078756.
[48] J.D. Murray. Mathematical Biology II: Spatial Models and Biomedical Applications. Springer, Berlin, third edition, 2003. doi: 10.1007/b98869.
[49] R. Hoyle. Pattern formation. An introduction to methods. Cambridge University Press, New York, first edition, 2007. doi: 10.1017/CBO9780511616051.
[50] A. Doelman. Pattern formation in reaction-diffusion systems - an explicit approach. In M.A. Peletier, R.A. van Santen, and E. Siteur, editors, Complexity Science, an Introduction, chapter 4, pages 129–182. World Scientific, 2018. doi: 10.1142/9789813239609 0004.
[51] The MathWorks Inc. Matlab version: 9.14.0 (r2023a), 2023. URL https://www. mathworks.com.
[52] F. Veerman and A. Doelman. Pulses in a Gierer-Meinhardt equation with a slow nonlinearity. SIAM Journal on Applied Dynamical Systems, 12(1):28–60, 2013. ISSN 1536-0040. doi: 10.1137/120878574. URL http://dx.doi.org/10.1137/120878574.
[53] A. Doelman and F. Veerman. An explicit theory for pulses in two component, singularly perturbed, reaction-diffusion equations. Journal of Dynamics and Differential Equations, 27(3):555–595, 2015. ISSN 1040-7294. doi: 10.1007/s10884-013-9325-2. URL http: //dx.doi.org/10.1007/s10884-013-9325-2.
[54] A. Doelman, R.A. Gardner, and T.J. Kaper. Large stable pulse solutions in reactiondiffusion equations. Indiana University Mathematics Journal, 50(1):443–507, 2001. ISSN 0022-2518. doi: 10.1512/iumj.2001.50.1873. URL http://dx.doi.org/10.1512/iumj. 2001.50.1873. [link].
[55] E.J. Doedel. AUTO: A program for the automatic bifurcation analysis of autonomous systems. Congressus Numererantium, 30(265-284):25–93, 1981.
[56] N. Fenichel. Geometric singular perturbation theory for ordinary differential equations. J. Differential Equations, 31(1):53–98, 1979. ISSN 0022-0396. doi: 10.1016/0022- 0396(79)90152-9.
[57] C.K.R.T. Jones, T.J. Kaper, and N. Kopell. Tracking invariant manifolds up to exponentially small errors. SIAM J. Math. Anal., 27(2):558–577, 1996. doi: 10.1137/s003614109325966x.
[58] C. Kuehn. Multiple Time Scale Dynamics. Springer International Publishing, 2015. doi: 10.1007/978-3-319-12316-5. URL https://doi.org/10.1007%2F978-3-319-12316-5.

grifofarequilibriumtravellingpulses2024.md 71 KB Историја Датотека