PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
 
J Biol Dyn. 2010 July; 4(4): 346–380.
Published online 2009 October 28. doi:  10.1080/17513750903301954
PMCID: PMC2903777

Nonlinear stability analyses of vegetative pattern formation in an arid environment

Abstract

The development of spontaneous stationary vegetative patterns in an arid isotropic homogeneous environment is investigated by means of various weakly nonlinear stability analyses applied to the appropriate governing equation for this phenomenon. In particular, that process can be represented by a fourth-order partial differential time-evolution logistic equation for the total plant biomass per unit area divided by the carrying capacity of its territory and defined on an unbounded flat spatial domain. Those patterns that consist of parallel stripes, labyrinth-like mazes, rhombic arrays of rectangular patches, and hexagonal distributions of spots or gaps are generated by the balance between the effects of short-range facilitation and long-range competition. Then those theoretical predictions are compared with both relevant observational evidence and existing numerical simulations as well as placed in the context of the results from some recent nonlinear pattern formation studies.

Keywords: pattern formation, weakly nonlinear analysis, vegetation distribution

1. Introduction

Rayleigh-Bénard buoyancy-driven convection has to date provided perhaps the best-studied example of nonlinear pattern selection (reviewed by Koschmieder [10]). One of the methods traditionally used to predict such pattern selection is a weakly nonlinear stability analysis that, although incorporating nonlinearities of the relevant model system, basically pivots a perturbation procedure about the critical point of linear stability theory (reviewed by Wollkind et al. [31]). The advantage of such an approach over strictly numerical procedures is that it allows one to deduce quantitative relationships between system parameters and stable patterns, which are valuable for experimental design and difficult to accomplish by using simulation alone. Recently, there has been considerable interest generated in pattern formation and selection during normal-incidence ion-sputtered erosion of metallic or semiconductor solid surfaces and during diffraction of radiation induced by the injection of a laser pump field into a ring cavity containing a purely absorptive two-level atomic sodium vapour medium. In order to predict the sequence of solid surface morphologies and nonlinear optical patterns actually observed during such ion-sputtered erosion and laser-induced diffraction, respectively, Pansuwan et al. [19] and Wollkind et al. [30] performed both rhombic- and hexagonal-planform weakly nonlinear stability analyses on the governing damped Kuramoto-Sivashinsky and modified Swift-Hohenberg equations, respectively, appropriate for modelling these phenomena.

We wish to continue our examination of nonlinear phenomena by investigating the development of spontaneous stationary equilibrium vegetative patterns in an arid isotropic homogeneous environment (reviewed by Rietkerk et al. [21]). Pattern formation in this phenomenon can be modelled for the isotropic case by a partial differential time-evolution logistic equation describing the total plant biomass per unit area divided by its carrying capacity that under a weak gradient-low density truncation is reduced to the fourth order in its spatial variables defined on an unbounded two-dimensional flat domain [13]. Under appropriate conditions, those patterns that consist of parallel stripes, labyrinth-like mazes, rhombic arrays of rectangular patches, and hexagonal distributions of spots or gaps can be generated in an initially homogeneous environment by the balance between the effects of short-range facilitation and long-range competition. We shall perform these same rhombic- and hexagonal-planform weakly nonlinear stability analyses on that truncated fourth-order partial differential time-evolution logistic equation, and then compare the theoretical predictions obtained with both relevant observational evidence and existing numerical simulations as well as place them in the context of some recent pattern formation studies. We begin below with a brief sketch of the reduction procedure required to obtain the governing truncated fourth-order partial-differential time-evolution logistic equation.

Lefever and Lejeune [12] and Lefever et al. [13] proposed a nondimensional propagator-inhibitor logistic equation of the form

equation image

in order to model spatio-temporal patterning for botanical ecosystems. Here ρ(x, y, t) is the total plant biomass per unit area divided by its carrying capacity where (x, y) [equivalent] a two-dimensional spatial coordinate system and t [equivalent] the time while F1 · F2 and F3 represent a nonlocal birth term and a local death term, respectively. In particular, F1, which represents the actual propagation distribution function and involves the growth effects of reproduction, seed dissemination, germination, and maturation of young plants, depends mainly on vegetative aerial structures like the canopy in the case of trees or bushes, while F2, which represents the inhibition distribution function and involves effects that inhibit the development of vegetation such as root systems and the uptake of scarce resources like water and nutrients, depends mainly on the availability of free soil for vegetative colonization and on soil vegetative resources of limited supply. Specifically, for the isotropic case employing Gaussian weighting functions in conjunction with a Taylor series expansion for ρ, these distribution functions take the following form [13]:

equation image

where the gradient operator [nabla] [equivalent] ([partial differential]/[partial differential]x, [partial differential]/[partial differential]y), [nabla]0 [equivalent] 1, and [nabla]2 [equivalent] [nabla] · [nabla] while L and Λ are the facilitation-to-competition range and interaction ratios, respectively. Finally,

equation image

where the mortality-to-growth rate ratio μ serves as a measure of the environment's aridity. We note that in Equation (1) time and space have been scaled with the inverse of the vegetation growth rate, which corresponds roughly to the period spent by the plants to achieve adult size and with the interplant competition range, which can be approximated by the radius of the superficial root system of the dominant plants, respectively.

It is now possible to truncate the partial derivatives appearing in Equation (1) by introducing scaling laws such that

equation image

equation image

and hence to third order

equation image

Thus upon substitution of Equations (1c) and (2c), Equation (1a) reduces to the truncated form

equation image

Observe that this development is consistent with the scaling laws

equation image

for which each term in Equation (3a) is of O(ε3). Here, the parameter ε can be identified with the mean biomass density α, to be defined below in Equation (8), which is relatively small for arid or semi-arid environments. These scaling laws differ from those introduced by Lefever et al. [13], and that discrepancy will be discussed in more detail in the last section of this paper. Therefore as observed by Lejeune and Tlidi [14], Couteron and Lejeune [2], Lejeune et al. [15], and Lejeune et al. [16], Equation (3) has been derived in a low density and weak gradient limit from a generalized logistic equation describing nonlocal interactions. This truncated equation applies to water-poor ecosystems where average plant phytomass density is low with respect to carrying capacity (close-packing density) of unstressed vegetation and for which pattern wavelengths are large in comparison with the average size of the dominant plant form but very small when compared with the territorial length-scale characteristic of that arid environment. Since under this latter condition, the actual territorial boundary of that arid environment does not significantly influence the patterns [8], it then seems reasonable as a first approximation for us to consider Equation (3a) on an unbounded two-dimensional flat spatial domain.

Rewriting that equation in the form

equation image

where

equation image

we first seek a uniform stationary solution, ρ = ρ0, of it satisfying

equation image

and find that

equation image

Note that ρ0 = 0 which represents bare ground always exists while when ρ±0 are real and positive, these other two solutions represent spatially homogeneous plant distributions. If 0 ≤ Λ, μ ≤ 1 then ρ+0 ≥ 0, ρ0 ≤ 0 and hence only ρ+0 is biologically meaningful. If Λ > 1, then ρ+0 > 0 for 0 ≤ μμ* = 1 + (Λ − 1)2/4 while ρ0 > 0 for 1 ≤ μμ*. Performing a linear stability analysis of these distributions to homogeneous disturbances by considering a solution of (4) of the form

equation image

yields the perturbation equation

equation image

where

equation image

Since

equation image

and we have linear stability in Equation (6) when f′ (ρ0) < 0 and instability when f′ (ρ0) > 0, it can be concluded from Equation (7) that ρ0 = 0 is linearly stable for μ > 1 and unstable for 0 ≤ μ < 1, ρ0 = ρ+0 is linearly stable everywhere it exists, and ρ0 = ρ0 is linearly unstable everywhere it exists.

