Hostname: page-component-cd9895bd7-dzt6s Total loading time: 0 Render date: 2024-12-23T07:32:03.792Z Has data issue: false hasContentIssue false

Analysis and numerical simulations of travelling waves due to plant–soil negative feedback

Published online by Cambridge University Press:  07 December 2023

Annalisa Iuorio*
Affiliation:
Faculty of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, Vienna, 1090, Austria Department of Engineering, Parthenope University of Naples, Centro Direzionale - Isola C4, Naples, 80143, Italy
Nicole Salvatori
Affiliation:
Department of Agri-Food, Animal and Environmental Sciences (DI4A), University of Udine, via delle Scienze 206, Udine, 33100, Italy Primo Principio Società Cooperativa, PortoConte Ricerche, viale Lincoln 5, Alghero, 07041, Italy
Gerardo Toraldo
Affiliation:
Department of Mathematics and Physics (DMF), University of Campania Luigi Vanvitelli, viale Lincoln 5, Caserta, 81100, Italy
Francesco Giannino
Affiliation:
Department of Agricultural Sciences, University of Naples Federico II, via Università 100, Portici, 80055, Italy
*
Corresponding author: Annalisa Iuorio; Email: [email protected]
Rights & Permissions [Opens in a new window]

Abstract

In this work, we carry out an analytical and numerical investigation of travelling waves representing arced vegetation patterns on sloped terrains. These patterns are reported to appear also in ecosystems which are not water deprived; therefore, we study the hypothesis that their appearance is due to plant–soil negative feedback, namely due to biomass-(auto)toxicity interactions.

To this aim, we introduce a reaction-diffusion-advection model describing the dynamics of vegetation biomass and toxicity which includes the effect of sloped terrains on the spatial distribution of these variables. Our analytical investigation shows the absence of Turing patterns, whereas travelling waves (moving uphill in the slope direction) emerge. Investigating the corresponding dispersion relation, we provide an analytic expression for the asymptotic speed of the wave. Numerical simulations not only just confirm this analytical quantity but also reveal the impact of toxicity on the structure of the emerging travelling pattern.

Our analysis represents a further step in understanding the mechanisms behind relevant plants‘ spatial distributions observed in real life. In particular, since vegetation patterns (both stationary and transient) are known to play a crucial role in determining the underlying ecosystems’ resilience, the framework presented here allows us to better understand the emergence of such structures to a larger variety of ecological scenarios and hence improve the relative strategies to ensure the ecosystems’ resilience.

Type
Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press

1. Introduction

Vegetation patterns have been widely studied in recent years due to their importance as ecological indicators [Reference Chen, Kolokolnikov, Tzou and Gai10, Reference Gilad, von Hardenberg, Provenzale, Shachak and Meron16, Reference Klausmeier19, Reference Lei, Zhang and Zhou20, Reference Rietkerk, Boerlijst, van Langevelde, HilleRisLambers, van de Koppel, Kumar, Prins and de Roos32, Reference Rietkerk, Dekker, de Ruiter and van de Koppel33]. Changes in their structure, in fact, provide a valuable insight into ecological changes and illustrate the flexibility and resilience of an ecosystem [Reference Scheffer, Carpenter, Foley, Folke and Walker37, Reference Yizhaq, Gilad and Meron40]; vegetation patterns can therefore be used as early-warning indicators for desertification [Reference Rietkerk, Dekker, de Ruiter and van de Koppel33]. These patterns can appear in different forms depending on the type of terrain considered: while spots, gaps, and labyrinths can be observed on flat grounds, stripes and arced bands typically appear on sloped terrains. These types of patterns were observed as early as 1950 [Reference Macfadyen23] under different slope conditions [Reference Valentin, d’Herbès and Poesen38] and have been studied in numerous works including [Reference Consolo, Currò and Valenti11, Reference Deblauwe, Barbier, Couteron, Lejeune and Bogaert12, Reference Gandhi, Werner, Iams, Gowda and Silber15]. As reported by recent experimental work, this type of pattern appears also in other biological systems, such as fungi [Reference Allegrezza, Bonanomi, Zotti, Carteni, Moreno, Olivieri, Garbarino, Tesei, Giannino and Mazzoleni1]. Partial differential equations (PDE) have been widely used with the aim to reproduce experimental observations by linking the occurrence of such patterns to space-dependent factors [Reference Bastiaansen, Chirilus-Bruckner and Doelman2, Reference Bastiaansen and Doelman3, Reference Gandhi, Werner, Iams, Gowda and Silber15, Reference Rietkerk and van de Koppel34, Reference Valentin, d’Herbès and Poesen38]. Most models consist of two/three reaction-diffusion-advection equations accounting only for biomass and water interactions. While this assumption may result reasonably when approaching arid and semi-arid environments, aerial photographs and data collection show that these patterns can emerge also when water is not a scarce resource [Reference Bonanomi, Incerti, Stinca, Cartenì, Giannino and Mazzoleni6, Reference Leprun21, Reference Marasco, Iuorio, Cartenì, Bonanomi, Tartakovsky, Mazzoleni and Giannino26]. In such cases, further biological mechanisms should hence be included in the model in order to justify their emergence.

