Received 26 March 2015; accepted 25 May 2015; published 28 May 2015

ABSTRACT

A rhombic planform nonlinear cross-diffusive instability analysis is applied to a particular interaction-diffusion plant-ground water model system in an arid flat environment. This model contains a plant root suction effect as a cross-diffusion term in the ground water equation. In addition a threshold-dependent paradigm that differs from the usually employed implicit zero-threshold methodology is introduced to interpret stable rhombic patterns. These patterns are driven by root suction since the plant equation does not yield the required positive feedback necessary for the generation of standard Turing-type self-diffusive instabilities. The results of that analysis can be represented by plots in a root suction coefficient versus rainfall rate dimensionless parameter space. From those plots regions corresponding to bare ground and vegetative patterns consisting of isolated patches, rhombic arrays of pseudo spots or gaps separated by an intermediate rectangular state, and homogeneous distributions from low to high density may be identified in this parameter space. Then, a morphological sequence of stable vegetative states is produced upon traversing an experimentally-determined root suction characteristic curve as a function of rainfall through these regions. Finally, that predicted sequence along a rainfall gradient is compared with observational evidence relevant to the occurrence of leopard bush, pearled bush, or labyrinthine tiger bush vegetative patterns, used to motivate an aridity classification scheme, and placed in the context of some recent biological nonlinear pattern formation studies.

Keywords:

Leopard Bush, Tiger Bush, Pearled Bush, Root Suction, Nonlinear Stability Analyses, Threshold-Dependent Pattern Formation

1. Introduction

von Hardenberg et al. [1] devised a plant-ground water (sometimes called soil water) interaction-diffusion system to model self-organized vegetative pattern formation in arid environments (reviewed by Rietkerk et al. [2] ). Here the positive feedback for an activator consumer (e.g., plants) in the plant equation and the self-diffusivity advantage for an inhibitory limiting resource (e.g., ground water) provided the necessary conditions for the onset of Turing [3] pattern formation. von Hardenberg et al.’s [1] model also included the effect of plant root suction by adding a cross-diffusion term in their ground water equation. Rietkerk et al. [4] performed numerical simulations using reflecting boundary conditions on a similar interaction-diffusion model system consisting of three partial differential equations describing the spatio-temporal behavior of plant density, soil water, and surface water, respectively, but excluding the root suction cross-diffusion term in their soil water equation. That model had been carefully developed by HilleRisLambers et al. [5] for a flat semi-arid grazing system.

We wish to formulate an interaction-diffusion model system for plant biomass density (gm/m²) and ground (soil) water content (mm of depth), where two-dimensional coordinate system (m, m) and time (d), based upon the interaction terms of Rietkerk et al. [4] and the diffusion terms of von Hardenberg et al. [1] , defined on a flat unbounded arid environment. Toward that end, we first introduce the auxiliary dependent variable surface water content (mm of depth) and the coupled interaction- diffusion model system given by

(1.1)

where ([1] [4] )

(1.2)

(1.3)

with and for. Here maximum specific

water uptake rate by the plants, conversion of water uptake by the plants to plant growth, specific loss rate of plant density due to mortality, half-saturation soil water constant relevant to specific plant growth and water uptake, specific loss rate of soil water due to evaporation and drainage, rainfall rate, maximum specific water infiltration rate, plant saturation shaping constant for water infiltration , fraction of maximum specific water infiltration rate in the absence of plants, dispersal coefficient for plants, diffusion coefficient for ground water, diffusion coefficient for surface water, coefficient of plant root suction, and.

More than forty years ago, Keller and Segel [6] proposed the initiation of slime mold aggregation viewed as an instability in a landmark paper with that title. They formulated a mathematical reaction-diffusion model system for the aggregation of the cellular slime mold Dictyostelium discoideum involving four dependent variables: Namely, the density of this amoeba; the concentrations of the acrasin and acrasinase produced by it which mediate its aggregation; and the concentration of an intermediate complex formed by these chemicals in a reversible reaction. Keller and Segel [6] then simplified that model to a two-component system involving just the density of the amoeba and the concentration of the acrasin by making Haldane’s assumption that the complex was in chemical equilibrium and the additional assumption that the total concentration of the acrasinase enzyme in both its free and bound state was a constant. This simplified model included flux terms deduced from Fickian self-diffusion of its two components and a chemotaxis term for the amoeba generated by the introduction of cross-diffusion involving the acrasin gradient.