We are interested in determining conditions under which nonzero homogeneous distributions of vegetation that are linearly stable in the absence of spatial effects can become unstable to heterogeneous perturbations. This is reminiscent of the diffusive instabilities that occur during chemical Turing pattern formation [32]. Thus we need only consider the solution ρ+0 which will now be designated by

equation image

In order to eliminate the possibility of bistability between that homogeneous solution and the zero state, in what follows after Lejeune et al. [16] we shall restrict our attention to

equation image

In this context, we adopt the far-field condition that

equation image

which is implicitly satisfied by ρ [equivalent] α. Observe that from Equations (5b) and (8) we can deduce the general relation

equation image

which for the special representative case for Λ = 1 reduces to

equation image

It is the weakly nonlinear stability of this homogeneous solution to one-dimensional and two-dimensional rhombic- or hexagonal-planform heterogeneous perturbations with which we shall be concerned in Sections 2–4, respectively.

2. One-dimensional patterns: a nonlinear stability analysis of the Stuart-Watson type

As a necessary prelude to our two-dimensional rhombic- and hexagonal-planform nonlinear stability analyses of the truncated fourth-order logistic equation (4) to be presented in the next two sections, we perform a weakly nonlinear stability analysis of its homogeneous state in this section by seeking a real one-dimensional Stuart [25]–Watson [28] type solution of it through third-order terms of the form [31]

equation image

where the amplitude function A1(t) satisfies the Landau equation

equation image

and q = 2π/λ, λ being the wavelength of the class of spatially periodic perturbations under investigation. Substituting the solution of Equation (12) into Equation (4), we obtain a sequence of problems, one for each pair of values m and n which corresponds to a term of the form Λn1 (t) cos(mqx) appearing in Equation (12a). Then the n = m = 0 problem is satisfied identically by virtue of the fact that

equation image

while the linear stability problem for n = m = 1 yields the secular equation

equation image

which is a parabola in the q2σ plane. Note that if β ≥ α, then σ < 0 for the parameter range of Equation (9) while if β < α then this parabola has a maximum value at its vertex (q2c, σc) with

equation image

where the discriminant of the quadratic in q2 of Equation (14a) is given by

equation image

If, in addition, D(α, β) > 0, 0 < σ < σc(α, β) for q21 < q2 < q22 where

equation image

This scenario is plotted in Figure 1 for the special case of Λ = 1 when α = 0.2 and β = 0.02. Hence, in what follows we shall equate the q and σ of Equation (12) to

equation image

Then

equation image

where βc = βc(α) is defined implicitly by D(α, βc) = 0 or explicitly by

equation image

Therefore, the homogeneous solution is linearly stable for β > βc(α), neutrally stable for β = βc(α), and unstable for 0 < β < βc(α). Thus β = βc(α) serves as the marginal stability curve in the α − β plane. This situation is plotted in Figure 2 for the special case of Λ = 1. In that instance

equation image

where

equation image

while, since βc′;(α) = 1 − 3α1/2, βc′(α0) = 0 implies

equation image

Here we are restricting our analysis to the critical wavenumber alone. In Section 5, we shall discuss the results of the so-called side-band instabilities analysis originally introduced contemporaneously by Segel [23] and Newell and Whitehead [18] that examines the consequences of investigating other wavenumbers in the side band centred about that critical wavenumber as well.

Figure 1.
A plot of Equation (14) in the q2σ plane for Λ = 1, α = 0.2, and β = 0.02.
Figure 2.

A plot in the α − β plane of the marginal stability curve β = βc (α) of Equation (16) for Λ = 1 on which the growth rate σc(α, β) = 0 denoted by the solid line. Here the dashed ...

Continuing our description of the results of this one-dimensional expansion procedure, the second-order problems corresponding to n = 2 and m = 0 or 2 can be solved in a straightforward manner to yield

equation image

and

equation image

Although there are also two third-order problems, it is permissible to concentrate our attention exclusively on the one corresponding to n = 3 and m = 1 which contains the Landau constant a1 for the Fredholm-type method of solvability we shall use. That problem may be represented as

equation image

where

equation image

Here we have only denoted the β dependence of the quantities in question for ease of exposition. Now taking the limit of Equation (18) as ββc = βc(α), employing Equations (14b), (14c), (15), and (17), and assuming the requisite continuity at β = βc, we obtain the relevant solvability condition represented symbolically as

equation image

or given explicitly by

equation image

A main feature of the weakly nonlinear stability theory implicit to the formulation of our expansion of Equation (12) with q [equivalent] qc = qc(α, β) is a phenomenological interpretation of the problem under examination based upon the long-time behaviour of the solution A1 (t) to the truncated amplitude Equation (12b). In general, the dynamics of such a Landau equation with real coefficients can be divided into the four qualitatively different cases represented by the possibility of σ, its growth rate, and a1, its Landau constant, being either positive or negative. These cases can be catalogued as follows:

  • Case (i). σ, a1 > 0. There exists a stable equilibrium solution A2e = σ/a1. Since σ < 0(β < βc), the linear theory would predict instability whereas our nonlinear analysis shows the existence of this finite-amplitude supercritically stable equilibrium state.
  • Case (ii). σ > 0, a1 < 0. The undisturbed state of A1 = 0, which corresponds to a uniform homogeneous distribution, is unstable in this case as well; however, now finite-amplitude effects tend to enhance such disturbance growth.
  • Case (iii). σ, a < 0. There exists an unstable equilibrium solution A2e = σ/a1. This is an instance of a subcritical instability in that the linear theory predicts stability of the uniform homogeneous distribution to infinitesimal disturbances since σ < 0(β > βc) whereas the nonlinear theory shows that it can be unstable to disturbances the amplitudes of which satisfy A21 > σ/α1. Such an occurrence is often termed a metastable state.
  • Case (iv). σ < 0, a1 > 0. The uniform homogeneous distribution is stable to both infinitesimal and finite-amplitude disturbances.

From the description of the four cases just summarized, we can see that the stability behaviour of the Landau equation is crucially dependent upon the sign of a1. Hence in order to determine this behaviour, it is necessary that we next examine the formula for a1 given by Equation (18d). Towards that end we plot a1 of Equation (18d) versus α in Figure 3 with Λ = 1. From this figure, we observe that a1 > 0 whenever

equation image

where

equation image
Figure 3.
A plot of a1 of Equation (18d) versus α for Λ = 1.

Having identified those regions of Figures 2 and and33 corresponding to the cases catalogued above, we now offer a morphological interpretation of each of these cases in the αβ plane.

  • Case (i). Since σ, a1 > 0 for 0 < β < βc(α) and α1 < α < α2, we can conclude that in this parameter range the instability under examination represents a periodic one-dimensional equilibrium pattern consisting of stationary parallel stripes having a characteristic wavelength of
    equation image
    Then designating this equilibrium amplitude by
    equation image
    we have
    equation image
    which is plotted in the x − ρ plane of Figure 4 for Λ = 1, α = 0.12, and β = 0.02. These supercritical stripes are represented in the contour plot of Figure 5, where after Lejeune et al. [16] regions of higher density (ρ > α) appear dark (green) and those of lower density (ρ < α), light (tan).
    Figure 4.
    A plot of the equilibrium density of Equation (21) versus x for Λ = 1, α = 0.12, and β = 0.02.
    Figure 5.
    A contour plot of the supercritical stripes of Figure 4.
  • Case (ii). For 0 < β < βc(α) and α < α1 or α > α2, the uniform homogeneous distribution is unstable since σ > 0 and a1 < 0.
  • Case (iii): For β > βc(α) and α < α1 or α > α2, since σ, a1 < 0, there is metastability between the uniform homogeneous and subcritical striped states.
  • Case (iv). For β > βc (α) and α1 < α < α2, the uniform homogeneous distribution is stable since σ < 0 and a1 > 0.

3. Two-dimensional patterns: a rhombic planform analysis

In order to investigate the possibility of occurrence for the truncated fourth-order logistic equation of the type of two-dimensional vegetative patterns mentioned earlier, we shall first consider a rhombic planform solution of Equation (4) of the form [19,30]

