Implications of tristability in pattern-forming ecosystems

Ecosystem functioning is often impaired by transitions from one stable system state to an alternative stable state induced by changes in environmental conditions or by disturbances.¹ Examples of such transitions (also called “regime shifts”) are loss of transparency and vegetation in shallow lakes subject to human-induced eutrophication,² climate-induced bleaching in coral reefs,³ and loss of biological productivity in drylands, or desertification,⁴ as a result of natural and human-related causes, such as droughts, deforestation, over-grazing, and others. Studying the dynamics of desertification is particularly challenging as one of the alternative stable states is often a periodic vegetation pattern.^5,6

The significance of bistability between uniform vegetation and a periodic vegetation pattern has been studied recently in the context of Namibian grasslands showing the spectacular phenomenon of “fairy circles.”^7–9 Fairy circles are circular bare-soil gaps (6 m in diameter on average) that self-organize to form nearly periodic hexagonal patterns on large, landscape scales. Model studies of this system reveal a range of bistability between these patterns and a uniform vegetation state with a subrange within which multiple hybrid states exist.¹⁰ These are stationary localized states consisting of fixed periodic-pattern domains embedded in the uniform vegetation state and various combinations thereof. This type of behavior has been found and studied in many other contexts and is commonly referred to as homoclinic snaking.^11–14 It is related to the capability of fronts separating uniform and patterned states to be pinned in place over a range of the control parameter (called the “snaking range”),¹⁵ in contrast to fronts separating two uniform states, which in general propagate and only become stationary at isolated values of the control parameter—the Maxwell point(s). While bistability of two vegetation states of different productivity may naively imply abrupt desertification transitions associated with crossings of tipping points, the spatial extent of the system and the existence of multiple hybrid states offer many additional response forms, including gradual state transitions and incomplete transitions.^16,17

Quite often, along with stable uniform and patterned vegetation states, the bare-soil state remains stable as well, giving rise to tristability ranges. Such situations can result from local feedbacks that stabilize the bare-soil state even at relatively high precipitation rates. An example of such a situation is fast soil-water evaporation in bare soil relative to water uptake by seedlings, and reduced evaporation by shading in vegetated soil. Another example is fast overland-water evaporation in bare soil relative to water infiltration into the soil. The tristability range of bare soil (BS), uniform vegetation (UV), and periodic patterns (PP) includes three bistability ranges, BS-UV, UV-PP, and PP-BS, where various forms of localized structures, including fronts, pulses, and hybrid states, may exist. To what extent do these bistability structures persist in tristable systems, where localized structures involving the third state exist, and what additional localized structures can be expected?

In this paper, we address a few aspects of these questions using a simple model for dryland vegetation. We study (i) the manner by which the snaking range associated with pinned UV-PP fronts is affected by the presence of UV-BS fronts, (ii) possible transitions between different front types, and (iii) new front solutions with more complex inner structures. We further discuss the significance of our results for regime shifts in dryland ecosystems subjected to rainfall variability.

Several continuum models have been used to investigate spatial patterns in dryland ecosystems.^18,19 Most of these models assume the form of partial differential equations (PDEs) for the areal densities of biomass and water variables that capture pattern-forming biomass-water feedbacks.^20,21 We focus here on one of the simpler models, the Klausmeier–Gray–Scott model,^22–24 which for flat terrain assumes the nondimensional form

where B is the vegetation biomass and W the water content. In these equations, vegetation growth in the presence of water is modeled by the nonlinear term WB² in both equations, vegetation dies at rate m, water is lost by evaporation at a rate normalized to unity, and precipitation increases the concentration of water uniformly across space at a constant rate a. Both the biomass and water diffuse laterally in space with the water diffusing d times faster. In the following, we assume that m = 0.5, focusing on the system in one spatial dimension with either Neumann or periodic boundary conditions, and explore its behavior when the remaining two parameters a and d vary.  Appendix A provides details on the numerical methods we use. We note that more realistic descriptions of water dynamics require three-variable models that distinguish between overland water and soil water.^25,26