For the sake of model analysis, HilleRisLambers et al. [5] introduced a quasi steady-state approximation in (1.1) by taking. We next, after Keller and Segel’s [6] employment of Haldane’s assumption, simplify our model system (1.1)-(1.2)-(1.3) by assuming that the surface water is in hydrological equilibrium and making the quasi stationary approximation that

(1.4)

Such an assumption is dependent upon being relatively much greater than either or and O having a much faster time scale than either W or N [4] .

Then, introducing (1.4) into (1.1)-(1.2)-(1.3) by replacing the rate of infiltration term in with R, we obtain the final formulation of our interaction-diffusion plant-ground water model system

(1.5)

where

(1.6)

Our main purpose in doing so is to devise a model system of this sort that demonstrates root suction alone can generate the two-dimensional vegetative patterns (e.g., leopard bush, pearled bush, and labyrinthine tiger bush) occurring in arid flat environments as described by Rietkerk et al. [2] ).

2. Equilibrium Points and Their Linear Stability

There exist two equilibrium points of model system (1.5)

(2.1)

satisfying

(2.2)

given by

(2.3)

and

(2.4)

Note that (2.3) which exists for all parameter values corresponds to a bare ground or no vegetation situation while (2.4) which only exists for or, equivalently,

(2.5)

corresponds to a community equilibrium point or a state exhibiting a nontrivial homogenous vegetative distribution. In this context, we adopt the far-field condition that

(2.6)

We next wish to examine the linear stability behavior of these critical points and shall proceed sequentially. That is we begin with (2.3) by considering a separation of variables solution to system (1.5) and boundary condition (2.6) of the form

(2.7)

(2.8)

(2.9)

and find that

(2.10)

Upon examination of (2.10) it follows that identically and provided

(2.11)

Hence we can conclude that (2.11) represents the linear stability criterion for this bare ground equilibrium point. Since the community equilibrium point of (2.4) does not exist for these parameter values by virtue of (2.5), an exchange of stabilities between the two critical points of system (1.5) occurs at.

The primary focus of our research is on the stability behavior of this community equilibrium point. Now introducing the nondimensional variables and parameters

(2.12)

(2.13)

we transform system (1.5) and far-field condition (2.6) into

(2.14)

where,

(2.15)

and

(2.16)

Observe that the equilibrium point in question corresponds to

(2.17)

in our dimensionless formulation since

(2.18)

Here we are concerned with the stability of this solution to initially infinitesimal one-dimensional perturbations. Toward that end, we consider a reduced form of our basic system with; seek a separation of variables solution to it satisfying

(2.19)

where and are the wavenumber and growth rate of the linear perturbation quantities (i.e.) with for the constants and; and obtain

(2.20)

where

(2.21)

which are tabulated below for the relevant values of p and s.

Note that in (2.20) we have implicitly made use of the fact that. Upon substitution of these expansion coefficients from Table 1 into (2.20), we obtain the explicit secular equation for

(2.22)

where for those satisfying (2.5). Thus, since quadratics of the form of (2.22) have roots with negative real parts provided their coefficients are positive, we can conclude that the community equilibrium point is stable to linear homogeneous perturbations for which Further observe that in the absence of plant root suction () this equilibrium point would be linearly stable to heterogeneous perturbations, for which, as well. When root suction is considered, to guarantee the onset of such a cross-diffusive instability we require that

(2.23)

For fixed, , , and, the curve in the first quadrant of the plane is marginal since it serves as a boundary between the linearly stable region where and the unstable region of (2.23). This marginal stability curve has a minimum point at given by

Thus, when there exists a band of squared wavenumbers centered about corresponding to growing disturbances for which where represents the most dangerous mode of (2.22) while when there exists no such band (see Figure 1).

The restriction of in the expansions relevant to the weakly nonlinear stability analyses of the next three sections, allows us to consider, for fixed, , and, such that when, and when (see Figure 1). Hence, the locus with defined by (2.24), is a marginal stability curve in the plane, with, , and fixed. We plot that locus in Figure 2 for, , , and, corresponding to the typical parameter values ([1] [4] )

(2.25)

(2.26)

3. One-Dimensional Analysis: Stuart [7] -Watson [8] Nonlinear Stability Results

In the previous section we deduced the critical conditions for the occurrence of cross-diffusive instabilities when. To ascertain both the long-time behavior and spatial pattern of such growing perturbations we must consider the nonlinear terms in our basic equations. Defining the vector quantities

(3.1)

it is standard operating procedure to examine the weakly nonlinear stability of our community equilibrium point by letting [9]

(3.2)

where, and satisfies the amplitude equation