equation image

where

equation image

and

equation image

Here we are employing the notation ρnjmk for the coefficient of each term in Equation (22a)(22b) of the form An1(t)Bj1(t) cos[qc(mx + kz)] where ρ1010 = ρ0101 = 1, implicitly. Then substituting this rhombic planform solution of Equation (22) into Equation (4), we again obtain a sequence of problems, each of which corresponds to one of these terms. Solving those problems we find that

equation image

as defined in the previous section while

equation image

and b1, in particular, satisfies

equation image

Now taking the limit of Equation (23c) as ββc(α) and making use of Equations (23a) and (23b) as well as Equations (14b), (14c), and (17a), we obtain the solvability condition

equation image

Having developed these formulae for its growth rate and Landau constants, we turn our attention to the rhombic planform amplitude Equation (23c), which possess the following equivalence classes of critical points (A0, B0):

equation image

Assuming that a1, a1 + b1 > 0 and investigating the stability of the critical points of Equation (24) by seeking a solution of our amplitude equations of the form

equation image

one finds that p has the associated roots [29]

equation image

which yield the stability criteria that I is stable for σ < 0; II, for σ > 0, b1 > a1; and V, for σ > 0, a1 > b1. Note that I and II, as in the one-dimensional analysis of the previous section, represent the uniform homogeneous and supercritical striped states, respectively, while V can be identified with a rhombic pattern possessing the characteristic angle ϕ (see below and [32]). We next use these criteria to refine our one-dimensional predictions relevant to the former states due to the presence of the latter. Towards that end, we examine the signs of a1 ± b1 for Λ = 1, α [sm epsilon] (α1, α2), and ϕ [sm epsilon] (0, π/2) with ϕ = π/2 (or equivalently 90°) representing a square planform. We first illustrate this procedure by defining the ratio of Landau constants [5]

equation image

and plotting that quantity versus α in Figure 6 for a fixed value of ϕ, namely ϕ = 90°. Here, there exist two α-intervals of stable square rhombic patterns, where a1 ± b1 > 0 or equivalently − 1 < γ < 1, given by

equation image

provided in addition that σ > 0 or 0 < β < βc(α). These intervals have been indicated by shading in this figure and observe that for α values between them γ > 1or b1 > a1 > 0. Hence for α [sm epsilon] (0.1176, 0.2094) and 0 < β < βc(α) there exist stable supercritical striped patterns. Now, repeating the process used to produce Figure 6 but for other values of ϕ, we find the same generic behaviour as for ϕ = π/2 in that there are two intervals of stable rhombic patterns for α [sm epsilon] (αm, αl) or α [sm epsilon] (αr, αM) with stable stripes between them for α [sm epsilon] (αl, αr) when 0 < β < βc(α) where

equation image

and summarize these results in Table 1.

Figure 6.

A plot of γ of Equation (27) versus α for ϕ = π/2 and Λ = 1.

Table 1.

α-range for stable rhombic patterns versus ϕ.

In addition, we can deduce from Equations (18d) and (23d) that

equation image

or

equation image

which implies that only stripes can be stable with respect to rhombic patterns of such small characteristic angle. Next, we plot our ratio of Landau constants versus ϕ [sm epsilon] [0, π] for the fixed values of Λ = 1 and α = 0.12 in Figure 7. Restricting ourselves to the interval of interest ϕ [sm epsilon] [0, π/2], we see that there are two bands of stable rhombic patterns flanking ϕ = π/3 for ϕ [sm epsilon] (ϕm, ϕl) or ϕ [sm epsilon] (ϕr, ϕM) when 0 < β < βc(α) where

equation image

which have been designated by shading, with no stable patterns between these bands for ϕ [sm epsilon] (ϕl, ϕr) and stable stripes outside of them for ϕ [sm epsilon] [0, ϕm) or ϕ [sm epsilon] (ϕM, π/2] when 0 < β < β c(α). This figure has been drawn for the extended interval ϕ [sm epsilon] [π/2, π] in order to demonstrate graphically the symmetry about ϕ = π/2 characteristic of rhombic patterns since

equation image

Here, properties (29) and (31) are a consequence of mode interference occurring at exactly ϕ = 0 and modal interchange, respectively [3]. Note in this context that

equation image

or

equation image

and thus, for α = αc = 0.16 and 0 < β < βc = 0.032, stripes are stable versus rhombic patterns. Again, repeating the process used to produce Figure 7 but for other ααc, we find the same generic behaviour as for α = 0.12 and summarize these results for selected values in Table 2.

Figure 7.
A plot of γ of Equation (27) versus ϕ for α = 0.12 and Λ = 1.
Table 2.

Angle range for stable rhombic patterns versus α.

We observe from Figure 7 and Table 2 that there exist no stable rhombic patterns of characteristic angle ϕ = π/3 [congruent with] 1.047 by virtue of the fact that for α [sm epsilon] (α1, α2)

equation image

The occurrence of this vertical asymptote is a consequence of b1 of Equation (29) possessing a negative numerator and a denominator having a double root at ϕ = π/3.

We close this section with a morphological instability interpretation of the potentially stable critical points II and V of our amplitude equations when σ > 0 or 0 < β < βc (α) relative to the vegetative patterns under investigation. Then, to the lowest order the deviation from its homogeneous configuration of the function associated with these critical points is given by

equation image

The contour plot in the x − y plane for this deviation function with A0 > 0 and B0 = 0 relevant to critical point II is identical to that in Figure 4, of the previous section. Here, as indicated earlier, regions of higher density appear dark and those of lower density, light. Clearly, such alternating dark and light parallel bands should be identified with a striped vegetative pattern as anticipated above. In order to make an analogous interpretation of critical point V, we consider the deviation function (34) with B0 = A0 > 0 and allow ϕ to take on the values 90° and 58°, sequentially. We represent the contour plot for that function with ϕ = 90° in Figure 8. From the checkerboard appearance of the array, it is equally clear that this critical point should be identified with a vegetative pattern of square planform. Finally, in a similar manner, we generate the contour plot of Figure 9 for 58°, which forms a family of rectangles. Here, ϕ = 58° (or equivalently 1.01) is a representative rhombic angle associated with stable patterns for α = 0.15 (see Table 2). To demonstrate that the same thing happens with any ϕ [sm epsilon] (0, π/2], we first put Equation (34) for critical point V in the form

equation image

where

equation image

or

equation image

and

equation image

From Equations (35) and (36), we can then deduce that the intersecting level curves in the associated contour plots are two families of straight lines possessing slopes of

equation image

respectively, which intersect at right angles since

equation image

Given that, in general, this is a rectangular planform we need to explain in what sense such a critical point can be identified with a rhombic pattern. To do so, we form the quadrilateral depicted by solid lines in Figure 9. Its sides are each composed of two half-diagonals, collectively contained in the four light rectangles surrounding a dark one. Thus, each side of the quadrilateral has the length of one of these diagonals

equation image

and hence, the quadrilateral is a rhombus of characteristic angle ϕ. Further, ϕ also plays a role in characterizing the family of rectangles. Each member of that family has aspect ratio

equation image

where w and l are its width and length, respectively, while ϕ/2 serves as its angle of inclination as well, an assertion most easily verified by the relation

equation image

Therefore, in what follows we shall refer to such a vegetative pattern as a rhombic array of rectangular patches.

Figure 8.
A contour plot of Equation (34) for critical point V of Equation (24) with Λ = 1, α = 0.1, β = 0.02, and ϕ = π/2.
Figure 9.
A contour plot of Equation (34) for critical point V of Equation (24) with Λ = 1, α = 0.15, β = 0.02, and ϕ = 1.01. Here the quadrilateral formed by the solid lines depicts the rhombic symmetry of the rectangular pattern. ...

4. Two-dimensional patterns: a hexagonal planform analysis