In this context, an ecological factor which has proven to play a crucial role in vegetation dynamics is the so-called autotoxicity, produced by plant litter decomposition [Reference Berg and McClaugherty4] and whose occurrence is related to the exposure to fragmented self-DNA [Reference Mazzoleni, Bonanomi, Giannino, Rietkerk, Dekker and Zucconi28, Reference Mazzoleni, Bonanomi, Incerti, Chiusano, Termolino, Mingo, Senatore, Giannino, Cartenì, Rietkerk and Lanzotti29]. In particular, the involvement of accumulated extracellular DNA in plant–soil negative feedback has been observed experimentally, providing further evidence of its importance in shaping plant communities [Reference Bonanomi, Zotti, Idbella, Termolino, De Micco and Mazzoleni7]. Mathematical models incorporating biomass-autotoxicity interactions have been increasingly studied in recent years in order to understand the impact of autotoxicity on several relevant ecological issues, such as species coexistence and biodiversity [Reference Bonanomi, Giannino and Mazzoleni5, Reference Mazzoleni, Bonanomi, Giannino, Incerti, Dekker and Rietkerk27, Reference Mazzoleni, Bonanomi, Giannino, Rietkerk, Dekker and Zucconi28]. Moreover, including explicit spatial components using reaction-diffusion-ODE models has allowed for an investigation of autotoxicity effects on plants’ spatial organisation, such as clonal rings [Reference Bonanomi, Incerti, Stinca, Cartenì, Giannino and Mazzoleni6, Reference Cartenì, Marasco, Bonanomi, Mazzoleni, Rietkerk and Giannino8], fairy rings [Reference Karst, Dralle and Thompson18, Reference Salvatori, Moreno, Zotti, Iuorio, Cartenì, Bonanomi, Mazzoleni and Giannino35] and vegetation patterns [Reference Marasco, Giannino and Iuorio25, Reference Marasco, Iuorio, Cartenì, Bonanomi, Tartakovsky, Mazzoleni and Giannino26]. Its influence on sloped patterns, however, has not been fully understood yet.

In this work, following the PDE model presented in [Reference Cartenì, Marasco, Bonanomi, Mazzoleni, Rietkerk and Giannino8], we consider a novel advection term with the aim to understand the effect of autotoxicity on the formation of arced bands on sloped terrains. Our analysis reveals the absence of Turing patterns; however, as suggested by our numerical investigation, the system supports the occurrence of travelling waves. In particular, focusing on one-dimensional sections of the emerging profiles, we observe that while travelling fronts arise when the effect of toxicity is neglected from the biomass dynamics, these fronts turn into pulses when the dynamics of these variables are coupled. These transient patterns often emerge in different kinds of biological and ecological systems and have therefore raised a considerable interest in the scientific community in the past few years (see, among others, [Reference Deng and Zhang13, Reference Doelman, Kaper and Zegeling14, Reference Iuorio and Veerman17, Reference Scheel36, Reference van Saarloos39]). This result is consistent with the current available literature results, since the presence of a slope induces a movement in the emergent pattern (see e.g. [Reference Carter and Doelman9, Reference Klausmeier19]). The structure of our model allows us to obtain an analytical result for the critical wave speed which is confirmed by numerical observations.

The paper is organised as follows: in Section 2, we introduce the model, which we subsequently analyse in Section 3. In particular, we investigate stability with respect to homogeneous and heterogeneous perturbations (which exclude the emergence of Turing patterns) and derive a formula for the asymptotic wave speed of travelling wave solutions. Numerical results aimed to visualise the system’s dynamics and support the existence of travelling waves as well as investigate the slope intensity on the height of the biomass peak are discussed in Section 4. We finally conclude our work with a discussion of the results and future research perspectives in Section 5.

2. The mathematical model

Our modelling framework consists of two partial differential equations (PDEs) describing the dynamics of biomass $B$ ( $\textrm{kg} \cdot \textrm{m}^2$ ) and toxicity $T$ ( $\textrm{kg} \cdot \textrm{m}^2$ ) within a smooth domain $\Omega \in \mathbb{R}^2$ . These variables are assumed to evolve as follows: biomass grows according to a growth rate $g$ inhibited by toxicity via a sensitivity coefficient $s$ , dies with a mortality rate $d$ and diffuses proportionally to the constant parameter $D_B$ . In particular, we assume that (even in the absence of toxicity) the biomass dynamics follows a logistic law, i.e. the biomass growth is limited by resource abundance via a carrying capacity term given by $g/d$ . On the other hand, toxicity increases proportionally to the biomass death via a coefficient $c$ , decays at rate $q$ , diffuses in the soil proportionally to $D_T$ and additionally spreads with a velocity $V_T$ in the downhill direction of the terrain whose slope is described by $m$ . The intensity of the total advection term is hence given by the product between these two quantities, which have been kept separate to emphasise their ecological derivation and to include an explicit space dependence linked with real topographical data in the future.

The above assumptions result in the following system of reaction-diffusion-advection equations:

(2.1a) \begin{align} \frac{\partial B}{\partial t} &= g\, B(1-s\,T)- d\,B^2 + D_B \, \Delta B, \end{align}
(2.1b) \begin{align} \frac{\partial T}{\partial t} &= c\, d\, B^2 - q\, T + D_T \, \Delta T + V_T \, m \, \partial _x T. \end{align}

A detailed description of the model parameters, together with their values and units, is illustrated in Table 1.

Table 1. Model parameters and their values for equations (2.1a), (2.1b)

3. Analysis

In order to investigate the emergence of patterned structures as well as travelling waves, we first reduce equations (2.1a), (2.1b) to a suitable nondimensional form. To this aim, we rescale the model variables as follows

(3.1) \begin{equation} \begin{gathered} \tilde{t} = g \, t, \quad \tilde{x} = \sqrt{\frac{g}{D_B}} x, \quad \tilde{y} = \sqrt{\frac{g}{D_B}} y \\[5pt] \tilde{B} = \frac{d}{g} B, \quad \tilde{T} = s \, T. \end{gathered} \end{equation}

Therefore, the corresponding nondimensional system is (we omit the tildes for the sake of readability)

(3.2a) \begin{align} \frac{\partial B}{\partial t} &= B(1-T)-B^2+ \Delta B, \end{align}
(3.2b) \begin{align} \frac{\partial T}{\partial t} &= \alpha B^2-\beta T + \gamma \, \partial _x T + \delta \Delta T, \end{align}

where we have introduced the following nondimensional parameters:

(3.3) \begin{equation} \alpha = \frac{c\, g\, s}{d}, \quad \beta =\frac{q}{g}, \quad \gamma =\frac{V_T \, m}{\sqrt{g\, D_B}}, \quad \delta = \frac{D_T}{D_B}. \end{equation}

