## Abstract

Pattern formation at the ecosystem level is a rapidly growing area of spatial ecology. The best studied example is vegetation stripes running along contours in semi-arid regions. Theoretical models are a widely used tool for studying these banded vegetation patterns, and one important model is the system of advection–diffusion equations proposed by Klausmeier. The present study is part of a series of papers whose objective is a comprehensive understanding of patterned solutions of the Klausmeier model. The author focuses on the region of parameter space in which the propagation speed of the patterns is close to its maximum possible value. Exploiting the large value of one of the model parameters, a leading order approximation is obtained for the maximum propagation speed, and the author undertakes a detailed investigation of the parameter region in which there are patterns with speeds close to this maximum.

## 1. Ecological background

Pattern formation at the level of whole ecosystems is a rapidly growing area of spatial ecology. Examples include regular isolated spots of trees and shrubs in savannah grasslands (Lejeune *et al.* 2002; Ben Wu & Archer 2005); patterns of open-water pools in peatlands (Belyea 2007; Eppinga *et al.* 2009); tussock patterns in freshwater marshes (van de Koppel & Crain 2006; Yu 2010); labyrinthine patterns in mussel beds (van de Koppel *et al.* 2005, 2008); striped patterns of tree lines (‘ribbon forests’) in the Rocky Mountains (Hiemstra *et al.* 2006; Bekker *et al.* 2009); and spots, stripes and labyrinths of vegetation in semi-arid environments (Galvin *et al.* 2007; Scanlon *et al.* 2007). Of these, by far the best studied are striped patterns of vegetation, separated by bare ground, running along contours on gentle slopes in semi-arid regions (for reviews, see Valentin *et al.* 1999; Rietkerk *et al.* 2004). These banded patterns are very widespread, occurring in Africa, Australia, North and South America and Asia (see Valentin *et al.* 1999, table 1 and fig. 3). Their study is motivated both by their intrinsic fascination, which is particularly clear from satellite images, and by the susceptibility of these fragile ecosystems to abrupt shifts to total desert (Kéfi *et al.* 2007).

There are no laboratory replicates of banded vegetation, so that empirical study is limited to observation of existing patterns. Because the time scale of pattern evolution is very slow (decades), such observational data are ineffective as a basis for assessing the implications of changes in environmental parameters such as rainfall. Therefore, theoretical models are an important and widely used tool for studying these patterns (Borgogno *et al.* 2009). One group of models is based on the ‘water redistribution hypothesis’. This proposes that rain water falling onto bare ground mostly runs off to nearby vegetated areas. There the infiltration capacity of the soil is higher because of the presence of organic matter and roots (Hills 1971; Callaway 1995; Valentin *et al.* 1999; Rietkerk *et al.* 2000), and thus the rain water can infiltrate. The resulting positive feedback in plant growth provides a potential mechanism for patterning, and this paper is concerned with one mathematical model for this process, due originally to Klausmeier (1999). However, it is important to emphasize that there are alternative hypothesized mechanisms for banded vegetation pattern formation which have also been studied using mathematical models. In particular, some authors argue that the non-locality of water uptake as a result of extended root systems is important. Some mathematical models consider this in addition to water redistribution (e.g. Gilad *et al.* 2007; Kletter *et al.* 2009; von Hardenberg *et al.* 2010), while others set it alongside local facilitation because of shading (e.g. Lefever & Lejeune 1997; Barbier *et al.* 2008; Lefever *et al.* 2009), providing a ‘short range activation, long range inhibition’ mechanism.