We continue our investigation of two-dimensional vegetative patterns by next seeking a hexagonal planform solution of Equation (4) that to lowest order satisfies [19,30]

equation image

where

equation image

equation image

equation image

The weakly nonlinear stability behaviour of the amplitude-phase equations (39c), (39d) to be described below depends only on the values of their growth rate and Landau constants. We can determine the solvability conditions for these Landau constants most easily by introducing the transformation

equation image

which reduces Equation (39) to

equation image

where

equation image

equation image

Hence, we consider solutions of Equation (4) of the form of Equation (39a) with [19,30]

equation image

Now we are employing the notation ρnjMk for the coefficient of each higher-order term in Equation (42) of the form Λn1(t)Bj1(t) cos(Mqcx/2) cos(k√3qcy2). Then substituting the solution of Equations (39)–(42) into Equation (4) and solving the resultant sequence of problems, we find specifically that

equation image

and, in particular, the other two Landau constants satisfy

equation image

equation image

Finally, taking the limit of Equation (44) as β → βc(α) and employing Equation (43) as well as Equations (14b), (14c), (15), and (17), we obtain the solvability conditions

equation image

equation image

Note from Equations (44a) and (44b) that the expression a2(α) of Equation (44d) does not contain the component limβ→βc(α) ρ1111 since the coefficient of that quantity, namely limβ→βc(α) (4α0 + Λ − 1 − 3α + q2c/2 − q4c/8), is identically equal to zero by virtue of Equation (44a). As pointed out by Wollkind et al. [31] such independence may be expected as a general rule, thus eliminating the necessity of determining ρ1111.

Having developed formulae for their growth rate and Landau constants, we return to the six-disturbance hexagonal-planform amplitude-phase equations (39c) and (39d). In cataloguing the critical points of these equations and summarizing their orbital stability behaviour, it is necessary to employ the quantities

equation image

Then, those critical points can be represented as equivalence classes of the form (A0, B0, B0, 0, 0, 0) corresponding to ϕ1 = ϕ2 = ϕ3 = 0 with

equation image

where it is assumed that a1, a1 + 4a2 > 0. In order to examine the stability of these critical points of Equation (46), we seek solutions of the amplitude-phase equations (39c) and (39d) of the form

equation image

equation image

Using similar but slightly more complicated considerations than those employed in the last section, one can deduce orbital stability criteria for these critical points that are posed in terms of σ [29]. Thus critical point I is stable in this sense for σ < 0 while the stability behaviour of II and III±, which depends on the signs of a0 and 2a2a1 as well, has been summarized in Table 3 under the assumption that a1 + a2 > 0.

Table 3.

Orbital stability behaviour of critical points II and III±.

In this parameter range, A+0 > 0 and A0 < 0. Finally, critical point IV, which reduces to II for σ = σ1 and to III± for σ = σ2 and is hence called a generalized cell [11], is not stable for any value of σ. Here, we use the term orbital stability of pattern formation to mean a family of solutions in the plane that may interchange with each other but do not grow or decay into a solution type from a different family. Such an interpretation depends upon the translational and rotational symmetries inherent to the amplitude-phase equations, these invariancies also limiting each equivalence class of critical points to a single member that must be considered explicitly. The same reasoning was implicit to our cataloguing of the critical points of the rhombic planform amplitude equations in the previous section.

We next offer a morphological interpretation of the potentially stable critical points catalogued above relative to the vegetative patterns under investigation. Then, critical points I and II represent the uniform homogeneous and supercritical striped states, respectively, described in both of the last two sections. Observing that to lowest order, the deviation function associated with critical points III± satisfies

equation image

where

equation image

and noting that the function g0(x, y) exhibits hexagonal symmetry [32], we can deduce that these critical points represent hexagonal lattices possessing individual hexagons with dark circular regions at their centres and light boundaries for III+ and with light circular regions at their centres and dark boundaries for III. The contour plots relevant to these critical points are depicted in Figures 10 and and1111 where Figure 10 is for III+ and Figure 11, for III. Hence, in what follows, we shall identify hexagonal close-packed distributions of vegetative spots or gaps with critical points III+ or III, respectively.

Figure 10.
A contour plot of Equation (48) for critical point III+ of Equation (46) with Λ = 1, α = 0.1, and β = 0.015.
Figure 11.
A contour plot of Equation (48) for critical point III of Equation (46) with Λ = 1, α = 0.2, and β = 0.01.

Having summarized those stability criteria and morphological identifications, we now return to our expressions for the Landau constants of Equations (39c) and (39d) given by Equations (18d), (44c), and (44d). First, we examine the signs of a1 + 4a2, a1 + a2, a0, and 2a2a1 versus α in Figures 12 and and1313 for Λ = 1. From the results of this examination, we observe that besides α1, α2, and αc defined in Equation (19) and Equation (28b) there exist the following other significant values of α

equation image

such that

equation image

equation image

equation image

equation image

These values of α are compiled in Table 4 for Λ = 0.825 and 1.

Figure 12.

Plots of a1 + 4a2 and a1 + a2 of Equations (18d) and (44d) versus α with Λ = 1 where the plot of a1 of Figure 3 is presented for the purpose of comparison.

Figure 13.

Plots of a0 and 2a2a1 of Equations (18d), (44c), and (44d) versus α with Λ = 1.

Table 4.

The α values of Equation (49).

In order to compare our theoretical predictions with numerical simulations and relevant observational evidence, we must represent the stability results of Table 3 graphically in the α − β plane. To do so, it is now necessary for us to generalize the approach which produced the marginal stability σc = 0 curve β = βc(α) of Equation (15) from Equation (14) in order that we may generate the analogous loci associated with σ = σj for j = −1, 1, and 2, respectively. That is, employing Equations (14), (18d), (44c), (44d), and (45), we solve

equation image

for β to obtain

equation image

respectively. Since all the quantities required for the identification of the vegetative patterns of Table 3 have been evaluated, we can represent graphically the regions corresponding to these patterns in the α − β plane of Figure 14 with Λ = 1, where the loci of Equations (15) and (50b) are denoted by σc and σj with j = −1, 1, and 2, respectively, in that figure which has been plotted for 0.0600 ≤ α ≤ 0.2108. Then we observe from Figure 14 in conjunction with Tables 3 and and44 that for a1 + 4a2 > 0 all (when 2a2a1 < 0) or part (when 2a2a1 > 0) of the region (σc, a1 > 0), where the one-dimensional analysis of Section 2 predicted striped vegetation patterns, is further divided into two subregions characterized by hexagonal patterns consisting of either vegetative spots (when a0 < 0) or gaps (when a0 > 0), respectively. In the overlap regions satisfying σ1 < σ < σ2 or

equation image

where stripes and spots (α1 < α < αc) or stripes and gaps (αc < α < α2) are predicted, the initial conditions determine which stable equilibrium pattern of each pair will be selected. In this context, we define

equation image

which from Figure 4 implies

equation image

There also exists a region of bistability corresponding to σ−1 < σ < 0 or

equation image

the homogeneous distribution being stable for σ < 0 or β > βc (α) and hexagons for σ−1 < σ < σ2. Given that σ−1 < 0 for a0 ≠ 0 or ααc = 0.16, the hexagons persisting in this overlap region would be subcritical in nature. Finally, to justify the truncation procedure inherent to the asymptotic representation of Equations (39c) and (39d), it is necessary that its coefficients satisfy the size constraint [31]

equation image

which, as can be seen in Figure 15, is valid in the parameter range of interest. Hence, we may conclude that such a truncation procedure is justified for our hexagonal planform weakly nonlinear stability analysis of Equation (4). Further, given the hexagonal close-packed nature of the distributions associated with III±, we have referred to them collectively as hexagons above.

Figure 14.
Stability diagram in the α − β plane for the propagator-inhibitor logistic equation (4) with Λ = 1 denoting the predicted vegetative patterns summarized in Table 3.
Figure 15.
Plots of (a1 + 4a2)2 and a0 of Figures 12 and and1313 versus α.

