1 Introduction

Vegetation patterns are a ubiquitous feature of ecosystems in semi-arid climate zones. Occurrences of such mosaics of plants and bare soil have been reported from all continents except Antarctica, including the African Sahel (Deblauwe et al. 2012) and the Horn of Africa (Gowda et al. 2018), Western Australia (Gandhi et al. 2018), northern Chile (Fernandez-Oto et al. 2019), Israel (Sheffer et al. 2013), the Chihuahuan Desert in North America (Deblauwe et al. 2012) and Southeastern Spain (Lesschen et al. 2008). A detailed understanding of the evolution of vegetation patterns is of considerable importance as they hold valuable information on the health of ecosystems. For example, changes to a pattern’s properties such as its wavelength, its recovery time from perturbations, or the area fraction covered by biomass can act as early warning signals of desertification (Corrado et al. 2014; Dakos et al. 2011; Gowda et al. 2016; Kéfi et al. 2007; Rietkerk et al. 2004; Saco et al. 2018; Zelnik et al. 2018). Desertification processes are a major threat to economies in semi-deserts as agriculture provides a significant contribution to GDP (United Nations Convention to Combat Desertification 2017). For example, the livestock sector, which depends in part on animals grazing on spatially patterned vegetation, accounts for 20% of GDP in Chad and involves 40% of its population (Dickovick 2014; United Nations Food and Agriculture Organization 2005).

A number of feedback mechanisms may be involved in the pattern formation process (see Meron 2018 for a review), but it is widely agreed that a central mechanism is the vegetation-infiltration feedback loop, which results in a redistribution of water towards areas of high biomass. On bare soil, the formation of physical and biological soil crusts inhibits water infiltration into the soil (Eldridge et al. 2000). Thus, water run-off towards existing vegetation patches occurs. The enhancement of environmental conditions in these sinks for the limiting resource drives further plant growth and thus closes the feedback loop (Thompson et al. 2010).

Dryland plants have developed a range of seed production and dispersal strategies to cope with the environmental stress in their habitats (Ellner and Shmida 1981; van Rheede van Oudtshoorn and van Rooyen 2013). One such mechanism, commonly observed in water-controlled ecosystems, is ombrohydrochory, the dispersal of seeds caused by an opening of the seed container due to contact with water (Parolin 2006; Navarro et al. 2009; van Rheede van Oudtshoorn and van Rooyen 2013). One particular form, exhibited by members of the Aizoaceae family in semi-arid regions of the Sahel, Australia and South America, is ballistic dispersal, which uses the kinetic energy of raindrops to expulse the plants’ seeds (Parolin 2006; Friedman et al. 1978). Some semi-arid environments such as those in the Mediterranean are characterised by seasonal fluctuations in their environmental conditions and in particular in their precipitation patterns (Noy-Meir 1973). In combination with processes that allow plants to store diaspores during periods of drought, ombrohydrochory yields a synchronisation of seed dispersal with the beginning of the wet season in such seasonal environments. This synchronisation has, for example, been reported in Mesembryanthemum crystallinum and Mesembryanthemum nodiflorum in Southeastern Spain (Navarro et al. 2009). If seed dispersal strategies different from ombrohydrochory are dominant, most species disperse their seeds during the dry season (Navarro et al. 2009; Shabana et al. 2018).

The seasonal synchronisation of seed dispersal splits the annual life-cycle of a plant population into two distinct stages. During the wet season, seeds germinate, new seedlings emerge and adult plants increase their biomass, but no spatial movement takes place. Seed dispersal only occurs during, or at the end of the dry season, while growth processes are dormant (Baudena and Provenzale 2008). By contrast, most mathematical models for dryland vegetation patterns consist of partial differential equations and thus assume that seed dispersal occurs continuously in time. A widely used approach to account for the temporal structure of the annual life cycle is the use of integrodifference equations. This splits the system into 2 distinct, non-overlapping phases, which are both described as discrete, instantaneous processes: a growth phase during which dispersal processes are either not present or negligible and a dispersal phase during which no growth occurs. The application of integrodifference equations to biological and ecological systems in which spatial dispersal plays a significant role was in part pioneered by Kot and Schaffer (1986), and has become a well-established tool in the description of biological and ecological systems since then [e.g. by Musgrave and Lutscher (2014a), Musgrave and Lutscher (2014b), Powell and Zimmermann (2004), Clark et al. (2003), Neubert et al. (1995)].

The spatial and temporal scales associated with the evolution of vegetation patterns do not allow their recreation in laboratory settings. Instead, a range of mathematical models have been proposed to address different aspects of the pattern dynamics (Borgogno et al. 2009; Zelnik et al. 2013). A significant amount of modelling work is based on systems of partial differential equations, most notably by Gilad et al. (2004), HilleRisLambers et al. (2001); Rietkerk et al. (2002) and Klausmeier (1999). The reaction-advection-diffusion Klausmeier model (Klausmeier 1999) is a deliberately basic description of dryland ecosystems based on the vegetation-infiltration feedback loop. Its relative simplicity provides a rich framework for model analyses and extensions [e.g. the work by Bennett and Sherratt (2018), Consolo et al. (2019), Eigentler and Sherratt (2018), Eigentler and Sherratt (2019), Siteur et al. (2014a), Ursino and Contarini (2006), Sherratt (2010), Sherratt (2011), Sherratt (2013a), Sherratt (2013b), Sherratt (2013c)]. The recent development of new remote sensing technology, using temporal sequences of satellite images, allows for comparisons between model predictions and field data (Bastiaansen et al. 2018; Gandhi et al. 2018).

In the Klausmeier model, seed dispersal is modelled by a diffusion term. In reality, the dispersal of seeds is affected by nonlocal processes, such as ballistic dispersal or long range dispersal (e.g. via mammals or wind) (Pueyo et al. 2008; Bullock et al. 2017). The Klausmeier model has been extended to account for such nonlocal processes (Eigentler and Sherratt 2018; Bennett and Sherratt 2018) and a similar approach has been applied to other models for dryland vegetation (Baudena and Rietkerk 2013; Pueyo et al. 2008, 2010). Integrodifference systems also provide a description of nonlocal dispersal effects through a convolution of the plant density with a kernel function. The kernel function is a probability density function describing the average distribution of seeds dispersed from a single plant. The dispersal kernel’s properties (in particular its shape and standard deviation) depend on both plant species and environmental conditions (Bullock et al. 2017).

In this paper we address the significance of seed dispersal synchronisation and its temporal separation from growth processes in seasonal dryland environments. To do so, we introduce an integrodifference model describing the plant–water dynamics in semi-arid ecosystems in Sect. 2. We base our model on the Klausmeier model, to compare our results to previous model analyses of models with no temporal structure. To aid comparisons to the PDE model, we review the most relevant results for the Klausmeier model in Sect. 2. Even though an integrodifference model cannot explicitly take into account the length of the plant growth stage, a consistency result (Proposition 1) yields a control parameter for the temporal separation of seed dispersal events through an appropriate parameter setting. In Sect. 3 we focus on this special case and perform a linear stability analysis to determine a condition for pattern onset in the model and investigate this condition under variations in the growth season length. The analytical derivation of this condition relies on a specific (but nevertheless biologically relevant) choice of the dispersal kernels. To relax this assumption we perform numerical simulations in Sect. 4 to determine the parameter region in which pattern onset occurs for other biologically relevant dispersal kernels. Finally, we discuss our results in Sect. 5.

2 The models

In this section we introduce the integrodifference model which we use to investigate the effects of seasonal synchronisation of seed dispersal on the onset of vegetation patterns in semi-arid environments. The model is based on the reaction-advection-diffusion model by Klausmeier (1999) and to facilitate the comparison of our results on the discrete model to that of the time-continuous model, we start by reviewing relevant results for the Klausmeier model. We relate the models through a convergence result that shows that the Klausmeier PDE model can be obtained from the integrodifference model in an appropriate limit.

2.1 Klausmeier model

One of the well-established models describing vegetation patterns in semi-arid environments is the Klausmeier model (Klausmeier 1999). It reduces the plant–water dynamics to a small set of basic processes (rainfall, plant mortality, evaporation/drainage, vegetation-infiltration feedback and spatial dispersal). The relative simplicity of this modelling approach provides a framework for a rich mathematical analysis [e.g. by Sherratt (2005), Sherratt and Lord (2007), Sherratt (2010), Sherratt (2011), Sherratt (2013a), Sherratt (2013b), Sherratt (2013c), Siteur et al. (2014a), Ursino and Contarini (2006)]. Suitably nondimensionalised (Klausmeier 1999; Sherratt 2005), the model is

$$\begin{aligned} \frac{\partial u}{\partial t}&= \overbrace{u^2w}^{\text {plant growth}} - \overbrace{Bu}^{\text {plant mortality}} + \overbrace{\frac{\partial ^2 u}{\partial x^2}}^{\text {plant dispersal}}, \end{aligned}$$
(1a)
$$\begin{aligned} \frac{\partial w}{\partial t}&= \underbrace{A}_{\text {rainfall}}-\underbrace{w}_{\begin{array}{c} \text {evaporation}\\ \text {and drainage} \end{array}}-\underbrace{u^2w}_{\begin{array}{c} \text {water consumption}\\ \text {by plants} \end{array}} + \underbrace{\nu \frac{\partial w}{\partial x}}_{\begin{array}{c} \text {water flow}\\ \text {downhill} \end{array}}+ \underbrace{d\frac{\partial ^2w}{\partial x^2}}_{\text {water diffusion}}. \end{aligned}$$
(1b)