The Klausmeier model was the first continuum model for patterning due to water redistribution. Many subsequent models incorporate more specific details of this process, typically involving separate variables for soil and surface water (e.g. HilleRisLambers *et al.* 2001; Rietkerk *et al.* 2002; Ursino 2007, 2009). Some authors have also incorporated features such as rainfall variability (Ursino & Contarini 2006; Guttal & Jayaprakash 2007) and a herbivore population (van de Koppel *et al.* 2002; for other recent extensions, see Kéfi *et al.* 2008; Pueyo *et al.* 2008). These various studies have made important contributions to the ongoing ecological debates on the causes of vegetation patterns, and the most effective management strategies. However, they are primarily simulation based; a recent exception to this is the work of Goto *et al.* (2011), who prove an existence theorem for the model of Gilad *et al.* (2007). In particular, there is only a limited mathematical understanding of even the relatively simple Klausmeier model. This paper is the second in a series whose objective is a comprehensive understanding of patterned solutions of the Klausmeier model, which can act as a springboard both for simulation-based studies of that model and for analysis of more detailed models.

## 2. Model equations and the form of patterns

When appropriately non-dimensionalized (see Klausmeier 1999; Sherratt 2005), the Klausmeier model is 2.1aand 2.1b

Here *u*(*x*,*t*) is plant density, *w*(*x*,*t*) is water density, *t* is time and *x* is space in a one-dimensional domain of constant slope, with the positive direction being uphill. It is important to emphasize that this model is deliberately simple and conceptual. In particular, it does not separate water into surface and sub-surface components. In this regard, it differs from the various more detailed models that include a more mechanistic representation of the water budget (e.g. HilleRisLambers *et al.* 2001; Rietkerk *et al.* 2002; Gilad *et al.* 2007; Ursino 2009). The model assumes constant rainfall, with water lost via both evaporation and uptake by plants. The assumption of proportionality between evaporation and water density is consistent with both more detailed modelling and field data (Rodriguez-Iturbe *et al.* 1999; Salvucci 2001). Per capita plant growth is assumed to be jointly proportional to water density and biomass. The dependence on biomass reflects the positive correlation between the infiltration capacity of soil and its vegetation level in semi-arid environments (Hills 1971; Callaway 1995; Valentin *et al.* 1999; Rietkerk *et al.* 2000). Plant loss is assumed to occur at a constant per capita rate, which may include herbivory. Finally, water flows downhill while plant dispersal is modelled by linear diffusion. The model contains three dimensionless parameters, *A*, *B* and *ν*, which reflect rainfall, plant loss and slope gradient, respectively. Ursino (2005) has performed a detailed calculation that relates *ν* to soil parameters; her work shows that, in general, there will also be a diffusive term for water, although this will be negligible when the mean capillary rise in the soil is small; this depends on soil type.

Pattern solutions of equations (2.1) move in the positive *x* direction (uphill) at a constant speed. This migration has been the subject of a long controversy in the ecological literature (see Tongway & Ludwig 2001, pp. 24–26 for a detailed discussion). A considerable number of field studies do report uphill migration (e.g. Valentin *et al.* 1999, table 5), and this is backed up by a new and very detailed study using satellite images (Deblauwe 2010, ch. 10) that confirms migration in three different geographical locations. The proposed cause of uphill migration is that moisture levels are higher on the uphill edge of the bands than on their downhill edge, leading to reduced plant death and greater seedling density (Montaña *et al.* 2001; Tongway & Ludwig 2001). However, other field data indicate stationary patterns (e.g. Dunkerley & Brown 2002). One suggested explanation for this is that seed dispersal may be preferentially oriented in the downslope direction because of transport in run-off (Saco *et al.* 2007; Thompson & Katul 2009). The model (2.1) neglects any such directed dispersal.

Mathematically, patterns moving with constant shape and speed can be studied via the ansatz *u*(*x*,*t*)=*U*(*z*) and *w*(*x*,*t*)=*W*(*z*), where *z*=*x*−*ct* with *c* being the migration speed. Substituting these solution forms into equations (2.1) gives the travelling wave equations as
2.2aand
2.2bPatterned solutions correspond to periodic solutions of equations (2.2). Previously, Gabriel Lord and I used numerical bifurcation analysis to study these periodic solutions (Sherratt & Lord 2007). We showed that, for a given value of the migration speed *c*, patterns occur for a range of rainfall parameter values *A*, and our results suggested that this range is bounded by a Hopf bifurcation point and a homoclinic solution; details of the relevant homogeneous equilibria are given in §2. A typical result is illustrated in figure 1, which shows the loci of the Hopf bifurcation point and the homoclinic solution in the *A*–*c* parameter plane, for fixed values of *B* and *ν*.