5. Comparisons, discussion, and conclusions

We are now ready to compare these theoretical predictions developed in Sections 3 and 4 with relevant numerical simulations and observational evidence as well as place them in the context of some recent pattern formation studies. We begin naturally by comparing our predictions with the numerical simulations of Lejeune and his co-workers: namely with those contained in Lejeune and Tlidi [14], Couteron and Lejeune [2], Lejeune et al. [15], and Lejeune et al. [16]. In particular, Lejeune et al. [16] numerically integrated Equation (3) on a square grid with periodic boundary conditions for Λ = 0.825 and L = 0.1 after constructing the corresponding bifurcation diagram of associated vegetative states as a function of μ to select appropriate values of this aridity control parameter favouring hexagonal pattern formation. These authors reported the following sequence of spatially periodic vegetation states obtained from that bifurcation diagram as aridity was increased: first, a pattern of spots of sparser vegetation or gaps appeared forming a hexagonal lattice, which then transformed into an alternation of stripes of sparser and thicker vegetation, and finally into a hexagonal pattern of vegetation spots separated by bare ground. Recalling that β = L2 and α is inversely related to μ by Equation (11b), we can see that the sequence of morphologies obtained from Figure 14 when the horizontal line

equation image

is traversed in the direction of increasing α through the patterned region

equation image

is in qualitative agreement with the sequence just catalogued. Here after Borckmans et al. [1], we employ the notation B/C to indicate bistability between patterns B and C. Next, Lejeune et al. [16] performed the numerical simulations for the fixed values of Λ and L as described above and generated vegetation patterns constituted of hexagonal gaps, alternating stripes, or hexagonal spots when μ was assigned the values of 0.965, 0.970, or 0.980, sequentially. Although Λ = 1 chosen in Sections 2–4 for ease of exposition is a special case since such a choice eliminates the quadratic term of f (ρ), it is a representative value of that parameter for our purposes. The qualitative behaviour of our stability analyses is the same for the other values of 0 < Λ < 1 defined by Equation (9) as it was for Λ = 1. This is demonstrated conclusively relevant to the simulation results of [16] by Table 4, which lists the special α values of Equation (49) for Λ = 0.825, which was the value employed during that simulation, as well as those for Λ = 1. To compare our theoretical predictions with these simulation results, we must first construct a stability diagram analogous to Figure 14 but for Λ = 0.825 rather than Λ = 1. Figure 16 represents such a diagram plotted in the α − β plane. Note in this context that αc as defined implicitly by Equation (49e) has the explicit formulation as a function of Λ given by

equation image

and hence relevant to Figures 14 and and1616

equation image

This critical value of density has been denoted by the vertical line α = αc in both these figures and separates the two types of vegetative hexagonal patterns from each other since spots can only occur for α < αc while gaps can only occur for α > αc. Then using Equation (11a) with Λ = 0.825, it is possible for us to convert Figure 16 into Figure 17, which represents the corresponding stability diagram plotted in the μ − β plane. Again note that the vertical line in this figure is given by

equation image

where the latter quantity corresponds to the αc of Equation (57) when Λ = 0.825. Finally, in order to compare the simulation results just described with our theoretical predictions, we now select β = L2 = 0.01 consistent with the value of L = 0.1 employed by Lejeune et al. [16]. Upon examination of Figure 17, we can deduce the predicted morphological sequence of vegetative patterns summarized in Table 5 as the horizontal transit curve β = 0.01 is traversed in the direction of increasing μ. Observe from Table 5 that the complete patterned region occurs for the interval μ [sm epsilon] (μc, μc+) where

equation image

are the two critical points of bifurcation as identified by Lejeune et al. [16] while the three μ values of 0.965, 0.970, or 0.980 chosen by the latter authors for their simulation purposes lie within the subintervals for which only gaps, only stripes, or only spots are stable, respectively. Hence our theoretical predictions are in quantitative agreement with these simulation results. Further, note from Figures 16 and and1717 that the term homogeneous appearing in the stable pattern column of the first two entries of Table 5 refers to uniform vegetative distributions of relatively high density while the same term appearing in the corresponding column of its last two entries refers to similar distributions of relatively low density. The results of our weakly nonlinear analyses are only strictly valid in the vicinity of the marginal stability curve. Often, however, in situations of this sort it is possible to extend these results to regions of the relevant parameter space substantially removed from the marginal curve. In order to determine the validity of such an extrapolation, it is standard operating procedure to perform numerical simulations of the original governing partial differential equation(s) in those extended regions and compare these simulation results with those theoretical predictions. This is precisely what was done by Firth and Scroggie [4] in their nonlinear optical pattern formation problem involving an atomic sodium vapour ring cavity and more recently by Golovin et al. [6] in their chemical Turing pattern formation problem of the Brusselator model with superdiffusion. In the latter instance, the theoretical predictions given in their Figure 3 were compared with the corresponding numerical simulations appearing in their Figure 5. Indeed Lejeune et al. [16] also compared their theoretical arid environmental vegetative pattern formation predictions and numerical simulations by means of their Table II which included the maximum and minimum phytomass density values for the simulated patterns of their Figure 5 and the corresponding values for the analytical bifurcation diagram given in their Figure 4 that was a plot of the stationary states of their amplitude equations versus μ. We reproduce those results in our Table 6. Considering the amplitude of the difference of these maximum and minimum values the table shows the relative deviation between the analytical and numerical values to be less than 10% in each of those cases, as pointed out by Lejeune et al. [16]. These cases, summarized earlier upon the introduction of Table 5, were for Λ = 0.825 and β = 0.01 with the aridity parameter μ = 0.965, 0.970, and 0.980 corresponding to gaps, stripes, and spots, respectively, all of which, as can be seen from Figure 17, lie relatively deep within the patterned region.

Figure 16.
Stability diagram in the α − β plane for the propagator-inhibitor logistic equation (4) with Λ = 0.825.
Figure 17.
Stability diagram in the μ − β plane corresponding to Figure 16.
Table 5.

Morphological stability predictions versus μ for β = 0.01 and Λ = 0.825.

Table 6.
Comparison of analytical and numerical ρ results for β = 0.01 and Λ = 0.825 from [16].

So far, we have been concerned with the bifurcation behaviour of Equation (4) when 0 < Λ ≤ 1 only for β in the range defined by Equation (55). We shall now examine what transpires should we consider Figure 14 in conjunction with the horizontal transit curve

equation image

Here βcrit = 0.015 corresponds to the local minimum value of the curve β = β1 (α) for α > αc. Then we obtain the sequence of morphologies

equation image

as α is increased along the transit curve of Equation (60). The simulation results of Lejeune and Tlidi [14] and Couteron and Lejeune [2] are in a parameter range that yields bifurcation behaviour of this generic sort. In particular, these authors numerically integrated Equation (3) on a square grid with periodic boundary conditions for Λ = 1.2 and L = 0.2 obtaining a hexagonal pattern of gaps for μ = 0.98, alternating stripes for μ = 1.0, and a hexagonal pattern of spots for μ = 1.01. To compare our theoretical predictions with these simulation results, we again construct stability diagrams analogous to Figures 16 and and1717 but for Λ = 1.2 in both the α − β plane of Figure 18 and the μ − β plane of Figure 19. In doing so, we first point out that this requires a relaxation of our original constraint of 0 ≤ Λ, μ ≤ 1 and second, that neither Lejeune and Tlidi [14] nor Couteron and Lejeune [2] constructed a bifurcation diagram relevant to their two-dimensional simulations although Lejeune et al. [15] did generate a so-called one-dimensional bifurcation diagram for this case (see below). Observe that Figures 16 and and1717 generated for Λ = 0.825 are in qualitative agreement with Figures 18 and and19,19, respectively, generated for Λ = 1.2. We also note that the vertical lines in those figures are given by

equation image