The parameter $\alpha$ , representing the nondimensional production of $T$ due to $B$ , collects all the parameters contributing to the toxicity production; in particular, it describes the ratio between the biomass growth and death rates ( $g$ and $d$ , respectively) times the factors of conversion and sensitivity to toxicity ( $c$ and $s$ , respectively). Moreover, we note that the parameter $\beta$ here coincides with $\gamma$ in [Reference Cartenì, Marasco, Bonanomi, Mazzoleni, Rietkerk and Giannino8] – coherently with the structure for the $T$ equation presented here – and represents the characteristic rate of toxicity dynamics relative to plant growth rate. In particular, for fixed values of $g$ this parameter can be interpreted as the nondimensional version of the toxicity decay rate. The parameter $\gamma$ represents the ratio between the lateral flow of the toxicity and the radial diffusion of the biomass, whereas $\delta$ corresponds to the ratio between the two diffusion coefficients. According to the ranges identified for the dimensional parameters in Table 1, we have that

(3.4) \begin{equation} \alpha \in \left [0, 30\right ], \quad \beta \in \left [ \frac{1}{30}, \frac{1}{3} \right ], \quad \gamma \in \left [ 0, \sqrt{30} \right ], \quad \delta \in \left [ \frac 85, 4 \right ]. \end{equation}

3.1. Equilibria and stability

Equations (3.2a), (3.2b) admit two equilibria, obtained by solving

(3.5) \begin{equation} \begin{aligned} 0 &= B(1-T)-B^2, \\[5pt] 0 &= \alpha \, B^2-\beta \, T, \end{aligned} \end{equation}

and discarding the one which does not satisfy the feasibility assumption $B \geq 0$ . These two equilibria are

(3.6) \begin{equation} \begin{aligned} (B_0,T_0) &= (0,0) \quad &\text{(bare soil)},\\[5pt] (B_c,T_c) &= \left ( \frac{-\beta + \sqrt{\beta (4\alpha +\beta )}}{2 \alpha }, \frac{2\alpha + \beta -\sqrt{\beta (4\alpha +\beta )}}{2 \alpha } \right ) \quad &\text{(coexistence)}. \end{aligned} \end{equation}

The steady-state $(B_c,T_c)$ can also be written as $(\Gamma, 1-\Gamma )$ , with

(3.7) \begin{equation} \Gamma = \frac{-\beta + \sqrt{\beta (4\alpha +\beta )}}{2 \alpha }. \end{equation}

We note that the coexistence steady-state is monotonic in $\beta$ ; in particular, $B_c$ is increasing and $T_c$ is decreasing with respect to the parameter $\beta$ . This is consistent with the idea that the faster the toxicity is removed from the soil, the higher the biomass density and the lower the toxicity density, respectively. Conversely, $B_c$ is decreasing and $T_c$ is increasing with respect to the parameter $\alpha$ : as expected, in fact, higher values of the factors related to biomass which increase toxicity lead to higher toxicity and lower biomass at steady state (see Figure 1).

Figure 1. $B$ (left panel) and $T$ (right panel) components of the spatially homogeneous coexistence equilibria $(B_c, T_C)$ as functions of $\beta$ , for $\alpha =6$ (dotted line), $\alpha =18$ (dashed line) and $\alpha =30$ (solid line). $B_c$ is increasing in $\beta$ and decreasing in $\alpha$ , whereas $T_c$ is decreasing in $\beta$ and increasing in $\alpha$ .

Figure 2. Bifurcation diagram associated with equations (3.2a), (3.2b) with respect to parameter $\beta$ for $\alpha =18$ . The blue and purple lines correspond to $B_0$ and $B_c$ , respectively. Solid (resp. dashed) lines represent stability (resp. instability) with respect to spatially homogeneous perturbations, respectively.

3.1.1. Stability to homogeneous perturbations

The Jacobian matrix associated with equation (3.5) is given by

(3.8) \begin{equation} J= \begin{pmatrix} 1-2\,B_\ast -T_\ast\;\;\;\;\; & -B_\ast \\[5pt] 2 \, \alpha \, B_\ast\;\;\;\;\; & -\beta \end{pmatrix}, \end{equation}

where $(B_\ast, T_\ast )$ represents a generic equilibrium. We therefore have that (see also Figure 2)

  • The bare soil equilibrium is a saddle, since $J$ evaluated at $(B_0,T_0)$ admits the eigenvalues $\lambda _1=1,$ $\lambda _2=-\beta$ .

  • The coexistence equilibrium is stable, since $J$ evaluated at $(B_c,T_c)$ satisfies

    \begin{equation*} \begin {aligned} \textrm {tr}\,J \vert _{(B_c,T_c)} &= \frac {\beta (1-2 \alpha ) -\sqrt {\beta(4 \alpha +\beta )}}{2 \alpha } \lt 0, \\[5pt] \textrm {det}\,J \vert _{(B_c,T_c)} &= \frac {\beta \!\left (4 \alpha +\beta -\sqrt {\beta (4 \alpha +\beta )} \right )}{2 \alpha } \gt 0, \end {aligned} \end{equation*}
    for any $\alpha, \, \beta \gt 0$ .

3.1.2. Stability to heterogeneous perturbations

In order to investigate the stability of the steady-states (3.6) with respect to heterogeneous perturbations, we introduce the following non-uniform infinitesimal perturbations:

(3.9) \begin{equation} \begin{aligned} B(t,x,y) &= B_\ast + \tilde{B}(0)\, e^{i (k x+ l y)+\lambda t}, \\[5pt] T(t,x,y) &= T_\ast + \tilde{T}(0)\, e^{i (k x+ l y)+\lambda t}. \\[5pt] \end{aligned} \end{equation}

In particular, the norm of the spatial perturbation is defined as $h=\sqrt{k^2+l^2}$ .

When we linearise equations (3.2a), (3.2b) around the abovementioned steady states, we obtain the following system for the perturbations $\tilde{B}$ , $\tilde{T}$ (again we omit the tildes for sake of simplicity):

(3.10) \begin{equation} \begin{aligned} \lambda \, B &= (1-2\,B_\ast -T_\ast )\,B-B_\ast \,T-h^2 B, \\[5pt] \lambda \, T &= 2\, \alpha \,B_\ast \, B-\beta \, T+ i\, k\, \gamma \,T-\delta h^2\, T. \end{aligned} \end{equation}