Intuitively, sufficiently high rainfall levels give rise to only homogeneous vegetation, whereas very low rainfall levels imply a complete desert. Intermediate rainfall levels are too low to maintain homogeneous plant cover, but are compatible with vegetation bands.

The parameters *A*, *B* and *ν* will vary between different ecosystems. Estimates of *A* (rainfall) and *B* (plant loss) lie in the ranges 0.1–3.0 and 0.05–2.0, respectively (Klausmeier 1999; Rietkerk *et al.* 2002). Throughout the paper, I make the assumption that *B*<2. This guarantees that equations (2.1) have simple local dynamics: phenomena such as limit cycles can occur in the kinetics for larger values of *B*, but these are not relevant in applications. In comparison to *A* and *B*, the slope parameter *ν* is much larger, with a value of 200 being typical. The non-dimensionalization (see Klausmeier 1999; Sherratt 2005) shows that *ν* depends on the ratio of the advection rate of water and the (square root of the) plant diffusion coefficient; hence, it is very large, even though the slopes on which banded vegetation occurs are quite gentle (a few per cent).^{1} Mathematically, the large value of *ν* suggests an investigation of the asymptotic form of periodic solutions of equations (2.2) for large *ν*. This is the approach taken in this paper, and my focus is the maximum value of the migration speed *c* for which patterns occur. I will show that this speed is *O*_{s}(*ν*) as , meaning that there are constants *k*_{1} and *k*_{2} such that the maximum speed satisfies *k*_{1}*ν*≤*c*≤*k*_{2}*ν* for *ν* sufficiently large.^{2} Further, I will show that, for fixed *B*, the scaling *c*=*O*_{s}(*ν*), *A*=*O*_{s}(1) incorporates an important part of the region of the *c*–*A* parameter plane giving patterns, including the crossing of the homoclinic solution and Hopf bifurcation loci that is shown in the inset in figure 1. It is important to note that *c*=*O*_{s}(*ν*) does not imply that *c* and *ν* have comparable numerical values, since the constants *k*_{1} and *k*_{2} could be either very large or very small. This is exemplified by the crossing point of the homoclinic solution and Hopf bifurcation loci in figure 1, which occurs at *c*≈20, while *ν*=182.5; nevertheless, I will show that this wave speed is *O*_{s}(*ν*) as . Figure 2 shows a typical pattern solution of equations (2.2) close to the crossing point.

In §3, I show that, for large values of the parameter *ν*, the third-order travelling wave equations (2.2) reduce to a second-order system, to leading order, and in §4 I present two main analytical results on periodic solutions of this system. In §5, I describe a numerical study of periodic solutions, which builds on the analysis. Section 6 contains the proof of a theorem stated in §4, and §7 summarizes my results and discusses them in the wider context of ecological travelling waves.

## 3. Leading order travelling wave equations

