Existence and influence of mixed states in a model of vegetation patterns

The Rietkerk vegetation model is a system of partial differential equations, which has been used to understand the formation and dynamics of spatial patterns in vegetation ecosystems, including desertification and biodiversity loss. Here, we provide an in-depth bifurcation analysis of the vegetation patterns produced by Rietkerk's model, based on a linear stability analysis of the homogeneous equilibrium of the system. Specifically, using a continuation method based on the Newton–Raphson algorithm, we obtain all the main heterogeneous equilibria for a given size of the domain. We confirm that inhomogeneous vegetated states can exist and be stable, even for a value of rainfall for which no vegetation exists in the non-spatialized system. In addition, we evidence the existence of a new type of equilibrium, which we call “mixed state”, in which the equilibria are always unstable and take the form of a mix of two equilibria from the main branches. Although these equilibria are unstable, they influence the dynamics of the transitions between distinct stable states by slowing down the evolution of the system when it passes close to it. Our approach proves to be a helpful way to assess the existence of tipping points in spatially extended systems and disentangle the fate of the system in the Busse balloon. Overall, our findings represent a significant step forward in understanding the behaviour of the Rietkerk model and the broader dynamics of vegetation patterns.

1 Introduction

In semi-arid regions, vegetation tends to be spatially organized around patterns (Barbier et al., 2006; Deblauwe et al., 2008, 2011, 2012). This phenomenon appears in various parts of the world where water is the limiting factor for plants' growth. Vegetation patterns can be modelled and explained with reaction–diffusion equations. Those types of equation exhibit the existence of a homogeneous stable equilibrium which is unstable to heterogeneous perturbation (Turing, 1952). In general terms, heterogeneous equilibria result from the joint effects of a short-range activation mechanism and long-range inhibition. For vegetation on ferruginous soil in semi-arid regions, the short-range activation effect is related to the positive feedback of vegetation on soil water availability, as vegetation limits water loss by runoff by enhancing water infiltration (Meron, 2015).

The rapid diffusion of surface water, however, acts against vegetation growth on bare soil by limiting water availability. The contrast between slow soil water diffusion and fast surface water diffusion produces the spatial, heterogeneous patterns. The mechanisms are described and captured in the Rietkerk model (Rietkerk et al., 2002). This model features three (prognostic) variables: biomass (B) (g m⁻²), soil water (W) (mm) and surface water (O) (mm); all are functions of time and space. The evolution of these quantities are governed by three equations:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial B}{\partial t}}& =c{g}_{\max }{\displaystyle \frac{WB}{W+k_{\mathrm{1}}}}-\mathrm{d}B+D_B\mathrm{\Delta }B,\\ {\displaystyle \frac{\partial W}{\partial t}}& =\mathit{\alpha }O{\displaystyle \frac{B+k_{\mathrm{2}}{w}_{\mathrm{0}}}{B+k_{\mathrm{2}}}}-g_{\max }{\displaystyle \frac{WB}{W+k_{\mathrm{1}}}}-r_{\mathrm{w}}W+D_W\mathrm{\Delta }W,\\ {\displaystyle \frac{\partial O}{\partial t}}& =R-\mathit{\alpha }O{\displaystyle \frac{B+k_{\mathrm{2}}{w}_{\mathrm{0}}}{B+k_{\mathrm{2}}}}+D_O\mathrm{\Delta }O,\end{array}\end{array}$$

where Δ is the Laplacian operator and R is the rainfall (mm d⁻¹). The rainfall is the external forcing of the system which we consider to be a spatially independent function. The first term in the biomass equation represents water uptake by the plant. The first term in the soil water equation is linked to the infiltration rate of water in the soil that is enhanced by the presence of biomass. The factors in front of the Laplacians (ΔB, ΔW and ΔO) are the diffusion constants of the different quantities. The diffusion constant for surface water is considered to be much bigger than those of biomass and soil water. This contrast is essential for pattern creation. In the following, the values of the parameters (see Table 1) are as in Rietkerk et al. (2002).

Table 1 Parameters for Rietkerk's model.

As explained above, the existence of reaction–diffusion processes enables stable equilibria in the form of patterns. Compared to a system without spatial dynamics, such pattern equilibria tend to broaden the range of rainfall compatible with the presence of vegetation.

The classical configuration for analysing such a system of equations uses periodic boundary conditions. A pattern is then defined as a spatially periodic equilibrium of the differential equations. Here, we will consider and analyse in depth a region of the parameter space called the Busse balloon (Busse, 1978). This is the region with heterogeneous equilibria, that is, the parameter space admitting at least one stable, spatially periodic equilibrium.

The motivating scientific question comes from the following observation. In a system without spatial dynamics – thus without diffusion – a catastrophic transition occurs between sustained vegetation and bare soil when rainfall decreases below a critical point. The transition point corresponds to a fold bifurcation and can be qualified as a tipping point (Lenton et al., 2008). Rietkerk et al. (2021) suggested that spatial dynamics and the existence of the Busse balloon smooth the transition between fully fledged vegetation and bare soil. In their terms, the Busse balloon “evades” the tipping point. This result implies that spatial dynamics effectively lower the precipitation threshold above which vegetation can be sustained.

Siteur et al. (2014) showed that a Busse balloon appears in the Klausmeier vegetation model (Klausmeier, 1999), where it occupies a region of the parameter space with lower rainfall than necessary to sustain homogeneous vegetation. However, the nature of this transition through the Busse balloon may be complex. For the non-spatial Rietkerk model, there is no such thing as a fold bifurcation. But we will see that adding a spatial component creates fold bifurcation.