Here u(x, t) denotes the plant density, w(x, t) the water density, \(x\in {\mathbb {R}}\) the space domain where x is increasing in the uphill direction and \(t>0\) the time. Originally, the model only focussed on a sloped spatial domain, but the addition of a water diffusion term to account for the possibility of a description on flat terrain is a well established addition (Kealy and Wollkind 2012; Siteur et al. 2014a; van der Stelt et al. 2013; Zelnik et al. 2013). To emphasise on the description of seed dispersal as a local process, we refer to this model as the “local Klausmeier model” throughout the paper. Water input to the system is assumed to occur at a constant rate, evaporation and drainage effects are proportional to the water density (Rodriguez-Iturbe et al. 1999; Salvucci 2001) and the plant mortality rate is density-independent. The nonlinearity in the description of water uptake and plant growth processes arises due to a soil modification by plants. The term is the product of the density of the consumer u and of the available resource uw, the amount of water that is able to infiltrate into soil layers where plant roots consume water. The dependence on the plant density u in the latter term occurs due to a positive correlation between the plant density and the soil surface’s permeability (Rietkerk et al. 2000; Valentin et al. 1999; Cornet et al. 1988). Finally, plant growth is assumed to be proportional to the amount of consumed water (Rodriguez-Iturbe et al. 1999; Salvucci 2001). The parameters A, B, \(\nu \) and d are combinations of different dimensional parameters but can be interpreted as rainfall, plant loss, the slope and water diffusion, respectively.

In a previous paper (Eigentler and Sherratt 2018) we have introduced nonlocal seed dispersal effects to the model by replacing the plant diffusion term by a convolution of a dispersal kernel (a probability density function) \(\phi \) and the plant density u. The resulting model is referred to as the “nonlocal Klausmeier model” and is

$$\begin{aligned} \frac{\partial u}{\partial t}&= u^2w - Bu + C\left( \phi (\cdot ) *u(\cdot ,t) - u(x,t)\right) , \end{aligned}$$
(2a)
$$\begin{aligned} \frac{\partial w}{\partial t}&= A-w-u^2w + \nu \frac{\partial w}{\partial x}+ d\frac{\partial ^2w}{\partial x^2}. \end{aligned}$$
(2b)

The additional parameters C and a represent the rate of plant dispersal and reciprocal width of the dispersal kernel, respectively. Note that the convolution \((\phi *u)(x,t)\) accounts for all plant biomass dispersed to the space point x, including the fraction of biomass that is not dispersed. The final term in (2a) ensures that the total biomass over the whole domain remains unchanged by the seed dispersal term. The nonlocal model (2) and the local model (1) are related through a convergence result. If the dispersal kernel \(\phi \) is decaying exponentially as \(|x|\rightarrow \infty \), then the local model (1) can be obtained from the nonlocal model (2) in the limit \(C\rightarrow \infty \) and \(\sigma \rightarrow 0\) with \(C=2/\sigma ^2\), where \(\sigma \) denotes the standard deviation of \(\phi \) (Eigentler and Sherratt 2018).

Linear stability analysis of both the local and the nonlocal Klausmeier model with the Laplace kernel

$$\begin{aligned} \phi (x) = \frac{a}{2}e^{-a|x|}, \quad a>0, x\in {\mathbb {R}}. \end{aligned}$$
(3)

provides analytically derived conditions for pattern onset to occur in the system. On flat ground, i.e. \(\nu = 0\), Turing-type patterns form due to a diffusion-driven instability, i.e. there exists a threshold \(d_c>0\) on the diffusion coefficient such that an instability occurs for all \(d>d_c\). In the local model (1), the threshold is

$$\begin{aligned} d_c(A,B) = \frac{8B\sqrt{-A^2+A\sqrt{A^2-4B^2}+4B^2}-2A^2+2A\sqrt{A^2-4B^2}+16B^2}{B\left( A-\sqrt{A^2-4B^2}\right) ^2}. \end{aligned}$$
(4)

A corresponding threshold \(\widetilde{d_c}(A,B,C,a)\) for the nonlocal model (2) with the Laplace kernel (3) can be derived explicitly, but it omitted due to its algebraic complexity.

On sloped ground (\(\nu \ne 0\)) pattern onset has been studied close to a Turing-Hopf bifurcation, which is characterised by an upper bound on the rainfall parameter A that has been derived analytically valid to leading order in \(\nu \) as \(\nu \rightarrow \infty \) for both models (Eigentler and Sherratt 2018; Sherratt 2013b). The calculation of this upper bound on the precipitation parameter for the nonlocal model with the Laplace kernel shows that long range dispersal of seeds inhibits the formation of patterns by decreasing the size of the parameter region that supports the onset of patterns. On flat ground an increase of the dispersal kernel’s standard deviation causes an increase in the threshold on the diffusion coefficient, while on sloped ground an increase in the dispersal kernel’s width inhibits the formation of patterns by decreasing the upper bound on the rainfall parameter.

The analytical derivation of pattern onset conditions in the nonlocal model is facilitated by the simple algebraic form of the Laplace kernel’s Fourier transform and the associated polynomial structure of the dispersion relation in the linear stability analysis. For other biologically relevant seed dispersal kernels, conditions for pattern onset are not analytically tractable. Numerical simulations, however, confirm the qualitative trends obtained for the model with the Laplace kernel. Simulations further suggest that the dispersal kernel’s decay at infinity has an influence on the rainfall threshold. For narrow dispersal kernels, those that account for more rare long-range dispersal events (algebraic decay rather than exponential) have an inhibitory effect on the formation of patterns, while for sufficiently wide kernels those that decay algebraically at infinity promote pattern formation compared to exponentially decaying kernels.

2.2 Integrodifference model

Integrodifference models are a common type of model widely used in the description of systems in which dispersal processes are temporally separated from other dynamics such as growth/birth and decay/death. To account for the separation of plant growth and seed dispersal stages in dryland ecosystems, we propose the integrodifference model

$$\begin{aligned} u_{n+1}(x)&= C \phi *f(u_n,w_n), \end{aligned}$$
(5a)
$$\begin{aligned} w_{n+1}(x)&= D \phi _1 *g(u_n,w_n), \end{aligned}$$
(5b)

where

$$\begin{aligned} f\left( u,w\right)&= u^2w-Bu + \frac{1}{C}u, \\ g(u,w)&= A - u^2w- w + \frac{1}{D}w. \end{aligned}$$

Here \(u_n(x)\) denotes the plant density, \(w_n(x)\) the water density after \(2n, n\in {\mathbb {N}}\) seasons and location \(x\in {\mathbb {R}}\), where x increases in the uphill direction. The formulation of the model splits the processes involved into two phases: a growth and evolution phase described by the functions f(u, w) and g(u, w) during which no dispersal occurs, and a dispersal phase modelled as a convolution of the evolved densities with dispersal kernels. As in the nonlocal Klausmeier model (2), the plant dispersal kernel \(\phi \) is symmetric and represents isotropic dispersal of plants. To model the flow of water downhill, the water dispersal kernel \(\phi _1\) is in general asymmetric with mean \(\mu _{\phi _1}\le 0\). The special case of a symmetric kernel \(\phi _1\) corresponds to the model on flat ground, which is the main aspect of the study in this paper. The model is based on the Klausmeier models (2) and (1) and thus the functions f(u, w) and g(u, w) consist of the terms describing the rate of change in the original model, appropriately scaled by the coefficients C and D to reflect the time between steps in the discrete model, added to the existing densities.

As the integrodifference model (5) arises directly from the local Klausmeier model (1), the two models can be linked through a consistency result in an appropriate limit which shows that the integrodifference model (5) tends to the local Klausmeier model (1) as \(T\rightarrow 0\). To show this, we consider the parameter setting

$$\begin{aligned} C=T, \quad \sigma _\phi ^2=2T, \quad D=T, \quad \mu _{\phi _1} = -\nu T, \quad {\widetilde{\sigma }}_{\phi _1}^2 = 2dT, \end{aligned}$$
(6)

where \(\mu \) and \(\sigma \) denote the mean and standard deviation of the respective kernels and \( {\widetilde{\sigma }}_{\phi _1}^2 = \int _{-\infty }^\infty \phi _1(x)x^2 {{\text {d}}}x\), the second raw moment of the kernel function \(\phi _1\). Further, we define operators \(P, P_T: C^\infty ({\mathbb {R}}\times [0,\infty ), [0,\infty )^2) \rightarrow C^\infty ({\mathbb {R}}\times [0,\infty ), [0,\infty )^2)\) by

$$\begin{aligned} P\varvec{v}(x,t) = \frac{\partial \varvec{v}}{\partial t} (x,t) - \Gamma \varvec{v}(x,t) - h_1(\varvec{v}(x,t)), \end{aligned}$$
(7)

for any function \(\varvec{v}(x,t) = (u(x,t),w(x,t)) \in C^\infty ({\mathbb {R}}\times [0,\infty ), [0,\infty )^2)\), where

$$\begin{aligned} \Gamma = {\text {diag}}\left( \frac{\partial ^2}{\partial x^2}, \nu \frac{\partial }{\partial x} + d\frac{\partial ^2}{\partial x^2} \right) , \quad h_1\left( \varvec{v}\right) = \left( \begin{array}{c} u^2w-Bu\\ A-u^2w-w\\ \end{array} \right) , \end{aligned}$$

and

$$\begin{aligned} P_T\varvec{v}(x,t) = \frac{1}{T} \left( \varvec{v}(x,t+T) - h_2(\varvec{v}(x,t))\right) , \end{aligned}$$
(8)

where

$$\begin{aligned} h_2\left( \varvec{v}(x,t)\right) = \left( \begin{array}{c} -C \phi (\cdot ) *f(u(\cdot ,t),w(\cdot ,t))\\ - D \phi _1(\cdot ) *g(u (\cdot ,t),w(\cdot ,t))\\ \end{array} \right) . \end{aligned}$$