The model (1) has two time-independent spatially uniform solutions, corresponding to bare soil (B = 0) and uniform vegetation (B > 0). The bare soil state (BS) exists and is stable for all parameter values, while the uniform vegetation state exists for $a>2m$ only. The uniform vegetation state is unstable to spatially nonuniform perturbations when the diffusion ratio d is high and the precipitation a is low (but still larger than 2m). This is a Turing instability (see  Appendix B) that results in stationary periodic patterns of different wavelengths. Previous studies of this model^23,27 have focused on the case of a supercritical instability, present for sufficiently high d, and it has been assumed that the Turing instability is always supercritical in this model.²⁸ However, as shown numerically in the context of both chemical reactions^29,30 and dryland vegetation,³¹ the instability may in fact be subcritical. The subcriticality leads in turn to bistability between the uniform vegetation state (UV) and a periodic state (PP), and hence to tristability between the BS, UV, and PP states. The extent of tristability can be determined analytically as shown in  Appendix C. Tristability of this kind has also been found in more realistic models of dryland vegetation.^10,31,32

We begin by addressing three types of localized structures that arise from pairwise bistability between the UV, PP, and BS states, and play essential roles in the dynamical behaviors to be described. They are summarized in Fig. 1 which shows partial bifurcation diagrams in terms of the L² norm of the vegetation cover, $||B||=L−1∫0LB2(x) dx$⁠, as a function of the precipitation rate a. These show the UV and BS states in black, the PP state in blue, and the localized structures in red. Sample solution profiles B(x) and W(x) of the biomass (green curves) and water (blue curves) distribution corresponding to the locations indicated by open circles are shown in the bottom row. We now discuss each case in turn.

Within the UV-BS bistability range (Fig. 1, left column), there are simple fronts connecting these two uniform states. Such fronts are stable solutions and move with a constant velocity that depends on the model parameters, with UV invading BS for values of a exceeding the Maxwell point and vice versa for values below that point. Figure 1 (left column) shows the stationary front at the Maxwell point. When periodic patterns are bistable with the bare-soil state, a pulse-like solution exists, describing a single patch of vegetation in otherwise bare soil (Fig. 1, right column). This solution represents the last vegetation form that survives at low values of a. In contrast, bistability between PP and UV leads to a different type of localized state, as shown in the middle column of Fig. 1. In this bistability range, we find homoclinic snaking and a multitude of hybrid states, describing localized groups of nearly bare-soil patches in otherwise uniform vegetation. Localized structures of this kind may affect the nature of the transition between the UV and the PP states.^10,13

To acquire a better understanding of the bifurcation structure and its consequences for the system dynamics, we now examine the behavior of the system over a large part of the (a, d) parameter plane. Figure 2 summarizes the regions of existence and stability of different states with the boundaries between them shown in different colors.

The Maxwell point between BS and UV is indicated by the magenta line so that for lower values of a or d the BS state invades the UV state, and vice versa for larger a or d. The magenta line, together with the two red lines, forms a triangle-like shape where hybrids of the UV and PP states are stable, i.e., the snaking region. Thus stable hybrid states are never present in the region where the BS state invades the UV state. The existence range of the UV state is bounded on the left by the vertical black line but this state is unstable between this line and the leftmost blue line denoting the location of a Turing bifurcation. The blue line to the right marks the location of a saddle-node bifurcation where the PP solution branch that extends to the highest precipitation value disappears. The stability range of the single-patch state is marked by the two orange lines, so that nonuniform solutions are confined to the region between the leftmost orange line and the rightmost blue line. Finally, traveling front solutions connecting BS and UV are only stable in part of the region of bistability between BS and UV: the two dashed gray lines, together with the leftmost red line that connects them, mark the stability boundary for such front states; for larger d or smaller a such fronts break up into a periodic pattern.

We now turn to look at the implications of the bifurcation structure described in Sec. III for the dynamics of the system. We choose several points in the parameter space, labeled by letters from A to F, along the vertical cut a = 1.01 through the (a, d) plane (Fig. 2, right panel), and by G to L, along the horizontal cut d = 8; note that D and J mark the same point. We perform simulations at each point starting with similar initial conditions and report the results in Fig. 3. Each simulation starts with a spatial mixture of UV, PP, and BS, each of which occupies one third of the domain, UV in the left third, PP in the middle, and BS in the right third. The UV and BS states are initialized with values corresponding to the uniform states at the corresponding parameter values, with the periodic state in the middle consisting of 10 periods of a sinusoidal oscillation oscillating between the two uniform states on either side. We use periodic boundary conditions, so that each of the three regions interacts with the other two, and run the simulation for a fixed length of time, t = 4000. The top row shows the results for A to F, i.e., for increasing values of the diffusion ratio d, while the bottom row shows results for G to L, i.e., for decreasing precipitation a.