The objective of the present study is, in the idealized context of a spatial domain of one dimension, to characterize the intermediate states that may emerge during the transition from full vegetation to bare soil, and to examine the dynamics that underlie potential transitions between these states.

Specifically, we provide an in-depth analysis of the Busse balloon, demonstrate the coexistence of multiple equilibria for a given rainfall intensity and foresee the circumstances which may trigger transitions between these different equilibria. In that sense, we propose an extension of the work by Zelnik et al. (2013), who compute partially the bifurcation diagram for Rietkerk's model. We also develop the method for finding equilibrium branches and characterize their stability. Finally, we highlight the existence of another type of equilibrium different than the regular patterns. We call them mixed state because of their shape and show how they can have an influence on the dynamics of the system.

2 Bifurcation diagram and stability

In this section, we present a method to construct the bifurcation diagram for Rietkerk's model and determine the form of the different equilibrium branches as a function of rainfall.

2.1 Linear analysis and Turing zone

The classical approach is to consider, as in Siero (2020), first, the static homogeneous equilibrium. To this end, we define the equilibrium $\overline{B}=B(x,t)$ where $\overline{B}$ is a constant in time and space – likewise with the other variables, which satisfy the relationships:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{0}=c{g}_{\max }{\displaystyle \frac{\overline{W}\overline{B}}{\overline{W}+k_{\mathrm{1}}}}-\mathrm{d}\overline{B},\text{(3)}& {\displaystyle }\mathrm{0}=\mathit{\alpha }O{\displaystyle \frac{\overline{B}+k_{\mathrm{2}}{w}_{\mathrm{0}}}{\overline{B}+k_{\mathrm{2}}}}-g_{\max }{\displaystyle \frac{\overline{W}\overline{B}}{\overline{W}+k_{\mathrm{1}}}}-r_{\mathrm{w}}\overline{W},\text{(4)}& {\displaystyle }\mathrm{0}=R-\mathit{\alpha }\overline{O}{\displaystyle \frac{\overline{B}+k_{\mathrm{2}}{w}_{\mathrm{0}}}{\overline{B}+k_{\mathrm{2}}}}.\end{array}$$

Two solutions exist depending on the value of R. One solution with vegetation is

$$\begin{array}{}\text{(5)}& {\displaystyle }\overline{B}=c\left({\displaystyle \frac{R}{d}}-{\displaystyle \frac{k_{\mathrm{1}}{r}_{\mathrm{w}}}{c{g}_{\max }-d}}\right),\text{(6)}& {\displaystyle }\overline{W}={\displaystyle \frac{\mathrm{d}{k}_{\mathrm{1}}}{c{g}_{\max }-d}},\text{(7)}& {\displaystyle }\overline{O}={\displaystyle \frac{R\left(\right(c{g}_{\max }-d\left)\right(cR+d{k}_{\mathrm{2}})-cd{k}_{\mathrm{1}}{r}_{\mathrm{w}})}{\mathit{\alpha }\left(\right(c{g}_{\max }-d\left)\right(cR+d{k}_{\mathrm{2}}{w}_{\mathrm{0}})-cd{k}_{\mathrm{1}}{r}_{\mathrm{w}})}}.\end{array}$$

Physical equilibria for positive parameters must be positive and hence satisfy the relations cg_max>d and $R>k_{\mathrm{1}}{r}_{\mathrm{w}}d/(c{g}_{\max }-d)$. The homogeneous equilibria without vegetation are, by definition, ${\overline{B}}_{\mathrm{0}}=\mathrm{0}$, which implies

$$\begin{array}{}\text{(8)}& {\displaystyle }{\overline{W}}_{\mathrm{0}}& {\displaystyle }={\displaystyle \frac{R}{r_{\mathrm{w}}}},\text{(9)}& {\displaystyle }{\overline{O}}_{\mathrm{0}}& {\displaystyle }={\displaystyle \frac{R}{w_{\mathrm{0}}\mathit{\alpha }}},\end{array}$$

which are, again, valid for positive parameters. For the parameters chosen here, the homogeneous vegetated equilibrium exists for R>1 and the non-vegetated equilibrium always exists.¹ Figure 1 shows these equilibria as a function of rainfall R.

Figure 1 Homogeneous equilibria of Rietkerk's model. From left to right are represented the biomass, the soil water and the surface water. The blue line is the equilibrium with vegetation and yellow is without.

We now consider the stability of these equilibria. Again, following standard practice we consider spatially periodic perturbations as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }B(x,t)& {\displaystyle }=\overline{B}+\mathit{ϵ}\mathit{\delta }B(x,t),\text{(11)}& {\displaystyle }W(x,t)& {\displaystyle }=\overline{W}+\mathit{ϵ}\mathit{\delta }W(x,t),\text{(12)}& {\displaystyle }O(x,t)& {\displaystyle }=\overline{O}+\mathit{ϵ}\mathit{\delta }O(x,t),\end{array}$$

with

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\delta }B(x,t)& {\displaystyle }=\mathit{\delta }B{e}^{\mathrm{\Omega }t+i\mathit{\kappa }x},\text{(14)}& {\displaystyle }\mathit{\delta }W(x,t)& {\displaystyle }=\mathit{\delta }W{e}^{\mathrm{\Omega }t+i\mathit{\kappa }x},\text{(15)}& {\displaystyle }\mathit{\delta }O(x,t)& {\displaystyle }=\mathit{\delta }O{e}^{\mathrm{\Omega }t+i\mathit{\kappa }x}.\end{array}$$