Equation (3.10) is equivalent to an eigenvalue problem $J_h\, U = \lambda \, U$ , where $U=(B,T)$ and

(3.11) \begin{equation} J_h= \begin{pmatrix} 1-2\,B_\ast -T_\ast - h^2 & -B_\ast \\[5pt] 2 \, \alpha \, B_\ast & -\beta + i\, k\, \gamma -\delta \,h^2 \end{pmatrix}. \end{equation}

The characteristic polynomial, necessary to investigate the stability of (3.6), is given by

(3.12) \begin{equation} \lambda ^2-\textrm{tr}\,J_h \,\lambda +\textrm{det}\,J_h=0. \end{equation}
  • The eigenvalues of $J_h$ evaluated at $(B_0,T_0)$ are $\lambda _1=1-h^2$ , $\lambda _2=-\beta + i\, k\, \gamma -\delta \,h^2$ . Since this steady state is unstable with respect to homogeneous perturbations (see Section 3.1.1), no diffusion-driven instability occurs, and therefore, the stability character of the bare soil equilibrium is not relevant in the investigation of vegetation patterns. However, we observe that for $h$ such that $h^2\lt 1$ the bare soil equilibrium retains its saddle structure, whereas for $h$ such that $h^2\gt 1$ it becomes stable (hence diffusion may induce a stabilising effect, as often expected).

  • As for the coexistence equilibrium, we have that the trace and the determinant of $J_h$ evaluated at $(B_c,T_c)$ are

    (3.13) \begin{equation} \begin{aligned} \textrm{tr}\,J_h &= -\beta -\Gamma -h^2(\delta +1)+ i\,k\,\gamma, \\[5pt] \textrm{det}\,J_h &= (h^2+\Gamma )(\beta +\delta \, h^2-i\,k\,\gamma )+2\,\alpha \,\Gamma ^2. \end{aligned} \end{equation}
    Equation (3.12) is solved by two roots given by
    (3.14) \begin{equation} \lambda _\pm =\frac 12 \left ( -\beta -\Gamma -h^2(\delta +1)+i\,k\,\gamma \right ) \pm \sqrt{X+i\,Y}, \end{equation}
    where
    (3.15) \begin{equation} \begin{aligned} X &= (\beta - \Gamma +h^2(\delta -1))^2 - k^2\,\gamma ^2-8\,\alpha \,\Gamma ^2, \\[5pt] Y &= -2\,k\,\gamma \left ( \beta -\Gamma +h^2(\delta -1) \right ). \end{aligned} \end{equation}
    The real and imaginary part of the eigenvalues (3.14) can be written as
    (3.16a) \begin{align} \textrm{Re}\,\lambda &= \frac 12 \left ( -\beta -\Gamma -h^2(\delta +1)+j\sqrt{\frac{X+\sqrt{X^2+Y^2}}{2}} \right ), \end{align}
    (3.16b) \begin{align} \textrm{Im}\,\lambda &= \frac 12 \left ( \,k\,\gamma +j\, \textrm{sgn}(Y) \sqrt{\frac{\sqrt{X^2+Y^2}-X}{2}} \right ), \end{align}

    where $j=\pm 1$ . In order for patterns to emerge, there must exist $k$ , $l$ such that $\textrm{Re}\,\lambda \gt 0$ [Reference Murray30]. It is already clear that this condition cannot be satisfied in the case $j=-1$ ; therefore, we focus our forthcoming analysis on the case $j=+1$ . Since the right-hand side of (3.16a) is given by the sum of a positive term (the square root) and a negative term (the rest), for this quantity to be positive it must be

    \begin{equation*} \sqrt {\frac {X+\sqrt {X^2+Y^2}}{2}} \gt \beta +\Gamma +h^2(\delta +1). \end{equation*}
    Using the ranges determined for the nondimensional parameters in (3.4), we have that $\textrm{Re}\,\lambda \lt 0$ for any $h$ and for any $\alpha$ , $\beta$ , $\delta$ in (3.4). Therefore, no Turing patterns can emerge.

3.2. Travelling waves

In this section, we consider the formation of travelling wave profiles, i.e. profiles which move through the spatial domain keeping constant shape and speed.

Extensive simulations of equations (2.1a), (2.1b) reveal that one-dimensional sections of compact initial data $B(0,x,y)$ and $T(0,x,y)$ numerically evolve into travelling waves for a wide parameter range (see Figure 3, Section 4). In particular, the unstable nature of the trivial steady-state $(B_0, T_0)$ suggests we are in the presence of pulled fronts, where the linear spreading of small perturbations pulls the front into the linearly unstable bare soil steady-state [Reference Lewis, Li and Weinberger22, Reference Maini, McElwain and Leavesley24, Reference van Saarloos39]. The (nondimensional) asymptotic travelling speed of such structures $\tilde{c}_\infty$ therefore coincides with the linear spreading speed $\tilde{c}_\infty$ . In order to compute it, we introduce the characteristic polynomial of the Jacobian matrix associated with equations (3.2a), (3.2b) evaluated at $(B_0,T_0)$ focusing on one-dimensional domains:

(3.17) \begin{equation} p_0(\lambda ) \;:\!=\; \lambda ^2+ \mathcal{A}_1 \, \lambda + \mathcal{A}_0, \end{equation}

where

(3.18) \begin{equation} \begin{aligned} \mathcal{A}_1 &\;:\!=\; k \, (k \, (1+\delta ) - i \, \gamma ) + \beta -1, \\[5pt] \mathcal{A}_0 &\;:\!=\; \left ( k^2-1 \right ) \left ( k \, (\delta \, k - i \, \gamma )+\beta \right ). \end{aligned} \end{equation}

The dispersion relation, corresponding to the solution to $p_0(\lambda )=0$ with the largest real part, is given by

(3.19) \begin{equation} \lambda (k) = 1-k^2. \end{equation}

The linear spreading speed $\tilde{c}_\infty$ is obtained by solving (see [Reference Murray30])