Note that the operator P arises from the local Klausmeier model (1), because \(P\varvec{v} = 0\) for any \(\varvec{v}\) that satisfies (1). Similarly, \(P_T\) represents the integrodifference model (5), because a sequence \(\varvec{v}_n(x) = \varvec{v}(x,nT)\) satisfies (5) if \(P_T\varvec{v}_n = 0\) for all \(n\in {\mathbb {N}}\). Utilising this reformulation of both models, it is possible to show the following result.

Proposition 1

Consider the parameter setting (6) and let the kernel functions \(\phi \) and \(\phi _1\) have finite moments of all orders and decay exponentially as \(|x| \rightarrow \infty \). Then the integrodifference model (5) is consistent with the local Klausmeier model (1), i.e.

$$\begin{aligned} P\varvec{v} - P_T\varvec{v} \rightarrow 0 \quad \text {as} \quad T\rightarrow 0^+, \end{aligned}$$

for any \(\varvec{v} \in C^\infty ({\mathbb {R}}\times [0,\infty ), [0,\infty )^2)\).

In other words, the model equations (5) converge to the model equations (1) as \(T\rightarrow 0^+\). The notion of consistency is widely used in the field of numerical analysis, and crucially it does not imply convergence of model solutions. While we are unable to construct an argument to prove convergence, numerical simulations suggest that solutions of the integrodifference model (5) converge to solutions of the local Klausmeier model (1) in the parameter setting (6) as \(T\rightarrow 0^+\) (Fig. 2).

On sloped ground Proposition 1 requires that \(\nu = o(T^{-1})\), so that \(T \nu \rightarrow 0\) as \(T\rightarrow 0^+\) and \(\nu \rightarrow \infty \), to facilitate any asymptotic analysis in \(\nu \) similar to that of the local Klausmeier model (Sherratt 2005, 2010, 2011, 2013a, b, c). On flat ground, \(\phi _1\) is symmetric and thus \(\mu _{\phi _1}=0\) and \({\widetilde{\sigma }}_{\phi _1}\) coincides with the kernel’s standard deviation \(\sigma _{\phi _1}\).

The parameter T can be interpreted as the time between separate dispersal events and the scalings (6) are thus the main focus of the model’s analysis in Sect. 3. While the time between two seed dispersal events in a seasonal environment is usually fixed, we are interested in variations of T as this parameter establishes a connection between the local Klausmeier model (1) and the integrodifference model (5). In particular, as \(T\rightarrow 0^+\) in the model, the length of each season tends to zero. As a consequence, this limit corresponds to the disappearance of any seasonality in the model and all processes are assumed to occur continuously in time, as, for example, in the Klausmeier model (1).

One kernel function satisfying the conditions in Proposition 1 is the Laplacian kernel (3). We define the corresponding asymmetric Laplace kernel by \(\phi _1(x) = Ne^{-a_2x}\) for \(x\ge 0\) and \(\phi _1(x) = Ne^{(a_2-a_1)x}\) for \(x < 0\), where \(N = (a_2-a_1)a_2/(2a_2-a_1)\) and \(a_2>a_1>0\). The parameter \(a_1\) controls the extent of the asymmetry of the kernel function and \(a_1=0\) yields the symmetric Laplace kernel (3). The model with this particular kernel function is studied in some detail in this paper as the Fourier transform of the symmetric Laplacian kernel \({\widehat{\phi }}(k) = a^2/(a^2+k^2)\) provides a significant simplification in the analysis of pattern onset.

Proof of Proposition 1

Firstly, we show that

$$\begin{aligned} P_T\varvec{v}(x,t) = \frac{\varvec{v}(x,t+T) - \varvec{v}(x,t)}{T} - \Gamma \varvec{v}(x,t) - h_1(\varvec{v}(x,t)) + O\left( T^2\right) . \end{aligned}$$
(9)

To this end, we define \(\phi (x) = \sigma _{\phi }^{-1}\varphi (\sigma _{\phi }^{-1}x)\) and \(\phi _1(x) = {\widetilde{\sigma }}_{\phi _1}^{-1}\varphi _1({\widetilde{\sigma }}_{\phi _1}^{-1}x)\) Under the changes of variables \(y=x-\sigma _{\phi }z\) and \(y=x-{\widetilde{\sigma }}_{\phi _1}z\), respectively, \(P_T \varvec{v} = ((P_T\varvec{v})_1,(P_T\varvec{v})_2)\) satisfies

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _1 = u(x,t+T) - C\int _{-\infty }^\infty \varphi (z) f\left( u\left( x-\sigma _{\phi }z,t \right) , w\left( x-\sigma _{\phi }z,t \right) \right) {{\text {d}}}z , \end{aligned}$$
(10a)
$$\begin{aligned}&\begin{aligned} T\left( P_T\varvec{v}\right) _2&= w(x,t+T) - D\int _{-\infty }^\infty \varphi _1(z) g\left( u\left( x-{\widetilde{\sigma }}_{\phi _1}z,t \right) , w\left( x-{\widetilde{\sigma }}_{\phi _1}z,t \right) \right) {{\text {d}}}z. \end{aligned} \end{aligned}$$
(10b)

Due to the parameter setting (6), small values of T correspond to small values of \(\sigma _{\phi }\) and \({\widetilde{\sigma }}_{\phi _1}\) Hence, to investigate the system’s behaviour for \(T\ll 1\), consider the Taylor expansions of \(u(x-\sigma _{\phi }z,t)\), \(w(x-\sigma _{\phi }z,t)\), \(u(x-{\widetilde{\sigma }}_{\phi _1}z,t)\) and \(w(x-{\widetilde{\sigma }}_{\phi _1}z,t)\) about x, which give

$$\begin{aligned}&f\left( u\left( x-\sigma _{\phi }z,t \right) , w\left( x-\sigma _{\phi }z,t \right) \right) \nonumber \\&\quad = u(x,t)^2w(x,t)+\left( \frac{1}{C}-B\right) u(x,t) \nonumber \\&\qquad -\,\sigma _{\phi }z \left( u(x,t)^2 w_x(x,t) + \left( \frac{1}{C}-B\right) u_x(x,t) + 2u(x,t)u_x(x,t)w(x,t)\right) \nonumber \\&\qquad +\,\sigma _{\phi }^2z^2 \left( \frac{1}{2}u(x,t)^2 w_{xx}(x,t) + u_x(x,t)^2w(x,t) + \frac{1}{2}\left( \frac{1}{C}-B\right) u_{xx}(x,t) \right. \nonumber \\&\qquad \left. +\,u(x,t)u_{xx}(x,t)w(x,t) + 2u(x,t)u_x(x,t)w_x(x,t) \right) + O\left( \sigma _{\phi }^3\right) , \end{aligned}$$
(11)

and similarly

$$\begin{aligned}&g\left( u\left( x-{\widetilde{\sigma }}_{\phi _1}z \right) , w\left( x-{\widetilde{\sigma }}_{\phi _1}z \right) \right) \nonumber \\&\quad = A - u(x,t)^2w(x,t) + \left( \frac{1}{D} -1\right) w(x,t) \nonumber \\&\qquad -\, {\widetilde{\sigma }}_{\phi _1}z \left( -u(x,t)^2w_x(x,t) + \left( \frac{1}{D} -1\right) w_x(x,t) - 2u(x,t)u_x(x,t)w(x,t)\right) \nonumber \\&\qquad +\, {\widetilde{\sigma }}_{\phi _1}^2z^2 \left( -u_x(x,t)^2w(x,t) - \frac{1}{2}u(x,t)^2w_{xx}(x,t) + \frac{1}{2}\left( \frac{1}{D} -1\right) w_{xx}(x,t) \right. \nonumber \\&\qquad \left. -\,u(x,t)u_{xx}(x,t)w(x,t) - 2u(x,t)u_x(x,t)w_x(x,t) \right) + O\left( {\widetilde{\sigma }}_{\phi _1}^3\right) , \end{aligned}$$
(12)

where the subscripts of u and w denote partial differentiation. Substitution of this into (10) and term-wise integration using Watson’s Lemma [e.g. (Miller 2006)] gives

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _1 = u(x,t+T) - C\left( u(x,t)^2w(x,t) + \left( \frac{1}{C}-B\right) u(x,t) \right. \\&\quad \left. +\, \left( u(x,t)^2w_{xx}(x,t) + 2(u_{x}(x,t))^2w(x,t) + \left( \frac{1}{C}-B\right) u_{xx}(x,t) \right. \right. \\&\quad \left. \left. +\,2u(x,t)u_{xx}(x,t)w(x,t) + 4u(x,t)u_{x}(x,t)w_{x}(x,t) \right) \sigma _{\phi }^2 \int _{-\infty }^\infty \varphi (z)z^2{{\text {d}}}z + O\left( \sigma _{\phi }^3\right) \right) , \end{aligned}$$

and

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _2 = w(x,t+T) \\&\quad -\, D\left( 2\left( A-u(x,t)^2w(x,t) +\left( \frac{1}{D}-1\right) w(x,t)\right) \right. \left. \int _{-\infty }^\infty \varphi _1(z){{\text {d}}}z \right. \\&\quad \left. +\, \left( u(x,t)^2w_{x}(x,t)-\left( \frac{1}{D} -1\right) w_{x}(x,t) +2u(x,t)u_{x}(x,t)w(x,t)\right) {\widetilde{\sigma }}_{\phi _1}\int _{-\infty }^\infty \varphi _1(z)z{{\text {d}}}z \right. \\&\quad \left. +\, \left( -2(u_{x}(x,t))^2w(x,t) - u(x,t)^2w_{xx}(x,t) +\left( \frac{1}{D}-1\right) w_{xx}(x,t) \right. \right. \\&\quad \left. \left. -\,2u(x,t)u_{xx}(x,t)w(x,t) -4u(x,t)u_{x}(x,t)w_{x}(x,t) \right) \frac{{\widetilde{\sigma }}_{\phi _1}^2}{2}\int _{-\infty }^\infty \varphi _1(z)z^2{{\text {d}}}z + O\left( {\widetilde{\sigma }}_{\phi _1}^3\right) \right) . \end{aligned}$$