We now introduce the perturbation into the original equations, develop a Taylor expansion for small ϵ and get rid of the 0th-order term by using the vegetated homogeneous equilibrium. The linearized equations obtained can be recast as an eigenvalue problem:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}A& =\left(\begin{array}{ccc}d-\frac{c{g}_{\max }\overline{W}}{k_{\mathrm{1}}+\overline{W}}+D_B{\mathit{\kappa }}^{\mathrm{2}}& -\frac{c{g}_{\max }{k}_{\mathrm{1}}\overline{B}}{(k_{\mathrm{1}}+\overline{W}{)}^{\mathrm{2}}}& \mathrm{0}\\ \frac{g_{\max }\overline{W}}{k_{\mathrm{1}}+\overline{W}}+\frac{k_{\mathrm{2}}\mathit{\alpha }(w_{\mathrm{0}}-\mathrm{1})\overline{O}}{(k_{\mathrm{2}}+\overline{B}{)}^{\mathrm{2}}}& r_{\mathrm{w}}+\frac{g_{\max }{k}_{\mathrm{1}}\overline{B}}{(k_{\mathrm{1}}+\overline{W}{)}^{\mathrm{2}}}+D_W{\mathit{\kappa }}^{\mathrm{2}}& -\frac{\mathit{\alpha }(k_{\mathrm{2}}{w}_{\mathrm{0}}+\overline{B})}{k_{\mathrm{2}}+\overline{B}}\\ -\frac{k_{\mathrm{2}}\mathit{\alpha }(w_{\mathrm{0}}-\mathrm{1})\overline{O}}{(k_{\mathrm{2}}+\overline{B}{)}^{\mathrm{2}}}& \mathrm{0}& \frac{\mathit{\alpha }(k_{\mathrm{2}}{w}_{\mathrm{0}}+\overline{B})}{k_{\mathrm{2}}+\overline{B}}+D_O{\mathit{\kappa }}^{\mathrm{2}}\end{array}\right)\\ & \cdot \left(\begin{array}{c}\mathit{\delta }B\\ \mathit{\delta }W\\ \mathit{\delta }O\end{array}\right)=-\mathrm{\Omega }\left(\begin{array}{c}\mathit{\delta }B\\ \mathit{\delta }W\\ \mathit{\delta }O\end{array}\right).\end{array}\end{array}$$

Every pair of values (R,κ) defines a new eigenvalue problem. It corresponds to the linear dynamics obtained by perturbing a homogeneous equilibrium that exists at a particular value of R with wavenumber κ. We are interested in finding out whether there are values of (R,κ) for which there exists an exponentially growing solution and thus positive Ω. It is then said that the homogeneous equilibrium for a specific value of R is linearly unstable to perturbations of wavenumber κ. For the eigenvalue problem to have a non-trivial solution, the determinant det(A+Ω) must be 0. This leads to a cubic equation for Ω:

$$\begin{array}{}\text{(17)}& a_{\mathrm{1}}(R,\mathit{\kappa }){\mathrm{\Omega }}^{\mathrm{3}}+a_{\mathrm{2}}(R,\mathit{\kappa }){\mathrm{\Omega }}^{\mathrm{2}}+a_{\mathrm{3}}(R,\mathit{\kappa })\mathrm{\Omega }+a_{\mathrm{4}}(R,\mathit{\kappa })=\mathrm{0},\end{array}$$

which we solve numerically on a physical domain defined by $\mathit{\{}R,\mathit{\kappa }\mathit{\}}\in \mathbb{R}$ with R≥0. This leads to three roots at each point in the domain, which we can refer to as Ω₁(R,κ), Ω₂(R,κ) and Ω₃(R,κ). Two of the roots, say ${\mathrm{\Omega }}_{(\mathrm{2},\mathrm{3})}(R,\mathit{\kappa })$, always have a real negative value and hence correspond to exponentially decaying modes. One of the roots, Ω₁(R,κ), which is real in all the domain, has a positive real part in a region of the (R,κ) plane. The contour of this area is shown in Fig. 2 by the black line. This defines the Turing zone. Hence, for those values of (R,κ) with a positive real part of Ω₁ (inside of the contour) the corresponding homogeneous equilibrium is linearly unstable to inhomogeneous perturbations of wavenumber κ.

Figure 2 Contour plot for Ω₁(R,κ). The thick black line is the contour Ω₁=0, and the region inside that contour corresponds to the parameter region (R,κ) in which the homogeneous equilibria are linearly unstable against inhomogeneous perturbations; i.e. Ω>0.

Zero modes correspond to marginally non-growing inhomogeneous solutions of the linearized equations and correspond to perturbations with Ω=0. Their existence usually signals the presence of an instability that can send the system towards an inhomogeneous time-invariant solution in the full nonlinear model. Specifically, Fig. 2 suggests the existence of two zero modes for any value of rainfall comprised between critical bounds 1.0 and 1.25 mm d⁻¹. Such zero modes can satisfy the boundary conditions for a (low) wavenumber that depends on R and therefore require a domain large enough to develop. If the domain is so small that its fundamental mode is in the stable region (κ≳0.6), it is always stable. In this section, we found zero modes along the homogeneous branch compatible with a specific domain. Those zero modes indicate the start of inhomogeneous branches of equilibria. In the following section we will rely on those zero modes to compute equilibria of the full nonlinear model.

2.2 Bifurcation diagram and continuation method