respectively. Next we select β = L2 = 0.04 consistent with the value of L = 0.2 employed by Lejeune and Tlidi [14] and Couteron and Lejeune [2]. Then upon examination of Figure 19, we can deduce the predicted morphological sequence of vegetative patterns summarized in Table 7 as the horizontal transit curve β = 0.04 is traversed in the direction of increasing μ.

Figure 18.
Stability diagram in the α − β plane for the propagator-inhibitor logistic equation (4) with Λ = 1.2.
Figure 19.

Stability diagram in the μ − β plane corresponding to Figure 18.

Table 7.

Morphological stability predictions versus μ for β = 0.04 and Λ = 1.2.

Observe from Table 7 that this morphological sequence is both in qualitative agreement with the sequence of Equation (61) and quantitatively compatible with the simulation results under examination once the μ values in that table have been rounded off to the two places of accuracy reported by Lejeune and Tlidi [14] and Couteron and Lejeune [2]. Upon examination of Figure 19, we can see that the parameter values of Table 7 lie even deeper within the patterned region of this figure than the parameter values of Table 6 did in Figure 17. In this context, Lejeune et al. [16] reported their simulation μ values to three places of accuracy. We note that for all the simulation results of Lejeune and his co-workers described in this section which generated striped patterns,

equation image

to the relevant degree of accuracy. In order for us to demonstrate that this occurrence was not coincidental, it is necessary that we consider our rhombic planform results of Section 3 in conjunction with the morphological stability predictions summarized in Table 3.

Towards this end, we begin by noting that, although these rhombic results of Section 3 were developed explicitly for the special case of Λ = 1, they are qualitatively valid for any other value of Λ satisfying the condition in Equation (3b) – e.g. specifically for Λ = 0.825 or 1.2. Consistent with Equation (52), we implicitly define α±2 by

equation image

and observe from Tables 3, ,5,5, and and77 and Figures 14, ,16,16, and and1818 that our hexagonal planform analysis predicted stripes as the only stable pattern for

equation image

From the results of Section 3, we can conclude that for any

equation image

all such stripes are unstable with respect to rhombic patterns having characteristic angle

equation image

In particular, for α in these intervals rhombic arrays of this sort are difficult to distinguish from hexagonal distributions since the allowable bands of their characteristic angles closely flank π/3 [congruent with] 1.047 (see Table 2 and Figure 9). Hence, if we substitute this type of rhombic pattern as the stable morphology for the intervals

equation image

then hexagonal or nearly hexagonal arrays can be anticipated for all ααc over the instability region. In this context, we note that all our pattern-formation figures have used the homogeneous solution value of α as the threshold to trigger the colour change from light to dark. Thus, all spatial regions characterized by ρα appear dark and those by ρ < α light. If in the α-intervals of Equation (64e) we adopt the protocol that αc represents this critical threshold instead, then in the high threshold interval where α < αc light regions predominate and hence the rhombic arrays in that interval would be characterized by dark rectangular patches with rounded corners completely surrounded by the light background or pseudo-‘spots’, whereas in the low threshold where α > αc dark regions predominate and hence those arrays in that interval would be characterized by light rectangular patches completely surrounded by the dark background or pseudo-‘gaps’. Hence, after Wollkind and Stephenson [32] using such a high threshold–low threshold type of argument, it is plausible to assume that these pseudo-‘spots’ are the preferred patterns for α2 < α < αc while pseudo-‘gaps’ are preferred for αc < α < α+2. Finally, from Figures 14 and and161619, in conjunction with the morphological stability results of Section 3, we can conclude that the occurrence of stripes requires α = αc, or equivalently μ = μc since such patterns are only stable for both planforms at this critical value of α or μ and when β < βc(αc). Here where our analysis predicts a stable parallel striped pattern the equivalence class for the critical point designated as II in the last section actually contains the three solutions [22]

equation image

These represent a family of stripes aligned parallel to the y axis as per our original identification plus two similar families of stripes making angles of ± 60° with them, for which stable coexistence with a member of either the original family or one another is impossible [22]. Then for α = αc (or equivalently μ = μc) and β < βc(αc) as initial conditions vary from point to point over a flat environment, such families of stripes can give rise to polygonal arcs the boundaries of which would appear quite random in orientation. Indeed, the simulation results classified as striped patterns by Lejeune and Tlidi [14] and Couteron and Lejeune [2] have the appearance of such curved elongated stripes or interconnected bicontinuous patterns in the relevant reproductions contained therein.

So far we have limited our discussion to analyses for which the wavenumber was restricted to the critical wavenumber alone. In order to investigate the consequences of considering other wavenumbers in the instability side band centred about this critical wavenumber, it would be necessary to convert the Landau-type amplitude ordinary differential equations in time to Ginzburg-Landau type partial differential equations by adding the appropriate spatial derivative terms to them. As reviewed in detail by Wollkind et al. [31], such an analysis, although beyond the scope of our present investigation always yields two additional instabilities besides those discussed above: namely, zig-zag and cross-band relevant to the interaction of oblique and perpendicular modes, respectively. Specifically, it should be pointed out that the stripes simulated by Lejeune et al. [16] for β = 0.01, Λ = 0.825, and μ = 0.97 were of the zig-zag variety.

Up until now, we have been concerned with a1 > 0 type behaviour for the simulation results of Lejeune and his co-workers when Λ = 1.2 and L = 0.2. Lejeune et al. [15] constructed a one-dimensional bifurcation diagram corresponding to this problem for Λ = 1.2 and L = 0.2. They found isolated vegetation patches as a simulation outcome occurring outside a patterned aridity interval when μ > μ* = 1.01 and interpreted such localized structures to be a spatial compromise between the periodic patchy vegetative and bare stable states. We plot a1 versus α in Figure 20 for Λ = 1.2 relevant to this situation and note that α = 0.1 is equivalent to μ = 1.01. Note that Figure 4 for Λ = 1 is in qualitative agreement with Figure 20 for Λ = 1.2. Here the latter case yields the values α1 = 0.1028 and α2 = 0.2821. Observe from Figure 19 that for the μ range included in Table 7 the transit curve β = 0.04 lies totally within the region where σ, a1 > 0. Should μ be extended just beyond the bounds of that table, however, a1 < 0 although σ > 0 (see Figures 1820). This behaviour is highly reminiscent of the morphological phase separation instabilities (σ > 0) in thin liquid films analysed by Tian and Wollkind [26]. The main result of their analysis was that a1 > 0 type patterns of this sort consisting of uniform distributions of nanodroplets, bicontinuous ridges, and polygonal cells separated by ultrathin flat films could occur sequentially as mean layer thickness was increased over a particular interval with rupture occurring outside of that interval where a1 < 0. This dewetting-type rupture occurred by isolated drop formation in relatively thin layers, but by isolated hole formation in thicker ones. Analogous behaviour was discovered by Golovin et al. [6] in their examination of chemical Turing pattern formation for the Brusselator problem with superdiffusion. Given the similarity of behaviour among all these phenomena, we conjecture that outside the bounds of Table 7 where σ > 0 and a1 < 0 there will be desertification characterized by the formation of isolated vegetation patches for μ > 1.0098 and by the formation of isolated patches of bare ground for μ < 0.9769.

Figure 20.

A plot of a1 of Equation (18d) versus α for Λ = 1.2.