Starting with the vertical cut, we look at point A. For a diffusion rate this low (d = 4.5), there are no stable periodic patterns and the periodic state therefore collapses to the BS state immediately. Since we are to the left of the Maxwell point, the BS state invades UV, and we end up with the bare soil state. For higher diffusion rates, as shown by point B, periodic patterns are stable. This means that the PP state no longer collapses, but instead slowly expands, while the UV state eventually becomes a stable single-patch state. Because no patch-splitting²⁹ occurs for these parameters, the peaks in the periodic state deepen and gradually drift apart until they are evenly spaced and form a long-wavelength periodic pattern that fills the whole domain (not shown). At point C (d = 7.2), we have crossed the Maxwell line, and now the UV state invades BS. Since the single-patch state is still stable, no patch-splitting occurs, and the expansion of the pattern into BS is much slower than that of UV into BS. Thus, we eventually end up with a localized patch of the periodic state embedded in a UV background (not shown). At point D, we still have localized solutions, but no stable single-patch state. Here the BS region is invaded from both sides, by PP from the left and UV from the right. The invasion of BS by UV from the right occurs via a front that travels at a constant speed since both states are uniform in space. In contrast, the invasion front whereby the PP state invades BS occurs through a process known as front depinning.³³ Fronts of this type do not travel at constant speed—their speed oscillates around a nonzero mean. Once the BS region is eliminated, the PP and UV states begin to interact, and the competition between them ultimately results in front pinning, leading to the formation of a localized hybrid state, but one with more gaps than in the initial condition. Note that no peaks or gaps are added along the UV/PP boundary, i.e., the expansion of the PP state is highly asymmetric. At point E, outside the snaking range, the PP state invades both UV and BS, i.e., both fronts bounding the PP state depin, but the UV state remains linearly stable. Of particular interest is the fact that the front between the two uniform states at the right of the figure now nucleates a periodic state so that the PP state invades BS from both sides. This occurs since there is no longer a stable UV-BS front solution. We thus end up with a periodic pattern in the whole domain. Finally, for high diffusion rates such as that at point F (d = 10.2), the UV state is no longer stable and quickly collapses to a PP state, whereas the invasion of BS by PP remains a gradual process. Thus increasing the diffusion ratio d gradually favors biomass growth, as water is able to redistribute itself throughout the domain more rapidly.

Along the horizontal cut, we begin at the point G, with a high value of a (a = 1.025). This point is to the right of the Maxwell point, so that the UV state invades BS, and again there are no stable periodic patterns (in fact, they do not exist in this range). Hence, the UV state takes over the PP region rapidly, while invading BS more slowly. For lower values of a, such as the point H (a = 1.021), the PP pattern is stable but the UV state still takes over. We thus have two separate progressive transitions, whereby the UV state takes over the regions on either side. For yet lower a, such as point I (a = 1.015), we enter the snaking range where localized states are stable. As a result, the initial periodic region remains essentially unchanged, except for a small change in wavelength, while the BS region shrinks in size, finally becoming a single localized gap or hole. Note that both here and at point H the front between PP and BS serves as a nucleus for generating a growing inclusion consisting of a UV state so that BS is in fact invaded by UV from both sides. This is not the case at point J (same point as D), where the PP domain expands faster and the BS domain is invaded by PP from the left and UV from the right. Thus the width of the localized state that remains at the end depends on the relative speed of the UV-BS front and the PP-BS front. Moving to lower values of a, we exit the snaking range. At point K (a = 1.005), the UV state is linearly stable but is invaded by the PP state. However, point K also lies within the range where single-patch solutions are stable and so no patch-splitting occurs. As a result the PP domain expands slowly into BS while maintaining the number of peaks, until they are evenly spaced and form a long-wavelength periodic pattern that fills the whole domain (not shown). Finally, at point L (a = 1.002) the UV state is no longer linearly stable and it collapses into a PP state which then slowly invades the BS state. Thus decreasing the precipitation rate a favors the PP state over the UV state.