(3.20a) \begin{align} i \, \tilde{c}_\infty + \lambda '(k_0) &= 0, \end{align}
(3.20b) \begin{align} \textrm{Re} \left [ i \, \tilde{c}_\infty \, k_0 + \lambda (k_0) \right ] &= 0. \end{align}

Here, $k_0$ is a constant related to the asymptotic evaluation of the Fourier transform of the perturbation solving the linearised system via the method of steepest descents for $t \to \infty$ . Equation (3.20b), on the other hand, represents the assumption that such perturbation neither grows nor decays at the leading edge of the front.

Figure 3. Travelling arc moving uphill obtained by numerically simulating equations (2.1a), (2.1b) at four different time steps for biomass ( $B$ ), toxicity ( $T$ ) (first and second row panels) and a one-dimensional cross-section obtained by cutting the profiles in the middle of the two-dimensional domain along the direction of the slope. The parameter values used for the simulation are taken from Table 1 together with $m=0.5$ , $q=0.01$ , $s=3$ , $D_T=1.2$ , and $V_T=5$ .

Since $\tilde{c}_\infty \in \mathbb{R}$ , equation (3.20a) implies $k_0 = \dfrac{i \, \tilde{c}_\infty }{2}$ . Plugging this into equation (3.20b) (and focusing our attention on right-moving fronts, i.e. fronts with positive speed), we obtain

(3.21) \begin{equation} \tilde{c}_\infty = 2. \end{equation}

(This, in turn, leads to $k_0=i$ ). We note that this coincides with the critical wave speed of a travelling wave solution to the Fisher-KPP equation (see e.g. [Reference Murray31]). Heuristically, this can be explained by the absence of toxicity at the leading edge of the front, which turns equation (3.2a) precisely into the Fisher-KPP equation.

In order to compare this result with numerical simulations of the dimensional model (2.1a), (2.1b) and to facilitate the ecological interpretation of this quantity, we convert $\tilde{c}_\infty$ into its dimensional version $c_\infty$ by recalling that, due to equation (3.1), we have $\frac{dx}{dt} = \sqrt{D_B \, g} \, \frac{d\tilde{x}}{d\tilde{t}}$ , therefore

(3.22) \begin{equation} c_\infty = 2\sqrt{D_B \, g}. \end{equation}

The right-moving front hence moves with constant shape and speed $\tilde{c}_\infty$ , which solely depends on the biomass growth rate and its diffusion coefficient. The presence of the slope, regulated by the parameter $\gamma$ , induces the formation of a travelling arc (in the uphill direction) for $\gamma \gt 0$ (see Figure 3), i.e. the ‘tail’ of the ring is disrupted by the downhill flow of toxicity. Turning off the effect induced by toxicity via the slope term by considering $\gamma =0$ –representing a flat ground – substantially changes the structure of the emerging pattern: the travelling arc reduces to a ring symmetrically expanding in all directions and, depending on the intensity of the toxicity parameters, secondary travelling rings may form inside the primary ring (see Figure 4). In both cases, we recall that the stability analysis illustrated in Section 3.1.2 implies that, as $t \to \infty$ , we have convergence to the uniform equilibrium $(B_c,T_c)$ .

Figure 4. Travelling ring obtained by numerically simulating equations (2.1a), (2.1b) at four different time steps for biomass ( $B$ ), toxicity ( $T$ ) (first and second row panels) and a one-dimensional cross-section obtained by cutting the profiles in the middle of the symmetric two-dimensional domain. The parameter values used for the simulation are taken from Table 1 together with $m=0$ , $q=0.01$ , $s=3$ , and $D_T=1.2$ .

4. Numerical simulations

In this section, we perform a series of numerical simulations of equations (2.1a), (2.1b) equipped with Neumann boundary conditions on a two-dimensional square domain $[0, L] \times [0, L]$ using the software Matlab. As for the initial conditions, we set

(4.1) \begin{equation} B(0,x,y) =B_i(x,y), \qquad T(0,x,y) =T_i(x,y) \end{equation}

where $B_i(x,y) = 0$ everywhere except for $B_i(L/2,L/2)=0.5$ ( $\textrm{kg} \cdot{m}^{-2}$ ) and $T_i(x,y)=0$ ( $\textrm{kg} \cdot{m}^{-2}$ ) for all $x, y$ . We note that from here on we focus on the dimensional model, as dimensional parameters and quantities are easier to interpret from the ecological viewpoint.

In Figure 5, we investigate the influence of toxicity on the structure of the emerging travelling patterns by comparing the simulated one-dimensional profiles in the cases $s \gt 0$ (biomass dynamics affected by toxicity) and $s=0$ (biomass and toxicity dynamics decoupled) on a sloped terrain at a fixed instant. We observe that, while in the first case a travelling pulse arises moving in the uphill direction, in the second case this pulse reduces to two fronts symmetrically travelling in an outward direction. The maximum peak reached by the biomass is also sensibly lower in the case $s\gt 0$ ; toxicity further influences this quantity via the slope parameter $m$ , as investigated below (see Figure 6).

Figure 5. Comparison between two one-dimensional cross sections of biomass profiles obtained by simulating equations (2.1a), (2.1b) with initial conditions as in (4.1) for $D_B=0.5$ , $g=0.6$ , $d=0.06$ , $c=0.6$ , $V=5$ , $D_T=1.2$ , $q=0.01$ , $m=0.5$ and $s=3$ (solid line) vs $s=0$ (dashed line). In the first case, a travelling pulse forms which moves along the slope direction, whereas in the second case, we have two travelling fronts symmetrically moving outwards.

Figure 6. Variations in the peak of biomass ( $B$ , left vertical axis, continuous line) and toxicity ( $T$ , right vertical axis, dashed line) due to changes in slope ( $m$ ). The parameter values used for the simulation are taken from Table 1 with $k=0.01$ , $s=3$ , $D_T=1.2$ , and $V_T=5$ .

Figure 7. Analytical representation (equation (3.22), continuous line) and numerical simulations (circles) of the asymptotic travelling speed $c_\infty$ with respect to the biomass growth parameter $g$ (left panel) and the biomass diffusion parameter $D_B$ (right panel). The other parameter values used for the simulation are taken from Table 1 together with $m=0.5$ , $k=0.01$ , $s=3$ , $D_T=1.2$ and $V_T=5$ .

