Model of pattern formation in marsh ecosystems with nonlocal interactions

TL;DR

This paper develops a mathematical model using reaction-diffusion equations with nonlocal interactions to explain the formation of spatial patterns in marsh shorelines driven by grass-sediment feedbacks.

Contribution

It introduces a novel reaction-diffusion framework with a Mexican-hat kernel to capture scale-dependent feedbacks in marsh ecosystems.

Findings

Pattern formation depends on the scale and strength of feedback.

Conditions for spatial patterns are derived mathematically.

Model explains characteristic marsh shoreline variability.

Abstract

Smooth cordgrass Spartina alterniflora is a grass species commonly found in tidal marshes. It is an ecosystem engineer, capable of modifying the structure of its surrounding environment through various feedbacks. The scale-dependent feedback between marsh grass and sediment volume is particularly of interest. Locally, the marsh vegetation attenuates hydrodynamic energy, enhancing sediment accretion and promoting further vegetation growth. In turn, the diverted water flow promotes the formation of erosion troughs over longer distances. This scale-dependent feedback may explain the characteristic spatially varying marsh shoreline, commonly observed in nature. We propose a mathematical framework to model grass-sediment dynamics as a system of reaction-diffusion equations with an additional nonlocal term quantifying the short-range positive and long-range negative grass-sediment interactions. We use a Mexican-hat kernel function to model this scale-dependent feedback. We perform a steady state biharmonic approximation of our system and derive conditions for the emergence of spatial patterns, corresponding to a spatially varying marsh shoreline. We find that the emergence of such patterns depends on the spatial scale and strength of the scale-dependent feedback, specified by the width and amplitude of the Mexican-hat kernel function.

Tables

Table 1. Table 1: Biologically realistic parameters for the original system before non-dimensionalization

Symbol	Meaning	Unit	Value	Source
$\widehat{D_G}$	cordgrass diffusion coefficient	$m^2$ yr-1	0.06 - 0.135	(Adams et al., 2012)
$\widehat{D_S}$	sediment diffusion coefficient	$m^2$ yr-1	0.876	(Liu et al., 2014)
$c$	self-limiting growth rate of cordgrass	$m^2$ shoots-1 yr-1	0.0057	(Yang et al., 2014)
$\psi$	minimum erosion rate	yr-1	0.002-0.3	(Hardaway Jr & Byrne, 1999; Rosen, 1980)
$k_s$	cordgrass density at which marsh erosion is half-maximal	shoots $m^{-2}$	30-50	estimated
$g$	erosion constant in the absence of cordgrass	non-dimensional	5	(Mariotti & Fagherazzi, 2010; Sheehan & Ellison, 2015)
$\eta$	sediment deposition rate	$m$ yr-1	0.002-0.006	(Stumpf, 1983; Goodman et al., 2007)
$p^{\ast }$	intrinsic growth rate of cordgrass	yr-1	1.5	(Yang et al., 2014)
$l_1$	sediment threshold for cordgrass persistence	$m$	0.02	estimated
$l_1^{\ast }$	sediment elevation at which cordgrass growth is half-maximal	$m$	0.06	estimated
$\hat{\lambda }$	strength of nonlocal cordgrass-sediment interactions	$m^2$ shoots-1yr-1	0.0004-0.3	(Bouma et al., 2007)
${\sigma }_1$	standard deviation of the excitatory feedback for cordgrass	$m$	0.43	(Bouma et al., 2007)
${\sigma }_2$	standard deviation of the inhibitory feedback for cordgrass	$m$	0.68	(Bouma et al., 2007)

Equations