Using that \(\varphi (x) = \sigma _{\phi }\phi (\sigma _{\phi }x)\), \(\varphi _1(x) = {\widetilde{\sigma }}_{\phi _1}\phi _1({\widetilde{\sigma }}_{\phi _1}x)\), and the definition of the moments of a probability distribution give

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _1 = u(x,t+T) - C \left( u(x,t)^2w(x,t)+ \left( \frac{1}{C}-B\right) u(x,t) + \frac{\sigma _{\phi }^2}{2C}u_{xx}(x,t) \right. \\&\quad \left. +\, \sigma _{\phi }^2 \left( \frac{1}{2}u(x,t)^2w_{xx}(x,t) + (u_{x}(x,t))^2w(x,t) - \frac{1}{2}Bu_{xx}(x,t) \right. \right. \\&\quad \left. \left. +\,u(x,t)u_{xx}(x,t)w(x,t) +2u(x,t)u_{x}(x,t)w_{x}(x,t)\right) + O\left( \sigma _{\phi }^3\right) \right) , \end{aligned}$$

and

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _2 = w(x,t+T) - D\left( A-u(x,t)^2w(x,t) +\left( \frac{1}{D}-1 \right) w(x,t) \right. \\&\quad \left. -\,\frac{\mu _{\phi _1}}{D} w_{x}(x,t) +\frac{{\widetilde{\sigma }}_{\phi _1}^2}{2D}w_{xx}(x,t) +\mu _{\phi _1}\left( u(x,t)^2w_{x}(x,t)+w_{x}(x,t) \right. \right. \\&\quad \left. \left. +\,2u(x,t)u_{x}(x,t)w(x,t) \right) + {\widetilde{\sigma }}_{\phi _1}^2\left( -(u_{x}(x,t))^2w(x,t)-\frac{1}{2}u(x,t)^2w_{xx}(x,t) \right. \right. \\&\quad \left. \left. -\,\frac{1}{2}w_{xx}(x,t) -u(x,t)u_{xx}(x,t)w(x,t) -2u(x,t)u_{x}(x,t)w_{x}(x,t) \right) + O\left( {\widetilde{\sigma }}_{\phi _1}^3\right) \right) . \end{aligned}$$

Applying (6) yields

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _1 = u(x,t+T) \\&\quad -\,\left( u(x,t) + T\left( u(x,t)^2w(x,t)-Bu(x,t)+u_{xx}(x,t)\right) \right) + O\left( T^2\right) , \end{aligned}$$

and

$$\begin{aligned}&T\left( P_T\varvec{v}\right) _2 = w(x,t+T) \\&\quad -\,\left( w(x,t) + T\left( A-u(x,t)^2w(x,t)-w(x,t)+\nu w_x(x,t)+dw_{xx}(x,t) \right) \right) + O\left( T^2\right) , \end{aligned}$$

which shows (9).

The Taylor expansions \(u(x,t+T) = u(x,t) + T u_t (x,t) +O(T^2)\) and \(w(x,t+T) = w(x,t) + T w_t (x,t)+O(T^2)\) yield

$$\begin{aligned} P_T\varvec{v}(x,t) = \frac{\partial \varvec{v}}{\partial t}(x,t) - \Gamma \varvec{v}(x,t) - h_1(\varvec{v}(x,t)) +O\left( T^2\right) , \end{aligned}$$

and thus

$$\begin{aligned} P\varvec{v} - P_T\varvec{v} = O\left( T^2\right) , \end{aligned}$$

which tends to zero as \(T\rightarrow 0\).

\(\square \)

3 Linear stability analysis

A common approach to study the onset of spatial patterns in a model is linear stability analysis. Spatial patterns occur if a steady state that is stable to spatially homogeneous perturbations becomes unstable if a spatially heterogeneous perturbation is introduced. In this section we show that such a linear stability analysis of the integrodifference model (5) on flat ground with the Laplacian kernels in the parameter setting (6) yields a condition for pattern onset that is equivalent to the corresponding condition for the local Klausmeier model (1). This implies that pattern onset is independent of the parameter T, the temporal separation of seed dispersal events.

The steady states of (5) are identical with those of the Klausmeier models (1) and (2), i.e.

$$\begin{aligned} ({\overline{u}}_1,{\overline{w}}_1)&= \left( 0,A\right) , \quad ({\overline{u}}_2,{\overline{w}}_2) = \left( \frac{2B}{A-\sqrt{A^2-4B^2}},\frac{A-\sqrt{A^2-4B^2}}{2}\right) , \\ ({\overline{u}}_3,{\overline{w}}_3)&= \left( \frac{2B}{A+\sqrt{A^2-4B^2}},\frac{A+\sqrt{A^2-4B^2}}{2}\right) . \end{aligned}$$

Existence of \(({\overline{u}}_2,{\overline{w}}_2)\) and \(({\overline{u}}_3,{\overline{w}}_3)\) requires \(A>A_{\min }:=2B\). The steady states are independent of C, D and the dispersal widths a, \(a_1\) and \(a_2\) and are thus independent of frequency changes to the temporal intermittency when using the scalings (6). For the Klausmeier models \(({\overline{u}}_1,{\overline{w}}_1)\) and \(({\overline{u}}_2,{\overline{w}}_2)\) are stable to spatially homogeneous perturbations, while \(({\overline{u}}_3,{\overline{w}}_3)\) is unstable to spatially homogeneous perturbations in the biologically relevant parameter region \(B<2\) (Klausmeier 1999; Sherratt 2005; Eigentler and Sherratt 2018). Preservation of this structure of the steady states in the integrodifference model (5) is only achieved in a certain parameter region.

Proposition 2

If

$$\begin{aligned} D= \ell {\overline{D}}, \ \ell< 1, \quad C=\frac{\ell _1{\overline{D}}}{B(m-\ell _1{\overline{D}})}, \ m > 2, \ell _1<1, \end{aligned}$$
(13)

where

$$\begin{aligned} {\overline{D}} = \frac{2\left( A^2-A\sqrt{A^2-4B^2} - 2B^2\right) }{A^2-A\sqrt{A^2-4B^2}}, \end{aligned}$$
(14)

then \(({\overline{u}}_1,{\overline{w}}_1)\) and \(({\overline{u}}_2,{\overline{w}}_2)\) are stable to spatially homogeneous perturbations, and \(({\overline{u}}_3,{\overline{w}}_3)\) is unstable to spatially homogeneous perturbations.

This condition is sufficient but not necessary. Outside this region further restrictions on the rainfall parameter A can be imposed to guarantee conservation of the steady state structure. In the limiting case (6) such a restriction on the rainfall parameter cannot be avoided. The following condition ensures that (13) holds in the limiting case (6).

Corollary 3

If

$$\begin{aligned} A^2< A_+^2:=\min \left\{ \frac{4B^2}{(2-T)T}, \frac{B(BT+1)^2}{T}\right\} , \quad T<\frac{1}{2}, \quad B<2, \end{aligned}$$
(15)

in (5) with \(C=D=T\), then \(({\overline{u}}_1,{\overline{w}}_1)\) and \(({\overline{u}}_2,{\overline{w}}_2)\) are stable to spatially homogeneous perturbations, and \(({\overline{u}}_3,{\overline{w}}_3)\) is unstable to spatially homogeneous perturbations.

In the limit \(T\rightarrow 0^+\) this becomes the whole A-B parameter region considered for the continuous-time Klausmeier models, providing a reasonable framework for a comparison of the two models. The upper bounds on T and A do, however, introduce a significant restriction on the model as no arbitrarily large time between dispersal events or large precipitation volumes A can be considered. In this, as well as the parameter region given by (13), the plant density \(u_n(x)\) and the water density \(w_n(x)\) remain positive for initial conditions close to the steady states. This is sufficient for the linear stability analysis and simulations that follow. In the parameter region in which \(({\overline{u}}_2,{\overline{w}}_2)\) is unstable, four different behaviours of the system’s solution can be observed; (i) convergence to the desert steady state, (ii) divergence, (iii) a chaotic solution or (iv) a periodic solution for which period doubling occurs as T is increased. However, these different behaviours can yield negative densities of the system’s quantities and are thus not considered further in this paper.

Spatial patterns of (5) arise if the steady state \(({\overline{u}}, {\overline{w}}):=({\overline{u}}_2, {\overline{w}}_2)\), which is stable to spatially homogeneous perturbations, becomes unstable if a spatially heterogeneous perturbation is introduced.

Proposition 4

The steady state \(({\overline{u}}, {\overline{w}})\) is stable to spatially heterogeneous perturbations if \(|\lambda (k)|<1\) for both eigenvalues of the Jacobian

$$\begin{aligned} J = \left( \begin{array}{cc} C{\widehat{\phi }}(k)\alpha &{} C{\widehat{\phi }}(k)\beta \\ D\widehat{\phi _1}(k)\gamma &{} D\widehat{\phi _1}(k)\delta \\ \end{array} \right) , \end{aligned}$$
(16)

for all \(k>0\), where

$$\begin{aligned} \begin{aligned} \alpha&= f_u({\overline{u}},{\overline{w}}) = \frac{BC+1}{C}, \\ \beta&= f_w({\overline{u}},{\overline{w}}) = \frac{4B^2}{\left( A-\sqrt{A^2-4B^2}\right) ^2}, \quad \gamma = g_u({\overline{u}},{\overline{w}}) = -2B, \\ \delta&= g_w({\overline{u}},{\overline{w}}) = -\frac{2\left( A^2D-AD\sqrt{A^2-4B^2} -A^2+A\sqrt{A^2-4B^2} +2B^2\right) }{D\left( A-\sqrt{A^2-4B^2}\right) ^2}. \end{aligned} \end{aligned}$$
(17)

