Rise and Fall of Periodic Patterns for a Generalized Klausmeier–Gray–Scott Model

In this paper we introduce a conceptual model for vegetation patterns that generalizes the Klausmeier model for semi-arid ecosystems on a sloped terrain (Klausmeier in Science 284:1826–1828, 1999). Our model not only incorporates downhill flow, but also linear or nonlinear diffusion for the water component.

To relate the model to observations and simulations in ecology, we first consider the onset of pattern formation through a Turing or a Turing–Hopf bifurcation. We perform a Ginzburg–Landau analysis to study the weakly nonlinear evolution of small amplitude patterns and we show that the Turing/Turing–Hopf bifurcation is supercritical under realistic circumstances.

In the second part we numerically construct Busse balloons to further follow the family of stable spatially periodic (vegetation) patterns. We find that destabilization (and thus desertification) can be caused by three different mechanisms: fold, Hopf and sideband instability, and show that the Hopf instability can no longer occur when the gradient of the domain is above a certain threshold. We encounter a number of intriguing phenomena, such as a ‘Hopf dance’ and a fine structure of sideband instabilities. Finally, we conclude that there exists no decisive qualitative difference between the Busse balloons for the model with standard diffusion and the Busse balloons for the model with nonlinear diffusion.

1. Notice however the extra b ² in the coefficient in front of the nonlinear term. In Morgan et al. (2000), b is scaled out of the matrix \(\mathcal{M}_{c}\) in formula (3.24). That is, the matrix \(b\mathcal{M}_{c}\) in Morgan et al. (2000) plays the role of our matrix \(\mathcal{M}_{\mathrm{i}\omega_{c}}(a_{c},0,\mathrm {i}k_{c})\). This is equivalent to scaling \(\mathcal{A}\to b\mathcal{A}\) in (2.47).