The analysis in the previous section showed that homogeneous vegetation is linearly unstable against spatially periodic perturbations in the rain range $\mathrm{1}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}<R<\mathrm{1.25}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$. This linear analysis suggested the existence of periodic inhomogeneous equilibria in the nonlinear system within that rainfall range. Now, our focus shifts towards explicitly computing these nonlinear equilibria. To achieve this, numerical methods are employed to solve the equations involved.

Figure 3 Zero modes present in a domain of L=100 m. Horizontal lines correspond to harmonics that fit inside the periodic L=100 m domain. The intersection of the horizontal lines with the Ω₁=0 contour provides the values of rain R at which new branches of equilibria might appear. Note that each of the nine horizontal lines that intersect the contour intersect it twice; hence, there are 18 zero modes.

For the time being we assume a periodic domain with a finite size of L=100 m. The choice of the size of the domain influences the number of zero modes exhibited on the homogeneous equilibrium branch. This choice together with periodic boundary conditions effectively discretizes the set of zero modes (Fig. 3). Indeed, setting $\mathit{\kappa }=n\frac{\mathrm{2}\mathit{\pi }}{L}$, with n the wavenumber, we find two pairs (κ,R) corresponding to zero modes, with the range $\mathit{\kappa }=[\mathrm{0}:\mathrm{0.6}]$ corresponding to perturbing the homogeneous equilibrium $\overline{S}\left(R\right)$ with a perturbation δS(x)=δScos (κx).

For each pair of (κ,R) in the set of zero modes, we identify the corresponding homogeneous equilibrium ($\overline{S}$) and perturb it with the corresponding perturbation δScos (κx). The periodic boundary conditions are enforced by setting $\mathit{\kappa }=n\frac{\mathrm{2}\mathit{\pi }}{L}$, with n the wavenumber. The perturbed homogeneous solution

$$\begin{array}{}\text{(18)}& {\displaystyle }{B}_{\mathrm{0}}(x,R)& {\displaystyle }=\overline{B}+\mathit{\delta }B\cos \left(\mathit{\kappa }x\right),\text{(19)}& {\displaystyle }{W}_{\mathrm{0}}(x,R)& {\displaystyle }=\overline{W}+\mathit{\delta }W\cos \left(\mathit{\kappa }x\right),\text{(20)}& {\displaystyle }{O}_{\mathrm{0}}(x,R)& {\displaystyle }=\overline{O}+\mathit{\delta }O\cos \left(\mathit{\kappa }x\right),\end{array}$$

is then taken as the first-guess input in a Newton–Raphson iteration used for finding the corresponding nonlinear inhomogeneous equilibrium.

At this point we need to discretize the spatial domain into, say, N points. We used N=100. Specifically, the discretized first guess reads $u_{\mathrm{0}}\left(R\right)=\left[B_{\mathrm{0}}\right(R),W_{\mathrm{0}}(R),O_{\mathrm{0}}(R\left)\right]\in {\mathbb{R}}^{\mathrm{3}N}$ and the algorithm

$$\begin{array}{}\text{(21)}& u_{i+\mathrm{1}}\left(R\right)=u_i\left(R\right)-\mathit{ϵ}{J}^{-\mathrm{1}}\left(u_i\right(R\left)\right)f\left(u_i\right(R\left)\right),\end{array}$$

converges towards the equilibrium u_S(R) with J(u_i(R)) being the Jacobian matrix of the system and f(u_i(R)) being the right-hand side as in Eq. (1). The Laplacian is discretized with a periodic pseudospectral method. This equilibrium, u_S(R), is then used as the first guess to solve for the equilibrium with rainfall R+δR,

$$\begin{array}{}\text{(22)}& u_{\mathrm{0}}(R+\mathit{\delta }R)=u_S\left(R\right)+\mathit{\delta }u,\end{array}$$

with δu≪ϵ. This iterative procedure leads to the construction of the full branch of equilibria.

A new code has been developed to implement this continuation method. This code is available online in the repository mentioned in the “Code availability” section.

Each branch obtained by the above technique is denoted by n between 1 and 9, with n the wavenumber associated with the perturbation. The periodic boundary conditions, together with a finite domain size, fixed the number of zero modes present within the Turing zone and, therefore, the number of inhomogeneous branches that exist. For the periodic L=100 m domain, only nine zero modes fit within the Turing zone (see Fig. 3). Higher-wavenumber modes correspond to linearly stable perturbations that decline exponentially.

Figure 4 Bifurcation diagram for the Rietkerk's model with L=100 m. Each branch is labelled by an integer corresponding to the order of the wavenumber associated with the zero mode. Solid lines correspond to stable states and dashed lines to unstable states.

All the branches obtained with this approach are shown in Fig. 4. We now see that the range of rainfall where vegetation is sustained is much wider than expected from the linear analysis. The linear analysis in the previous section predicted the growth of vegetation patterns in the Turing zone: between 1.0 and 1.25 mm d⁻¹. The analysis of the full nonlinear system, presented in Fig. 4, actually reveals time-invariant pattern equilibria on branches attached to zero modes in the range 0.5–1.3 mm d⁻¹. Although the nonlinear equilibrium may differ in shape from the linear perturbation, it is found that wavenumber n used to perturb the homogeneous equilibrium describes reasonably well the shape of the pattern along the branch obtained from that perturbation. One such equilibrium is displayed for illustration in the top panels of Fig. 5. For this particular equilibrium, with wavenumber n=2, the biomass accumulates in two places, and soil water peaks there due to enhanced infiltration rate. Surface water accumulates around those areas.