Due to the asymmetry of \(\phi _1\) some of the entries of the Jacobian (16) are complex-valued. A significant simplification can therefore be achieved by considering the integrodifference model (5) on flat ground. This corresponds to \(a_1 =0\) in \(\phi _1\). As a consequence, the Jury conditions [see e.g. Murray (1989)] can be used to determine the steady state’s stability to spatially heterogeneous perturbations. To study this in more detail, and in particular to show that the model does not provide information on effects the temporal separation of seed dispersal events, we focus on the limiting case (6) and the Laplacian kernel (3).

Proposition 5

The steady state \(({\overline{u}}, {\overline{w}})\) of the integrodifference model (5) under the scalings (6) on flat ground with the Laplacian kernels (3) is unstable to spatially heterogeneous perturbations if

$$\begin{aligned} 1+\det (J)- |{{\text {tr}}}(J)| <0, \quad \text {for some } k>0, \end{aligned}$$
(18)

where J is the Jacobian given in Proposition 4 with \(a_1=0\).

In other words, Proposition 5 provides a sufficient condition for spatial patterns to occur. The following proposition shows that (18) is equivalent to the stability condition (4) of \(({\overline{u}},{\overline{w}})\) in the local Klausmeier model. In other words, a diffusion driven instability causes the occurrence of spatial patterns in the integrodifference model, i.e. given a level of rainfall A, an instability occurs for \(d>d_c(A,B)\), where \(d_c(A,B)\) is given in (4).

Proposition 6

The steady state \(({\overline{u}}, {\overline{w}})\) of the integrodifference model (5) under the scalings (6) on flat ground with the Laplacian kernels (3) is unstable to spatially heterogeneous perturbations if \(d>d_c(A,B)\), where the threshold \(d_c\) is identical with the corresponding threshold (4) for the local Klausmeier model.

The condition’s independence of T yields that the integrodifference model does not provide any information on the effects of the temporal separation of seed dispersal events on the onset of spatial patterns. The equivalence of the condition to that of the local Klausmeier model follows directly from the condition’s independence of T and Proposition 1, which shows that the integrodifference model converges to the local Klausmeier model as \(T\rightarrow 0^+\). Thus for sufficiently small values of T, Proposition 6 does indeed provide the exact same information as the diffusion threshold obtained for the local Klausmeier model. For larger T the model does not provide any information on the transition between uniform and patterned vegetation as the decrease in the upper bound \(A_+\) on the rainfall parameter reduces the size of the rainfall interval for which the derivation of \(d_c\) is valid.

Proof of Proposition 2

Stability of a steady state \(({\overline{u}},{\overline{w}})\) is determined by the Jury conditions applied to the Jacobian

$$\begin{aligned} J = \left( \begin{array}{cc} C(2{\overline{u}}{\overline{w}}-B)+1 &{} C{\overline{u}}^2 \\ -2D{\overline{u}}{\overline{w}} &{} -D({\overline{u}}^2+1)+1 \\ \end{array} \right) . \end{aligned}$$

The steady state \(({\overline{u}}_3,{\overline{w}}_3)\) is unstable in the whole parameter region, because

$$\begin{aligned} 1+\det (J)-|{{\text {tr}}}(J)| = -\frac{2BCD\left( A^2+A\sqrt{A^2-4B^2}-4B^2 \right) }{\left( A+\sqrt{A^2-4B^2}\right) ^2} <0. \end{aligned}$$

The desert steady state \(({\overline{u}}_1,{\overline{w}}_1)\) is monotonically stable if \(C<B^{-1}\) and \(D<1\). If \(1<D<2\) or \(B^{-1}<C<2B^{-1}\) it is still stable but solutions are oscillating about (0, A), which is biologically impossible. Finally, the Jury conditions yield that \(({\overline{u}}_2,{\overline{w}}_2)\) is stable to spatially homogeneous perturbations if \(\min \{\overline{C_2}, \overline{C_3}\}< C < \overline{C_1}\), where

$$\begin{aligned} \overline{C_1}&= \frac{AD\left( A-\sqrt{A^2-4B^2}\right) }{B\left( (D-1)A\left( \sqrt{A^2-4B^2}-A\right) +2B^2(2D-1) \right) }, \\ \overline{C_2}&= \frac{2\left( (D-2)\left( A\sqrt{A^2-4B^2}-A^2\right) -4B^2\right) }{B\left( (D-2)\left( A^2-A\sqrt{A^2-4B^2} \right) -4B^2(D-1) \right) }, \\ \overline{C_3}&= \frac{(D-2)\left( A^2-A\sqrt{A^2-4B^2}\right) +4B^2}{B\left( A^2-A\sqrt{A^2-4B^2}-2B^2\right) }. \end{aligned}$$

Combined, this gives that the steady state structure of the continuous time model is preserved if

$$\begin{aligned} D<1 \quad \text {and} \quad \max \left\{ 0,\min \left\{ \overline{C_2},\overline{C_3} \right\} \right\}< C < \min \left\{ \frac{1}{B}, \overline{C_1} \right\} . \end{aligned}$$
(19)

If \(D > 1/2\), then \(\min \left\{ 1/B, \overline{C_1} \right\} = 1/B\), because

$$\begin{aligned} \overline{C_1} - \frac{1}{B} = -\frac{2\left( D-\frac{1}{2}\right) \left( A^2-A\sqrt{A^2-4B^2}-2B^2\right) }{B\left( \left( D-\frac{1}{2} \right) \left( A^2-A\sqrt{A^2-4B^2}-4B^2\right) -\left( A^2-A\sqrt{A^2-4B^2}\right) \right) }>0, \end{aligned}$$

since \(A^2-A\sqrt{A^2-4B^2}-2B^2>0\) and \(A^2-A\sqrt{A^2-4B^2}-4B^2<0\). Similarly, if \(D < 1/2\), then \(\min \left\{ 1/B, \overline{C_1} \right\} = \overline{C_1}\). Further, if \(D < {\overline{D}}\) (defined in (14)), then \(\max \{0,\min \{\overline{C_2},\overline{C_3} \} \} = 0\) and similarly, if \(D > {\overline{D}}\), then \(\max \{0,\min \{\overline{C_2},\overline{C_3} \} \} = \overline{C_2}\).

Hence, (19) can be simplified by splitting it into different parameter regions. It becomes (i) \(C<\overline{C_1}\) if \(D<1/2\) and \(D<{\overline{D}}\), (ii) \(\overline{C_2}<C<\overline{C_1}\) if \(D<1/2\) and \({\overline{D}}<D<1\), (iii) \(C<1/B\) if \(1/2<D<1\) and \(D<{\overline{D}}\) and (iv) \(\overline{C_2}<C<1/B\) if \(1/2<D<1\) and \({\overline{D}}<D<1\).

This classification is used below to show that if C and D are defined as in (13), then (19) is satisfied in the whole parameter plane that is considered in the continuous-time PDE models (\(A>2B\), \(B<2\)). To show this it is sufficient to show that (i) and (iii) are satisfied because \(\ell <1\). For case (iii) note that

$$\begin{aligned} C=\frac{\ell _1D}{B(m-\ell _1D)}<\frac{1}{B} \iff \ell _1D<\frac{m}{2}, \end{aligned}$$

which is satisfied since \(\ell _1D<1\) and \(m>2\). For case (i) note that

$$\begin{aligned} C=\frac{\ell _1D}{B(m-\ell _1D)}< \overline{C_1} \iff \ell _1D < \frac{2B^2+(m-1)\left( A^2-A\sqrt{A^2+4B^2}\right) }{4B^2} :=\overline{{\overline{D}}}. \end{aligned}$$

This is always satisfied because \(\ell _1D<1\) and

$$\begin{aligned} \overline{{\overline{D}}}>1 \iff m>\frac{2B^2}{A^2-A\sqrt{A^2+4B^2}} +1 :={\overline{m}}, \end{aligned}$$

which holds true since \(m>2\) and

$$\begin{aligned} {\overline{m}} <2 \iff A^2-A\sqrt{A^2-2B^2}-2B^2 > 0, \end{aligned}$$

which is clearly satisfied.

\(\square \)

Proof of Proposition 4

Linearisation of the model (5) about the steady state \(({\overline{u}}, {\overline{w}})\) gives \(u_{n+1}(x) = C \phi (\cdot )*(\alpha u_n(\cdot ) + \beta w_n(\cdot ))\) and \(w_{n+1}(x) = D \phi _1(\cdot )*(\gamma u_n(\cdot ) + \delta w_n(\cdot ))\). Taking the Fourier transform of both equations yields \({\widehat{u}}_{n+1}(k) = C {\widehat{\phi }}(k) (\alpha {\widehat{u}}_n(k) + \beta {\widehat{w}}_n(k))\) and \({\widehat{w}}_{n+1}(k) = D \widehat{\phi _1}(k)(\gamma {\widehat{u}}_n(k) + \delta {\widehat{w}}_n(k))\), where \({\widehat{\phi }}\) and \(\widehat{\phi _1}\) denote the Fourier transforms of the kernels \(\phi \), and \(\phi _1\), respectively. Under the assumption that \({\widehat{u}}_n(k)\) and \({\widehat{w}}_n(k)\) are proportional to \(\lambda ^n{\tilde{u}}(k)\) and \(\lambda ^n{\tilde{w}}(k)\), respectively, where \(\lambda \in {\mathbb {C}}\) denotes the growth rate, the system becomes \(\lambda {\tilde{u}}(k) = C {\widehat{\phi }}(k) (\alpha {\tilde{u}}(k) + \beta {\tilde{w}}(k))\) and \(\lambda {\tilde{w}}(k) = D \widehat{\phi _1}(k)(\gamma {\tilde{u}}(k) + \delta {\tilde{w}}(k))\), i.e. \(\lambda \) is an eigenvalue of the Jacobian J. \(\square \)