Despite some significant differences in the dynamics along the vertical and horizontal cuts, some interesting parallels can be drawn. For high values of d or low values of a (points E, F, K, L), the PP state quickly takes over the whole system, while for low values of d or high values of a (points A, B, G, H), the system typically ends up in a spatially uniform state (except for point B, which is also taken over by PP, albeit after a very long time). In the middle (points C, I, D = J), there is a hybrid situation in which localized states survive but with different dynamics taking them there. It is worth noting that in many situations the third of the available states emerges from a front between the two other states, such as a pattern nucleated from a front between the two uniform states (points E, F, K, L) or a UV state nucleated from a front between PP and BS states (points H and I).

One main consequence of tristability, as can be gleaned from Fig. 2, is that the snaking range is limited by the Maxwell point of fronts between UV and BS. Indeed, as can be seen by comparing panels in Fig. 4, the snaking range shrinks when we decrease d. More precisely, while for d = 8 the Maxwell point and the snaking range are quite far apart (Fig. 4, right panel), for lower values of d the Maxwell point moves to higher values of a, so that around d = 7.2 (middle panel) it meets the snaking range. Once this occurs, the snaking range shrinks (Fig. 4, left panel) as the Maxwell point continues to move to higher values of a, a phenomenon known from earlier studies of the Swift-Hohenberg equation.³³ 

We can understand what happens when the Maxwell point and snaking range collide by looking at the different states before and after the collision, as shown in Fig. 5. For d = 8 (Fig. 5, right panels), we focus on two distinct families of solutions, localized state (hybrid state) made up of several gaps in a UV background, and single-patch solutions in a BS background. While these two solution families do not meet when d = 8, they do so when d = 7 (Fig. 5, left panels). As the diffusion ratio decreases from d = 8, the regions of stability of these two solutions types of solutions move closer together and eventually overlap. This overlap results in the merger of the localized and single-patch states in the vicinity of the Maxwell point between the UV and BS states.

The merging of the single-patch state with the first rung of the snaking branch occurs near d = 7.2 (middle panel of Fig. 4). The periodic patterns begin to develop increasingly larger gaps of bare soil as one moves down the snaking branch on the left side of the snaking region,³¹ as shown in the middle panels of Fig. 5. We also speculate that a complex zoo of localized states exists in this region where bare soil gaps may open anywhere within any localized state in the snaking structure, but this remains to be verified.

The vertical portion of the branch that begins to form near d = 7.2 and is fully connected by d = 7 is located at the Maxwell point at which a stationary front between BS and UV is present. If we start from a single gap in a UV state and move down the bifurcation diagram, this gap opens up until it transforms the solution profile into two effectively separated and opposing fronts between BS and UV, as can be seen in the middle of the left panels in Fig. 5. To understand the formation of this state, it is important to recall that it may be translated by $L/2$⁠, where L is the domain length. As one continues down the solution branch, the BS region broadens, until a single patch of biomass remains. A similar process occurs when one starts from a UV state with a different number of gaps. In this case, each gap becomes two separate fronts facing in opposite directions, and as one proceeds down the branch each gap grows to become a BS domain, until finally the system is made up of several isolated pulses of vegetation. In effect, we may say that the bifurcation structure we see in the left panel of Fig. 4 is a consequence of the collision of the Maxwell point between the UV and BS states and the localized-state snaking range.

The collision between the Maxwell point of UV-BS fronts and the snaking range of localized states is significant for state transitions from the biologically productive UV state to the less productive PP state induced by droughts. The presence of the snaking range and the associated multitude of stable hybrid states¹⁰ provide more response forms to droughts or to local disturbances (grazing, clearcutting, fires). Such natural and human drivers of change may induce partial or incomplete transitions, culminating in hybrid states, rather than global desertification transitions to the less productive PP state in the absence of a snaking range. The vanishing snaking range as d decreases may therefore result in global rather than local desertification. These outcomes may be relevant to grasslands showing uniform-vegetation states as well as nearly hexagonal gap patterns, such as the fairy circles of Namibia, or the recently observed fairy circles in north-western Australia.³² The Namibian grassland ecosystem appears particularly relevant because of the strongly nonlinear soil-water diffusion (hydraulic conductivity) associated with its sandy soil. The sharp drop in the hydraulic conductivity with decreasing rainfall may result in the elimination of the snaking range even with very moderate rainfall decrease.