It remains for us to compare our theoretical predictions with observational evidence. Striking periodic vegetation patterns covering widespread areas of arid or semi-arid regions of Africa, Australia, the Americas, and the Near East became noticeable through aerial photography nearly 70 years ago [2,13,21]. In this instance, an arid region refers to an environment that is characterized by an extended dry season where yearly potential evaporation exceeds yearly rainfall and water availability is a limiting factor of plant growth [16]. The patterns reported consisted of gaps (‘pearled or spotted bush’), labyrinths, stripes (‘tiger bush’), and spots (‘leopard bush’). These included bushy vegetation punctuated by bare spots in Niger, labyrinths of perennial grasses in Israel, striped patterns of bushy vegetation in Niger, and spots of trees or shrubs in Ivory Coast and French Guiana [21]. From our one-dimensional morphological interpretation of the equivalence class for critical point II of Equation (65) it follows that the patterns usually described by ecologists as tiger bush, lanes, groves, and bands should be associated with what we have termed parallel stripes while those described as labyrinths, mazes, and arcs should be associated with what we have termed bicontinuous patterns. Similarly leopard or pearled bush should be associated with critical points or III+ or III, respectively. We wish to re-examine further these hexagonal patterns in order to determine sufficient conditions under which they can give rise to regions characterized by bare ground in agreement with relevant observational evidence. Towards that end, we return to our hexagonal symmetry function g0 of Equation (48) and consider it in more detail. In particular, focusing our attention upon a single hexagonal cell we note that each such individual g0 cell has an elevated central region with a maximum elevation of 3 at its centre which is bounded by a level curve of zero elevation appearing as the nearly circular loci found in Figures 10 and and11.11. The peripheral portion of each cell exterior to that central region is depressed, with the hexagonal cellular boundary of variable depth ranging from − 3/2 at its vertices to − 1 at the midpoint of its edges. Specifically, at the centre of these cells, the density function characteristic of III has a value

equation image

while on its boundaries the density function characteristic of III+ ranges from

equation image

Even though (see Figures 14, ,16,16, and and1818)

equation image

we see that there is about the same tendency for these density values of Equation (66) to be equal to, or less than, zero, and hence represent bare ground in agreement with the actual appearance of leopard or pearled bush, respectively. Finally, we shall investigate the predicted wavelength of these vegetation patterns. From Equations (14b) and (20), we can deduce that

equation image

Observe that, for β [equivalent] β0, λc is a decreasing function of α. Then evaluating this function of Equation (67a) on the marginal stability curve β = βc(α) of Equation (15b) we obtain

equation image

which for the special case of Λ = 1 reduces to

equation image

these critical wavelengths being generated by the interplay between short-range facilitative and long-range competitive plant interactions. Specifically for α = α±0 of Equation (66) we can see from Equation (67c) that

equation image

and hence conclude that the vegetative distributions in leopard bush have a tendency to be more widely spaced than the bare patches which regularly punctuate the vegetation cover in pearled bush. Further the characteristic wavelength of tiger bush which corresponds to α = αc = 0.16 for this case satisfies

equation image

and is intermediate in value between the distances λ±c of Equation (67d). After Lejeune and Tlidi [15] employing a length scale of 5 m in this situation yields the associated dimensional wavelength

equation image

consistent with field observations.

We close by placing our theoretical predictions in the context of some recent vegetative pattern formation studies. Rietkerk et al. [20] developed an interaction-diffusion model system consisting of a set of three partial differential equations describing the spatio-temporal behaviour of plant density, soil water, and surface water, respectively. They performed numerical simulations on their model system with reflecting boundary conditions. Typical spatial patterns on flat ground generated in this manner were represented over a two-dimensional domain by Rietkerk et al. [20] for different amounts of rainfall R measured in mm/day. For R = 0.75 a hexagonal spotted pattern formed, changing into labyrinths for R = 1; while a hexagonal gap pattern was generated for R = 1.25. Von Hardenberg et al. [27] proposed a similar interaction–diffusion model system that only considered plant biomass and ground water but introduced more complicated plant growth and ground water flux terms. They studied this model by integrating their system numerically at different precipitation values ρ. For a plain (flat) landscape, von Hardenberg et al. [27] found that vegetation patterns of various forms evolved for p in a particular interval with hexagonal spots occurring at a relatively low p value, holes in a uniform coverage at a relatively high p value, and stripes (or labyrinths) at a p value intermediate to these. In both instances, the average plant biomass density increased with increasing R or p, respectively, although the scale of the patterns was in the range of metres or tens of metres for the Rietkerk et al. [20] model but only on the order of centimeters for the von Hardenberg et al. [27] model due to the fact that surface water diffusion which served as the driving mechanism for the former is much greater than ground water diffusion, the driving mechanism for the latter [21]. Klausmeier [9] devised a related interaction–advection–diffusion model for vegetation and surface water that included a spatial gradient effect for its water rather than a diffusive one. Numerical integration of that system yielded no periodic vegetative pattern formation for flat ground and stripes when the downhill flow of rainwater was allowed for hillsides. Sherratt [24] presented a detailed theoretical analysis of pattern formation in the Klausmeier [9] model, derived formulae for the wavelength and migration speed of the predicted patterns, and systematically investigated how these depended on model parameters. Observe that the morphological sequence results described above that were produced by numerical simulation of the Rietkerk et al. [29] and von Hardenberg et al. [27] model systems are in qualitative agreement with our theoretical predictions for the propagator-inhibitor Lefever et al. [13] model equation once one realizes that the aridity parameter of the latter is inversely related to the rainfall and precipitation parameters of those systems. Also, we note that our systematic nonlinear stability analyses of vegetative pattern formation for this model equation is presented in the same spirit as was Sherratt's [24] detailed theoretical analysis of the Klausmeier [9] model system.

We conclude this discussion with some additional comments about our analyses in relation to those of Lejeune and his co-workers. As pointed out in the sketch of the development of the propagator-inhibitor model equation (3a) presented in Section 1, Lefever et al. [13] introduced scaling laws at variance to those catalogued in Equation (3b). Specifically, they took

equation image

while retaining our other scales of Equation (3b). Then, the terms

equation image

while all the other terms in Equation (3a) remain of O(ε3). Within the framework of our nonlinear stability analyses, let us examine the validity of the scaling laws of Equation (3b) where they differ from Equation (68a): namely

equation image

The first condition of Equation (69a) is consistent with all asymptotic scalings for the Stuart–Watson method of weakly nonlinear stability theory [17] while, using the marginal stability curve of Equation (15b) since our analysis is expected to be asymptotically valid as β → βc(α), we find that

equation image

which is consistent with the second condition of Equation (69a). Observe as noted earlier that with the scaling laws inherent to Equation (69) all the terms included in Equation (3a) are now of O(ε3). Indeed, in their seminal work Lefever and Lejeune [12] actually retained some terms of higher order appearing in the intermediate steps of Equation (2) which is the reason we did not attempt to compare our theoretical predictions with the simulation results contained therein. Also, Lejeune and Tlidi [14] claimed that they never observed rhombic patterns in their numerical simulations because such patterns are unstable with respect to stripes and hexagons. Our rhombic planform results of Section 3 show this assertion to be in error with respect to stripes except when α = αc. Actually, their simulation results were presented for both small- and large-scale windows. For the extended domains the striped patterns in the small windows were converted to zig-zag patterns while the hexagonal patterns become much less regular. All of these developments are consistent with our earlier discussion concerning the potential bistability of various types of stripes for α = αc and between rhombic and hexagonal patterns for α ≠ αc. In fact, Golovin et al. [7] found bistability between squares and hexagons in certain regions of their gravity number-capillary number parameter space for a Bénard–Marangoni surface tension-driven convection problem with poorly conducting boundaries.