2. Notice that this rescaling of the Klausmeier system can be acquired from Klausmeier’s original nondimensional system (see Klausmeier 1999),

   $$ \everymath{\displaystyle} \left\{ \begin{array}{@{}rcl} u_{T} & = & \nu u_X + a-u - uv^2; \\[1mm] v_{T} & = & \delta^{2\sigma}v_{XX} - mv + u v^2, \end{array} \right. $$

   (2.52)

   by rescaling with \(x=\frac{a^{2}}{\nu}X\), t=a ² T, v=aV, u=aU, \(A=\frac{1}{a^{2}}\), \(B=\frac{m}{a^{2}}\) and by further introducing \(0<\delta^{\sigma}:=\frac{a}{\nu}\ll1\). (From the estimates for a and ν in Klausmeier (1999), it can be deduced that indeed \(0<\frac{a}{\nu}\ll1\).)

3. Note that we silently switched back to the original parameters A, B, and C in (1.5). We will comment on the relation between A, B, and C on the one hand and a, b, and c on the other hand in Sect. 3.2.

4. With slight abuse of terminology, in the context of perturbations of periodic patterns we abbreviate the Turing–Hopf instability to ‘Hopf instability’ in the rest of this paper.

The authors thank Max Rietkerk for sharing his insights and the stimulating discussions.

Authors and Affiliations

1. Korteweg-de Vries Instituut, Universiteit van Amsterdam, P.O. Box 94248, 1090 GE, Amsterdam, The Netherlands

   Sjors van der Stelt & Geertje Hek

2. Centrum Wiskunde en Informatica (CWI), Science Park 123, 1098 XG, Amsterdam, The Netherlands

   Sjors van der Stelt & Jens D. M. Rademacher

3. Mathematisch Instituut, Universiteit Leiden, P.O. Box 9512, 2300 RA, Leiden, The Netherlands

   Arjen Doelman

4. Institut International de Lancy, Avenue Eugène Lance 24, 1218, Geneva, Switzerland

   Geertje Hek

Corresponding author

Correspondence to Sjors van der Stelt.

Communicated by Alan R. Champneys.

Appendix: Derivation of the Ginzburg–Landau Equation

In this appendix, we outline the derivation of the Ginzburg–Landau equation for the amplitude \(\mathcal{A}\) of the pattern that appears at the Turing–Hopf bifurcation. Each of the four different cases of Fig. 1 can be derived from the expressions given in this appendix by considering either γ=1 or c=0, or both.

In Appendix A.1 we derive the Ginzburg–Landau equation for the special case that γ=1 and \(c=\sqrt{\frac{2}{3}b}\) that was presented in Sect. 2.5.3.

The Ginzburg–Landau Ansatz for patterns that emerge at the Turing–Hopf instability can, for the rescaled GKGS-system (2.38), be written as

(A.1)

where \(\mathcal{A}\) and X _ij are functions of ξ=εx and τ=ε ²(x−c _g t) and c _g the group velocity defined by (2.26). Substituting this expansion in the GKGS-system (2.38) and collecting terms of equal powers of ε and the Fourier modes \(\mathrm{e}^{\mathrm{i}(k_{*} x + \omega_{*} t)}\), we derive expressions for X _02,12,22,13 and Y _02,12,22,13 subsequently. Notice that the scaling in (2.38) has the advantage that the terms of order ε ² only play a role in the equations for X ₁₃ and Y ₁₃.

As mentioned in paragraph 2.4, the solvability condition can be applied to solve an equation of the form

The equations for X _1j ,Y _1j , j=2,3 are of this form, with

(A.3)

We briefly point out the construction of the set of solutions for (A.2). The construction for c=0 differs from that for c≠0.

If c≠0, the matrix in (A.3) has two eigenvalues, λ ₊=0 and

If \(y\in\operatorname{Rg} \mathcal{M}_{\mathrm{i}\omega_{*}}(a_{*},k_{*},c)\) and \(c\not =0\), the set of solutions to (A.2) is given by

On the other hand, if c=0, we know from Proposition 1 that

It is then straightforward to show that

(A.5)

and that λ ₋=0 if γ=1. It is also straightforward to show that both columns of \(\mathcal {M}_{\mathrm{i}\omega _{*}}(a_{*},k_{*},0)\) span the range. We call the second column v ₁. So if y from (A.2) lies in the range, there exists an α∈ℝ such that y=αv ₁. Hence, if c=0, the set of solutions to (A.2) can be presented as

(A.6)

By plugging in the expansion (A.1) into (2.38) one obtains at order \(\mathcal{O}(\varepsilon)\) an equation for (X ₀₂,Y ₀₂)^T,

We use shorthands x ₀₂,y ₀₂ for \(X_{02}=x_{02}|\mathcal{A}|^{2}\), \(Y_{02}=y_{02}|\mathcal{A}|^{2}\). The values for x ₀₂ and y ₀₂ can be read from (A.7). At order \(\mathcal{O}(\varepsilon E)\), we find equations of the form

which can be solved if \((x_{1j},y_{1j})^{\mathrm{T}}\in\operatorname{Rg} \mathcal {M}_{\mathrm{i}\omega _{*}}(a_{*},k_{*},c)\). We find

It can be checked that, indeed, \((x_{1j},y_{1j})^{\mathrm{T}}\in\operatorname{Rg} \mathcal{M}_{\mathrm{i}\omega_{*}}(a_{*},k_{*},c)\). At order \(\mathcal{O}(\varepsilon E^{2})\) we find

with

At order ε ² E, we obtain equations for X ₁₃ and Y ₁₃. These equations can be written as

The right-hand sides I ₁,I ₂ are built up by several terms. The nonlinear terms from the reaction kinetics generate

We define L _tot as the sum of these expressions:

The nonlinear terms that appear from working out the nonlinear diffusion terms generate

Based on this we define

We then obtain for the right-hand side of the system

To derive the Ginzburg–Landau equation, we impose the solvability condition (2.42) to (A.10):

and obtain

1.1 A.1 The Special Case that γ=1 and \(c=\sqrt{\frac {2}{3}b}\)

In this section we present the expressions for x _ij ,y _ij , ij=02,12,22 for the special case that γ=1 and \(c=\sqrt{\frac {2}{3}b}\). In Proposition 2 we computed that for γ=1 and \(c=\sqrt{\frac{2}{3}b}\) one has

Also, one computes the second component of a basis vector of the kernel of \(\mathcal{M}_{\omega_{*}}(a_{*},k_{*},c)\) and the nonzero eigenvalue of \(\mathcal{M}_{\omega_{*}}(a_{*},k_{*},c)\) as

These values are used to derive

The nonlinear terms from the reaction kinetics are

The sum of these expressions is

The nonlinear diffusion terms are, of course, zero, so \(L_{\mathrm {NLD}}=L_{\mathcal{A},\mathrm{NLD}}=0\). We get for the right hand components as in (A.10):

These give the Ginzburg–Landau equation in (2.51):

1.2 A.2 Derivation of the Ginzburg–Landau Equation for the GKGS-System with c=0

In this section we present the expressions for x _ij ,y _ij , ij=02,12,22 for the special case that c=0. In Proposition 1 we computed that for c=0 one has

Also, one computes the second component of a basis vector of the kernel of \(\mathcal{M}_{\omega_{*}}(a_{c},k_{c},c)\) and the nonzero eigenvalue of \(\mathcal{M}_{\omega_{*}}(a_{c},k_{c},c)\) as

These values are used to derive

The sum of the nonlinear terms from the kinetics is

Also we have

This gives the Ginzburg–Landau equation in (2.49):

with

1.3 A.3 Derivation of the Ginzburg–Landau Equation for the Case c≫1: The Klausmeier Model and the GKGS-Model for c≫1

This appendix to Sect. 2.6 deals with an elaborate account on the scalings introduced in Sect. 2.6 that were used to derive the Klausmeier system (2.53) from the GKGS-system. Also, we derive the GLE for the Klausmeier system (2.53).

Scaling Analysis for the Klausmeier System as a Limit Case of the GKGS System

We remark that the equilibria for both systems (1.5) and (2.53) are the same. Patterns close to the equilibrium (U ₊,V ₊) can be described as

Substitution of these expansions in (1.5) gives the leading order formulation (2.38). We are interested in the behavior of the GKGS-model for \(0<1/\sqrt{c}\ll\varepsilon\ll 1\). Since we know from Proposition 2 that \(a_{c}=\mathcal {O}(\sqrt{c})\), we put \(a_{*} = \bar{a}_{*} \sqrt{c}\) and obtain for (2.38),

In order to derive the Klausmeier model, we must scale the components U and V such that the diffusion coefficient in the first component of (2.38) is of higher order in \(1/\sqrt{c}\) than the other terms in the equations. The other terms must balance at the same, highest order. In order to obtain this, we scale U, V, and r such that

With these scalings we obtain for (A.15), by neglecting higher orders of δ and \(1/\sqrt{c}\),

We can now scale out b by putting

and obtain, to leading order in ε and neglecting higher order terms of δ and \(\frac{1}{\sqrt{c}}\),

This system is the one presented in (2.55).

Derivation of the GLE for the Klausmeier System: The Regime \(0<1/\sqrt{c}\ll\varepsilon\ll1\)

From (A.19) one derives the dispersion relation associated to the linearization about the background state (U ₊,V ₊) in the Klausmeier model, neglecting higher orders of ε,

We apply conditions (2.7a), (2.7b) to derive critical parameters. Working out the dispersion relation (A.20) using condition (2.7a),

(A.21a)

(A.21b)

From these relations one derives

and by solving equation (A.21a) for \(\tilde{\omega}\) we get

and thus

Differentiating (A.24) with respect to \(\tilde{k}\) and substituting equation (A.21b) into the result, we get

Solving (A.22) and (A.25) for \(\tilde{a}\) and \(\tilde{k}\) then gives

which are the expressions for large c that we had derived in Proposition 2. From these expressions for \(\tilde{a}_{*}\) and \(\tilde{k}_{*}\) we further derive the critical frequency and the group speed

From (A.20) it follows that the kernel and range of the linearization about the equilibrium \((\tilde{U}_{+},\tilde{V}_{+})\) equal

and

From the expression for the range of \(\mathcal{M}_{\mathrm{i}\tilde{\omega}_{*}}(\tilde{a}_{*}, \tilde{k}_{*})\) we derive that the equations

can be solved for x if and only if \(f\in\operatorname{Rg} \mathcal {M}_{\mathrm{i}\tilde{\omega}_{*}}(\tilde{a}_{*}, \tilde{k}_{*})\), that is, if f fulfills the solvability condition

Since \(\mathrm{Det}\,\mathcal{M}_{\mathrm{i}\tilde{\omega}_{*}}(\tilde{a}_{*}, \tilde{k}_{*}) = 0\), it follows that the unique solution to the equation \(\mathcal{M}_{\mathrm{i}\tilde{\omega}_{*}}(\tilde{a}_{*}, \tilde{k}_{*}) x = f\) is

with λ ₋ the nonzero eigenvalue of \(\mathcal{M}_{\mathrm{i}\tilde{\omega}_{*}}(\tilde{a}_{*}, \tilde{k}_{*})\),

By using (A.28), the expansion (U,V) that describes the pattern near its onset (i.e., for \(\tilde{a} = \tilde{a}_{*} - \tilde{r}\varepsilon^{2}\)) can be written out as

with the shorthand

As pointed out in at the begin of the Appendix, in the Ginzburg–Landau formalism one subsequently derives equations of the form \(X_{02}=x_{02}\mathcal{A}^{2}, Y_{02}=y_{02}\mathcal{A}^{2}, X_{12}=x_{12}\mathcal{A}_{\xi}, Y_{12}=y_{12}\mathcal{A}_{\xi}, X_{22}=x_{22}|\mathcal{A} |^{2}\), and \(Y_{22}=y_{22}|\mathcal{A}|^{2}\). The formulas for x _ij and y _ij are derived by substituting the expansion (A.31) in the leading order system (A.19) and collecting terms of order ε ^j−1 E ⁱ (with shorthand \(E=\mathrm{e}^{\mathrm{i}(\tilde{k}_{c} x + \tilde{\omega}_{c} t)}\)) and solving them for X _ij and Y _ij :

The Ginzburg–Landau equation that describes the onset of patterns in the Klausmeier system (A.19) reads

(A.33)

with

If one works out the parameter values for \(\tilde{a}_{*}^{2} = \sqrt{2}-1\) and \(\tilde{k}_{*}^{2} = \sqrt{2}-1\), the leading order system (A.19) becomes

The matrix that describes the linear leading order part of this system is then

Working out the levels for the different expressions (A.32), we get

and

The Ginzburg–Landau equation then becomes

or, equivalently,

which is the equation presented in (2.56).

The GLE for the GKGS Model for c≫1: The Regime \(0<\varepsilon \ll1/\sqrt{c}\ll1\)

As in the previous section, we scale the leading order system of the GKGS model according to (2.54). We then obtain the leading order system (2.57). To first order in ε, the leading order system (2.57) reads

First, we remark that for γ≥1 the linear part of the nonlinear diffusion in the GKGS-model is in leading order \({\leq}\mathcal{O}(c^{-1/2})\). Secondly, we remark that to leading order in c,

with \(\mathcal{M}_{\lambda}(\tilde{a}_{*},\mathrm{i}\tilde{k})\) as defined in (A.20). Therefore, the critical \(\tilde{k}_{*}\), \(\tilde{a}_{*}\), \(\tilde{\omega}_{*}\), and \(\tilde{c}_{g}\) are to leading order in c as in (A.26) and (A.27).

The solvability condition is as in (A.29). Using the leading order system (2.57), we compute

In (A.38) it is understood that all l _i , i=0,…,4 do not depend on c and that θ _i ≥0 for i=0,…,4. We have not computed the l _i and θ _i explicitly. This gives for the GLE of the GKGS-model in general form

(A.39)

with

We have not bothered about calculating the constant denoted by “const”, since for asymptotically large c the associated expressions only play a role at higher order. That is, for asymptotically large c≫1 it is immediate that (A.39) reduces to (A.33). Using the expressions (A.38) and inserting \(\tilde{k}\), \(\tilde{a}\), \(\tilde{\omega}\), and c _g one obtains the GLE for the Klausmeier system (2.56).

Reprints and permissions