Figure 5 On the upper panels, the equilibrium for a given R (rainfall) (0.9 mm d⁻¹) on the n=2 branch is shown. The black dot on the left-hand-side panel shows the position of the equilibrium on the bifurcation diagram. The three other panels show, from left to right, B (biomass) (g m⁻²), W (soil water) (mm) and O (surface water) (mm). In the lower left-hand-side panel, eigenvalues of the Jacobian matrix associated with the equilibrium are shown. The real part is on the x axis and the imaginary part on the y axis. The larger eigenvalue is represented by a red dot. The black vertical line marks the 0 position. The three other lower panels show, from left to right, the eigenvector associated with the larger eigenvalue projected on B (biomass), W (soil water) and O (surface water), and a dashed black line indicates the first spatial derivative of the equilibrium rescaled.

2.3 Stability analysis

Characterizing the stability of equilibria will further help us to understand the dynamics of the system inside the Busse balloon. The Busse balloon is the region in the parameter space (κ,R) where at least one stable pattern equilibrium exists. A linear stability analysis gives partial but valuable information about a given solution and its possible evolutions. As we will see in the following, while stable equilibria tend to be the end state of time evolutions, unstable equilibria may still be relevant for the dynamics. The stability of equilibria is classically estimated based on the Jacobian matrix evaluated at each equilibrium, positive real parts of the eigenvalues signalling instability.

As a first example, the equilibrium for rainfall equal to 0.9 mm d⁻¹ of the branch n=2, the associated eigenvalues and the first eigenvector are displayed in Fig. 5. The majority of eigenvalues have a negative real part. One eigenvalue has a very small, positive real part ($\mathrm{3}\times {\mathrm{10}}^{-\mathrm{10}}$), which suggests quasi-neutral stability. Further inspection of the associated eigenvector shows that this quasi-neutrally stable mode is the spatial derivative of the equilibrium pattern. Therefore, it corresponds to a translation mode, which is indeed expected with periodic boundary conditions. Hence, we have a stable pattern, which may however translate consistently with the periodic character of the boundary conditions. We therefore attribute the small real part to a numerical artefact (it should be zero) and consider equilibria with positive real part of the highest eigenvalue of the order of 10⁻¹⁰ as stable.

This stability analysis is repeated for all the equilibria shown in Fig. 4. Stable sections of the branches are drawn as solid lines and unstable ones as dashed lines. This gives us a more precise idea about what is happening inside the Busse balloon. We now know what the shape of those stable states is and how and at what rainfall value they lose stability. The existence of multiple stable states for a given value of rainfall is also consistent with previous work based on models and real systems (Bel et al., 2012; Bastiaansen et al., 2018).

Some branches of equilibria, such as n=5, 6, 7, 8, 9, are unstable across all their existence domain. Other branches, such as n=1, 2, 3, 4, change stability along the branch. These branches of equilibria start as unstable from the low-vegetation zero modes on the homogeneous equilibrium. Then, they become stable and eventually become unstable again shortly before joining the highly vegetated homogeneous equilibrium at the corresponding zero mode. The change in stability may take place at the extreme rain values for which the equilibrium exists, as in the n=1 branch, or at a rainfall value in the middle of a branch. As we know, when an equilibrium loses stability, the real part of one of its eigenvalues changes from negative to positive, which indicates that the equilibrium becomes linearly unstable in the direction of the associated eigenvector. In all of the cases here the change of sign of the eigenvalue happens through zero (as opposed to, through infinity), meaning that there are zero modes at the intersection of the stable and unstable branch sections. As we show next, this indicates the branching of another branch of the equilibrium of the full nonlinear system. A simple way to identify those zero modes is proposed in Appendix B.

2.4 Mixed states

In the previous section we studied the stability of the “main” branches of equilibria shown in Fig. 4. These branches of equilibria originated from the zero modes present on the homogeneous vegetated branch. We also saw that some of these branches of equilibria had zero modes at the intersection of their stable and unstable parts. These zero modes act as bifurcation points from which new branches of equilibria emerge. The latter can be found by starting from the equilibrium at the bifurcation point, perturbed in the direction of the eigenvector associated with the newly positive eigenvalue (zero mode) corresponding to a new unstable direction. Again, a Newton–Raphson iteration allows us to find the new equilibrium. From there, the continuation method explained in the previous section allows us to trace the full new branch.

Figure 6 Panel (a): bifurcation diagram with the addition of mixed-state branches. The mixed-state branches originating the zero modes along the n=4 branch are called n=4 bis₁, n=4 bis₂ and n=4 bis₃. The bifurcation point is marked with a black dot. Panel (b): an enhancement of the area around the bifurcation point. Panels (c)–(f): the biomass profile of the n=4 and the mixed states for a particular value of rainfall. The position on these equilibria on the bifurcation diagram is marked by a circle in panel (b).

The result of this routine applied to the bifurcation point of the branch n=4 is shown in Fig. 6a and b. At the transition between stability and instability, there are three zero modes. Hence three branches of equilibria start at this transition, and they are shown in Fig. 6. They all begin at the transition between the stable and unstable part of the n=4 branch, and they reconnect with it for higher rainfall. We call those branches n=4 bis₁, n=4 bis₂ and n=4 bis₃. These three branches are close (with regard to the mean biomass) to that of the n=4 branch, but if we look closely at their profile (Fig. 6c–f) we see how they differ. For a given value of rainfall, the mixed-state equilibria look like a modulation of the n=4 equilibrium by another wavenumber.