(3.3)

Then, the substitution of this solution into (2.14)-(2.15)-(2.16) and the expansion of its interaction terms in a Maclaurin series, yields six vector problems: The O(1) problem is satisfied identically given that are the components of the uniform homogeneous solution; the problem is the same as the linear one of

the previous section with, , , and; the solutions of the

two problems are catalogued in the Appendix; of the two problems, we are only concerned with the one proportional to containing the Landau constant. Using the usual Fredholm solvability condition on that problem, yields the following formula for the Landau constant as a function of, , , and

(3.4)

where the components of

(3.5)

are catalogued in the Appendix.

The stability behavior of the Landau equation (3.3) is dependent upon the sign of. Thus, to predict this behavior, we must analyze the formula for given by (3.4) with its interaction coefficients as listed in Table 1 and the components of (3.5) as defined in the Appendix. Hence we plot in Figure 3 for the same -domain and choice of parameter values as used in Figure 2 (note corresponds to).

From Figure 3 we observe that has a zero at such that

(3.6)

Given these conditions, we may conclude that there exist two supercritical cases when: Namely, for, the community equilibrium point is stable, resulting in a uniform homogeneous vegetative distribution; while for, there occurs a re-equilibration producing stationary parallel vegetative stripes

(3.7)

of amplitude and characteristic wavelength

(3.8)

in dimensionless and dimensional variables, respectively. These supercritical stripes are plotted in Figure 4(a) where regions of higher density appear dark and those of lower density appear light.

When that bifurcation is subcritical and there also exist two cases: Namely, again, for and, respectively. In general such subcritical behavior requires us to take into account higher order terms in our expansions which can cause the development of a uniform homogeneous vegetative distribution and isolated local vegetative patches (see the last section for a detailed discussion of this topic).

Finally we synthesize the one-dimensional pattern formation results of this and the previous section in the plane of Figure 4(b). We plot the cross-diffusive instability boundary curve of Figure 3 and the vertical lines, , in that figure. Then the regions, ,; and,; can be identified with bare ground, uniform homogeneous distributions of vegetation, and stationary striped vegetative patterns, respectively, in that parameter space. In this context, observe that the line serves as a horizontal asymptote for while

(3.9)

4. Two-Dimensional Analysis: Rhombic-Planform Nonlinear Stability Results

Wishing to refine our one-dimensional predictions summarized in Figure 4(b) and to investigate the possibility

of occurrence of the two-dimensional vegetative patterns mentioned earlier, we next consider a rhombic-plan- form solution of system (2.14) of the form [10]

(4.1)

where

(4.2)

with an analogous expansion for, such that

(4.3)

(4.4)

Here we are employing the notation for the coefficient of each term in (4.1) of the form

. Then substituting this rhombic-planform solution of (4.1) into system (2.14),

we again obtain a sequence of problems, each of which corresponds to one of these terms. Solving those problems we find that

(4.5)

while applying the same method of analysis, as employed for deducing (3.4), to the, , system yields the Fredholm-type solvability condition for the second rhombic-planform third-order Landau constant

(4.6)

where the components of

(4.7)

as well as the solutions for the relevant second-order systems are catalogued in the Appendix.

The rhombic-planform amplitude Equations (4.3)-(4.4) possess the following equivalence classes of critical points

(4.8)

Assuming that and investigating the stability of these critical points one finds that [10] :

(4.9)

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 characteristic angle [10] . In the next section we shall use these criteria to refine our one-dimensional predictions of Figure 4(b) relevant to the former states due to the presence of the latter. Toward that end, we examine the sign of. We first illustrate this procedure by defining the ratio of Landau constants [11]

(4.10)

and plotting that quantity versus in Figure 5 for a fixed value of, namely, and with the other parameters taking on their values of Figure 3. Here, there exists an interval of stable rhombic patterns, where or equivalently, given by

(4.11)

provided in addition that or.

Now, repeating the process used to produce Figure 5 but for other values of we find the same generic behavior in that there exists an interval of stable rhombic patterns for where when and summarize these results in Table 2. Then, we plot the ratio of Landau constants of (4.10) versus in Figure 6 for and the other parameters taking on their values of Figure 3.

Restricting ourselves to the interval of interest, we see from Figure 6 that there exists a band of stable rhombic patterns for where when. Observe from Figure 6 the intercept and symmetry properties

(4.12)

Here, these properties of (4.12) are a consequence of mode interference occurring exactly at and modal interchange, respectively [12] . Again, repeating the process used to produce Figure 6, but for other, we find the same generic behavior as for and summarize these results for selected values in Table 3.