The next step in our numerical investigation focuses on the influence of the toxicity on the arising travelling patterns via the slope parameter $m$ . First, we analyse the scenario of a sloped terrain represented by $m\gt 0$ in Figure 3, which shows the evolution of the system’s dynamics at increasing time steps corresponding to 10, 50, 100 and 150 months. In this simulation, the slope coefficient $m$ is set to $0.5$ , with the uphill direction corresponding to the right-hand side of the domain. The first two rows represent the accumulation of biomass and toxicity in the square grid, respectively. As expected, biomass starts growing radially until it is inhibited by the presence of toxicity which, due to the advection term in equations (2.1a), (2.1b), diffuses prominently in the downhill direction. This induces the emergence of a travelling arced shape of biomass. This behaviour is displayed also in the third row of Figure 3 via a one-dimensional horizontal section of the square numerical domain showing the evolution in time of both biomass and toxicity.

An analogous temporal evolution of the system’s dynamics in the case of flat grounds (represented by the case $m=0$ ) is shown in Figure 4. Here, the radial growth of the biomass is symmetrically affected by toxicity and hence not disrupted, leading to rings expanding in the two-dimensional domain. With the parameter values adopted in this simulation, secondary rings form within the primary ring as a result of the fact that, after some time, the toxicity removed from the centre of the domain creates the conditions for new biomass growth. This result is consistent with the findings presented in [Reference Cartenì, Marasco, Bonanomi, Mazzoleni, Rietkerk and Giannino8].

Due to the importance of the slope intensity in the toxicity advection term, we evaluate numerically the dependence between the slope parameter $m$ and the height of the leading travelling pulse both for biomass and toxicity. As expected by the mathematical formulation, the higher the slope, the higher the advection which reduces, even if only slightly, the toxicity and its maximum peak (Figure 6). This decrease in the toxicity then induces a corresponding increase in the biomass peak of the main travelling pulse.

Finally, as the asymptotic travelling speed $c_\infty$ was shown to depend only on $D_B$ and $g$ (see equation (3.22)), we compare this expression with the numerically calculated velocity by varying both of these parameters. The numerical speed was calculated as the ratio of the number of pixels covered by the travelling peak to the number of timesteps (converted to m/month) with parameter values as in Table 1 and $k=0.01$ , $s=3$ , $D_T=1.2$ , $V_T=5$ . In both cases, numerical results are in good agreement with the analytical predictions (see Figure 7).

5. Discussion and conclusion

In this work, we introduce and analyse a reaction-diffusion-advection model describing vegetation spatial distributions (in the form of arcs on slopes) observed in nature and also in environments where water is not limited. Our investigation allows us to confirm the ecologically feasible condition of the transient nature of such patterns. Moreover, thanks to the ecological constraints reflected by the model, we are able to compute the asymptotic speed of the corresponding emerging travelling waves. This expression does not depend on the slope coefficient $m$ , nor does the emergence of such a transient pattern.

Despite the fact that the slope (represented by the parameter $m$ ) does not affect the biomass dynamics directly but only via the toxicity flow, the main effect produced by $m$ on the long-term dynamics consists of the height of the main travelling pulse of $B$ (as shown in Figure 6). In fact, as $m$ increases, the downhill flow of $T$ accordingly increases, therefore reducing the maximum of $T$ reached in the proximity of the $B$ peak. Consequently, the maximum biomass in the travelling pulse can reach a higher value – net of other relevant ecological factors not considered in our model, such as the ability of the root system to successfully establish.

Due to the ecological focus of our model, a natural next step consists of the comparison and validation with real-life data. In order to further increase potential connections with relevant ecological scenarios, a future extension we plan to consider is given by a space-dependent slope coefficient in our equations accounting for topographical differences in existing sloped profiles.

Financial support

AI acknowledges support from an FWF Hertha Firnberg Research Fellowship (T 1199-N) and is a member of Gruppo Nazionale per la Fisica Matematica, Istituto Nazionale di Alta Matematica (INdAM). FG and GT are members of Gruppo Nazionale per il Calcolo Scientifico, Istituto Nazionale di Alta Matematica (INdAM). GT was supported by the PRIN grant (Bando 2022) Numerical Optimization with Adaptive Accuracy and Applications to Machine Learning (Prot. 2022N3ZNAX).

Competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