Equilibria branching out of zero modes in the homogeneous equilibrium tend to exhibit a single perturbation mode; see Fig. 5. By contrast, equilibria branching out of zero modes in those inhomogeneous equilibria exhibit a mixture of modes (Fig. 6c–f), hence the name “mixed states”. The equilibria along the mixed-state branch n=4 bis₃ are unstable, with positive eigenvalues. The branches n=4 bis₂ and n=4 bis₃ reconnect around 1.14 mm d⁻¹. There are numerous mixed-state branches, obtained from the different zero modes at bifurcation points found along the main branches. These can also be found to emerge from zero modes present on the unstable sections of the main branches. For example, in the low vegetated unstable part of the n=3 branch (see Fig. A2), we find two zero modes at the same value of rainfall, from which two mixed-state branches emerge: one labelled n=2 bis because it connects to the n=2 branch and one labelled n=3 loc, which connects to the n=4 branch. The n=3 loc refers to the fact that the equilibrium of that branch is localized in space. The n=1 main branch is special in two ways. First, as Zelnik et al. (2013) stated, it is the only stable localized state. Second, it starts from a zero mode on the homogeneous equilibrium, but it ends up connecting to the n=2 branch as if it was a mixed state; see Fig. A1. The evolution of the shape along this branch goes from one localized vegetation peak to two localized vegetation peaks in the unstable part close to the connection with the n=2 branch.

All mixed-state branches share similar characteristics, with states appearing as a modulation of a main branch. Even though they are unstable, we expect mixed-state equilibria to influence the dynamics of nearby trajectories, depending on the value of the positive eigenvalue.

3 Numerical simulations of trajectories

To assess the relevance of the bifurcation diagram for understanding transient dynamics, we performed two series of numerical simulations.

Figure 7 Trajectory of the Rietkerk model under a decreasing rainfall scenario. Panel (a): the trajectory in the bifurcation diagram plotted as the blue dashed–dot line. Panel (b): the rainfall with respect to time; panel (c): the time evolution of the biomass.

3.1 Trajectories in a changing environment

In order to assess the relevance of the various branches for the dynamics, we propose the following numerical experiment. The Rietkerk model is run with rainfall changing over time at a rate of $\frac{\mathrm{d}R}{\mathrm{d}t}=-\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}\text{mm}{\text{d}}^{-\mathrm{2}}$. The starting point is the vegetated homogeneous equilibrium with rainfall $R=\mathrm{1.4}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$. The final rainfall is 0.5 mm d⁻¹. The result is shown in Fig. 7. Until $R\sim \mathrm{1.2}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, the system stays in the homogeneous equilibrium. At this point, the system undergoes a Turing bifurcation and jumps to a heterogeneous equilibrium. The patterned branch on which the system lands is the n=2 branch. We can understand this transition by looking at Fig. 3. Indeed we see there that the first mode to lose stability in a decreasing rainfall scenario is the n=2. So this simple linear analysis can give us a first idea about the way the system will be destabilized in a slow-change scenario. After that, the system tracks the n=2 branch until this equilibrium loses its stability. Finally, the system switches to the only other stable equilibrium, the n=1. For this experiment, the mixed states do not have an influence on the trajectory. But in the following we will see that even if the mixed states are unstable, they can still play a role in the dynamics.

3.2 Effect of a mixed state on trajectories

Here we address the question of whether unstable mixed-state equilibria are also able to influence system-transient trajectories. In this case the rainfall is set to rainfall $R=\mathrm{1.05}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, for which a mixed-state equilibrium exists.

Figure 8 Panel (a): trajectories in a summary phase space from an initial condition close to the n=5 equilibrium with a fixed rainfall $R=\mathrm{1.05}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$. The relevant equilibria are shown; the stable equilibria are represented by a circle and the unstable equilibria by a cross. In panel (b), we show the time evolution of the biomass. In panels (c)–(e), the dynamical solution at three different times is represented (in green) along with the shape (in black) of the mixed state denoted by n=4 bis₃.

Figure 9 Bifurcation diagram for the Rietkerk's model with L=200 m. Each branch is labelled by an integer corresponding to the order of the wavenumber associated with the zero mode. Solid lines correspond to stable states and dashed lines to unstable states.

The initial condition is an unstable equilibrium (n=5) with a small perturbation along the direction of the first eigenvector. Figure 8a and b summarize the evolution of the biomass pattern over time and the corresponding trajectory in a summary phase space. This summary phase space consists of a two-dimensional phase space where the two dimensions are the mean biomass and the maximum biomass. With that in mind, we observe that the system leaves the n=5 unstable equilibrium by first reorganizing itself: mean biomass remains constant while maximum biomass increases. This means that one or more vegetation bumps are growing at the expense of others. After that initial reorganization, an excursion on higher mean and maximum biomass takes place. It is during this out-of-equilibrium excursion that one of the vegetation bumps is lost. At this point in time the vegetation profile, containing four vegetation bumps, passes close to the three mixed states originating from zero modes on the n=4 branch explained identified in the previous section. Then, the system undergoes another transition to the n=3 equilibrium. To have a better representation of this trajectory, we also show the time evolution in Fig. 8b. We see that the system quite slowly leaves the n=5 unstable equilibria. The state is indeed unstable, but the associated positive eigenvalues are small, such that the dynamics around the equilibrium are slow. After abruptly departing from the n=5 equilibrium, one biomass bump disappears, leading to a rearranged state with four bumps.