Periodic solutions of equations (2.2) can be calculated numerically via continuation methods (Sherratt & Lord 2007). Such numerical solutions suggest that for values of *c* close to the maximum, periodic solutions of equations (2.2) have a period that scales with *ν* as . Therefore, I rewrite equations (2.2) using the rescaled coordinate *ζ*=*z*/*ν*, giving
and
Clearly, *c*=*O*_{s}(*ν*) is a distinguished scaling for these equations, which implies the leading order equations as follows:
3.1aand
3.1bSubstituting
3.2simplifies equations (3.1) to
3.3aand
3.3bAlmost all of the remainder of the paper is concerned with periodic solutions of equations (3.3).

For all values of *α* and *β*, equations (3.3) have a steady state (0,*α*); ecologically, this corresponds to a total desert, and it is always an unstable node in equations (3.3). When *α*>2*β*, there are two other steady states;
3.4The eigenvalues at satisfy
3.5Since *α*>2*β*, and . It follows from equation (3.5) that is always a saddle point, with the nature of depending on parameters. I write ; these are the values of at which the discriminant of equation (3.5) is zero. Then, laborious but straightforward investigation of the roots of equation (3.5) shows that the nature of is as follows:
For completeness, I comment that, in the special case *α*=2*β*, the unique positive steady state (1,*β*) has eigenvalues 2−*β* and 0, with the centre manifold associated with the latter being unstable.

A Hopf bifurcation occurs in equations (3.3), initiating a branch of pattern (periodic) solutions, when with *β*>2. Using equation (3.4), this simplifies to
3.6Rewriting these in terms of the original model parameters, using equation (3.2), gives the conditions for Hopf bifurcation as
3.7(Recall that I am assuming *B*<2.) The first of these is the leading order form (for large *ν*) of the locus of Hopf bifurcation points, which will be important in determining the range of rainfall levels giving patterns. The second part of equation (3.7) gives the maximum migration speed for patterns, again to leading order as . This confirms the numerical observation that the maximum speed is proportional to the slope parameter *ν*, and shows that it is an increasing function of the plant loss parameter *B*.

## 4. Analytical results on pattern solutions

I have been unable to obtain a complete analytical picture of the form of the periodic solution branch emanating from the Hopf bifurcation point (3.7), but the following two results give valuable partial information.

### Proposition 4.1

*The Hopf bifurcation point (3.6) of equations (3.3) is subcritical if β*>4 *and supercritical if β*<4.

Using equation (3.2), this implies that for given values of *B*, *ν* and *c*=*O*_{s}(*ν*), the periodic solution branch leaves *A*=*A*_{hb} (defined in equations (3.7)) in the direction of decreasing *A* when *c*<*νB*/(4−*B*), and in the direction of increasing *A* when *c*∈(*νB*/(4−*B*),*νB*/(2−*B*)). From an ecological viewpoint, the proposition implies that patterns with low amplitude occur for values of the rainfall parameter *A* below/above *A*_{hb} when the migration speed *c* is below/above the critical value *νB*/(4−*B*).

### Proof.

This is a standard normal form calculation for the solution amplitude near a Hopf bifurcation point; see for example Guckenheimer & Holmes (1983, §3.4). The website www.ma.hw.ac.uk/~jas/supplements/kl2/ contains a Maple worksheet that performs this calculation. ■

### Theorem 4.2

*For any β>2, there is a value of α>2β at which equations (3.3) have a solution that is homoclinic to the steady state* *, defined in equations (3.4).*

To avoid interrupting the flow of the manuscript, I present the proof of this result in §6. Figure 3 shows phase portraits of equations (2.2) for *α* below, at, and above the value giving a homoclinic solution, for one value of *β*.

From an ecological viewpoint, the theorem concerns patterns for given values of the slope parameter *ν* and the plant loss parameter *B* (<2), and for a fixed migration speed *c* that is large enough that *c*=*O*_{s}(*ν*), but is below the upper limit for patterns (see equations (3.7)). Then, the theorem implies that as the rainfall parameter *A* is varied, there is no upper limit on pattern wavelength. In practice, this translates into a marked sensitivity of pattern wavelength to changes in the rainfall, and this is discussed further in §5.

An important caveat to the theorem is that the proof in §6 gives no information on the relative locations, in parameter space, of the Hopf bifurcation point and the homoclinic solution. To investigate this important issue, I have had to rely on numerical study.

## 5. Numerical investigation of pattern solutions