Finally, we present a morphological interpretation of the stable rhombic patterns that can be associated with critical point V for the values of the characteristic angle relevant to Table 3. Then, to lowest order, the equilibrium vegetative pattern corresponding to that critical point satisfies

(4.13)

(4.14)

The three parts of Figure 7 are threshold contour plots of (4.14) for with threshold values of 1, 0, and, respectively. Hence from right to left the parts of this figure can be identified with what Wollkind and Stevenson [10] , Boonkorkuea et al. [13] , and Cangelosi et al. [14] called upper, zero, and lower threshold patterns, respectively. In this context note that. Traditionally, most pattern formation analyses of this type have used the dimensional homogeneous vegetative solution value of (2.4) as the threshold to trigger the color change from dark to light (see Figure 4(a)). Thus all spatial regions characterized by appear dark and those characterized by, light, where again dark regions correspond to high plant biomass density and light ones to low plant biomass density or bare ground. This is equivalent to our zero threshold case of Figure 7. Note from (2.4) and (2.13) that

(4.15)

For fixed values of the other parameters and satisfying (2.11) we may consider and to be increasing straight line functions of R alone given by and, respectively, where

(4.16)

is a dimensionless measure of the rainfall rate R. We now wish to select a particular and adopt the protocol that

(4.17)

represents this threshold instead. Then, when or, an upper threshold pattern of the type depicted in Figure 7 would occur while, when or, a lower threshold type would occur.

Given their appearance in Figure 7 we label these upper, zero, and lower threshold type rhombic vegetative arrays as pseudo spots, rectangles, and pseudo gaps and denote them by, , and, respectively, in what follows. We shall defer the specific choice for, our rationale for making that selection, and its morphological interpretations until the comparison of these results with some recent vegetative pattern formation studies included in the next section.

5. Morphological Interpretations and Comparisons.

As a prelude to the morphological interpretations to be developed in this section, we first demonstrate that our model system does not generate any hexagonal patterns. We do so by considering a hexagonal-planform expansion for of (2.12)-(2.13) with terms to first order given by [10]

(5.1)

where, for even permutations of (1, 2, 3),

(5.2)

(5.3)

with a similar expansion for. Segel [15] , who developed this six-disturbance methodology to analyze the Rayleigh-Bénard model of buoyancy-driven hexagonal-cellular convection, showed that the simplest way to deduce the additional Landau constants and contained in (5.2)-(5.3) was to consider its two-distur- bance reduction [16] obtained by setting

(5.4)

which reduces (5.1) and (5.2)-(5.3) to

(5.5)

with a similar expansion for where

(5.6)

(5.7)

See equation (4.4) of Wollkind and Stephenson [10] for the explicit higher-order terms in these expansions of

(5.6)-(5.7) proportional to with proportionality constants and.

Then upon substitution of (5.6)-(5.7) into (2.14) we find that (4.5) holds again while

(5.8)

for, and

(5.9)

as defined in Section 3. If we proceed in the same manner as we did with the rhombic-planform expansions of the last section, the Fredholm-type solvability conditions for and, respectively, yield

(5.10)

(5.11)

where the components of

(5.12)

for and 1 as well as the solutions for the other relevant second-order systems are catalogued in the Appendix. Observe, as pointed out by Wollkind and Stephenson [10] , that the expression for of (5.11) does not contain the component since its coefficient vanishes identically in this limit by virtue of (5.10), and

hence is often referred to as a free mode which is why that component is not catalogued in the Appendix.

The six-disturbance hexagonal-planform amplitude phase-Equations (5.2)-(5.3) have equivalence classes of critical points given by and

I:;

II:;

:;

IV: with.

Critical points I and II represent the uniform homogeneous and supercritical striped states described in the previous two sections; hexagonal close-packed arrays of spots and gaps, respectively; and IV, a genera-

lized cell that reduces to II for and to for [17] . Since the existence

and orbital stability of these critical points depends, in part, upon the signs of various combinations of the Landau constants of (5.2)-(5.3), we examine the signs of, , and by plotting those quantities as well as versus for, , and in Figure 8 and observe that they are all identically negative.

Although I is stable for, II is not stable for and does not exist for while IV is not stable for any set of parameter values [18] . Hence, we have demonstrated that this hexagonal- planform analysis does not yield any additional stable stationary heterogeneous vegetative patterns for our model system.

Having established this fact, let us return to the subject with which we ended the last section: Namely, the selection of the proper value to be assigned for and hence,. In order to motivate our specific choice, we