Associated with the tristability precipitation range are three types of fronts: fronts separating UV and BS domains, UV and PP domains, and BS and PP domains [see Fig. 3(d), for example]. We refer to these fronts as “simple” because of their simple spatial structure. In each case, these fronts connect a pair of stable states, and thus represent so-called pushed fronts, i.e., fronts whose speed is determined by nonlinear processes. Outside the tristability range another class of simple fronts exists, describing stable states invading unstable states,³⁴ e.g., a stable PP state invading an unstable UV state [Figs. 3(f) and 3(l)]. Fronts of this type are called pulled because their speed is generally determined by the linear instability of the unstable state ahead of the front.³⁴ 

Simple fronts connecting two stable states can still give rise to complex dynamics as a result of longitudinal and transverse front instabilities.^11,35–37 While such instabilities may be found in models of dryland vegetation and may have significant ecological implications, we focus here on another type of instability associated directly with the presence of tristability. As shown in Fig. 6, fronts between UV and BS are not stable throughout the bistability range of UV and BS, and when they are unstable the instability nucleates a periodic structure that invades both the region ahead of the front and that behind it [Figs. 3(e) and 3(k)]. Figure 6 shows solution branches that describe a UV domain invading a BS domain (solid black line) and a PP domain invading a BS domain (green line). The former front disappears in a saddle-node front bifurcation as a is decreased or d increased. The disappearance of the UV-BS front results in a front transition to the remaining stable front, PP invading BS, as indicated by the red arrows, explaining the behaviors found in Figs. 3(e) and 3(k).

Besides the effect of tristability on simple-front dynamics, an additional repercussion is the existence of multiple different complex fronts in the parameter range typified by point C in Fig. 3. At this location, the UV state invades BS but does so slowly. Figure 7 shows that multiple different types of fronts between the UV and BS states are present at this location, characterized by different numbers of oscillations or peaks in between the two homogeneous states. Each of these different fronts moves in a constant speed, which is different for each type of front, so that the more peaks there are, the slower is the speed with which the UV state invades the BS region. Figure 8 shows that the numerically computed inverse front velocity $vf−1$ associated with the front motion scales linearly with the number n of peaks according to

A careful look at the profiles in Fig. 7 reveals that the spacing between the spatial oscillations is not uniform across the front, as expected of localized structures with broken left-right symmetry:³⁸ the oscillations are more spread out near the leading edge than at the trailing edge. Figure 9 shows that the spacing between the peaks decreases linearly as a function of log(a–a_M), where a_M is the value of a at the Maxwell point. This behavior may be rationalized by assuming that the front extrapolates between a periodic structure near the back of the front and a larger wavelength structure near the leading edge resembling isolated peaks separated by larger distances. This differentiation becomes more pronounced for fronts with more oscillations in them, which also affects the parameter range over which they exist.

As shown in Fig. 10, all the fronts terminate at the Maxwell point for low values of a, while their limit for high values of a depends on the number n of oscillations within the front. The Maxwell point is a limit for these fronts since beyond it the BS state invades UV: beyond this point, the gap between the UV state and the peak just in front of it expands in time, filling with the BS state, and no fronts of constant form are found. For higher values of a, the range of the complex front solutions is largely determined by a nontrivial relation between the behavior of stationary single-patch solutions and the speed of the complex fronts. Once past the range where single-patch solutions are stable, such solutions tend to break apart and split into more patches. This happens to the complex front as well whenever the more widely separated peaks at its leading edge can acquire more water by moving away from the nearest neighbors. Thus, the parameter range within which complex fronts are found is reduced as the number n of peaks grows (see Fig. 9). This reduction in range is a consequence of an instability of the complex front that is inherited from the instability of a stationary single-patch state that it resembles at the leading edge (shown in a dashed black vertical line in Fig. 10). When a system with a complex front is taken outside its stability range, the result is a simple front taking over the entire domain instead of the complex front, as shown in Fig. 11.

The homoclinic-snaking behavior in the bistability range of uniform vegetation (UV) and periodic patterns (PP) has interesting implications for the ecosystem response to rainfall variability, as the multitude of associated hybrid states provides additional response forms.¹⁷ For example, a drought that temporarily pushes the ecosystem out of the UV stability range may result in a hybrid state rather than a global transition to the PP state which has lower vegetation coverage and is thus less productive. This flexibility can be lost in the presence of a stable BS state, as a result of a possible collision of the snaking range with the Maxwell point of a UV-BS front. In other words, conditions that stabilize the bare soil state at the relative high precipitation rates, where the UV state is stable, may result in more severe forms of desertification, characterized by lower vegetation coverage, as a result of the disappearance of hybrid states.