Using the analytical results of §4 as reference points, I undertook a detailed numerical study of periodic solutions of equations (3.3), using the continuation software AUTO (Doedel 1981; Doedel *et al.* 1991, 2006). It is most instructive to consider behaviour as *α* varies, with *β* fixed. Then, for *β*∈(2,4), the periodic solution branch proceeds monotonically as *α* increases from *α*_{hb} (defined in equation (3.6)), terminating at the homoclinic solution, at *α*=*α*_{hc} say (illustrated in figure 4*a*). For *β*>4, the Hopf bifurcation is subcritical (see proposition 4.1 in §4), and the periodic solution branch leaves the Hopf bifurcation point *α*=*α*_{hb} with *α* decreasing (figure 4*b*,*c*). The branch then folds, at *α*=*α*_{fold} say, after which *α* increases along the branch until it terminates at the homoclinic solution; again I denote the corresponding value of *α* by *α*_{hc}. Figure 4*d* shows the loci of *α*_{hb}, *α*_{hc} and *α*_{fold} in the *α*–*β* plane, showing that *α*_{hc}=*α*_{hb} at a critical value *β*=*β**. Setting *α*−*α*_{hb}≡*α*−*β*^{2}/(*β*−1)^{1/2} as an (overspecified) parameter in AUTO shows that *β**≈4.4036.

The proximity of *α*_{hc} and *α*_{hb} implies that the wavelength of pattern solutions is quite sensitive to the value of the rainfall parameter *A*, which is proportional to *α*. This is illustrated in figure 5, which shows the variation in wavelength with *α* for the values of *β* used in figure 4*a*–*c*.

In Sherratt & Lord (2007), we performed a numerical bifurcation study of the full travelling wave equations (2.2), and concluded that, for all values of *c*, the range of rainfall levels giving patterns is bounded by the loci of the Hopf bifurcation point on the one side and the homoclinic solution on the other side (see figure 1). The more detailed calculations in figure 4 show that this conclusion was incorrect, at least to leading order as . In fact, the solution region giving patterns is slightly wider for some values of the migration speed *c*, being bounded below by the locus of *α*_{fold}. As an illustration of this, and to confirm that the same behaviour does occur for finite values of *ν*, I repeated the numerical continuations used for figure 1, paying careful attention to the detail of behaviour when *c* is large. This confirmed that there is a fold in the periodic solution branch, whose locus is shown in figure 6. One particular implication of this is that there is a region of parameter space, in between the fold and Hopf bifurcation loci, in which there are two different pattern solutions.

As , *α*_{hb}, *α*_{hc} and *α*_{fold} all also. However, detailed numerical investigation suggests that the difference as . The limit corresponds to *c*=*o*(*ν*), suggesting that the Hopf bifurcation remains subcritical when *c*≪*ν*, but that the periodic solution branch does not fold. This is consistent with results that I have derived previously for the case *ν*^{1/2}≪*c*≪*ν* (Sherratt 2010). This suggests that the fold in the periodic solution branch, and hence the occurrence of two different patterns for a given set of parameters, is a phenomenon restricted to the case *c*=*O*_{s}(*ν*).

## 6. Proof of the theorem

### (a) Rescaling the equations