$$\left\{\begin{aligned} \partial_{\tau}\hat{G}&=\hat{D}_{\hat{G}}\partial^{2}_{x}\hat{G}+\underbrace{\hat{G}\Big{(}\hat{F}(\hat{S})-c\hat{G}\Big{)}}_{\text{Logistic growth}},&x\in\mathbb{R},\tau>0,\\ \partial_{\tau}\hat{S}&=\hat{D}_{\hat{S}}\partial^{2}_{x}\hat{S}+\underbrace{\eta}_{\text{Deposition}}-\underbrace{\hat{S}\hat{L}(\hat{G})}_{\text{Erosion}}+\underbrace{\hat{\lambda}\hat{S}\int_{-\infty}^{\infty}P(x^{\prime})\hat{G}(x-x^{\prime})dx^{\prime}}_{\text{Nonlocal deposition/erosion}},&x\in\mathbb{R},\tau>0,\\ \hat{G}(x,0)&=\hat{G_{0}}(x,0)\geq 0,\,\,\,\hat{S}(x,0)=\hat{S_{0}}(x,0)\geq 0,&x\in\mathbb{R}.\end{aligned}\right.$$

$$\hat{F}(\hat{S})=\frac{p^{\ast }(\hat{S}-l_1)}{\hat{S}+l_1^{\ast }},\hat{L}(\hat{G})=\frac{\psi (\hat{G}+k_{s}g)}{\hat{G}+k_s},$$

$$p^{\ast },c,l_1,l_1^{\ast },\psi ,k_s,g,\eta ,\hat{\lambda }\ge 0.$$

$$\displaystyle P(x)=\frac{1}{\sqrt{2\pi}}\Big{[}\frac{1}{\sigma_{1}}\exp\Big{(}-\frac{x^{2}}{2\sigma_{1}^{2}}\Big{)}-\frac{1}{\sigma_{2}}\exp\Big{(}-\frac{x^{2}}{2\sigma_{2}^{2}}\Big{)}\Big{]},\,\,\sigma_{1}<\sigma_{2}.$$

$${\int }_{-\infty }^{\infty }P(x)dx=0.$$

$$t$$

$$\left\{\begin{aligned} \partial_{t}G&=D_{G}\partial^{2}_{x}G+G\Big{(}F(S)-G\Big{)},&x\in\mathbb{R},t>0,\\ \partial_{t}S&=D_{S}\partial^{2}_{x}S+\phi\Big{(}-L(G)S+1\Big{)}+\lambda S\int_{-\infty}^{\infty}P(x^{\prime})G(x-x^{\prime})dx^{\prime},&x\in\mathbb{R},t>0,\\ G(x,0)&=G_{0}(x,0)\geq 0,\,\,\,S(x,0)=S_{0}(x,0)\geq 0,&x\in\mathbb{R}.\end{aligned}\right.$$

$$F(S)$$

$$\begin{array}{ccccc}{e}_1=\frac{\psi g{l}_1}{\eta },& p_1=\frac{\psi g{l}_1^{\ast }}{\eta },& e_3=\frac{k_{s}c}{p^{\ast }},& \delta =\frac{1}{g}& \\ \varphi =\frac{\psi g}{p^{\ast }},& D_G=\frac{{\hat{D}}_G}{p^{\ast }},& D_S=\frac{{\hat{D}}_S}{p^{\ast }},& \lambda =\frac{\hat{\lambda }}{c}.& \end{array}$$

$$\left\{\begin{aligned} \frac{dG}{dt}&=G\Big{(}F(S)-G\Big{)},&t>0,\\ \frac{dS}{dt}&=\phi(-L(G)S+1),&t>0.\end{aligned}\right.$$

$${\mathbf{J}}_{E_S}=\left[\begin{array}{cc}F(1)& 0\\ -\varphi \frac{dL}{dG}& -\varphi L(0)\end{array}\right].$$

$$A=1+p_1\delta ,B=-1+e_3+p_1{e}_3.$$

$$e_1^{\ast }=\frac{2e_{3}A({\delta }^2+B\delta -A{e}_3)+(2A{e}_3-B\delta )\sqrt{A^2{e}_3^2-A{e}_3\delta (B+\delta )}}{{\delta }^2\sqrt{A^2{e}_3^2-A{e}_3\delta (B+\delta )}}.$$

$$G$$

$$⟹e_1=\frac{1-G}{L(G)}-p_{1}G$$

$$⟹e_1=\frac{e_3-BG-A{G}^2}{G\delta +e_3}:=K(G),$$

$$G_K^{\pm }=\frac{-B\pm \sqrt{B^2+4e_{3}A}}{2A}.$$

$$K^{\prime }(G)$$

$$=\frac{-L(G)-L^{\prime }(G)(1-G)}{L(G{)}^2}-p_1,$$

$$K(0)=1,K^{\prime }(0)=\frac{-(B+\delta )}{e_3},$$

$${\tilde{G}}_{\pm }=\frac{-e_3\pm \sqrt{e_3^2-\frac{e_3\delta (B+\delta )}{A}}}{\delta }.$$

$$e_1^{\ast }=K({\tilde{G}}_{+})=\frac{2e_{3}A({\delta }^2+B\delta -A{e}_3)+(2A{e}_3-B\delta )\sqrt{A^2{e}_3^2-A{e}_3\delta (B+\delta )}}{{\delta }^2\sqrt{A^2{e}_3^2-A{e}_3\delta (B+\delta )}}.$$

$$G_{\pm }^{\ast }$$

$$S_{\pm }^{\ast }$$

$$\mathbf{J}(E_{GS})=\left[\begin{array}{cc}-G^{\ast }& G^{\ast }{F}^{\prime }(\frac{1}{L(G^{\ast })})\\ -\varphi \frac{1}{L(G^{\ast })}{L}^{\prime }(G^{\ast })& -\varphi L(G^{\ast })\end{array}\right].$$

$$\mathbf{J}=\left[\begin{array}{cc}-& +\\ +& -\end{array}\right].$$

$$\mathrm{Tr}\mathbf{J}$$

$$\mathrm{Det}\mathbf{J}$$

$$\displaystyle=\phi G^{*}\Big{(}L(G^{*})+\frac{1}{L(G^{*})}L^{\prime}(G^{*})\frac{p_{1}+e_{1}}{(\frac{1}{L(G^{*})}+p_{1})^{2}}\Big{)}$$

$$\displaystyle=\phi G^{*}\Big{(}L(G^{*})+\frac{1}{L(G^{*})}L^{\prime}(G^{*})\frac{p_{1}+\frac{1-G^{*}}{L(G^{*})}-p_{1}G^{*}}{(\frac{1}{L(G^{*})}+p_{1})^{2}}\Big{)}$$

$$\displaystyle=\phi G^{*}\Big{(}L(G^{*})+\frac{L^{\prime}(G^{*})(1-G^{*})}{1+L(G^{*})p_{1}}\Big{)}$$

$$\displaystyle=\phi G^{*}\Big{(}\frac{-K^{\prime}(G^{*})L^{2}(G^{*})}{1+L(G^{*})p_{1}}\Big{)}.$$

$$-e_1\delta -B$$

$$(e_1\delta +B{)}^2$$

Full text

Keywords: Pattern formation; nonlocal interactions; marsh ecosystem; reaction diffusion; steady state; cooperation.

MSC (2010): 92D40, 92D25, 35K57; 35B36

1 Introduction

Tidal marshes are among the richest and most productive ecosystems, supporting a variety of wildlife, serving as storm and erosion buffers, and playing an important role in improving water quality (Perry & Atkinson, 2009; Fagherazzi et al., 2013; Fagherazzi, 2014). The global loss of these ecosystems in the recent decades has motivated much research to understand their dynamics and aid in their restoration and management (Deegan et al., 2012; Priestas et al., 2015). Marsh evolution is dynamic and complex, combining various biological and morphological processes happening not only in the marsh itself, but also in the tidal flat that borders it. As a result of these forces and interactions, a sharp scarp separating the marsh and tidal flat becomes a characteristic feature. The processes that take place on this scarp (i.e., marsh edge) influence whether the marsh recedes or expands (Tonelli et al., 2010). Various configurations of the marsh edge can be observed in nature, ranging from a uniform to a more jaggedy, sinusoidal shoreline (Figure 1). While previous ecogeomorphic models have carefully considered the effects of sea-level rise, marsh vegetation colonization, wave activity, sediment fluxes, and underlying hydrodynamics (Mariotti & Fagherazzi, 2010; Tonelli et al., 2010; Fagherazzi et al., 2012; Schile, 2014), most of these have been numerical, computationally intensive models. We propose a simpler, phenomenological model to describe the large-scale evolution of the marsh in the horizontal direction in terms of two-way interactions between marsh vegetation and sedimentation. In particular, we are interested in the scale-dependent feedback present between marsh vegetation and sedimentation and whether this scale-dependent feedback can explain the spatially varying marsh shoreline, observed in nature.

Scale-dependent feedbacks are characterized by the presence of positive and negative interactions that happen at different spatial scales. In particular, scale-dependent feedbacks involving long-range negative interactions and short-range positive interactions are thought to be crucial for pattern development (Gierer & Meinhardt, 1972; Green & Sharpe, 2015; Hiscock & Megason, 2015), explaining spatially varying patterns in chemical (Castets et al., 1990; Rovinsky & Menzinger, 1993), biological (Nakamasu et al., 2009; Raspopovic et al., 2014) and ecological systems (Rietkerk & Van de Koppel, 2008). In ecological systems, such scale-dependent feedbacks are thought to explain patterns in arid savannas, mussel beds, coral and oyster reefs, mudflats and other ecosystems (Rietkerk & Van de Koppel, 2008; van der Heide et al., 2012; Dibner et al., 2015; de Jager et al., 2017; Pringle & Tarnita, 2017; Barbier et al., 2008). Previously, we proposed a mathematical framework to investigate the evolution of the marsh edge as a result of scale-dependent interactions between sedimentation dynamics and two common marsh species, ribbed mussels (Geukensia demissa) and smooth cordgrass (Spartina alterniflora) (Zaytseva et al., 2018), whose facilitatory nature and positive feedbacks have a significant effect on marsh development and proliferation (Bertness, 1984; Bertness & Grosholz, 1985; Watt et al., 2010; Altieri et al., 2007). While mussels are commonly found in tidal marshes, that is not always the case. Since we are interested in the marsh edge dynamics in the absence of mussels, in this paper we focus on the related model without the mussel population. The goal is to understand which conditions lead to a spatially varying marsh shoreline (Figure 1a) versus a spatially uniform marsh shoreline (Figure 1b) and what may be the implications of this spatial heterogeneity.

In general, there are three classes of deterministic models that explain pattern formation as a result of scale-dependent feedbacks: Turing-style activator inhibitor models, kernel-based models and differential flow models (Borgogno et al., 2009). In Turing-style activator inhibitor models, the pattern formation arises as a result of differences in the diffusion coefficients of the activator and inhibitor species (Turing, 1952; White, 1998; Parshad et al., 2014). Kernel-based models are typically integro-differential equations where the pattern formation arises from the spatial interactions modeled using a kernel function, describing the nature of the short-range and long-range interactions (Britton, 1990; Gourley et al., 2001; Murray, 2001; Billingham, 2003; Ninomiya et al., 2017). This kernel-based approach is a common feature in neural models (Amari, 1977), and has also been used in models of vegetation patterns in arid and semi-arid climates (D’Odorico et al., 2006; Borgogno et al., 2009; Merchant & Nagata, 2011; Martínez-García et al., 2013, 2014; Martínez-García & Lopez, 2018). Finally, differential flow models are similar to Turing models but the pattern formation now arises not just from the differences in diffusion coefficients, but also from the differences in the flow rates of the species, reflected in the additional advection terms (Rovinsky & Menzinger, 1993; Siero et al., 2015; Klausmeier, 1999).

The model we propose here combines elements of both the Turing-model and the kernel-based model and includes both diffusion terms and a kernel function that describes the short-range and long-range interactions between marsh grass and sediment volume. On a local scale, marsh grass attenuates hydrodynamic energy, enhancing sediment accretion and promoting further vegetation growth while the diverted water flow promotes formation of erosion troughs over longer distances (Bouma et al., 2007; Balke et al., 2012; Schwarz et al., 2015). We model this scale-dependent feedback using a Mexican-hat kernel function that quantifies the strength of positive and negative feedbacks neighboring individuals exert on each other (Fuentes et al., 2003; D’Odorico et al., 2006; Borgogno et al., 2009; Siebert & Schöll, 2015). Similar kernel-based approaches have been used to model nonlocal interactions in the context of predator-prey and competition dynamics (Merchant & Nagata, 2011; Bayliss & Volpert, 2015; Banerjee & Volpert, 2016). The interactions in our system are mostly cooperative and the impact of nonlocal interactions in such systems have not been studied in depth. Given the importance of facilitation in ecosystem dynamics (Bertness & Callaway, 1994; Halpern et al., 2007; Silliman et al., 2015; He et al., 2013), it therefore becomes imperative to study nonlocal interactions in cooperative systems. In addition, cooperative systems are likely to display bistable dynamics and the phenomenon of hysteresis (Kéfi et al., 2016; van de Koppel et al., 2001). This makes such systems especially prone to collapsing to an irreversible state as environmental conditions gradually worsen and a tipping point is reached (Dakos et al., 2011; Kéfi et al., 2014, 2016). Pattern formation has previously been suggested as a possible coping mechanism, allowing such systems to escape degradation past their tipping point (Chen et al., 2015). Due to the reported degradation of tidal marsh habitats around the world, the study of pattern formation in these systems becomes particularly important and can provide more insight into the possible pattern forming mechanism and its implication for the system’s resilience and adaptation to environmental changes.

Our paper is organized as follows: In Section 2, we introduce the nonlocal reaction-diffusion model and the background from ecological literature. Section 3 includes analysis and simulation results. By approximating our model using a steady state biharmonic approximation, we are able to derive conditions for the emergence of spatial patterns in our system. We then use numerical simulations to confirm our theoretical findings. Some concluding remarks are made in Section 4.

2 Model

We consider the two-way interactions between marsh grass and sediment (Figure 2). Marsh grass binds sediment, stabilizes the marsh edge and attenuates wave energy, helping to mitigate effects of erosion (Gleason et al., 1979; Gedan et al., 2011; Ysebaert et al., 2011; Möller et al., 2014). As a consequence of reduced erosion, the increased sediment levels promote vegetation growth by decreasing tidal currents (Nyman et al., 1993; van de Koppel, van der Wal, Bakker & Herman, 2005). Along with these local interactions, there is a nonlocal interaction that occurs between marsh vegetation and sediment (van Wesenbeeck et al., 2008; Schwarz et al., 2015; Bouma et al., 2009; Van Hulzen et al., 2007). Over short distances, marsh vegetation enhances sediment accretion through the attenuation of hydrodynamic energy, contributing to short-range activation. However, as the water gets diverted to the surrounding areas, those areas erode more quickly, contributing to long-range inhibition (Bouma et al., 2007; Balke et al., 2012; Bouma et al., 2013; Fagherazzi et al., 2013; Fagherazzi, 2014).

Incorporating all the above mentioned interactions, we obtain the following system:

[TABLE]

where

[TABLE]

with

[TABLE]

We consider the change in grass shoot density $\hat{G}(x,t)$ ($shoots/m^{2}$) and sediment height $\hat{S}(x,t)$ ($meters$) on an infinite domain with $x\in\mathbb{R}$, which represents the one-dimensional horizontal cross-section of the marsh edge (see Figure 3a). We assume logistic growth for the grass density and make an adjustment for the obligatory nature of grass-sediment interactions where below some minimum sediment height $l_{1}$, grass cannot persist. For the sediment equation, we include the baseline sediment deposition $\eta$ (van de Koppel, Rietkerk, Dankers & Herman, 2005; Liu et al., 2012, 2014). The erosion term is a decreasing function of grass density with $g>1$ where $\psi g$ corresponds to the minimum erosion rate in the total absence of grass (Mariotti & Fagherazzi, 2010; Silliman et al., 2012). In addition, each equation also includes a diffusion term to quantify spread along the shoreline with diffusion coefficients $\hat{D}_{\hat{G}}$ and $\hat{D}_{\hat{S}}$. To model the scale-dependent interactions, we use a convolution term with a Mexican-hat kernel function $P(x)$:

[TABLE]

The choice of the kernel function is appropriate given the nature of the scale-dependent feedback with short-range positive interactions and long-range negative interactions (Figure 3b). There are three main parameters that control the shape of the kernel: $\hat{\lambda}$, which modulates the amplitude and variances $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$, which specify the scale of the excitatory and inhibitory interactions, respectively. Further, the kernel function $P(x)$ is symmetric and satisfies the following property:

[TABLE]

For mathematical simplification, we non-dimensionalize system (2.1) by using the following rescaling:

[TABLE]

Then the original system (2.1) becomes:

[TABLE]

with

[TABLE]

and $P(x)$ still defined as before. The new parameters are all positive with the following rescaling:

[TABLE]

Not only does this rescaling simplify the notation, but it also allows for an easier interpretation of the functional forms of $F(S)$ and $L(G)$ (See Figure 10 in the Appendix). The scaled intrinsic growth rate of grass is now between [math] and $1$, and we can think of the threshold $e_{1}$ as the minimum amount of sediment necessary for the persistence of grass. Similarly, the erosion term given by $L(G)$ is scaled to be between $\delta$ and $1$ for ease of interpretation.

3 Results

In classic Turing models, spatially patterned solutions result from symmetry-breaking instability in which an otherwise stable spatially uniform steady state can become destabilized by the addition of diffusion and lead to the emergence of spatial patterns. The condition for the emergence of spatial patterns is contingent on the idea that the species in the model diffuse at significantly different rates, with the activator species diffusing much more slowly than the inhibitor species. The conditions for such a Turing-instability can be derived by performing a linear stability analysis around the positive steady state to obtain conditions under which the addition of diffusion acts to destabilize the system. Our model differs from the classic Turing model in that it lacks the classic activator-inhibitor dynamics and includes an additional kernel function term that models the scale-dependent feedback between grass and sediment volume. Assuming that the kernel function has a limited effect at relatively large distances, we can perform a biharmonic approximation of our system and decompose the integral term into two terms involving just partial derivatives, corresponding to short-range positive interactions and long-range negative interactions. We can then perform a linear stability analysis around the positive steady state and derive conditions under which this state is destabilized and leads to the emergence of a spatially periodic solution. Therefore, we first consider the spatially independent dynamics of our model and derive conditions under which the positive steady state is stable in the corresponding system of ODEs and then use these results to understand the spatial dynamics of the full model.

3.1 Spatially homogeneous model

Let’s assume that $G$ and $S$ do not vary and are spatially constant. Then, we can use the property in (2.3) and drop both the diffusion and integral terms. In this way, we are left with the following spatially independent system:

[TABLE]

We look for spatially uniform steady states $(G^{*},S^{*})$ of (3.1) which satisfy $\frac{dG}{dt}=0$ and $\frac{dS}{dt}=0$. There are two such types of steady states: the degraded (grass-free) state $E_{S}=(0,1)$ and the coexistence state $E_{GS}=(G^{*},\frac{1}{L(G^{*})})$ with both grass and sediment present, where $G^{*}$ satisfies $G=F(\frac{1}{L(G)})$. Since we are interested in physically realistic positive steady states, the coexistence state $E_{GS}$ exists if and only if $\frac{1}{L(G^{*})}>e_{1}$.

We first consider the degraded state $E_{S}=(0,1)$ and its stability. This result is summarized below.

Proposition 3.1**.**

The degraded steady state $E_{S}=(0,1)$ is locally asymptotically stable with respect to (3.1) if $e_{1}>1$ and is unstable with respect to (3.1) if $e_{1}<1$.

Proof.

The Jacobian matrix $\mathbf{J}$ of (3.1) evaluated at $E_{S}=(0,1)$ is given by:

[TABLE]

The two corresponding eigenvalues are $\lambda_{1}=F(1)=\frac{1-e_{1}}{1+p_{1}}$ and $\lambda_{2}=-\phi L(0)=-\phi$. It is clear that $\lambda_{2}$ is always negative. Further, $\lambda_{1}=\frac{1-e_{1}}{p_{1}+1}$ is negative for $e_{1}>1$. Therefore, the steady state $E_{S}$ is locally asymptotically stable for $e_{1}>1$ and unstable for $e_{1}<1$. ∎

The parameter $e_{1}$ is the minimal steady state sediment elevation needed for the persistence of grass. For the trivial steady state $E_{S}=(0,1)$, as long as $e_{1}>1$, its value will exceed the steady state value of sediment, leading to a negative growth rate for grass and a stable trivial state.

We now consider the positive coexistence steady state $E_{GS}=(G^{*},\frac{1}{L(G^{*})})$ where $G^{*}$ satisfies $G=F(\frac{1}{L(G)})$ and obtain the following theorem:

Theorem 3.2**.**

Suppose that $p_{1},e_{1},e_{3},\phi>0$ and $0<\delta<1$. Let

[TABLE]

(Case I)*  If $B+\delta<0$, then there exists a saddle-node bifurcation point $e_{1}=e_{1}^{*}>1$ such that (3.1) has one positive steady state $(G^{*}_{+},\frac{1}{L(G^{*}_{+})})$ for $0<e_{1}\leq 1$ and $e_{1}=e_{1}^{*}$, two positive steady states $(G^{*}_{\pm},\frac{1}{L(G_{\pm})})$ for $1<e_{1}<e_{1}^{*}$, and no positive steady state for $e_{1}>e_{1}^{*}$. The bifurcation point $e_{1}^{*}$ is defined as follows:*

[TABLE] 2. 2.

(Case II)*  If $B+\delta\geq 0$, then there exists a unique positive steady state $(G^{*}_{+},\frac{1}{L(G^{*}_{+})})$ for all $0<e_{1}<1$, and no positive steady state for $e_{1}\geq 1$.*

Proof.

We can rewrite $G=F(\frac{1}{L(G)})$ as

[TABLE]

where $A$ and $B$ are defined as in (3.2). The function $K(G)$ from (3.4) crosses the horizontal axis at

[TABLE]

Since $A\geq 0$, the roots in (3.5) have to be of opposite sign. Therefore, the graph of $K(G)$ has one positive and one negative root. Note that the vertical asymptote of $K(G)$ is irrelevant as it is located where $G$ is negative and outside of the physically realistic range.

Differentiating $K(G)$ in (3.4) with respect to $G$ yields:

[TABLE]

and

[TABLE]

Further, we can set $K^{\prime}(G)=0$ to obtain the maximum and minimum points of the function:

[TABLE]

We then have two cases arising depending on the sign of $K^{\prime}(0)$ in (3.7) (Figure 4).

Case I. The first case corresponds to $B+\delta<0$ and a more physically realistic parameter regime where grass is more effective at attenuating erosion. In this parameter regime, $\delta$ is smaller and therefore, the erosion rate decays faster as a function of grass. From (3.8), it is clear that $\tilde{G}_{-}<0<\tilde{G}_{+}$ and there exists only one peak for positive values of G, given by the value of $\tilde{G}_{+}$. Further, since $K(0)=1$, this means that for $e_{1}<1$, there exists only one positive steady state and for $e_{1}>1$, there exist two positive steady states (Figure 4a). The two positive steady states collide and annihilate each other at the saddle-node bifurcation point $e_{1}^{*}$, given by:

[TABLE]

Case II Case II corresponds to $B+\delta\geq 0$, a parameter regime in which cordgrass is less effective at attenuating sediment erosion. From (3.8), it is clear that $\tilde{G}_{-}<\tilde{G}_{+}<0$ and there exist no peaks for positive values of G. Therefore, since $K(0)=1$, for $e_{1}<1$, we have one positive steady state, while for $e_{1}>1$ there is no positive steady state (Figure 4b). ∎

Now that we know how many positive steady states can be expected, we evaluate their stability and obtain the following theorem:

Theorem 3.3**.**

Suppose that $p_{1},e_{1},e_{3},\phi>0$, $0<\delta<1$ and let $A$ and $B$ be defined as in (3.2). For the positive steady states $E_{GS}=(G^{*}_{\pm},S^{*}_{\pm})$ defined as:

[TABLE]

we have the following cases:

(Case I)*  Let $B+\delta<0$. For $1<e_{1}<e_{1}^{*}$ where $e_{1}^{*}$ is defined as in (3.3), the high density positive steady state $(G^{*+},\frac{1}{L(G_{+}^{*})})$ is locally asymptotically stable and the low density positive steady state $(G_{-}^{*},\frac{1}{L(G_{-}^{*})})$ is unstable. For $0<e_{1}\leq 1$, there is only one positive steady state $(G^{*+},\frac{1}{L(G_{+}^{*})})$ which is locally asymptotically stable.* 2. 2.

(Case II)*  Let $B+\delta\geq 0$. Then, for all $0<e_{1}<1$, the unique positive steady state $(G_{+}^{*},\frac{1}{L(G_{+}^{*})})$ is locally asymptotically stable.*

Proof.

We first evaluate the Jacobian matrix $\mathbf{J}$ of system (3.1) at the positive steady state $E_{GS}$. This is given by:

[TABLE]

This is just the general form of the Jacobian evaluated at the positive steady state type. From Theorem 3.2, we can have either two such positive states (high and low) or just one, depending on the parameter regime. We will consider both cases in this proof. We note the special form of the Jacobian matrix, reflecting the cooperative nature of our system:

[TABLE]

From $\mathbf{J}$, we can define the trace and determinant as follows:

[TABLE]

In order for $E_{GS}$ to be locally asymptotically stable, we need $\operatorname{Tr}{\mathbf{J}}<0$ and $\operatorname{Det}{\mathbf{J}}>0$. Since $L(G^{*})\geq 0$ and $G^{*}$ is a positive quantity, the trace of $\mathbf{J}$ is always negative. Note that since $\operatorname{Tr}{\mathbf{J}}<0$, a Hopf bifurcation cannot occur from the positive steady state. Therefore, to assess stability, we need to determine the sign of $\operatorname{Det}{\mathbf{J}}$. Using (3.4) and (LABEL:eq:slope), we can rewrite $\operatorname{Det}{\mathbf{J}}$ in terms of $K^{\prime}(G)$ to obtain:

[TABLE]

From equation (3.4), we can solve the steady states explicitly as in (3.9). We now consider two cases from Theorem 3.2. For Case I, both low and high steady states $(G_{+}^{*},S_{+}^{*})$ and $(G_{-}^{*},S_{-}^{*})$ are positive, while for Case II, only the high positive steady state $(G_{+}^{*},S_{+}^{*})$ is positive. These are the steady states we consider and assess their stability.

Case I For Case I ($B+\delta<0$), we have the following scenarios:

• •

(i) For $1<e_{1}<e_{1}^{*}$, there are two positive steady states $(G_{+}^{*},S_{+}^{*})$ and $(G_{-}^{*},S_{-}^{*})$. From the definitions of the steady states in (3.9), it follows that

[TABLE]

Evaluating $K^{\prime}(G)$ from equation (LABEL:eq:slope) at $G=G_{+}^{*}$ and $G=G_{-}^{*}$ yields:

[TABLE]

with

[TABLE]

Using conditions from (3.10), we can show:

[TABLE]

and

[TABLE]

Therefore, since $K^{\prime}(G)<0$ on the $(G_{+}^{*},S_{+}^{*})$ branch, $\operatorname{Det}{\mathbf{J}}$ evaluated at $(G_{+}^{*},S_{+}^{*})$ is positive and $(G_{+}^{*},S_{+}^{*})$ is locally asymptotically stable. Similarly, since $K^{\prime}(G)>0$ on the $(G_{-}^{*},S_{-}^{*})$ branch, $\operatorname{Det}{\mathbf{J}}$ evaluated at $(G_{-}^{*},S_{-}^{*})$ is negative and $(G_{-}^{*},S_{-}^{*})$ is unstable.

• •

(ii) For $0<e_{1}\leq 1$, there is only one positive steady state branch corresponding to $(G_{+}^{*},\frac{1}{L(G_{+}^{*})})$. Further, we can show that

[TABLE]

From (3.11), it then follows that $K^{\prime}(G_{+}^{*})<0$. Since $\operatorname{Det}{\mathbf{J}}$ evaluated at $(G_{+}^{*},S_{+}^{*})$ is positive, $(G_{+}^{*},S_{+}^{*})$ is locally asymptotically stable.

Case II For Case II ($B+\delta\geq 0$) , there is a unique positive steady state branch corresponding to $(G_{+}^{*},S_{+}^{*})$. From equation (LABEL:eq:slope) it is clear that $K^{\prime}(G)<0$ for all positive values of $G$. Therefore, since $\operatorname{Det}{\mathbf{J}}$ evaluated at $(G_{+}^{*},S_{+}^{*})$ is positive, the steady state $(G_{+}^{*},S_{+}^{*})$ is locally asymptotically stable. ∎

The results from Proposition 3.1, Theorem 3.2 and Theorem 3.3 are summarized in Figure 5. In Case I, the system displays bistability for values $1<e_{1}<e_{1}^{*}$, where both the high positive steady state and the trivial steady state are stable, separated by an unstable positive steady state branch. The two positive steady states then merge in a saddle-node bifurcation at $e_{1}^{*}$, after which only the stable trivial steady state $E_{S}$ remains. Bistability is not surprising given the highly cooperative nature of this system and large role that the grass plays in erosion mitigation. In Case II, which corresponds to the scenario where grass is less effective at attenuating erosion, the unique stable positive state gradually decreases and eventually undergoes a transcritical bifurcation at $e_{1}=1$ at which it exchanges stability with the trivial steady state $E_{S}$. Note that the positive steady state in Case II ceases to exist for smaller values of $e1$ than in Case I. This is intuitive as Case I corresponds to a more cooperative parameter regime that makes population persistence more possible.

3.2 Generalized Cooperative System with Nonlocal Interactions

We now consider the spatially extended system to investigate the emergence of a patterned solution. Given the complexity of the spatially extended system (2.4), we carry out a steady state biharmonic approximation of this system, allowing us to perform linear stability analysis on the approximated system and gain insight into the dynamics of the original system (2.4) (Murray, 2001; D’Odorico et al., 2006; Borgogno et al., 2009).

Let’s consider the following general form of our system (2.4):

[TABLE]

where $P(x)$ is defined the same as in (2.2), and $f,g$ are general smooth functions. Following standard procedure, we assume the kernel has a limited effect at relatively large distances and perform a Taylor’s expansion of the integral term around $x^{\prime}=x$ (Murray, 2001, pages 482-489):

[TABLE]

This is a reasonable assumption in the context of our model as the scale-dependent grass-sediment feedback is thought to occur on a relatively small spatial scale ($1-4$ meters). We can then define the moments $P_{m}$ in the following way:

[TABLE]

Given the symmetry of the kernel $P(x)$, the odd-power moments vanish, as does $P_{0}$ since $\int_{-\infty}^{\infty}P(x^{\prime})dx^{\prime}=0$. From the specific form of the Mexican-hat kernel in (2.2), we can obtain exact expressions for $P_{2}$ and $P_{4}$ in term of the variances $\sigma_{1}$ and $\sigma_{2}$ of the excitatory and inhibitory effects, respectively:

[TABLE]

Truncating the expansion at the fourth partial derivative, the original system (3.14) can now be approximated by the following biharmonic system (Bates & Ren, 1996, 1997; Couteron & Lejeune, 2001):

[TABLE]

In this way, the evolution of $u$ and $v$ now depends not only on their own diffusion as in the classic reaction-diffusion system, but also on the additional short-range cross-diffusion $\partial^{2}_{x}u$ and long-range cross-diffusion $\partial^{4}_{x}u$ terms. Here, $\lambda P_{2}$ and $\lambda P_{4}$ represent the corresponding cross-diffusion coefficients. We are interested in the conditions that lead to the emergence of a spatially patterned solution in such a system. In general, spatial patterns arise in such systems through Turing instability, a symmetry breaking mechanism in which an otherwise stable spatially uniform steady state is destabilized by the addition of diffusion and cross-diffusion terms. To derive conditions for such an instability, we perform a classic Turing type linear stability analysis on the approximated system (3.16).

We expand our system (3.16) about a spatially uniform positive steady state $(u^{*},v^{*})$ with $u^{*}>0$ and $v^{*}>0$. Substituting

[TABLE]

into (3.16) and dropping any nonlinear terms, the resulting linearized system about $(u^{*},v^{*})$ becomes:

[TABLE]

with

[TABLE]

where

[TABLE]

Here, we consider a cooperative form of $\mathbf{J}$ with $f_{u},g_{v}<0$ and $f_{v},g_{u}>0$:

[TABLE]

Note that this is different from the classic Turing model activator-inhibitor form of $\mathbf{J}$ where $f_{u}$ and $g_{v}$ are of opposite sign.

Following standard convention, we let

[TABLE]

Here, $a$ and $b$ are constants, and $k$ is the corresponding wavenumber, with $1/k$ being proportional to the wavelength of the emergent patterns. Since $e^{ikx}$ is periodic and bounded, the sign of $\alpha$ plays an important role in determining whether these small perturbations away from the steady state will grow or decay.

Substituting (3.21) into (3.17) and looking for a nontrivial solution, we require

[TABLE]

This yields the following dispersion relation:

[TABLE]

where

[TABLE]

Using this dispersion relation, we can then derive conditions for Turing type instability, summarized in the following theorem:

Theorem 3.4**.**

Let $(u^{*},v^{*})$ be a constant steady state solution of (3.16) with $\mathbf{D}$ defined as in (3.18) with $d_{11},d_{22}>0$ and $P_{2}<0$ and $P_{4}<0$ defined in (3.15). Also, let $\mathbf{J}$ be defined as in (3.20) with $f_{u},g_{v}<0$ and $f_{v},g_{u}>0$. If

[TABLE]

*then $(u^{*},v^{*})$ is locally asymptotically stable with respect to the corresponding ODE system, but is unstable with respect to system (3.16). *

Proof.

The solution to (3.23) yields:

[TABLE]

Note that $k^{2}=0$ corresponds to the spatially homogeneous case. For Turing instability, we require the spatially homogeneous state $(u^{*},v^{*})$ to be stable in the absence of spatial variation ($k^{2}=0$). Therefore, for $k^{2}=0$, the eigenvalues given in (3.26) have to be negative. This occurs when the trace of $\mathbf{J}$ is negative and the determinant of $\mathbf{J}$ is positive. From the special form of our matrix $\mathbf{J}$ in (3.20), it is clear that the trace of $\mathbf{J}$ is always negative, and from the first assumption in (3.25) know that the determinant of $\mathbf{J}$ is positive. So $(u^{*},v^{*})$ is locally asymptotically stable with respect to the corresponding ODE system.

For the emergence of a non-constant spatially periodic solution, we further require that for some $k^{2}\neq 0$, Re$(\alpha_{+}(k^{2}))>0$, guaranteeing that the perturbation will grow. Since $\operatorname{Tr}{\mathbf{D}}>0$ and $\operatorname{Tr}{\mathbf{J}}<0$, a necessary but not sufficient condition is that

[TABLE]

This happens as long as the following condition is satisfied:

[TABLE]

Under the condition (3.27), the minimum of $c(k^{2})$ is achieved at some $k^{2}=k_{m}^{2}>0$. Minimizing $c(k^{2})$ with respect to $k^{2}$ yields:

[TABLE]

Guaranteeing that $c_{min}<0$, we then have the following final condition:

[TABLE]

∎

Now, the same range of wavenumbers $k$ that makes $c(k^{2})<0$ in (3.24) also guarantees that Re$(\alpha(k^{2}))>0$. We can further calculate the relevant range of wavenumbers $k_{-}^{2}<k^{2}<k_{+}^{2}$ by computing the zeros of the function $c(k^{2})$ such that $c(k_{-}^{2})=c(k_{+}^{2})=0$. Then

[TABLE]

where

[TABLE]

The spatial patterns that emerge have a corresponding wavelength $\omega$, defined as

[TABLE]

with $k_{m}$ defined in (3.28) in the interval (3.30) and the one for which the positive eigenvalue $\alpha_{+}(k^{2})$ from (3.26) achieves a maximum, corresponding to the most unstable and fastest growing mode.

We can gain further insight into the result from Theorem 3.4 by visualizing the instability conditions in the $P_{2}P_{4}$-plane (Figure 6). Letting

[TABLE]

the second stability condition from Theorem 3.4 is equivalent to

[TABLE]

Further rearrangement of (3.33) leads to the following condition for the instability of the uniform solution:

[TABLE]

Additionally, we have

[TABLE]

Note that the first stability condition from (3.25) in Theorem 3.4 is independent of $P_{2},P_{4}$ and $\lambda$. Therefore, given that this first condition holds, rearranging the other instability condition in Theorem 3.4 in the form of (3.34) and using (3.35) as well as the fact that $M<0$ and $N>0$, it is clear that the instability region corresponds to the area to the left of the downward facing parabola defined on the right hand side of (3.34) in the fourth quadrant of the $P_{2}P_{4}$-plane (Figure 6). Without the activator-inhibitor dynamics, a cooperative system cannot be destabilized by diffusion alone. Theorem 3.4 makes it clear that the additional cross-diffusion terms given by $P_{2}$ and $P_{4}$ make spatial heterogeneous patterns possible in this cooperative system. In particular, the cross-diffusion term with $P_{2}$ plays a crucial role in the pattern forming mechanism since in its absence ($P_{2}=0$), the term $c(k^{2})$ from the dispersion relation in (3.23) can never be negative for any $k^{2}$. Since the biharmonic parameter $P_{4}$ acts as a stabilizing force, we also note that as its absolute value increases, the window for spatial patterns decreases (Figure 6). Increasing the value of the strength parameter $\lambda$ offsets the effect of the biharmonic parameter $P_{4}$ and increases the size of the window in which spatial patterns are possible. We note that in the absence of the biharmonic long-range cross-diffusion term ($P_{4}=0$), the conditions for Turing instability can still be satisfied. In this case, the system is reduced to a special case of the reaction-diffusion model with cross-diffusion, for which Turing instability conditions have been previously derived (Madzvamuse et al., 2015).

Previous results in this section took into consideration the system on an infinite domain $\mathbb{R}$. In such a system, we will always find an unstable mode $k^{2}$ in the interval (3.30) if the conditions in Theorem 3.4 are satisfied. Numerical simulations require the choice of a finite domain with specific boundary conditions. Therefore, we now consider the scenario on a bounded domain $T=(-l,l)$ with periodic boundary conditions, where the size of $T$ also affects the pattern formation. This is a more restrictive situation than the infinite domain scenario as the wavenumbers $k$ are now discrete and depend on the size of the domain. In this case, we shall understand that the solution $(u,v)$ on $T$ are periodically extended to $\mathbb{R}$ so the integral terms in the original system is still integrated on $\mathbb{R}$.

The result for the bounded domain case is summarized in the following corollary:

Corollary 3.5**.**

Consider (3.16) on a finite domain $T=\{x\in\mathbb{R}:-l<x<l\}$ and the following periodic boundary conditions:

[TABLE]

Let $(u^{*},v^{*})$ be a constant steady state solution of (3.16) with $\mathbf{D}$ defined as in (3.18) with $d_{11},d_{22}>0$ and $P_{2}<0$ and $P_{4}<0$ defined in (3.15). Also, let $\mathbf{J}$ be defined as in (3.20) with $f_{u},g_{v}<0$ and $f_{v},g_{u}>0$. If

[TABLE]

then $(u^{*},v^{*})$ is locally asymptotically stable with respect to the corresponding ODE system, but is unstable with respect to system (3.16) on $T$ with boundary condition (3.36). Moreover the most unstable mode is given by $n\in\mathbb{N}$ such that

[TABLE]

*where $\alpha_{+}(k^{2})$ is defined in (3.26), and the corresponding wavelength is $\omega=2\pi/k_{m}=2l/n$. *

Proof.

The non-constant eigenfunctions that satisfy the corresponding eigenvalue problem $\phi^{\prime\prime}+\lambda\phi=0$ on the domain $(-l,l)$ with periodic boundary conditions are of the following form

[TABLE]

and the corresponding eigenvalues are $k_{i}=(i\pi/l)^{2}$ for $i\in\mathbb{N}$. Now, when $0<k_{-}^{2}<k_{i}^{2}=(i\pi/l)^{2}<k_{+}^{2}$ for some $i\in$ $\mathbb{N}$, where $k_{-}$ and $k_{+}$ are defined in (3.30), the eigenvalue $\alpha_{+}(k_{i}^{2})$ defined in (3.23) is positive for this $i$.

We then note that the discrete wavenumber $k$ increases by $\pi/l$ with each $i$. Therefore, to guarantee that we have at least one $k^{2}=(i\pi/l)^{2}$ in the interval given by $(k_{-}^{2},k_{+}^{2})$, it is sufficient that the length of the interval $(k_{-},k_{+})$ is larger than $\pi/l$ (Shi et al., 2011). Using

[TABLE]

and the expressions of $k_{-}$ and $k_{+}$ in (3.30), we obtain the second instability condition in (3.37). Now, for an interval of length $2l$, if the instability conditions (3.37) are satisfied, then a spatially patterned solution will emerge with the corresponding wavenumber $k=i\pi/l$ such that $k^{2}\in(k_{-}^{2},k_{+}^{2})$. The most unstable wavenumber $k_{m}$ is the one that maximizes $\alpha_{+}(k^{2})$ in (3.26). ∎

The second instability condition in (3.37) also defines a minimal length $l_{m}$ for the emergence of the spatial patterns:

[TABLE]

This implies that in numerical simulations, if one chooses $l<l_{m}$, then no spatial patterns can be observed. On the other hand, when the length $l$ is large, then the interval $(k_{-},k_{+})$ may contain multiple unstable wavenumbers $k=i\pi/l$, and the spatial patterns with all these wavenumebrs are possible but the one with most unstable wavenumber $k_{m}$ is the one most likely to be observed.

3.3 Grass-Sediment Cooperative System with Nonlocal Interactions

We now apply these results to our Grass-Sediment system (2.4). The biharmonic approximation yields the following approximated system:

[TABLE]

with $F(S),L(G),P_{2}$ and $P_{4}$ defined previously in Sections 2 and 3.2. From Section 3.1, at the stable uniform positive steady state $(G_{+}^{*},S_{+}^{*})$ defined in (3.9), we have:

[TABLE]

Note that the cooperative form of this system with $f_{v},g_{u}>0$ and $f_{u},g_{v}<0$. For numerical simulations, we consider this system on a finite domain $T=(-l,l)$ and the following periodic boundary conditions:

[TABLE]

Using the results from Theorem 3.4 and Corollary 3.5, we have the following condition necessary for Turing instability on $T$:

[TABLE]

where

[TABLE]

For Case I parameter regime from Section 3.1, we choose: $e_{1}=0.7,p_{1}=3.5,f=0.14,e_{3}=0.1140,\delta=1/7$. We then choose $D_{G}=0.04,D_{S}=0.6,\lambda=40$ for our nonlocal parameter values to satisfy the instability conditions (3.40) and numerically integrate the biharmonic system (3.39) on $T=(-l,l)=(-7\pi,7\pi)$. We find that a spatially patterned solution emerges, as predicted (Figure 7) and these simulations are also consistent with numerical simulations of the original system (2.4), suggesting that the theoretical results derived from the biharmonic system can be applied to the original system to give insight regarding under what conditions a spatially patterned solution emerges. The eigenvalue $\alpha_{+}(k^{2})$ is given by:

[TABLE]

Furthermore, on the domain $T=(-7\pi,7\pi)$, the range of wavenumbers for which the corresponding eigenvalue $\alpha_{+}(k^{2})$ is positive is given by:

[TABLE]

It can be calculated that for $2\leq i\leq 12$, (3.42) is satisfied, and when $i=7$, $\alpha_{+}(k^{2})$ is maximized. Hence the characteristic wavelength of the emerging patterns is $\omega=2l/7=2\pi$. This is consistent with simulation results which show $7$ peaks (Figure 7). Similar results are obtained for Case II parameter regime (see Figure 11 in the Appendix.)

Previously, we used biologically realistic parameters from Table 1 (see the Appendix) to perform all numerical simulations, including realistic parameters for the scale-dependent feedback ($P_{2}$, $P_{4}$, $\lambda$). Now, we are interested in how varying these scale-dependent feedback parameters may affect the nature of the spatial patterns in system (3.39). As predicted in Section 3.2, since the biharmonic term $P_{4}$ acts as a stabilizing force, as its value gets larger, the window for spatial patterns decreases and a larger value of $\lambda$ is necessary to offset its effect and allow spatial patterns to emerge (Figure 8a). In addition, choosing a larger value of $P_{4}$ results in an overall increase in the pattern wavelength (Figure 8a). This result can also be interpreted in the context of how the coefficients $P_{2}$ and $P_{4}$ are related to the shape of the kernel in (2.2) in the original system (2.4) (Figure 8b). The coefficient $P_{2}$ measures the difference of the variances $\sigma_{1}$ and $\sigma_{2}$ of the excitatory and inhibitory interactions, respectively. The coefficient $P_{4}$ is related to kurtosis and controls the weight of the kernel’s tails while $\lambda$ modulates the amplitude of the Mexican-hat kernel. Since for a fixed value of $P_{2}$, an increase in $P_{4}$ results in a wider, flatter kernel shape, the wider the range of the long-range effects given by $P_{4}$, the stronger these interactions need to be (given by $\lambda$) to have a significant effect and lead to the formation of spatial patterns. This makes biological sense, since the intensity of scale-dependent interactions tend to dissipate over larger distances and therefore need to be amplified to have any effect on spatial heterogeneity over longer ranges. In addition, we see that wider kernels result in wider spatial patterns characterized by longer wavelengths. Again, this makes biological sense as one would expect the scale of the spatial interactions to influence the resulting spatial patterns.

Finally, we compare our analytic results with numerical simulations of the approximated biharmonic system (3.39) and the original system (2.4) (Figure 9). The analytic results from (3.40) are consistent with numerical simulations of the biharmonic system and the original system. However, we note that the onset of patterns in the original system occurs sooner than in the biharmonic system. Nonetheless, this result suggests that using the biharmonic system can help find the relevant parameter regime in which spatial patterns are possible in the original system and gain understanding into how the nature of the scale-dependent feedbacks affects the development of spatial patterns.

4 Discussion

We propose a phenomenological model to describe the dynamics of the marsh edge in terms of two-way interactions between marsh grass Spartina alternifora and sedimentation. In nature, the marsh edge can frequently be observed in a number of configurations ranging from a spatially uniform to a more wave-like shoreline. The interest of this paper lies in understanding whether the well-known scale-dependent (nonlocal) feedback between marsh vegetation and sedimentation can lead to spatially variable shoreline configurations. Marsh grass promotes sediment accretion in its immediate surroundings by slowing down current acts as a facilitation mechanism. In turn, the diverted water flow contributes to increased erosion further away and acts as an inhibitory mechanism. We propose a system of reaction-diffusion equations with an additional integral term with a Mexican-hat kernel function that describes the nature of this scale-dependent feedback. Our system is highly cooperative; as cooperative systems often lack the classic activator-inhibitor mechanism necessary for pattern formation, it becomes of interest how and under what conditions spatial patterns may develop.

We perform a biharmonic approximation of our system and carry out analysis on the simpler biharmonic system that expresses the kernel function as separate short-range and long-range diffusion terms. Using the more mathematically tractable biharmonic system, we are then able to derive general condition for the formation of spatial patterns in a cooperative system such as ours. Further, using numerical simulations, we confirm that the biharmonic model, while an approximation, is consistent with the original model, and therefore we can apply the theoretical results from the biharmonic system to help gain insight into the formation of patterns in the original system. We parameterize the kernel function using a set of reasonable parameters from literature and find that spatial patterns can develop, given that the scale-dependent interactions between marsh vegetation and sediment dynamics are strong enough. The model thus provides further evidence that the presence of scale-dependent interactions is essential for pattern formation and that heterogeneous patterns cannot occur in the presence of weak scale-dependent interactions. Not surprisingly, we find that the choice of wider kernels tend to produce wider spatial patterns (characterized by longer wavelengths) and vice versa. The nature and strength of the grass-sediment scale-dependent interactions depends on many factors such as the underlying hydrodynamics and sediment composition, the exact spatial scale (corresponding to the widths of the Mexican-hat kernel) and relative strength of the scale-dependent feedback are difficult to estimate in the field and can vary substantially. We use one possible set of biologically realistic parameters for the kernel function (Table 1) and find that the patterns that emerge in simulations occur on a spatial scale consistent with what can observed in nature (4-10 meters between peaks) (Vandenbruwaene et al., 2011).

Furthermore, we find that there are two possible parameter regimes in the system. The first regime is especially of interest as it corresponds to a more realistic scenario where marsh vegetation is effective at attenuating erosion through the binding of sediment and decreasing the effect of wave erosion. Given the strong facilitatory nature of the grass-sediment interactions, bistability takes place in this parameter regime. In general, bistable dynamics makes a system especially prone to collapsing to an irreversible state as environmental conditions gradually worsen and a tipping point is reached (Dakos et al., 2011; Kéfi et al., 2014, 2016) through the phenomenon of hysteresis. Pattern formation has previously been suggested as a possible coping mechanism for systems close to degradation (Chen et al., 2015). The analysis in this paper gives more insight into this previously reported phenomenon as we also find this to be the case in our model (Zaytseva et al., 2018) where pattern formation allows the marsh edge to cope with harsher erosion through spatial variation. One limitation of our model is the lack of multiple spatial dimensions as only the dynamics on a one-dimensional cross-section of the marsh edge were considered. Hence, we were not able to observe the geometry of the protrusions. In addition, the model is meant to be phenomenological in nature, omitting processes such as the effect and variation of hydrodynamics and wave action, modeled in more detail previously (Fagherazzi et al., 2012). Despite the relatively simple dynamics of our one-dimensional model, it is able to capture the pattern formation on the marsh edge as a result of scale-dependent feedbacks between vegetation and sediment accumulation. The agreement between the model simulations and field observations suggests that important pattern-generating processes have been captured in the model and non-local interactions between plants and sedimentation can drive the formation of shoreline patterns. In addition, the results in this paper can be generalized to any cooperative system with scale-dependent feedbacks in the form of short-range activation and long-range inhibition, described using a Mexican-hat kernel function.

5 Appendix

Figure 10 shows the plot of the functional forms of $F(S)$ and $L(G)$.

Figure (11) shows numerical simulations of both the biharmonic system (3.39) and the original system (2.4) for the parameter regime in Case II from Section 3.1. We see that a spatially patterned solution emerges if the instability conditions in (3.40) are satisfied.

Table 1 shows the biologically realistic parameters for the original system and their sources. We use the parameter values from Table 1 to obtain the new re-scaled parameters from Section 2 to use in all numerical simulations performed in this paper.

All numerical simulations in this paper are performed using MATLAB. We use an implicit finite differencing scheme to numerically integrate the original equation. Although this scheme is more computationally intensive, it is chosen because it is always numerically stable and convergent. Because domain size plays an important role in the system’s ability to form patterns, a large enough domain has to be chosen to be able to fit patterns with their characteristic wavelength. We evaluate all integrals using the trapz function in MATLAB, which performs numerical integration using the trapezoidal rule. For the convolution term, we evaluate the integral of the product of the kernel and the periodically extended solution on the interval $(-3l,3l)$, to make sure an adequate number of kernels are considered in calculating the net effect. To numerically integrate the biharmonic system, we use an explicit finite differencing scheme in MATLAB. This scheme is less computationally intensive, and is easier to implement, given the extra biharmonic term. For both models, the numerical simulations are performed on a spatial domain $(-l,l)$ with $l=7\pi$ with periodic boundary conditions. We apply Turing’s idea of diffusion driven instability and use a spatially periodic perturbation of the stable steady state of the corresponding system of ODEs as the initial condition for our simulations.