Proof of Proposition 5

For an instability to occur, at least on of the Jury conditions \(\det (J) <1\) and \(1+\det (J) - |{{\text {tr}}}(J)|>0\) needs to be violated for some wavenumber \(k>0\). The former condition is satisfied for all \(k>0\). To show this, note that that \(\max \{\det (J) -1\}\) is at \(k=0\) because

$$\begin{aligned} \det (J) - 1 = \frac{\alpha _4 k^4 + \alpha _2 k^2 + \alpha _0}{ \left( d T k^2+1\right) \left( A-\sqrt{A^2-4B^2}\right) ^2\left( T k^2 +1\right) }, \end{aligned}$$
(20)

where

$$\begin{aligned} \alpha _4&= 2d T^2 \left( -A^2+A\sqrt{A^2-4B^2}+2B^2\right) , \\ \alpha _2&= -2T\left( A^2-A\sqrt{A^2-4B^2}-2B^2\right) (d+1), \\ \alpha _0&= 2T\left( \left( \frac{1}{2}B-1\right) \left( A^2-A\sqrt{A^2-4B^2}\right) \right. \\&\left. + \left( \frac{1}{2}B-TB\right) \left( A^2-A\sqrt{A^2-4B^2}-4B^2\right) \right) . \end{aligned}$$

The denominator of (20) is clearly positive and increasing for \(k>0\). Since further \(\alpha _4 <0\) and \(\alpha _2<0\), the numerator and thus the whole of (20) is decreasing for \(k>0\) and it attains its maximum at \(k=0\). The negativity of (20) then follows from that of \(\alpha _0\) which follows from \(B<2\) and \(T<1/2\).

\(\square \)

Proof of Proposition 6

Firstly, we note that \(\partial d_c / \partial A \ge 0\) for all \(A\ge 2B\). Hence, \(d_c\) attains its minimum on \(A=2B\), on which it simplifies to \(d_c = 2/B\). Since \(B<2\), \(d_c>1\). Next, we show that \({{\text {tr}}}(J)>0\). To do this, note that

$$\begin{aligned} {{\text {tr}}}(J) = \frac{\beta _2k^2+\beta _0}{ \left( d T k^2+1\right) \left( A-\sqrt{A^2-4B^2}\right) ^2\left( T k^2 +1\right) } >0, \end{aligned}$$

for all \(k>0\), where

$$\begin{aligned} \beta _2&= 2\left( A^2-A\sqrt{A^2-4B^2}-2B^2\right) \left( BT^2d+T+Td\right) -2T^2\left( A^2-A\sqrt{A^2-4B^2}\right) , \\ \beta _0&= 2\left( BT-T+2\right) \left( A^2-A\sqrt{A^2-4B^2}-2B^2\right) . \end{aligned}$$

The denominator is clearly positive and thus the condition for positivity of \({{\text {tr}}}(J)\) is \(\beta _2k^2+\beta _0>0\). The left hand side of this is decreasing in A since \(A^2-A\sqrt{A^2-4B^2}\) is decreasing in A and the assumptions on B and d, and thus obtains its minimum at \(A=A^+\), where \(A^+\) is given in (15). If \(B<1/(2-T)\), then \(A^+ = 4B^2/((2-T)T)\) and

$$\begin{aligned} {{\text {tr}}}(J)\left( \sqrt{A^+}\right)>0 \Longleftrightarrow k^2 > \frac{B}{1-d-BTd}, \end{aligned}$$

since \(d>1\). The right hand side is negative and thus \(\min ({{\text {tr}}}(J)) >0\) for \(B<1/(2-T)\). If \(B>1/(2-T)\), then \(A^+ = (BT+1)^2B/T\) and

$$\begin{aligned} {{\text {tr}}}(J)\left( \sqrt{A^+}\right)>0 \Longleftrightarrow k^2 > -\frac{TB^2+(2-T)B-1}{B^2T^2d+\left( (d+1)T-T^2\right) B-T}, \end{aligned}$$

since \(d>1\). Negativity of the right hand side follows from the lower bound on B and thus \(\min ({{\text {tr}}}(J)) >0\) for all \(B<2\). This shows that \({{\text {tr}}}(J)>0\). The stability condition (18) thus becomes \(1+\det (J)-{{\text {tr}}}(J)< 0 \Longleftrightarrow \gamma _4k^4+\gamma _2k^2+\gamma _0 <0\), where \(\gamma _4 = d(A^2-A\sqrt{A^2-4B^2}-2B^2)\), \(\gamma _2 = (A^2-A\sqrt{A^2-4B^2})(1-Bd)+2B^3d\) and \(\gamma _0 = B(A^2-4B^2)\). This condition and thus its minimum \(-\gamma _2^2/(4\gamma _4)+\gamma _0\) is independent of T. Determining the locus at which the minimum changes sign gives the threshold \(d_c(A,B)\). \(\square \)

4 Simulations

The preceding linear stability analysis relies on the use of the Laplace kernel. For other kernel functions whose Fourier transforms do not provide such a simplification numerical simulations of the model are considered to investigate the onset of patterns. In particular, this allows us to make comparisons between different dispersal kernels, similar to the analysis performed for the nonlocal model in Eigentler and Sherratt (2018). These show that both wide plant dispersal kernels and narrow water dispersal kernels inhibit the formation of patterns. Finally in this section, we show that as for the nonlocal Klausmeier model, the kind of decay of the plant dispersal kernel at infinity is also important.

Simulations are performed on the space domain \([-x_{\max }, x_{\max }]\) centred at \(x=0\). This domain is discretised into M equidistant points \(x_1,\dots ,x_M\) with \(-x_{\max } = x_1< x_2< \dots < x_M = x_{\max }\) such that \(\Delta x = x_2-x_1 = \dots =x_M - x_{M-1}\). On flat ground (5) then becomes

$$\begin{aligned} u_{n+1}(x_k)&= C \Delta x \left( \phi *f_n \right) _k, \end{aligned}$$
(21a)
$$\begin{aligned} w_{n+1}(x_k)&= D\Delta x \left( \phi _1 *g_n \right) _k, \end{aligned}$$
(21b)

where \(\phi , \phi _1\) denote the vectors consisting of the elements obtained by evaluating the corresponding function at each mesh point, \(f_n, g_n\) denote the vectors consisting of the elements obtained by evaluating the corresponding function at each \((u_n(x_k), w_n(x_k))\) and \(z_1 *z_2\) denotes the discrete convolution of two vectors \(z_1\) and \(z_2\). The convolution terms in (21a) and (21b) are obtained by using the convolution theorem and the fast Fourier transform, providing a significant simplification as this reduces the number of operations required to obtain the convolution from \(O(M^2)\) to \(O(M\log (M))\) [see e.g. Cooley et al. (1969)].

To mimic the infinite domain used for the linear stability analysis (Sect. 3), we define the initial condition of the system as follows; on a subdomain \([-x_{{\text {sub}}}, x_{{\text {sub}}} ]\) centred at \(x=0\) of the domain \([-x_{\max }, x_{\max }]\) considered in the simulation the initial condition is a random perturbation of the steady state \(({\overline{u}}, {\overline{w}})\), while on the rest of the domain the densities are initially set to equal the densities of the steady state \(({\overline{u}}, {\overline{w}})\). In other words, \(u_0(x_k) = {\overline{u}} + \delta (x_k)\) and \(w_0(x_k) = {\overline{w}} + \varepsilon (x_k)\) for \(x_k \in [-x_{{\text {sub}}}, x_{{\text {sub}}} ]\), where \(\Vert \delta \Vert _\infty < 0.1{\overline{u}} \) and \(\Vert \varepsilon \Vert _\infty < 0.1{\overline{w}}\) and \(u_0(x_k) = {\overline{u}}\) and \(w_0(x_k) = {\overline{w}}\) for \(x_k \notin [-x_{{\text {sub}}}, x_{{\text {sub}}} ]\). The size of the outer domain is chosen large enough so that any boundary conditions (which are set to be periodic) that are imposed on \([-x_{\max },x_{\max }]\) do not affect the solution in the subdomain in the finite time that is considered in the simulation. Figure 1 shows a typical patterned solution obtained by these simulations.

Based on the amplitude of the oscillation relative to the steady state of the solutions obtained by the simulations we set up a scheme to determine the critical rainfall level \(A_{\max }\) below which pattern onset occurs. Doing this allows us to investigate how certain changes of parameters and kernel functions affect the onset of patterns. Due to the random perturbation of the initial state of the system, all simulation results shown below are the averages taken over 100 simulations. For the symmetric dispersal kernels \(\phi \) and \(\phi _1\) we consider the Laplacian (3), the Gaussian

$$\begin{aligned} \phi _g(x) = \frac{a_g}{\sqrt{\pi }}e^{-a_g^2x^2}, \quad a>0, x\in {\mathbb {R}}, \end{aligned}$$
(22)

and the power law distribution

$$\begin{aligned} \phi _p(x) = \frac{(b-1)a_p}{2\left( 1+a_p|x|\right) ^b}, \quad a>0, b >3, x\in {\mathbb {R}}. \end{aligned}$$
(23)
Fig. 1
figure 1