Allegrezza, M., Bonanomi, G., Zotti, M., Carteni, F., Moreno, M., Olivieri, L.Garbarino, M.Tesei, G.Giannino, F. & Mazzoleni, S. (2022) Biogeography and shape of fungal fairy rings in the apennine mountains, italy. J. Biogeogr 49(2), 353363.10.1111/jbi.14306CrossRefGoogle Scholar
Bastiaansen, R., Chirilus-Bruckner, M. & Doelman, A. (2020) Pulse solutions for an extended Klausmeier model with spatially varying coefficients. SIAM J. Appl. Dyn. Syst. 19(1), 157.10.1137/19M1255665CrossRefGoogle Scholar
Bastiaansen, R. & Doelman, A. (2019) The dynamics of disappearing pulses in a singularly perturbed reaction–diffusion system with parameters that vary in time and space. Phys. D Nonlinear Phenom. 388, 4572.10.1016/j.physd.2018.09.003CrossRefGoogle Scholar
Berg, B. & McClaugherty, C. (2014) Does Humus Accumulate and Where? What Factors May Influence? In Plant Litter: Decomposition, Humus Formation, Carbon Sequestration, Springer, Berlin, Heidelberg, pp. 215234.10.1007/978-3-642-38821-7_11CrossRefGoogle Scholar
Bonanomi, G., Giannino, F. & Mazzoleni, S. (2005) Negative plant-soil feedback and species coexistence. Oikos 111(2), 311321.10.1111/j.0030-1299.2005.13975.xCrossRefGoogle Scholar
Bonanomi, G., Incerti, G., Stinca, A., Cartenì, F., Giannino, F. & Mazzoleni, S. (2014) Ring formation in clonal plants. Community Ecol. 15(1), 7786.10.1556/ComEc.15.2014.1.8CrossRefGoogle Scholar
Bonanomi, G., Zotti, M., Idbella, M., Termolino, P., De Micco, V. & Mazzoleni, S. (2022) Field evidence for litter and self-dna inhibitory effects on alnus glutinosa roots. New Phytol. 236(2), 399412.10.1111/nph.18391CrossRefGoogle ScholarPubMed
Cartenì, F., Marasco, A., Bonanomi, G., Mazzoleni, S., Rietkerk, M. & Giannino, F. (2012) Negative plant soil feedback explaining ring formation in clonal plants. J. Theor. Biol. 313, 153161.10.1016/j.jtbi.2012.08.008CrossRefGoogle ScholarPubMed
Carter, P. & Doelman, A. (2018) Traveling stripes in the Klausmeier model of vegetation pattern formation. SIAM J. Appl. Math. 78(6), 32133237.10.1137/18M1196996CrossRefGoogle Scholar
Chen, Y., Kolokolnikov, T., Tzou, J. & Gai, C. (2015) Patterned vegetation, tipping points, and the rate of climate change. Eur. J. Appl. Math. 26(6), 945958.10.1017/S0956792515000261CrossRefGoogle Scholar
Consolo, G., Currò, C. & Valenti, G. (2017) Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments. Appl. Math. Model. 43, 372392.10.1016/j.apm.2016.11.031CrossRefGoogle Scholar
Deblauwe, V., Barbier, N., Couteron, P., Lejeune, O. & Bogaert, J. (2008) The global biogeography of semi-arid periodic vegetation patterns. Global Ecol. Biogeogr. 17(6), 715723.10.1111/j.1466-8238.2008.00413.xCrossRefGoogle Scholar
Deng, D. & Zhang, D. (jan 2019) Existence of travelling waves with the critical speed for an influenza model with treatment. Eur. J. Appl. Math. 31(2), 232245.10.1017/S0956792518000700CrossRefGoogle Scholar
Doelman, A., Kaper, T. J. & Zegeling, P. A. (1997) Pattern formation in the one-dimensional gray - scott model. Nonlinearity 10(2), 523563.10.1088/0951-7715/10/2/013CrossRefGoogle Scholar
Gandhi, P., Werner, L., Iams, S., Gowda, K. & Silber, M. (2018) A topographic mechanism for arcing of dryland vegetation bands. J. R. Soc. Interface 15(147), 20180508.10.1098/rsif.2018.0508CrossRefGoogle ScholarPubMed
Gilad, E., von Hardenberg, J., Provenzale, A., Shachak, M. & Meron, E. (2004) Ecosystem engineers: from pattern formation to habitat creation. Phys. Rev. Lett. 93(9), 098105.10.1103/PhysRevLett.93.098105CrossRefGoogle ScholarPubMed
Iuorio, A. & Veerman, F. (2021) The influence of autotoxicity on the dynamics of vegetation spots. Phys. D Nonlinear Phenom. 427, 133015.10.1016/j.physd.2021.133015CrossRefGoogle Scholar
Karst, N., Dralle, D. & Thompson, S. (2016) Spiral and rotor patterns produced by fairy ring fungi. PLOS ONE 11(3), e0149254.10.1371/journal.pone.0149254CrossRefGoogle ScholarPubMed
Klausmeier, C. A. (1999) Regular and irregular patterns in semiarid vegetation. Science 284(5421), 18261828.10.1126/science.284.5421.1826CrossRefGoogle ScholarPubMed
Lei, C., Zhang, G. & Zhou, J. (2022) Pattern formation of a biomass–water reaction–diffusion model. Appl. Math. Lett. 123, 107605.10.1016/j.aml.2021.107605CrossRefGoogle Scholar
Leprun, J. C. (1999) The influences of ecological factors on tiger bush and dotted bush patterns along a gradient from mali to northern Burkina Faso. CATENA 37(1-2), 2544.10.1016/S0341-8162(98)00054-XCrossRefGoogle Scholar
Lewis, M. A., Li, B. & Weinberger, H. F. (2002) Spreading speed and linear determinacy for two-species competition models. J. Math. Biol. 45(3), 219233.10.1007/s002850200144CrossRefGoogle ScholarPubMed
Macfadyen, W. A. (1950) Vegetation patterns in the semi-desert plains of British Somaliland. Geogr. J. 4(6), 199.10.2307/1789384CrossRefGoogle Scholar
Maini, P. K., McElwain, D. L. S. & Leavesley, D. (2004) Travelling waves in a wound healing assay. Appl. Math. Lett. 17(5), 575580.10.1016/S0893-9659(04)90128-0CrossRefGoogle Scholar
Marasco, A., Giannino, F. & Iuorio, A. (2020) Modelling competitive interactions and plant–soil feedback in vegetation dynamics. Ric. Mat. 69(2), 553577.10.1007/s11587-020-00497-6CrossRefGoogle Scholar
Marasco, A., Iuorio, A., Cartenì, F., Bonanomi, G.Tartakovsky, D. M.Mazzoleni, S. & Giannino, F. (2014) Vegetation pattern formation due to interactions between water availability and toxicity in plant–soil feedback. Bull. Math. Biol. 76(11), 28662883.10.1007/s11538-014-0036-6CrossRefGoogle ScholarPubMed
Mazzoleni, S., Bonanomi, G., Giannino, F., Incerti, G., Dekker, S. C. & Rietkerk, M. (2010) Modelling the effects of litter decomposition on tree diversity patterns. Ecol. Model. 221(23), 27842792.10.1016/j.ecolmodel.2010.08.007CrossRefGoogle Scholar
Mazzoleni, S., Bonanomi, G., Giannino, F., Rietkerk, M., Dekker, S. & Zucconi, F. (2007) Is plant biodiversity driven by decomposition processes? an emerging new theory on plant diversity. Community Ecol. 8(1), 103109.10.1556/ComEc.8.2007.1.12CrossRefGoogle Scholar
Mazzoleni, S., Bonanomi, G., Incerti, G., Chiusano, M. L., Termolino, P., Mingo, A., Senatore, M., Giannino, F., Cartenì, F., Rietkerk, M. & Lanzotti, V. (2014) Inhibitory and toxic effects of extracellular self-DNA in litter: a mechanism for negative plant-soil feedbacks? New Phytol. 205(3), 11951210.10.1111/nph.13121CrossRefGoogle ScholarPubMed
Murray, J. D. (2001) Mathematical Biology II: Spatial Models and Biomedical Applications, Springer, New York, Vol. 3.Google Scholar
Murray, J. D. (2002) Mathematical Biology I: An Introduction. Heidelberg: Springer.10.1007/b98868CrossRefGoogle Scholar
Rietkerk, M., Boerlijst, M. C., van Langevelde, F., HilleRisLambers, R.van de Koppel, J.Kumar, L.Prins, H. H. T. & de Roos, A. M. (2002) Self-organization of vegetation in arid ecosystems. Am. Nat. 160(4), 524530.10.1086/342078CrossRefGoogle ScholarPubMed
Rietkerk, M., Dekker, S. C., de Ruiter, P. C. & van de Koppel, J. (2004) Self-organized patchiness and catastrophic shifts in ecosystems. Science 305(5692), 19261929.10.1126/science.1101867CrossRefGoogle ScholarPubMed
Rietkerk, M. & van de Koppel, J. (2008) Regular pattern formation in real ecosystems. Trends Ecol. Evol. 23(3), 169175.10.1016/j.tree.2007.10.013CrossRefGoogle ScholarPubMed
Salvatori, N., Moreno, M., Zotti, M., Iuorio, A., Cartenì, F., Bonanomi, G.Mazzoleni, S. & Giannino, F. (2022) Process based modelling of plants-fungus interactions explains fairy ring types and dynamics. Scientific Reports. 13(1), 19918.10.1038/s41598-023-46006-1CrossRefGoogle Scholar
Scheel, A. (2003). Radially Symmetric Patterns of Reaction-Diffusion Systems, American Mathematical Society, Vol. 165, No. 786.10.1090/memo/0786CrossRefGoogle Scholar
Scheffer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. (2001) Catastrophic shifts in ecosystems. Nature 413(6856), 591596.10.1038/35098000CrossRefGoogle ScholarPubMed
Valentin, C., d’Herbès, J. M. & Poesen, J. (1999) Soil and water components of banded vegetation patterns. CATENA 37(1-2), 124.10.1016/S0341-8162(99)00053-3CrossRefGoogle Scholar
van Saarloos, W. (2003) Front propagation into unstable states. Phys. Rep. 386(2-6), 29222.10.1016/j.physrep.2003.08.001CrossRefGoogle Scholar
Yizhaq, H., Gilad, E. & Meron, E. (2005) Banded vegetation: biological productivity and resilience. Phys. A Stat. Mech. Appl. 356(1), 139144.10.1016/j.physa.2005.05.026CrossRefGoogle Scholar
Figure 0