We see that even though the summary phase space may let us think that the system passes close to the mixed-state n=4 bis₂, we have to keep in mind that this two-dimensional space is obtained by the projection of the infinite-dimensional phase space on a two-dimension summary space. The biomass profiles at three points in time are represented in Fig. 8c–e, and their corresponding positions are also shown in panel a. In each of panels c–f, the shape of the n=4 bis₃ mixed state is also shown. The shape of the dynamical solution is actually close to the n=4 bis₃. That equilibrium is unstable, but the system lingers in its vicinity for a considerable period of time before another bump of biomass vanishes, propelling the system towards the stable equilibrium of n=3.

Switching from one state to another by losing one or more vegetation bumps is a known feature of the vegetation pattern model (Bastiaansen and Doelman, 2019; Bastiaansen et al., 2020). This mechanism is called sideband instability (Doelman et al., 2012; Siteur et al., 2014). We used the same setup with different rainfall values along the n=5 branch. The fact that the system passes close to a mixed state is consistent along the n=5 branch only for values of rainfall for which mixed-state equilibria exist. For rainfall lower than $\sim \mathrm{0.91}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, the system jumps directly to the n=3 equilibrium. Now, if we consider the same rainfall but we start from another unstable branch like n=6 or n=7, the dynamics are different. For n=6, we observe a different type of destabilization called period doubling (Doelman et al., 2012; Siteur et al., 2014). With that mechanism the system transitions to the n=3 equilibrium without passing through a mixed state, and the system spends less time around the unstable equilibrium of n=6 compared to the n=5 case. This might be explained by the fact that the positive eigenvalue on the n=6 solution is larger than the one on the n=5 solution. For n=7, the transition is even more rapid, as expected for an even larger positive eigenvalue, but the landing point is still n=3. In that case, we observe a destruction of four bumps similar to a sideband instability.

4 Discussion

We uncovered the existence of mixed states in Rietkerk's model and showed their importance for transient dynamics. As far as we are aware, such mixed states have never been identified nor described in models for vegetation patterns.

Mixed states emerge at the transition between the unstable and stable states along a branch of equilibria and have a functional form that appears as the combination of two equilibria from the main branches. We found that while these equilibria are unstable, they may still influence the system's dynamics by slowing down its evolution when it passes near them.

The influence of unstable modes on dynamics has been studied in ecological models. For example Sherratt et al. (2009) showed how spatio-temporal chaos appears in the wake of an unstable wave train for the complex Ginzburg–Landeau equation. In the more general context of ecology, Hastings et al. (2018) and more recently Morozov et al. (2020) proposed a classification of transient phenomena in ecological models based on a dynamical system approach. In their definition, a dynamical regime is transient if it persists for a sufficiently long time (quasi-stable) and if the transition between two regimes occurs on a much shorter timescale than the time of existence of the quasi-stable regime. According to that definition, the behaviour observed in Sect. 3.2 can be seen as a transient phenomenon. Indeed, the system spends hundred to thousands of days around unstable states (n=5, n=4 bis₁ and n=4 bis₃) and then takes only a few days to jump to the final stable equilibrium. By following the typology presented in those two papers, the transient observed here can be qualified as a “crawl by” transient. Indeed, the transient observed is linked to the existence of a saddle-type invariant set: the mixed state. Transient dynamics have also been studied in the context of coexistence in vegetation patterns by Eigentler and Sherratt (2019) and are characterized by the small size of a positive eigenvalue. Long transients can also be observed in very slow front invasion dynamics (Van De Leemput et al., 2015) The mechanism behind the transient is different but the effect is the same: the system spends a long time in a region of the phase space without an attractor.

As we see in Sect. 3.2, the history and the initial conditions of the system are important to see whether the mixed state will or will not appear in the dynamics Other vegetation pattern models exhibit this type of sensitivity to history and initial conditions (Sherratt, 2013; Adams et al., 2003; Alberti et al., 2015). Further research is needed to determine whether mixed states are a common feature of reaction–diffusion systems.

We adopted periodic boundary conditions, as in previous studies (Rietkerk et al., 2002; Dekker et al., 2007). As stated by Dijkstra (2011), the periodic boundary conditions allow for the existence of unstable equilibria with a high wavenumber. However, this boundary condition, widely used in reaction–diffusion models due to its simplicity and convenience, may not reflect real-world scenarios accurately. More generally, working with a periodic domain discretizes the set of admissible equilibria. Hence, the significance of mixed states for describing the dynamics of non-periodic infinite-size systems needs be further assessed.

As a step towards this objective, we computed the bifurcation diagram for a larger domain size of L=200 m instead of L=100 m. The result is shown in Fig. 9. As expected, we observed more branches in this case, but, remarkably, the branches and the stability of the branches shared between the two domains are the same as for L=100 m. Of course the even n number branches for the L=200 m domain are the same as for the L=100 m domain. The new branches are all designated with an odd n number as a consequence of the extended domain. The range of rainfall that can support a stable vegetation pattern also increased slightly with the larger domain size, with the n=1 branch extending this range.

These findings suggest that the existence of mixed states and their stability properties is robust across domain sizes, even though how they would manifest themselves in continuous large domains is still to be established.

5 Conclusions