Simulation of the integrodifference model. This figure shows a patterned solution obtained by simulating the integrodifference model on flat ground. The kernels used in these simulations are the symmetric Laplacian kernels, respectively. The parameter setting (15) with \(T=0.1\) is used in the simulation. The other parameters are \(A=0.9\), \(B=0.45\) and \(d=500\)

Full size image
Fig. 2
figure 2

Convergence of solutions. This figure visualises the convergence of solutions to the local PDE model (1) as \(T\rightarrow 0^+\), to complement the consistency result presented in Proposition 1. Solutions of the integrodifference model (5) are shown for \(T=0.3\), \(T=0.2\) and \(T=0.1\) and are compared with the solution of the local Klausmeier PDE model (1). Note that unlike in Fig. 1, the spatial domain is chosen to be small to impose the same wavelength restrictions on both models to aid the visualisation of the convergence

Full size image

We base our comparison on the kernels’ standard deviations, which are given by \(\sigma _\phi = \sqrt{2}/a\) for the Laplacian kernel (3), \(\sigma _{\phi _g} = 1/(\sqrt{2}\,a_g)\) for the Gaussian kernel (22) and \(\sigma _{\phi _p} = \sqrt{2}/(\sqrt{b^2-5b+6}\,a_p)\) for the power law kernel (23) provided \(b>3\). It is perfectly reasonable to perform simulations with kernels of infinite standard deviation (e.g. \(b<3\) in the power law kernel) but in the interest of comparing results for the kernels based on their standard deviation we consider only \(b=3.1\) and \(b=4\).

To investigate the model’s behaviour under changes to the dispersal kernels \(\phi \) and \(\phi _1\), we start by considering simultaneous changes in the kernel functions \(\phi \), and \(\phi _1\). The comparison between the kernel functions is based on the standard deviation of the plant dispersal kernel \(\phi \) and the width of the water dispersal kernel \(\phi _1\) is set to \(a_2 = 0.1a\) to obtain a ratio similar to that of the standard deviations under the scalings (6), which corresponds to the large value of the diffusion parameter d in the PDE and integro-PDE models. Figure 3 visualises the simulation results, which show that for small standard deviations, the rainfall threshold \(A_{\max }\) is close to its lower bound, before an increase in the kernel width causes it to peak before slowly decreasing as the kernel widths are further increased. For very narrow dispersal kernels very little spatial interaction takes place. In particular, as \(\sigma \rightarrow 0\), the kernel functions tend to the delta function \(\delta (x)\) centred at 0 and therefore the integrodifference system (5) becomes

$$\begin{aligned} u_{n+1}(x)&= u_n(x)+C\left( u_n(x)^2w_n(x)-Bu_n(x)\right) , \\ w_{n+1}(x)&= w_n(x)+D\left( A-u_n(x)^2w_n(x)-w_n(x)\right) . \end{aligned}$$

For this system, the steady state \((\overline{u_2},\overline{w_2})\), which was randomly perturbed to set the initial condition of the system in the simulation, is always stable. Therefore, no patterns exist and \(A_{\max } = 2B\) is the minimum value of the rainfall parameter for which vegetation is growing uniformly, recalling that for \(A<2B\), the steady state \((\overline{u_2},\overline{w_2})\) does not exist. Further, away from \(\sigma = 0\), a change in kernel width only has very little effect on \(A_{\max }\), an indication that an increase to the width of the plant dispersal kernel has the opposite effect on the tendency to form patterns as an increase to the width of the water dispersal kernel.

Fig. 3
figure 3

The maximum rainfall parameter \(A_{\max }\) under simultaneous changes of the dispersal kernels. This figure visualises variations of \(A_{\max }\) against simultaneous variations of both kernel functions. The standard deviation on the abscissa refers to the plant dispersal kernel \(\phi \), the width of the water dispersal kernel \(\phi _1\) is set to \(a_2=0.1a\). The rainfall threshold is determined up to an interval of length \(10^{-4}\) for \(\sigma _\phi = \{0.01, 0.02, \dots , 0.05, 0.1, 0.2, \dots , 2\}\). The parameter values used for this simulation are \(B=0.45\), \(\ell =\ell _1=0.5\), \(m=5\)

Full size image

To test this hypothesis, we investigate changes in the system’s behaviour as individual kernel functions are changed. First, we consider how the critical rainfall parameter \(A_{\max }\) is affected by a change of the shape of the dispersal kernel \(\phi \) in the plant equation (21a). The result (see Fig. 4a) is consistent with results of the integro-PDE model (Eigentler and Sherratt 2018) on sloped ground. Firstly, an increase in the width of the plant dispersal kernels reduces the size of the parameter region supporting pattern onset, where changes for larger values of the standard deviation \(\sigma _\phi \) are much smaller than close to \(\sigma _\phi =0\). Identical to the nonlocal Klausmeier model, a trend that for small standard deviations those kernel functions that decay algebraically at infinity predict a lower value of \(A_{\max }\) than those decaying exponentially, and vice versa for larger kernel widths, is also observed in these simulations.

Fig. 4
figure 4

The maximum rainfall parameter \(A_{\max }\) under separate variations of the dispersal kernels. a \(A_{\max }\) up to an interval of length \(10^{-4}\) with varying width (\(\sigma _\phi = \{0.05, 0.1, 0.2, \dots , 2 \}\)) and shape of the plant dispersal kernel \(\phi \), while b visualises the effects of changes in the water dispersal kernel \(\phi _1\). The latter was simulated for a larger range of the kernel’s standard deviation \(\sigma _{\phi _1}\), specifically \(\sigma _{\phi _1} = \{1, 2, \dots , 20 \}\), to account for the choice of \(a_2=0.1a\) in the previous simulation. Also in (b) \(A_{\max }\) is determined up to an interval of length \(10^{-4}\). The widths of the fixed kernels are set to \(a_2 = 0.1\) (a) and \(a=1\) (b), respectively. The other parameter values used both simulations are \(B=0.45\), \(\ell =\ell _1=0.5\), \(m=5\)

Full size image

Next, we perform a similar analysis for the symmetric water dispersal kernel \(\phi _1\). To be consistent with the setting \(a_2 = 0.1a\) in the simulation for the simultaneous change of the kernel functions, we consider a larger range of \(\sigma _{\phi _1}\) for this simulation. The results (Fig. 4b) show that for narrow kernels, \(A_{\max }\) is close to its minimum \(A=2B\), i.e. the rainfall interval supporting pattern formation is very small. In particular, as \(\sigma _{\phi _1}\rightarrow 0\), \(A_{\max } \rightarrow 2B\) and no patterns can occur. For the Laplace kernel, this can also be shown using linear stability analysis. If \(\sigma _{\phi _1} = 0\), then \(\widehat{\phi _1} \equiv 1\) and thus the Jacobian (16) becomes

$$\begin{aligned} J = \left( \begin{array}{cc} C{\widehat{\phi }}(k)\alpha &{} C{\widehat{\phi }}(k)\beta \\ D\gamma &{} D\delta \\ \end{array} \right) . \end{aligned}$$

Further, the stability condition is

$$\begin{aligned} k^2> \frac{BCa^2\left( A^2-A\sqrt{A^2-4B^2} -4B^2 \right) }{A^2-A\sqrt{A^2-4B^2}}. \end{aligned}$$

The right hand side is negative and thus the steady state is always stable to spatially heterogeneous perturbations. An increase of the kernel width then causes an increase in the rainfall threshold \(A_{\max }\), where those kernels that decay exponentially at infinity, yield a larger increase than those decaying algebraically.

The results above confirm that the plant dispersal kernel \(\phi \) and the water dispersal kernel \(\phi _1\) have opposite effects on the rainfall threshold \(A_{\max }\). While an increase in the width of the plant dispersal inhibits the onset of patterns, an increase in the standard deviation of the water dispersal kernel increases the tendency to form patterns. This explains the nearly constant value of \(A_{\max }\) in the simulations in which both kernel functions are varied simultaneously. Consequently, these results suggest that it is the ratio of plant dispersal to water dispersal, i.e. the ratio \(\sigma _{\phi }/\sigma _{\phi _1}\) that controls the tendency to form patterns. An increase in the ratio inhibits the onset of patterns, while a decrease has the opposite effect.

5 Discussion

The deliberately basic description of the plant–water dynamics in semi-arid environments by the Klausmeier model provides a rich framework for model extensions to address a range of different features of dryland ecosystems and their effects on vegetation patterns. Extensions include cross advection due to decreased surface water run-off resulting from an increase in infiltration in biomass patches (Wang and Zhang 2019); terrain curvature (Gandhi et al. 2018); nonlocal dispersal of seeds (Eigentler and Sherratt 2018; Bennett and Sherratt 2018); secondary seed dispersal due to overland water flow (Consolo and Valenti 2019); nonlocal grazing effects (Siero et al. 2019; Siero 2018); explicit modelling of a population of grazers (Fernandez-Oto et al. 2019); local competition between plants (Wang and Zhang 2018); the inclusion of autotoxicity (Marasco et al. 2014); multispecies plant communities (Eigentler and Sherratt 2019; Ursino and Callegaro 2016; Callegaro and Ursino 2018) and seasonality and intermittency in precipitation (Ursino and Contarini 2006; Eigentler and Sherratt 2020). One aspect that has not yet been considered in this context is the seasonal separation of plant growth and seed dispersal. In this paper we have considered the synchronised and seasonal occurrence of nonlocal seed dispersal through a system of integrodifference equations based on the Klausmeier reaction-advection-diffusion system.

While an integrodifference system cannot explicitly quantify the temporal separation of seed dispersal occurrences, the model’s derivation and an associated convergence result (Proposition 1) yield a parameter setting in which the length of the growth phase between dispersal stages can be accounted for. However, the main result of the linear stability analysis of the integrodifference model in this paper (Proposition 6) shows that conditions for pattern onset in the integrodifference model (5) are independent of the temporal separation of seed dispersal from plant growth. Moreover, due to the model’s derivation form the Klausmeier model (1), the pattern onset conditions for both models are equivalent.