first summarize the simulation results of Rietkerk et al. [4] . Their two-dimensional numerical simulations of model system (1.1), with and its other parameters set at values consistent with (2.25)-(2.26), yielded close-packed vegetative patterns of spots or gaps depending upon whether R was less than or greater than 1 mm/d, respectively. Motivated by the desire to replicate this behavior and given the similarity in appearance between these two types of patterns and the right- or left-hand parts of Figure 7, respectively, as well as the fact that for the Rietkerk et al. [4] model is equivalent to (2.4) of our model, we then select

(5.13)

Thus, from our rhombic planform analysis of the previous section, we can make the prediction that patterns will occur for or and ones, for or when. Influenced by that resemblance in appearance just cited, we have referred to these periodic rhombic arrays of as vegetative pseudo spots or pseudo gaps, respectively. We now incorporate these two-dimensional rhombic- planform morphological stability results for, , and in the plane of Figure 9 and identify regions corresponding to the predicted vegetative patterning. Note in this context that (4.15)-(4.16), (4.17), and (5.13) imply

(5.14)

Figure 9 represents the two-dimensional refinement of our one-dimensional predictions of Figure 4(b). Observe that the occurrence of striped and rhombic patterns is mutually exclusive by virtue of stability criteria (4.9). Hence since rhombic patterns occur for all in the patterned region this precludes the occurrence of any striped patterns there. Note that this is consistent with our hexagonal planform prediction of critical point II being identically unstable. Thus our major two-dimensional refinement of those one-dimensional predictions is the replacement of the whole region of striped vegetative patterns () appearing in Figure 4(b) with a rhombic (V) one instead. Specifically, these are identified in Figure 9 where as rhombic arrays of pseudo spots () for and of pseudo gaps () for in accordance with our morphological threshold introduced in (5.13). Further, note that the region of Figure 9 which can be identified with a uniform homogeneous vegetative pattern varies from a relatively sparse distribution () for to a relatively dense one () for and hence have been designated by, respectively.

So far, with the implicit exception of the above-mentioned, we have been concerned with the morphological stability behavior of our model system for. We next consider its behavior for. Wollkind et al. [9] showed a particular partial differential evolution equation containing fourth-order spatial derivatives could be used to mimic pattern formation in reaction-diffusion systems. A comparison of the simulation results of Lejeune et al. [19] with the weakly nonlinear stability ones of Boonkorkuea et al. [13] for a strongly related evolu-

tion equation describing vegetative pattern formation in arid isotropic environments led to the conjecture that when localized structures would occur where characterized by isolated patches of vegetation at low densities that were a spatial compromise between the periodic patchy vegetation and bare ground stable states. Chen and Ward [20] found local structures occurring in conjunction with such subcriticality for the Gray- Scott reaction-diffusion chemical model system. Cangelosi et al. [14] employed the same argument to identify a region of their relevant parameter space with isolated clusters for a mussel-algae interaction-diffusion model system. The resulting morphological sequence deduced from that identification provided close agreement with mussel bed patterning observations both in the field and laboratory [21] [22] . Given the similarity of behavior among all these phenomena we conjecture with some confidence that isolated patches of vegetation would occur for where and identify that region graphically in Figure 9 using the designation.

Taking into account the next higher-order term in the expansions and the amplitude equation of (3.1)-(3.2)- (3.3) when by, in particular, considering

and calculating, we might enhance our understanding of the morphological stability of this system in the subcritical regime. Should that equation will have three equilibrium points: Namely, 0 and

. Then 0 is stable for;, for; and in the overlap region where both can be stable depending on initial conditions 0 is stable for and, for, while only exists in that bistability region but is not stable there [23] . Here

the potentially stable equilibrium points 0 and would correspond to and, respectively.

Determining the generalized marginal curve for our problem analogous to (2.24) but with we would find that [see Kealy and Wollkind [24] for the derivation procedure involved]

where. Plotting the marginal curve associated with,

in Figure 9, we could make the proper morphological identifications. These would differ from those already appearing in this figure only due to the presence of the overlap region where as well as patterns could now occur. Should, as was the case for a similar situation in Kealy and Wollkind [24] involving hexagonal patterns, the marginal curves would coincide or causing the overlap region virtually to disappear and resulting in the exact same identifications as the ones appearing in Figure 9.