We present an in-depth analysis of the Rietkerk model's behaviour that we believe sheds a new light on the dynamics of vegetation patterns. Our analysis is based on the bifurcation diagram of the vegetation pattern model, which we constructed using a variety of techniques. First, we performed a linear stability analysis of the homogeneous equilibrium of the system, which allowed us to delineate the so-called Turing zone. This zone is characterized by the instability of the homogeneous equilibrium to small heterogeneous perturbations, and it serves as a starting point for our analysis. Next, we used perturbed states at the edge of the Turing zone, the so-called zero modes, to construct whole branches of equilibria associated with the system. To obtain these branches, we used a continuation method based on the Newton–Raphson algorithm. By doing so, we were able to obtain the heterogeneous equilibria for a given size of the domain. As expected, due to the system's nonlinearity, we discovered that the system allows for multiple equilibria for a given value of a parameter, specifically rainfall. We went a step further by assessing the stability of those equilibria by computing the eigenvalues and associated eigenvectors of the Jacobian matrix. This stability analysis showed that inhomogeneous vegetated states can exist and be stable, even for values of rainfall for which no vegetation exists in the non-spatialized system. Yet, we also found that the main branches of equilibrium originating from a zero mode are not the only ones present in the system. At the transition between unstable and stable states along one branch, a new type of equilibrium appears. These latter types of equilibria, which we called mixed states, are always unstable equilibria, and they appear as a mix of two equilibria from the main branches. Although these equilibria are unstable, they affect the dynamics by slowing down the evolution of the system when it passes close to it. This slowing effect has been described in ecology, but to the best of our knowledge, the mixed-state influence on the dynamics had not been previously shown for vegetation patterns. Overall, our approach allowed us to construct a bifurcation diagram that gives us valuable insights into the behaviour of the system. This approach is helpful to disentangle the fate of the system in the Busse balloon and could be used to assess the existence or non-existence of mixed states in spatially extended systems.

Appendix A: Random initial conditions with a fixed rainfall

Here, our objective is to verify that the stable pattern branches act as attractors in the dynamics of the system. To this end we specify a value of rainfall for which multiple stable equilibria exist. Starting with different initial conditions and evolving the system, we expect that different runs end in some of the different available stable states for the chosen rain value. To avoid favouring the attraction to a particular branch, we opt for bare-soil initial conditions with random noise. Rainfall is set to $R=\mathrm{1.1}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, for which the system has several stable equilibria (n=2, 3 and 4). The model is run for 5000 d with 500 distinct random initial conditions. These are random from a uniform distribution (positive) noise on top of a homogeneous bare-soil equilibrium. Trajectories, projected onto a summary space (mean and maximum biomass) are shown in Fig. A4. As this is a summary space – a projection of the full space – trajectories may appear to cross, but they do not in the full space. Out of the 500 runs, 465 land on the n=2 equilibrium, and the others land on the n=3 equilibrium. This indicates that the equilibria that we identify as stable are indeed those which attract trajectories. Having no trajectories on the n=4 equilibrium suggests that the basin of attraction of the latter is narrowly close to the homogeneous bare-soil equilibrium, with few chances for a trajectory starting from random initial conditions to be in the basin of attraction of that equilibrium or, in more informal terms, to evolve such as to pass close enough to the n=4 equilibrium and land on it. The n=4 equilibrium can however be reached with carefully designed initial conditions, such as a cosine with a wavenumber equal to $n\frac{\mathrm{2}\mathit{\pi }}{L}$.

Figure A1 A focus in the bifurcation diagram on the n=1 and n=2 branches; the stable (unstable) states are noted with a continuous (dashed) line. Upper right-hand-side panel: a zoom-in on the area of the bifurcation diagram where the n=1 branch connects with the n=2 branch. This connection is represented by a star. In panels (a–e), equilibria along the n=1 branch are shown. Each corresponds to a black dot in the upper left-hand-side panel.

Figure A2 A focus on the bifurcation diagram of the n=2, n=2 bis, n=3, n=3 loc and n=4 branches; the stable (unstable) states are noted with a continuous (dashed) line. In the upper right-hand-side panel we show a zoom-in of the bottom part of the bifurcation diagram. The n=3 branch exhibits a double-zero mode, indicated by the green star. Two equilibria emerge from that double-zero mode, the n=2 bis and the n=3 loc. The n=3 loc branch is characterized by a three-bump solution, which is localized. The first equilibrium connects with the n=2 branch, while the second one connects with the n=4 branch at the zero mode indicated by the orange star. In panels (a–e), equilibria along the n=3 loc branch are shown. Each corresponds to a black dot in the upper left-hand-side panel.

Figure A3 Bifurcation diagram for which the stability is computed based on the four highest eigenvalues. The solid line means that the eigenvalue considered is negative. The dashed line means that the eigenvalue is positive.

Figure A4 Trajectories on a summary phase space from random initial conditions with a fixed rainfall $R=\mathrm{1.1}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$. The relevant equilibria are shown. The stable equilibria are represented by a circle and the unstable equilibria by a cross.

Appendix B: Identification of zero modes along branches

In Sect. 2.4, we showed that new type of equilibrium (mixed state) can appear at the transition between unstable and stable part of a branch. This transition is linked to the transition from positive to negative of the eigenvalue. For example at the transition on the n=4 branch, we have three zero modes leading to the appearance of three mixed-state branches. In Fig. A3, we propose a way to identify and represent those zero modes. Each of the panels is linked to one of the four highest eigenvalues. If this eigenvalue is positive the line is dashed, and if the eigenvalue is negative, the line is solid. With that in mind, let us focus on the n=4 branch. We see that around $R=\mathrm{0.91}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$, the three first eigenvalues become negative. This means that there are three zero modes. We also observe that for the n=3 branch, two eigenvalues switch from positive to negative around $R=\mathrm{0.73}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$. Again, this means that two mixed-state branches emanate from the zero mode at the transition.