Since theorem (4.2) concerns only the case *α*>2*β*, it is possible to rescale (3.3) via
This is convenient because it maps the steady state to (1,1), with and (0,*α*) mapping to (1/*k*,*k*) and (0,1+*k*), respectively. The rescaled equations are
6.1aand
6.1bNote that there is a 1–1 correspondence between and *k*∈(0,1); however, at various stages in the proof, I will also consider *k*=0, which corresponds to the limit .

Calculation of the eigenvalues of the Jacobian matrix of equations (6.1) shows that, for all *β*>2 and *k*∈[0,1), is a saddle point. I denote by the unique trajectories leaving/entering this steady state into/from the region (figure 7).

### (b) Definition of the functions and

The nullclines are and , while the nullcline is (figure 7). I denote by the open region of the – plane with and (figure 7*a*), and by the open region with , , and with above both and the nullcline (figure 7*b*). Then, comparison of the eigenvectors of the Jacobian matrix of equations (6.1) at (1,1) with the slopes of the nullclines shows that, for all *k*∈[0,1), leaves (1,1) into , while for all *k*∈(0,1), approaches (1,1) from . Since throughout , must leave this region through either the axis or the boundary, except for the intermediate case in which and as along . I define by
6.2

In the region , as at each fixed , and thus cannot enter from infinity unless and as along . Moreover, on the part of the boundary of composed of the nullcline, , implying that all trajectories crossing this part of the boundary do so out of . Therefore, must enter either from or through . I define by
6.3Note that is defined at *k*=0, but is not: this is because the slope of as , so that in the limiting case *k*=0, approaches (1,1) along the boundary of the open region , rather than from the interior.

The points at which leaves and enters vary continuously with *k* (Kopell & Howard 1975), and thus and are continuous functions of *k*. Exploiting this continuity, I will show that there is a value of *k*∈(0,1) at which , corresponding to a homoclinic solution.

### (c) for *k* sufficiently small and positive

When *k*=0, equations (6.1) can be solved exactly,
and
Here, *Γ*_{1} and *Γ*_{2} are constants of integration. Therefore, is given by , implying that . Since is continuous, for *k*>0 sufficiently small, whereas for all *k*∈(0,1), by construction. Therefore, for *k*>0 sufficiently small.

### (d) for *k* sufficiently close to 1

Throughout this subsection, I write *k*=1−*ϵ*.

#### (i) as

Expanding the stable eigenvector of the Jacobian matrix of equations (6.1) at (1,1) as a power series in *ϵ* shows that the slope of at (1,1) is −2/*β*+*O*(*ϵ*) as ; a corresponding expansion of the equation of the nullcline shows that its slope at (1,1) is −1+*O*(*ϵ*). Since these slopes differ by an amount that is *O*_{s}(1) as , one expects intuitively that will remain above the nullcline until its distance from (1,1) is ≫*ϵ*.

Formally, I argue as follows. I define to be the half-line , , with *m*_{1}∈(2/*β*,1) fixed. I then define to be the open region enclosed by and the half-line , (figure 8*a*). Consider now a point on given by , , where 0<*ξ*=*O*(1) as . Then, power series expansion shows that , which is greater than −*m*_{1} for sufficiently small *ϵ*. Therefore, if a trajectory crosses at a point whose distance from (1,1) is *O*(*ϵ*) as , then it must be directed out of . Hence, the point at which (first) crosses must have a distance from (1,1) that is ≫*ϵ* as , i.e. as .

#### (ii) as

Power series expansion of the eigenvectors of the Jacobian matrix of equations (6.1) at (1,1) shows that the slope of at (1,1) is −1−*ϵ*/(*β*−2)+*O*(*ϵ*^{2}) as . I define to be the line in the – plane through (1,1) with slope −1−*ϵm*_{2}, where *m*_{2}>1/(*β*−2) (figure 8*b*), and I define to be the open region enclosed by , the nullcline , and the line (figure 8*b*). Then, enters , by construction. Now consider a point on given by , , where 0<*ξ*=*O*(1) as . Then, power series expansion shows that
for sufficiently small *ϵ*, since *m*_{2}>1/(*β*−2)>0. Therefore, for sufficiently small *ϵ*, all trajectories crossing the part of the boundary of are directed into . Also, on the part of the boundary of , and , so that again all trajectories are directed into . Therefore, must leave through the line (figure 8*b*). The point at which this occurs has
Recall that the only constraint on *m*_{2} is *m*_{2}>1/(*β*−2). Since the point is independent of the choice of *m*_{2}, it follows that
This estimate is actually more precise than needed: I will only use the fact that as .

The remainder of this section is the most delicate part of the proof. Recall from §3 that the steady state (1/(1−*ϵ*),1−*ϵ*) is a focus. The basic argument is that after the point , rotates around (1/(1−*ϵ*),1−*ϵ*), and reaches the nullcline at a point whose distance from (1/(1−*ϵ*),1−*ϵ*) is still *O*(*ϵ*^{2}). The complication is that, as *ϵ* decreases, not only does approach (1/(1−*ϵ*),1−*ϵ*), but also equations (6.1) change, and in fact is a singular limit. It is most convenient to work in polar coordinates centred at (1/(1−*ϵ*),1−*ϵ*),
and
It is straightforward to reformulate equations (6.1) in terms of *r* and *θ*. I then substitute *r*=*ϵ*^{2}*ρ*, where *ρ*=*O*(1) as , and expand the equations as power series in *ϵ*, giving
6.4and
6.5where with *ϕ*∈(*π*/4,*π*/2). At the point , *θ*=*π* ⇒ for sufficiently small *ϵ*, and *θ* continues to increase with provided 7*π*/4−*θ*≫*ϵ*. Along this part of , as , implying that there is an *O*(1) increase in and hence an *O*(1) increase in *ρ*, using equation (6.4). Moreover, I have shown that *r*=*O*(*ϵ*^{2}) at the point , and thus *r* remains *O*(*ϵ*^{2}) throughout this part of .

I must now consider how *θ* changes along in the vicinity of 7*π*/4, specifically when 7*π*/4−*θ*=*O*(*ϵ*) as . This is an important region because power series expansion shows that close to (1/(1−*ϵ*),1−*ϵ*), and for , the polar coordinate equation of the nullcline is *θ*=7*π*/4+*ϵ*+*O*(*ϵ*^{2}). To consider the solution for in this region, it is necessary to include the *O*_{s}(*ϵ*) contribution to as well as the *O*_{s}(1) terms given in equation (6.5). I omit the details for brevity; writing *δ*=*θ*−7*π*/4, the result is
Therefore, whenever *δ*≤*ϵ*, for sufficiently small *ϵ*, implying that continues to rotate around (1/(1−*ϵ*),1−*ϵ*) until it crosses the nullcline. Moreover, whenever *θ*−7*π*/4=*O*(*ϵ*). Therefore, during the part of in which *θ*−7*π*/4=*O*(*ϵ*), there is an *O*(1) increase in and hence an *O*(1) increase in *ρ*, using equation (6.4). Therefore, *r*=*O*(*ϵ*^{2}) throughout this part of , also.

Taken together, these results imply that the (first) point at which intersects is located at a distance from (1/(1−*ϵ*),1−*ϵ*) that is *O*(*ϵ*^{2}) as , implying that .

#### (iii) Comparison of and

The results in §6*d*(i), (ii) together imply that for *ϵ* sufficiently small, i.e. .

### (e) Completion of the proof

Since the functions and are continuous (see §6*b*), and recalling that *ϵ*=1−*k*, §6*c*,*d* together imply that there is a value of *k*∈(0,1) for which . Then, the trajectories and are the same, and constitute a solution of equations (6.1) that is homoclinic to .

## 7. Discussion

The Klausmeier model for banded vegetation in semi-arid environments is the oldest and simplest of a number of continuous models for patterning due to water redistribution. In this paper, I have focussed on patterns with migration speeds close to the maximum possible value, for large values of the slope parameter. I have shown that the mathematical origin of patterns is a Hopf bifurcation in the travelling wave equations that can be either subcritical or supercritical, and I have shown that the branch of pattern solutions terminates at a homoclinic solution, whose existence I have proved. From an ecological viewpoint, the most important implications of my results concern the way in which pattern solutions vary with the rainfall parameter. I have shown that, for a fixed migration speed, pattern wavelength depends sensitively on rainfall (see figure 5). In fact, field data (e.g. Valentin & d'Herbès 1999) and simulations of equations (2.1) (Sherratt & Lord 2007) both suggest that pattern wavelength remains constant when rainfall varies. My results then imply a significant variation in migration speed; moreover, this variation is very insensitive to the constant value of the wavelength (figure 9). The available field data on migration speeds are not sufficiently detailed to enable this to be tested. Moreover, the practicalities of measuring band migration (see Deblauwe 2010, ch. 10) mean that it can only be assessed over period of at least 5–10 years, so that measured speeds can only reflect multi-year averages of rainfall. Nevertheless, my prediction may be a possible explanation for the long-standing and on-going debate on the extent of migration in banded vegetation patterns (e.g. Tongway & Ludwig 2001).

Pattern solutions of the Klausmeier model (2.1) are periodic travelling waves, also known as wavetrains. Over the last decade, it has become increasingly clear that periodic travelling waves play an important role in spatial ecology. They have been studied almost exclusively for cyclic populations, meaning that even in the absence of spatial variation, the populations would exhibit multi-year oscillations in abundance. Field data on a number of cyclic populations reveal periodic travelling waves (see Sherratt & Smith 2008 for review). Moreover, periodic travelling waves are predicted by various types of mathematical models for cyclic ecological systems: reaction–diffusion equations (Sherratt *et al.* 1995; Petrovskii & Malchow 1999; Smith *et al.* 2008); integro-differential equations (Gourley & Britton 1993; Ashwin *et al.* 2002); integro-difference equations (Kot 1992); and cellular automata (Sherratt 1996). In the case of oscillatory reaction–diffusion equations, mathematical theory has enabled a detailed understanding of how a particular periodic travelling wave solution is selected from the wave family by initial and boundary conditions (Petrovskii *et al.* 1998; Sherratt 2001; Sherratt *et al.* 2003; Garvie 2007). A corresponding investigation for the Klausmeier model (2.1) is an important objective for future research. Preliminary results (Sherratt & Lord 2007) indicate that wave selection is dependent not only on initial and boundary conditions, but also on pattern history, suggesting that this will be a complex problem requiring careful study.

Compared to the substantial literature on periodic travelling waves in oscillatory reaction–diffusion systems, research on these solution forms in advection–reaction–diffusion equations such as equations (2.1) is very much in its infancy. The largest body of work concerns ‘differential flow-induced chemical instability’, where the advection term is due simply to fluid flow in an aqueous system of reactants. For a number of specific chemical systems, conditions for patterning have been derived, and in some cases tested against experiment (Merkin *et al.* 2000; Nekhamkina & Sheintuch 2003; Taylor 2003). There has also been some work on pattern formation in advection–reaction–diffusion equations applied to developmental biology (Perumpanani *et al.* 1995; Bernasconi & Boissonade 1997; Satnoianu & Menzinger 2002); in that context, a key issue is the phase difference between two interacting morphogens, which depends on the advection coefficients. In oscillatory reaction–diffusion equations, periodic travelling waves are important both in their own right, and because they play a key role in the transition to spatio-temporal chaos (Petrovskii & Malchow 2001; Sherratt *et al.* 2009). This suggests that a fuller understanding of periodic travelling waves in ecological models of advection–reaction–diffusion type has the potential not only to reveal details of periodic patterning, but also to provide new insights into other, more complex types of spatio-temporal dynamics in ecological systems.

## Acknowledgements

This work was supported in part by a Leverhulme Trust Research Fellowship. I am grateful to Jack Carr (Heriot–Watt University) for helpful discussions.

## Footnotes

↵1 On steeper slopes, rainwater does not flow downhill in a sheet; rather, gullies form. Then, the model (2.1) is invalid, and different phenomena occur.

↵2 Throughout the paper, I also use the standard notations

*O*(⋅) and*o*(⋅). For example,*f*=*O*(*ν*) as ⇔ there is a constant*C*_{1}>0 such that |*f*|≤*C*_{1}*ν*for all sufficiently large*ν*; intuitively, the order of magnitude of*f*is the same as, or smaller than, that of*ν*. For*o*(⋅),*f*=*o*(*ν*) as ⇔ as ; intuitively, the order of magnitude of*f*is smaller than that of*ν*. Thus*f*=*O*_{s}(*ν*) as ⇔*f*=*O*(*ν*) and*f*≠*o*(*ν*); intuitively, the order of magnitude of*f*is the same as that of*ν*.

- Received March 25, 2011.
- Accepted June 2, 2011.

- This journal is © 2011 The Royal Society