Some semi-arid environments in which vegetation patterning is a common phenomenon are characterised by large temporal and in particular seasonal fluctuations in their environmental conditions (Noy-Meir 1973; Chesson et al. 2004). For example, observed patterns in Spain, Israel and North America are all located in Mediterranean climate zones (Peel et al. 2007), in which precipitation mainly occurs during winter, while during the summer months little or no rainfall occurs. By contrast, most mathematical models describing these ecosystems employ partial differential equations. While PDE models provide a rich framework for mathematical model analysis, their use is based on the simplifying assumption that all processes occur continuously in time. The results presented in this paper emphasise the importance and significance of results obtained from such models. In the context of seed dispersal, the biologically more realistic temporal separation of plant growth and seed dispersal has no effect on the conditions for pattern onset to occur. We thus conclude that the results obtained for the Klausmeier PDE model are robust to changes in the temporal properties of seed dispersal processes and that the assumption of continuous seed dispersal provides a sufficiently accurate description.

The parameter setting used to establish a connection between the Klausmeier model (1) and the integrodifference model (5) couples the scale parameter a of the seed dispersal kernel to other model parameters. If, however, a more general parameter setting is considered, then the effects of changes to the average seed dispersal distance and the shape of the seed dispersal kernel can be analysed numerically. Our results, which are in full agreement with an earlier investigation of the nonlocal Klausmeier model (2) (Eigentler and Sherratt 2018), show that seed dispersal over longer distances inhibits the formation of patterns (Fig. 4a). Indeed, the threshold \(A_{\max }\) on the rainfall parameter above which no pattern onset occurs, tends to \(A_{\min }\), the minimum rainfall level required for the existence of a nontrival spatially uniform equilibrium, as dispersal distances become sufficiently large. Nevertheless, many plant species in semi-arid ecosystems have developed antitelechoric mechanisms which inhibit long range seed dispersal (Ellner and Shmida 1981; van Rheede van Oudtshoorn and van Rooyen 2013). While in the context of this paper this may appear as an evolutionary disadvantage, the development of narrow seed dispersal kernels is a side effect of other adaptations such as the development of seed containers as a protection to predation (Ellner and Shmida 1981). This suggests the existence of an evolutionary trade-off between seed dispersal distance and plant mortality. A numerical study of the threshold \(A_{\max }\) in the \(\sigma _\phi \)-B parameter plane (Fig. 5) gives some useful insight into this. The trade-off would restrict parameters to some increasing curve in the \(\sigma _\phi \)-B parameter plane. Depending on the exact functional form of such a trade-off, a decrease in the seed dispersal distance \(\sigma _\phi \) may cause a reduction in the precipitation threshold \(A_{\max }\), if the trade-off implies a sufficiently large simultaneous decrease in the plant mortality rate B. A lower \(A_{\max }\) value corresponds to an inhibition of pattern onset. We thus conclude that our model can capture the evolutionary advantage associated with the development of protective antitelechoric mechanisms if the trade-off between seed dispersal distance \(\sigma _\phi \) and plant mortality B is chosen appropriately, but emphasise that we are not aware of any data that provides quantitative information on the exact form of this trade-off.

Fig. 5
figure 5

The threshold \(A_{\max }\) in the \(\sigma _\phi \)-B parameter plane. The numerically obtained rainfall threshold \(A_{\max }\) is shown in the \(\sigma _\phi \)-B parameter plane as a contour plot, where \(\sigma _\phi \) denotes the standard deviation of the plant dispersal kernel \(\phi \). It was obtained on the spatial grid \(\{0.05,0.1,\dots ,1.95,2\}\times \{0.05,0.1,\dots ,1.95,2\}\) for the Laplace kernel (3) and \(a_2 = 0.1\), \(\ell = \ell _1 = 0.5\), \(m=5\). We speculate that there may be an evolutionary trade-off between dispersal distance and resistance to predation, which would restrict parameters to an increasing curve in the \(\sigma _\phi \)-B plane

Full size image

Our results further indicate that the shape of the seed dispersal kernel, and in particular its decay at infinity, has a significant effect on the onset of patterns. Fat-tailed kernels, for example, that account for a higher proportion of long-range dispersal events, yield a lower level of \(A_{\max }\) than kernel functions with exponential decay at infinity for a sufficiently small fixed standard deviation. This highlights the importance of obtaining knowledge of seed dispersal behaviour of plant species, a property that depends on both species and the environment (e.g. seed dispersal agent) (Bullock et al. 2017).

In our integrodifference model (5), we model the redistribution of water through a convolution similar to the modelling of the seed dispersal process. This nonlocal description can account for overland water flow from bare ground to biomass patches across larger distances during precipitation events. It does, however, rely on the assumption that the soil’s properties enhance overland water flow in regions of low biomass. Some ecosystems in semi-arid environments are characterised by soil conditions and soil types (e.g. sand) for which this assumption is invalid (Valentin et al. 1999). The formation of vegetation patterns under such environmental conditions can, however, be explained by other mechanisms, such as laterally extended root networks (Meron 2018). The integrodifference model presented in this paper is based on the assumption that little or no water infiltration occurs in regions of low biomass, and that the overland water flow towards regions of high biomass induced by this soil property is the main mechanism causing the self-organisation into patterns. In this context, our results show that water redistribution over longer distances yields the onset of patterns at higher precipitation levels (Fig. 4b). This is due to the enhancement of the pattern-inducing vegetation-infiltration feedback. Existing biomass patches deplete the water density locally, while regions of bare soil retain a higher water levels. Hence, any redistribution of water has a homogenising effect on the water density which yields to a redistribution of the limiting resource from areas of low biomass to areas of high biomass. An increase in the spatial range of the water redistribution kernel thus strengthens the pattern-inducing feedback and causes pattern onset under larger precipitation volumes.

The work in this paper shows that the description of seed dispersal as a synchronised event during a phase in which no plant growth occurs does not affect the condition for pattern onset compared to the continuous description of seed dispersal in the Klausmeier model (1). The stability of spatial patterns is equally important. A natural area of future work would therefore be an analysis of pattern stability in the integrodifference model (5) comparing results with stability results for the local Klausmeier model (Sherratt and Lord 2007) and the nonlocal Klausmeier model (Bennett and Sherratt 2018). For PDE models, the stability of spatial patterns can be determined through a calculation of their spectra. For this, a method based on numerical continuation has been developed by Rademacher et al. (2007) [for details see Rademacher et al. (2007); Sherratt (2012)]. For integrodifference equations, however, we are not aware of any methods that allow the determination of the stability of a patterned solution.

The integrodifference model (5) not only splits the dynamics of the plant population into separate growth and dispersal stages, but also that of the water dynamics into a water consumption stage and a water redistribution stage. In the model, spatial redistribution of water is synchronised with seed dispersal. This can provide an adequate description for species such as Mesembryanthemum crystallinum and Mesembryanthemum nodiflorum, which synchronise their seed dispersal with the beginning of the rain season (Navarro et al. 2009), but cannot provide a description of seed dispersal during drought periods or of water flow at any other time during the rain season. While a description of the water flow dynamics during precipitation events in the context of a vegetation model has been proposed by Siteur et al. (2014b), the exact dynamics on flat ground are the subject of ongoing research (e.g. Rossi and Ares 2017; Thompson et al. 2011; Wang et al. 2015) and could be utilised in a future extension of the integrodifference model (5). The description of the water density as one single variable would, however, be prohibitive for such an approach. Instead a distinction between surface water and soil moisture, such as in the Rietkerk et al. model (HilleRisLambers et al. 2001; Rietkerk et al. 2002) or the Gilad et al. model (Gilad et al. 2004), needs to be made to distinguish between surface water flow processes and water uptake processes that take place in the soil.

The integrodifference model (5) and its analysis presented in this paper is restricted to a one-dimensional space domain, motivated by the original formulation of the Klaumeier model and its mathematical accessibility (Klausmeier 1999). However, the consideration of a second space dimension is expected to give more insights into the ecohydroglogical dynamics, in particular on pattern existence and stability. For example, in related PDE models on two-dimensional space domains, different types of patterned solutions exist (gap patterns, labyrinth patterns, striped patterns and gap patterns) and phase transitions along the precipitation gradient can be investigated (Meron 2012). Moreover, even on sloped terrain, the impact of the consideration of a two-dimensional domain is significant, as the analysis on a one-dimensional domain may overestimate the size of the patterns’ stability regions (Siero et al. 2015). The analysis of the integrodifference model (5) on a two-dimensional domain presents a considerable challenge, in particular if one would want to obtain a wavenumber-independent results analogous to Proposition 6. Nevertheless, this would be a natural area of potential future work to further disentangle the complex ecosystem dynamics.

Finally, we remark that the integrodifference model (5) describes the discrete structure of plant growth mechanisms caused by the seasonality in precipitation. However, it does not capture the dynamics specific to drought periods between rainfall events and is thus only able to provide an insight into effects of accumulated rainfall volume rather than the temporal separation of precipitation seasons. In separate work, we account for a combination of rainfall, plant growth and seed dispersal pulses with the continuous nature of plant loss and water evaporation and drainage, using an impulsive model (Eigentler and Sherratt 2020). Such models combine partial differential equations with integrodifference equations (see for example Wang and Lutscher (2018) for an impulsive model in the context of predator-prey dynamics with synchronised predator reproduction). The impulsive model has its own limitations as it can only take into account a periodic separation of precipitation events, but not any seasonal patterns. A potential area of future work therefore consists of a combination of these approaches to describe both the seasonal and intermittent nature of rainfall in semi-arid climate zones.