Traditionally, morphological sequences of the sort referred to above have been generated from stability diagrams such as Figure 9 by traversing appropriate horizontal or vertical lines in that two-component parameter space (e.g., Cangelosi et al. [14] ; Kealy and Wollkind [24] ). A procedure of this sort is inherently dependent upon the implicit assumption that these two components are independent of each other. In the case of Figure 9, however, and, being nondimensional measures of rainfall and the coefficient of plant root suction, respectively, are actually related. To obtain the proper morphological sequence of vegetative states along a rainfall gradient predicted from Figure 9, it is first necessary for us to deduce that relationship. Toward this end, employing our basic definitions and the parameter values of (2.25)-(2.26), we find that

(5.15)

Note in (5.15) the units for as indicated below, are mm m⁴/(gm×d) consistent with being a dimensionless parameter. Adopting the root suction characteristic of Roose and Fowler [25] , we take

(5.16)

where

(5.17)

Here the relative water saturation in the soil while the parameters and m are determined from experimental data for different soils. To complete our formulation we let

(5.18)

where, specifically,

(5.19)

Upon recalling that

(5.20)

and substituting (5.20) into (5.18)-(5.19), yields

(5.21)

where, specifically,

(5.22)

Hence

(5.23)

Finally, selecting after Roose and Fowler [25] , and incorporating (5.23) and (5.16)-(5.17) into (5.15), we obtain the one-parameter family of root suction characteristic curves

(5.24)

where

(5.25)

We plot this curve with in Figure 9, where the assignment of that parameter has been made both for the purpose of definiteness and to be consistent with our silt loam soil choice for m [25] . Representing the stability boundary in that figure by, for ease of exposition, we observe that there exist two points of intersection between it and satisfying

(5.26)

where, specifically,

(5.27)

Here, from (5.20), these -values of (5.27) correspond to

(5.28)

respectively. In this context, we define

(5.29)

where. The morphological sequence of predicted stable vegetative states along a rainfall gradient obtained upon traversing the curve in the plane of Figure 9 is tabulated below. Note, in general, that

(5.30)

Thus, isolated patches are only predicted for transit curves of the form of (5.24)-(5.25) provided.

We next compare these theoretical predictions with relevant observational evidence. The relevant reported vegetative patterns [2] consist of spots (leopard bush) and gaps (pearled or spotted bush). After Boonkorkuea et al. [1] , we now associate our rhombic arrays of pseudo spots () or pseudo gaps () with these leopard or pearled bush vegetative patterns, respectively, and then investigate the predicted wavelength of those vegetative patterns. From (2.24) and (3.8), we can deduce that

(5.31)

Designating the’s associated with by, respectively, it follows from Table 4 that

(5.32)

Then, we can see from (5.31) and (5.32) that

(5.33)

and hence conclude that the vegetative distributions of spots in leopard bush have a tendency to be more widely spaced than the bare patches which regularly punctuate the vegetation cover in pearled bush [26] . Employing the length scale of (3.8) and the parameter values of (2.25)-(2.26), yields the associated dimensional wavelength relationships

(5.34)

consistent with field observations and in agreement with Boonkorkuea et al. [13] , who interestingly enough found an identical power law relationship to (5.31) between their pattern wavelength and plant biomass for the evolution equation of Lejeune et al. [19] .