In summary, we have performed systematic rhombic- and hexagonal-planform weakly nonlinear stability analyses on a truncated fourth-order partial differential time-evolution logistic equation describing the total plant biomass per unit area scaled with its carrying capacity over an unbounded flat spatial domain where vegetative pattern formation is induced by the intrinsic aridity of the region. Our main result was the production of plots of closed-form conditions in the relevant parameter space under which different vegetative patterns could be found. These patterns consisted of rhombic arrays of rectangular patches, parallel stripes or labyrinth-like mazes, and hexagonal distributions of spots or gaps. For a fixed value of the facilitation-to-competition interaction ratio, Λ, the stability behaviour of those vegetative patterns depended on two positive parameters α and β, which represented the mean plant biomass density and the square of the facilitation-to-competition range ratio, respectively. We showed that if β was greater than a certain threshold value only a uniform homogeneous distribution was possible, while should that not be the case, then the vegetative pattern selected depended on the size of α as well. In particular, one-dimensional parallel striped, labyrinth, or bicontinuous patterns were only possible if α = αc and β < βc; spots if α < αc and β2 < β < β−1; and gaps if α > αc and β2 < β < β−1; while if α2 < α < αc or αc < α < α+2 and β < βc, rhombic patterns were stable for characteristic angles ϕ in two narrow bands flanking π/3. In addition, the same value of α, namely αc, playing a central role with respect to stripe formation in both rhombic- and hexagonal-planform anlayses served as a partial but independent check on these analyses. Further, for Λ held constant, α was determined by μ, the measure of aridity, with α a decreasing function of μ within the patterned range. Such vegetative patterning marked the transition between homogeneous Savannas or deserts when βcrit < β < βc and between isolated patches of bare ground or vegetation when β < βcrit. Hence this propagation–inhibition model was unified in the sense that it could be used to generate a wide variety of vegetative patterns which are observable in arid or semi-arid environments. We then compared these theoretical predictions with relevant numerical simulations and that observational evidence, finding consistency in the former case and very good agreement when parameter values were chosen appropriately in the latter case. Finally, we placed those results in the context of some recent nonlinear vegetation pattern formation studies. In conclusion, this aridity-driven vegetative distributions generation problem, involving a single evolution propagator–inhibitor logistic equation defined over a flat environment, is compatible with our long-range aim of employing the simplest reasonable natural science models that preserve the essential features of pattern formation and are still consistent with observation.

Acknowledgements

We would like to acknowledge Oliver Lejeune (as well as an anonymous reviewer) for his constructive comments that greatly contributed to this revised version of our manuscript.

References

[1] Borckmans P., Dewel G., De wit A., Walgraef D. In: Turing bifurcation and pattern selection, in Chemical Waves and Patterns. Kapral R., Showalter K., editors. Dordrecht, The Netherlands: Kluwer; 1995. pp. 323–363.
[2] Couteron P., Lejeune O. Periodic spotted patterns in semi-arid vegetation explained by a propagation-inhibiton model. J. Ecol. 2001;89:616–628.
[3] Cross M.C., Hohenberg P.C. Pattern formation outside of equilibrium. Rev. Mod. Phys. 1993;65:851–1112.
[4] Firth W.J., Scroggie A.J. Spontaneous pattern formation in passive optical systems: Two-level atoms. Europhys. Lett. 1994;26:521–526.
[5] Geddes J.B., Indik R.A., Moloney J.V., Firth W.J. Hexagons and squares in a passive nonlinear optical system. Phys. Rev. A. 1994;50:3471–3485. [PubMed]
[6] Golovin A.A., Matkowsky B.J., Volpert V.A. Turing pattern formation in the Brusselator model with superdiffusion. SIAM J. Appl. Math. 2008;69:251–272.
[7] Golovin A.A., Nepomnyashchy A.A., Pismen L.M. Pattern formation in large-scale Marangoni convection with deformable interface. Physica D. 1995;81:117–147.
[8] Graham M.D., Kevrekidis I.G., Asakura K., Lauterbach J., Krischer K., Rotermund H.-H., Ertl G. Effects of boundaries on pattern formation: Catalytic oxidation of CO on platinum. Science. 1994;264:80–82. [PubMed]
[9] Klausmeier C.A. Regular and irregular patterns in semiarid vegetation. Science. 1999;284:1826–1828. [PubMed]
[10] Koschmeider E.L. Bénard Cells and Taylor Vortices. Cambridge: Cambridge University Press; 1993.
[11] Kuske R., Matkowsky B.J. On roll, square, and hexagonal cellular flames. Eur. J. Appl. Math. 1994;5:65–93.
[12] Lefever R., Lejeune O. On the origin of tiger bush. Bull. Math. Biol. 1997;59:263–294.
[13] Lefever R., Lejeune O., Couteron P. A case study of tiger bush in sub-Saharan Sahel. In: Maini P.K., Othmer H.G., editors. Mathematical Models for Biological Pattern Formation. New York: Springer-Verlag; 2001. pp. 83–112.
[14] Lejeune O., Tlidi M. A model for the explanation of vegetation stripes (tiger bush) J. Veg. Sci. 1999;10:201–208.
[15] Lejeune O., Tlidi M., Couteron P. Localized vegetation patches: A self-organized response to resource scarcity. Phys. Rev. E. 2002;66:010901-1–010901-4. [PubMed]
[16] Lejeune O., Tlidi M., Lefever R. Vegetation spots and stripes: Dissipative structures in arid landscapes. Int. J. Quantum Chem. 2004;98:261–271.
[17] Matkowsky B.J. Nonlinear dynamic stability: A formal theory. SIAM J. Appl. Math. 1970;18:872–883.
[18] Newell A.C., Whitehead J.A. Finite bandwidth, finite amplitude convection. J. Fluid Mech. 1969;38:279–303.
[19] Pansuwan A., Rattanakul C., Lenbury Y., Wollkind D.J., Harrison L., Rajapakse I., Cooper K. Nonlinear stability analyses of pattern formation on solid surfaces during ion-sputtered erosion. Math. Comput. Model. 2005;41:939–964.
[20] Rietkerk M., Boerlijst M.C., van Langevelde F., HilleRisLambers R., van de Koppel J., Kumar L., Prins H.H.T., de Roos A.M. Self-organization of vegetaion in arid ecosystems. Am. Nat. 2002;160:524–530. [PubMed]
[21] Rietkerk M., Dekker S.C., de Ruiter P.C., van de Koppel J. Self-organized patchiness and catastrophic shift in ecosystems. Science. 2004;305:1926–1929. [PubMed]
[22] Segel L.A. The nonlinear interaction of a finite number of disturbances in a layer of fluid heated from below. J. Fluid Mech. 1965;21:359–384.
[23] Segel L.A. Distant side-walls cause slow amplitude modulation of cellular convection. J. Fluid Mech. 1969;38:203–224.
[24] Sherratt J.A. An analysis of vegetation stripe formation in semi-arid landscape. J. Math. Biol. 2005;51:183–197. [PubMed]
[25] Stuart J.T. On the nonlinear mechanics of wave disturbances in stable and unstable flows. 1. The basic behaviour in plane Poiseuille flow. J. Fluid Mech. 1960;9:353–370.
[26] Tian E.M., Wollkind D.J. A nonlinear stability analysis of pattern formation in thin liquid films. Interface. Free Bound. 2003;5:1–25.
[27] von Hardenberg J., Meron E., Shachak M., Zarmi Y. Diversity of vegetation patterns and desertification. Phys. Rev. Letts. 2001;87:198101-1–198101-4. [PubMed]
[28] Watson J. On the nonlinear mechanics of wave disturbances in stable and unstable flows. 2. The development of a solution for plane Poiseuille flow and for plane Couette flow. J. Fluid Mech. 1960;9:371–389.
[29] Wollkind D.J. Rhombic and hexagonal weakly nonlinear stability analyses: Theory and application. In: Debnath L., editor. Nonlinear Instability Analyses II. Southampton, England: WIP Press; 2001. pp. 221–272.
[30] Wollkind D.J., Alvarado F.J., Edmeade D.E. Nonlinear stability analyses of optical pattern formation in an atomic sodium vapor ring cavity. IMA J. Appl. Math. 2008;73:902–935.
[31] Wollkind D.J., Manoranjan V.S., Zhang L. Weakly nonlinear stability analyses of prototype reacton-diffusion model equations. SIAM Rev. 1994;36:176–214.
[32] Wollkind D.J., Stephenson L.E. Chemical Turing pattern formation analyses: Comparison of theory with experiment. SIAM J. Appl. Math. 2000;61:387–431.

Articles from Taylor & Francis Open Select are provided here courtesy of Taylor & Francis