Table 1. Model parameters and their values for equations (2.1a), (2.1b)

Figure 1

Figure 1. $B$ (left panel) and $T$ (right panel) components of the spatially homogeneous coexistence equilibria $(B_c, T_C)$ as functions of $\beta$, for $\alpha =6$ (dotted line), $\alpha =18$ (dashed line) and $\alpha =30$ (solid line). $B_c$ is increasing in $\beta$ and decreasing in $\alpha$, whereas $T_c$ is decreasing in $\beta$ and increasing in $\alpha$.

Figure 2

Figure 2. Bifurcation diagram associated with equations (3.2a), (3.2b) with respect to parameter $\beta$ for $\alpha =18$. The blue and purple lines correspond to $B_0$ and $B_c$, respectively. Solid (resp. dashed) lines represent stability (resp. instability) with respect to spatially homogeneous perturbations, respectively.

Figure 3

Figure 3. Travelling arc moving uphill obtained by numerically simulating equations (2.1a), (2.1b) at four different time steps for biomass ($B$), toxicity ($T$) (first and second row panels) and a one-dimensional cross-section obtained by cutting the profiles in the middle of the two-dimensional domain along the direction of the slope. The parameter values used for the simulation are taken from Table 1 together with $m=0.5$, $q=0.01$, $s=3$, $D_T=1.2$, and $V_T=5$.

Figure 4

Figure 4. Travelling ring obtained by numerically simulating equations (2.1a), (2.1b) at four different time steps for biomass ($B$), toxicity ($T$) (first and second row panels) and a one-dimensional cross-section obtained by cutting the profiles in the middle of the symmetric two-dimensional domain. The parameter values used for the simulation are taken from Table 1 together with $m=0$, $q=0.01$, $s=3$, and $D_T=1.2$.

Figure 5

Figure 5. Comparison between two one-dimensional cross sections of biomass profiles obtained by simulating equations (2.1a), (2.1b) with initial conditions as in (4.1) for $D_B=0.5$, $g=0.6$, $d=0.06$, $c=0.6$, $V=5$, $D_T=1.2$, $q=0.01$, $m=0.5$ and $s=3$ (solid line) vs $s=0$ (dashed line). In the first case, a travelling pulse forms which moves along the slope direction, whereas in the second case, we have two travelling fronts symmetrically moving outwards.

Figure 6

Figure 6. Variations in the peak of biomass ($B$, left vertical axis, continuous line) and toxicity ($T$, right vertical axis, dashed line) due to changes in slope ($m$). The parameter values used for the simulation are taken from Table 1 with $k=0.01$, $s=3$, $D_T=1.2$, and $V_T=5$.

Figure 7

Figure 7. Analytical representation (equation (3.22), continuous line) and numerical simulations (circles) of the asymptotic travelling speed $c_\infty$ with respect to the biomass growth parameter $g$ (left panel) and the biomass diffusion parameter $D_B$ (right panel). The other parameter values used for the simulation are taken from Table 1 together with $m=0.5$, $k=0.01$, $s=3$, $D_T=1.2$ and $V_T=5$.