It remains for us to compare our results with those from some recent biological nonlinear pattern formation studies. We begin with the work of Gowda et al. [27] . These authors examined the standard sequence of patterned states (gaps ® labyrinth ® spots) generated in a general activator-inhibitor reaction-diffusion system as a bifurcation parameter was varied and then applied their results to the particular von Hardenberg et al. [1] plant-ground water model as its precipitation parameter was decreased. They employed both numerical simulation and analytical weakly nonlinear hexagonal-planform bifurcation methods. Gowda et al. [27] found, for the default set of parameter values von Hardenberg et al. [1] used in their numerical integration, that, although the simulation method reproduced the latter’s standard sequence, the hexagonal planform analysis as in our problem failed to predict vegetative patterns. These calculations, in accordance with (2.25)-(2.26), were performed for. When those calculations were repeated for, they found that the same standard sequence of vegetative patterns was produced for both the simulation and weakly nonlinear stability methods with the transition between the two hexagonal states occurring exactly where the second-order Landau constant changed sign [13] . Gowda et al. [27] concluded that weakly nonlinear stability theory failed to produce the correct results in the first instance because the simulated morphological sequence of vegetative patterns occurred for large amplitudes as the precipitation parameter was decreased. Given our results, we wish to suggest another possible explanation for this discrepancy: Namely, that a rhombic-planform weakly nonlinear stability analysis might yield a predicted morphological sequence involving pseudo gaps ® rectangles ® pseudo spots as the precipitation parameter was decreased in this case. As supporting evidence for such a conjecture we offer the following observation. The gap- and spot-type simulated patterns for appearing in both Gowda et al. [27] and von Hardenberg et al. [1] were much less regular in nature than were the corresponding simulated hexagonal patterns for appearing in Gowda et al. [27] . The simulated transition states between these two types of patterns were also different consisting of labyrinths in the former instance but of parallel stripes in the latter case. Since for each value of α we predict multistable rhombic states with an interval of characteristic angles (see Figure 6 and Table 3) and as initial conditions vary point by point over a flat environment these states can be selected quite randomly, it is possible to generate simulated patterns resembling those appearing in von Hardenberg et al. [1] from families of pseudo gaps and pseudo spots including labyrinths from families of rectangles. Incidentally these labyrinthine patterns have been associated with certain types of so-called tiger bush or banded thicket vegetative distributions found in arid or semiarid flat environments [28] . In this context, the numerically simulated two-dimensional patterns of Rietkerk et al. [4] used earlier to motivate our choice for also bore a strong resemblance to those of von Hardenberg et al. [1] including a labyrinthine transition state at. Unlike the von Hardenberg et al. [1] model ours is extremely robust to variations in. We performed additional rhombic and hexagonal planform nonlinear stability analyses on system (2.14) and found identical qualitative behavior for all.

We end this phase of our discussion by restating von Hardenberg et al.’s [1] claim that the power of model systems such as ours of (1.5) is their predicted sequence of stable states along a rainfall gradient such as the one summarized in Table 4 can be used to motivate an aridity classification scheme which is characterized by the three rainfall thresholds

(5.35)

Here we are employing the notation of von Hardenberg et al. [1] for these dimensionless rainfall (precipitation) rate thresholds and use them to introduce the following possible aridity classes based upon the inherent vegetative states of our system:

・ Dry-subhumid (): The only vegetative state the system supports corresponds to a dense homogeneous distribution.

・ Semiarid (): The only vegetative state the system supports corresponds to pseudo gaps of low threshold type.

・ Arid (): The only vegetative states the system supports correspond to either pseudo spots of high threshold type or isolated patches.

・ Hyperarid (): The only possible stable states the system supports correspond to either a sparse homogeneous vegetative distribution or bare ground.

As pointed out by von Hardenberg et al. [1] the advantage of the proposed aridity classification scheme pertains to the information it contains about dynamical aspects of drylands. Regions whose aridity classes imply the occurrence of upper threshold vegetative patterns, isolated patches, or a sparse homogeneous distribution are vulnerable to desertification. The mere knowledge of that threat, however, allows land managers to reverse this process for those regions by implementing crust disturbance, seed augmentation, or irrigation strategies. Meron et al. [29] suggested a cycling mechanism between plants and water to account for the formation of bare patches characteristic of vegetative patterning along such a precipitation gradient. Note that a process of this sort occurs in all directions for two-dimensional vegetative patterns (pseudo spots or gaps, rectangles, and isolated patches) but only in two directions for one-dimensional ones (stripes).

We close with a more detailed commentary on the role played by cross-diffusion in generating pattern formation instabilities for our two-component model system. Given that, our system violates the activator positive feedback necessary condition for the occurrence of a Turing self-diffusive instability which requires. Hence the cross-diffusive effect of plant root suction on ground water generates our instability since as noted earlier if its community equilibrium point would be identically linearly stable. Indeed, the other requirement of for a Turing self-diffusive instability to occur might also be violated should, and a cross-diffusive instability of this type could still be generated although, in our actual parameter range, it is not violated. Recently, Stancevic et al. [30] considered a reaction-chemotaxis-diffusion three-component in-host viral dynamics model system for the concentrations of uninfected or infected cells and the virus. They found that the cross-diffusive effect of chemotaxis toward the infected cells by the uninfected ones generated their pattern formation instability in a similar manner as for our two-component system. Since the community equilibrium point of their system was linearly stable in the absence of diffusion and chemotaxis, Stancevic et al. [30] referred to this as a Turing instability. To distinguish between these two cases, we shall refer to ours as a Turing cross-diffusive instability instead. The von Hardenberg et al. [1] two-component nondimensional model system also included the cross-diffusive effect of plant root suction on ground water. Since that system’s interaction terms satisfied the activator positive feedback condition for its community equilibrium point while in their default set of parameter values, the presence of this cross-diffusive effect mediated rather than generated their Turing self-diffusive instability. In particular, von Hardenberg et al. [1] took the coefficient of that term in this default set. Upon inspection of (2.14) we can see that this coefficient is related to our parameters by