The multiplicity of stable front solutions opens up new venues for controlling regime shifts based on shifting existing fronts to alternative stable fronts by local manipulations at the front zones.³⁹ Two types of front multiplicity have been studied here, multiplicity of simple fronts, examples of which are shown in Fig. 3, and multiplicity of complex fronts as Fig. 7 shows. Multiplicity of simple fronts is further demonstrated in the front bifurcation diagrams shown in Fig. 6. The higher-speed front represents a gradual regime shift of bare soil to uniform vegetation, while the lower-speed front represents a gradual regime shift of bare soil to periodic pattern. This result suggests the feasibility of shifting the low-speed front to the high-speed front by local vegetation plantation. The significance of such a manipulation is twofold: speeding up the recovery to a productive vegetation state and achieving a more productive state—uniform vegetation. Similar ideas can be applied to the case of complex fronts when encountered in the field.

Further analysis is needed to better understand the collision of the snaking range with a Maxwell point and the consequent disappearance of the snaking range. In fact, the snaking range can disappear for other reasons as well, such as the presence of a Belyakov-Devaney (BD) point.^40,41 Past this point (Fig. 12), all spatial eigenvalues are real and the basic mechanism responsible for the presence and ultimately stability of localized states, the overlap of oscillatory tails of the profiles of adjacent localized states, is absent.⁴² The possibility that the snaking range collides with the BD point thus provides a second scenario leading to the destruction of localized states.⁴³ 

In this brief report, we have avoided all mention of localized states that bifurcate from the folds of the UV branch and from the folds on the PP branch, which greatly complicate the possible transition scenarios. These states and their properties are the subject of a separate publication.⁴⁴ The extent to which the results reported here apply to more realistic models of dryland vegetation is to be addressed in future studies. Finally, although the specific context we consider here is dryland ecosystems, the results obtained herein are likely to be relevant to other naturally occurring systems exhibiting tristability between both uniform and patterned states.

This work was supported in part by the National Science Foundation Grant DMS-1440386 to the Mathematical Biosciences Institute (PG) and Grant No. DMS-1613132 to the University of California, Berkeley (EK). The work has also received funding from the Israel Science Foundation under Grant No. 305/13.

The space-time plots in Figs. 3, 7, and 11 were obtained using time-integration employing a semi-implicit Euler scheme to march forward in time together with a centered finite-difference scheme for spatial derivatives. Either periodic or Neumann boundary conditions were used. The integrations used a typical time step size of 0.1 and spatial discretization size of 0.5 and were also used to determine the speed and stability properties of the moving fronts (Figs. 2 and 6). Numerical continuation using AUTO⁴⁵ was used to find the bifurcation structure of stationary solutions in Figs. 1, 2, 4, 5, and 12 and the properties of fronts in Figs. 6, 8, 9, and 10. Stability information shown in Figs. 1, 4, 5, 6, and 10 was found by taking solutions from AUTO and calculating their linear stability numerically using a finite-difference scheme, and confirming the results using time-integration.

The separation lines in the parameter space of Fig. 2 were found using two-variable continuation in the variables a and d for the Maxwell line (magenta), the snaking range (red), and the single-patch range (orange). For the bistability range (blue), one-variable continuation for multiple parameter values was used to estimate the boundary at high a. The breakup of the front (dashed gray) was found with a binary search using time integration. In the bifurcation diagram for the front speed (Fig. 6), the curve of the simple front (black) was found using two-variable continuation (a and front velocity), while the speed of the pattern invading bare soil (green) was estimated from time-integration. The properties of simple and complex fronts (Figs. 8, 9, and 10) were determined using two-variable continuation (a and front velocity).

The extended Klausmeier–Gray–Scott model (1a)-(1b) admits a bare soil solution $(B,W)=(0,a)$ that is stable for all parameter values. In addition, the constant solutions $B±=g±g2−1, W±=m/B±=m(g∓g2−1)$ with $g=a/2m$ exist for g > 1. These two constant solution branches are created from a fold bifurcation at g = 1.