(5.36)

Thus that assignment would yield the root suction characteristic curve

(5.37)

which, as a decreasing function of is in qualitative accord with our formulation of (5.24). The three-compo- nent model systems of Rietkerk et al. [4] and HilleRisLambers et al. [5] by explicitly including surface water were able to generate Turing instabilities where none would have occurred for our simplified two-component version of that model without root suction. Finally, Wang et al. [31] conducted a definitive analysis of Turing instabilities for their predator-prey model system by including both self- and cross-diffusion terms in the prey and predator equations and performing weakly nonlinear hexagonal-planform bifurcations and numerical simulations on its community equilibrium point. Since the Allee positive feedback effect for the activator prey and the self-diffusivity advantage for the inhibitory predator were satisfied for their specific model, this was an investigation of cross-diffusion mediated rather than generated instabilities.

6. Conclusions

In summary, we formulated an interaction-diffusion plant-ground water model system in an arid flat environment. That system basically was formed from two existing models by coupling a simplified version of the interaction terms of one of those systems with the exact diffusion terms of the other. These terms included a cross- diffusion effect in the ground water equation due to plant root suction and a nonautocatalytic effect in the plant equation that precluded the occurrence of Turing self-diffusive instabilities. We then performed a rhombic planform nonlinear cross-diffusion instability analysis on that system and found an interval of characteristic angles for which stable rhombic patterns occurred. Defining a critical plant biomass threshold to interpret such rhombic arrays, those patterns were of an upper, zero, and lower threshold type which we labeled as pseudo spots, rectangles, and pseudo gaps, respectively. Our main result could be plotted in a coefficient of root suction versus a rate of rainfall parameter space. Traversing that parameter space along an experimentally determined root suction characteristic curve as a function of rainfall, we produced a predicted morphological sequence of vegetative patterns consisting of bare ground, isolated patches, these rhombic arrays, and uniform homogeneous distributions from low to high density. Then that predicted sequence of stable states along a rainfall gradient was shown to be in good qualitative and quantitative agreement with observations involving the occurrence of leopard bush, labyrinthine tiger bush, and pearled bush in arid flat environments including the wavelength behavior of such patterns. Finally, we introduced an aridity classification scheme, with classes based upon the inherent vegetative patterns included in that predicted morphological sequence along a rainfall gradient, which could be used to anticipate desertification and subsequently to implement land management reversibility strategies. Implicit to our continuum formulation were the assumptions that the pattern wavelength was relatively large when compared with the mean coverage diameter of an individual plant but quite small when compared with the territorial length scale characteristic of the arid environment which allowed us to have considered our interaction-diffusion equations on an unbounded spatial domain [32] .

We conclude by noting that although these results of our weakly nonlinear stability analyses are asymptotically valid only in the neighborhood of the marginal stability curve and the rhombic vegetative arrays along our specific coefficient of root suction characteristic curve occurred in this region, numerical simulations of pattern formation for several reaction-diffusion systems and evolution model equations have shown that such theoretical predictions can be extended to regions of the relevant parameter space relatively far from the marginal curve [13] [33] . Given that a weakly nonlinear hexagonal planform analysis of our system did not predict any stable patterns, we also offered an alternative explanation involving these stable rhombic patterns for why numerical simulations analyses of similar systems have in the past yielded periodic pattern formation over a parameter range where theoretical weakly nonlinear hexagonal ones did not. We finish by reiterating for the purpose of emphasis the fact that all of our pattern formation results for this model have been generated by the cross-diffusion process of root suction as opposed to the mediating effect it has often had on self-diffusion generated Turing patterns.

References

Appendix

In this appendix we catalogue the nonhomogeneous terms and second-order solutions relevant to our expansions of (3.1)-(3.2)-(3.3), (4.1)-(4.2)-(4.3)-(4.4), and (5.4)-(5.5)-(5.6)-(5.7).

For (3.1)-(3.2)-(3.3) we have:

(5.38)

where

(5.39)

with

(5.40)

and

(5.41)

For (4.1)-(4.2)-(4.3)-(4.4) we have:

(5.42)

(5.43)

(5.44)

where

(5.45)

for; and

(5.46)

For (5.4)-(5.5)-(5.6)-(5.7) we have:

(5.47)

(5.48)

(5.49)

(5.50)

where

(5.51)

(5.52)

and

(5.53)