In order to study the temporal stability of the spatially uniform states and the emergence of periodic patterns, we linearize the equations about the constant solution, $(B,W)=(B±+b,W±+w)$⁠, and obtain

After substitution for $W±$⁠, the dispersion relation for the growth rate σ of the perturbations (b, w) is

where $k2=k·k$ is the square of the spatial wavenumber.

We can rewrite this relation in the form

The lower branch $B−$ is always unstable while the upper branch gains stability through a Turing bifurcation in which a patterned state emerges when the fastest growing mode (⁠$∂kσ=0, ∂k2σ<0$⁠) is marginally stable (σ = 0). Implicit differentiation of Eq. (B3) yields

which reduces to

at the Turing bifurcation assuming that $kT≠0$⁠. We will also find the following alternate expression for the Turing wavenumber useful:

which has been obtained by combining Eq. (B3) with Eq. (B5). Combining expressions (B5) and (B6) provides the following condition that must be satisfied at the Turing point

which, upon solving for the biomass and substituting into Eq. (B5), gives

The saddle-node bifurcation of the constant solution occurs at B_sn = 1 (⁠$asn=2m$⁠) and the Turing bifurcation occurs at

Since $BT≥1$⁠, we have that $aT≥2m$⁠. Furthermore, $kT2≥0$ requires that $dm≥2$⁠. The temporal stability assignments are summarized in Fig. 12.

We use weakly nonlinear analysis to generate an amplitude equation for the spatially periodic solution near the Turing bifurcation. Earlier studies^23,46 have shown that the Turing bifurcation in the extended Klausmeier/Gray–Scott model (1a)-(1b) is supercritical under the assumption that $d≫1$⁠. We show that subcritical bifurcations are possible when this assumption is relaxed. Our analysis near the bifurcation along the upper constant branch $B+$ supposes that $a=aT+ϵ2α$⁠, where $ϵ≪1$ and $α=O(1)$⁠, and assumes that the envelope of the periodic state depends on a large spatial scale $X=ϵx$ and evolves on a long time scale $T=ϵ2t$⁠. We write the solution as an asymptotic series

Expanding Eqs. (1a) and (1b) in powers of ϵ, we find that the leading order equations satisfy

which are solved by the constant solution

as expected by construction. The order $O(ϵ)$ equations generate the linearization of the system about $v0$

where

The solution to Eq. (C5) can be written in terms of the eigenvector v of $M$ corresponding to the Turing wavenumber solution (viz., $Mv=−kT2v$⁠) as

where A₁₁ depends on the slow time T and long spatial scale X, and is determined from solvability conditions at higher order. We find that in terms of $b0=BT$ at a = a_T

We also find it useful to define

with the property that $v†[−kT2I−M]w=0$ for any w and $v†·v=0$⁠.

The $O(ϵ2)$ equations in the expansion of Eqs. (1a) and (1b) are

Taking the form of the solution to be

results in the three conditions

The condition (C12) requires that

where

Equation (C13) results in

where D₂₁ is an arbitrary function of X and T.

Finally, Eq. (C14) can be solved to give

Moving on to $O(ϵ3)$⁠, we find

The solution takes the form

Substitution of this ansatz generates, for terms proportional to $eikTx$

where the coefficients of the linear and cubic terms are given by

The solvability condition obtained by taking the inner product of the above expression with $v†$ produces a Ginzburg–Landau equation for the leading order amplitude A₁₁ of the patterned state that emerges from the uniform state with wavenumber k_T:

We can rewrite this equation as

where the coefficients are given by

Written in this way, the coefficients r, q, and g only depend on the product dm of the diffusion ratio d and the mortality rate m.

We assume that d > 1 so that the coefficient $m(d−1)>0$⁠. In this case, the constants r and q are positive for all $dm>2$⁠, but g changes sign from positive to negative as dm increases past $dm≈10.0698$⁠, i.e., within the assumed regime d > 1, m > 0, and $dm>2$⁠. Thus a subcritical pattern-forming bifurcation takes place when $dm ≲ 10.0698$⁠, i.e., in a regime in which the diffusion ratio is not too big. In contrast, in the Turing limit $d≫1$ the bifurcation is supercritical. These predictions have been confirmed in our numerical study of the model (1a) and (1b).