Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments

Accepted Manuscript

Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments Giancarlo Consolo, Carmela Curro, ` Giovanna Valenti PII: DOI: Reference:

S0307-904X(16)30636-9 10.1016/j.apm.2016.11.031 APM 11457

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

9 February 2016 22 August 2016 30 November 2016

Please cite this article as: Giancarlo Consolo, Carmela Curro, ` Giovanna Valenti, Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.11.031

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • vegetation patterns in semiarid areas are analyzed via a one-dimensional hyperbolic model • the hyperbolicity of the model brings some mathematical and physical advantages

CR IP T

• pattern features (wavelength and migration speed) can be modulated almost independently

• the characteristic velocity offers an additional degree of freedom to model complex vegetation dynamics

AC

CE

PT

ED

M

AN US

• the hyperbolicity of the model allows the search of exact solutions

1

ACCEPTED MANUSCRIPT

CR IP T

Pattern formation and modulation in a hyperbolic vegetation model for semiarid environments

Giancarlo Consolo1∗ , Carmela Curr`o1 and Giovanna Valenti2 Department of Mathematical, Computer, Physical and Earth Sciences, University of Messina, V.le F. Stagno D’Alcontres 31, I-98166 Messina, Italy. 2 Department of Engineering, University of Messina, C. Di Dio, I-98166 Messina, Italy. ∗ corresponding author Email: [email protected] Phone: +39.090.3977556 Fax: +39.090.3974170

AN US

1

M

Abstract

AC

CE

PT

ED

This paper presents analytical and numerical investigations of pattern formation and modulation in a one-dimensional hyperbolic extension of the Klausmeier reaction-diffusion-advection vegetation model for semiarid environments. The hyperbolicity of the model is shown to provide an additional degree of freedom which can be used to modify the dynamical response of the system away from the parabolic limit as well as to reproduce qualitatively the experimental data. Results of linear stability analysis reveal the possibility to modulate, almost separately, the main pattern features (wavelength and uphill migration speed) by acting on two independent parameters: the slope steepness and the characteristic velocity. Moreover, a possible justification of the mechanism underlying the uphill migration of the patterns is suggested in terms of the phase shift between the spatial distributions of vegetation biomass and water. The hyperbolicity of the system has also facilitated the search for a class of exact solutions, via the differential constraint method, for a hyperbolic vegetation model whose kinetic terms are more general than Klausmeier’s ones and whose applications are potentially broad. Keywords: hyperbolic model, patterned solutions, travelling waves, Hopf and Turing bifurcations, Wave bifurcation, differential constraints method.

1

Introduction

Highly organized vegetation patterns are typical of semi-arid environments and manifest themselves as bounded vegetated regions separated by bare ground areas. Multiple reasons justify the efforts made to describe the processes underlying the formation of such structures. First of all, vegetation patterns can be found in several landscapes around the world and exhibit either fascinating regular spatial configurations (such as stripes) or appear in irregular forms (labirinthic, random spots, etc) [1]-[7]. Moreover, it is certainly of relevant interest to investigate how a landscape manages its vital resources, especially when their availability is limited. It is likewise important to analyze the susceptibility of a given ecosystem to undergo sudden shifts towards the desert state as a result, for instance, of anthropogenic disturbances, climate changes, drought or pollution [8]. Studies on vegetation bands have been also stimulated by the observations concerning the process of upslope migration. However, despite this phenomenon has been detected and measured by several authors, its origin remains still controversial. Indeed, while the mechanism of band movement at fine scale is commonly associated with a non-uniform distribution of water and seeds, it is not fully clear

2

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

whether such a process arises from different local slopes or moisture levels observed at the uphill and the downhill sides of the pattern [13]-[19]. All these motivations have lead to the development of several mathematical models which aim to describe the complex competition between vegetation for the limited availability of water (e.g. see [9]-[12], to cite a few). Among these models, the one proposed by Klausmeier [9] has been widely and efficiently used in literature as it predicts the formation and the uphill migration of vegetation patterns along slopes. Despite several generalized versions of the Klausmeier model have been proposed in literature, here we focus our attention on its original formulation [9],[13]. In this manuscript we consider a hyperbolic reaction-diffusion-advection vegetation model for semi-arid environments which has been developed within the framework of Extended Thermodynamics [20],[21]. In particular, we aim to show that the hyperbolicity of the model, apart from removing the unphysical feature of infinite propagation speed in diffusion [22]-[29], entails some mathematical and physical consequences. First, it provides an additional degree of freedom which can be used to capture the complex vegetation dynamics occurring close to and away from the parabolic limit. In fact, while the results of previous papers [20],[24]-[26],[28]-[29] were obtained through the assumption of very small relaxation times, i.e. close to the parabolic limit, none of these works has tried to estimate that limit and has explored the dynamical response of the system beyond it. The current investigations are an attempt to fill this gap. In particular, through a linear stability analysis, we try to estimate a threshold value for the relaxation time that represents the parabolic limit and we show that the dynamical behavior of the hyperbolic system significantly deviates from the parabolic one [13]-[18] when the relaxation time is progressively increased from that threshold. Second, it allows to modulate, almost separately, the main pattern features, wavelength and uphill migration speed, by acting on two independent parameters without altering the equilibrium states. This feature makes the hyperbolic model more suitable to reproduce the behavior observed experimentally. Third, it also facilitates the search for a class of exact solutions, via the differential constraint method, for a more general hyperbolic vegetation model which could potentially capture a large variety of vegetation dynamics, not necessarily predicted by the Klausmeier model. Bearing in mind these goals, the paper has been organized as follows. In Section 2 we recall the one-dimensional hyperbolic reaction-diffusion-advection model proposed in [20] which predicts the formation and the uphill migration of vegetation stripes. Furthermore, we introduce the ODE system which describes the dynamics of travelling waves admitted by the PDE model in point. Finally, we classify the equilibrium states as a function of the rainfall value and introduce the quantities which will be used in the subsequent sections. In Section 3 we investigate the dynamical response of the steady states under small homogeneous and non-homogeneous perturbations. We characterize the bifurcation diagrams over the parameter space and put particular emphasis on the occurrence of diffusion-advection-driven (oscillatory Turing) instabilities giving rise to patterned solutions. In the same Section we also carry out numerical simulations in the nonlinear regime and test the capability of our model to reproduce experimental data extracted from literature. In Section 4, further information on pattern dynamics are extracted through the inspection of periodic travelling wave solutions arising from Hopf bifurcations for the system of ODEs. Finally, in Section 5 the search for a class of exact solutions is addressed via the differential constraints method. In particular, after a brief overview of the technique, we look for those response functions and exact solutions which bring an ecological meaning and that mimic, at qualitative level, the dynamics observed in the Klausmeier model.

2

Model description

Let us consider a generalized version of the Klausmeier model where the competition between surface water and vegetation in a semiarid environment is described through a one-dimensional hyperbolic model [9],[20].

3

ACCEPTED MANUSCRIPT

The governing system of PDEs in dimensionless form reads    ut + Jx = f (u, w)  wt − νwx = g(u, w)    Jt + γ 0 (u)ux = −γ 0 (u)J

which, in the vector form, becomes

Ut + M(U)Ux = N(U) 

 u U =  w , J



0 M= 0 γ0

0 −ν 0

 1 0 , 0



 f (u, w) N =  g(u, w)  −γ 0 J

(2)

(3)

AN US

being

CR IP T

(1)

M

where u(x, t) and w(x, t) represent, respectively, the vegetation biomass and water at time t and position x (the positive direction being uphill), J(x, t) is the dissipative flux, ν is the slope gradient (steepness), dγ >0 the subscript stands for the partial derivative with respect to the indicated variable and γ 0 (u) = du 1 is a response function related to the relaxation time τ through τ = 0 . We note that, in the limit case γ τ → 0, our hyperbolic model reduces to the corresponding parabolic one [13]-[18] where the dissipative flux obeys a Fick’s law J = −ux and, consequently, the vegetation biomass follows the evolution equation ut = f (u, w) + uxx . The source terms appearing in (3) can be conveniently expressed as

ED

f (u, w) = f1 (u)w + f2 (u), (4)

g (u, w) = g1 (u)w + g2 (u)

PT

and recalls those of the Klausmeier model [9] for g1 (u) = − (1 + f1 (u)) ,

CE

with

f1 (u) = u2 ,

g2 (u) = A,



f1 (u) u

0

>0

f2 (u) = −Bu.

(5)

(6)

AC

In (5),(6), the parameters A and B are representative of rainfall and plant loss, respectively. In particular, in order to use realistic values of parameters for plant growth in semi-arid environments, in accordance with previous works [15],[16],[18], the range of the plant loss is hereafter restricted to B < 2. From the mathematical point of view, the system here investigated is strictly hyperbolic since the characteristic velocities λ are real and distinct √ λ1 = −ν, λ2,3 = ∓ γ 0 . (7)

The left and right eigenvectors corresponding to these eigenvalues are respectively given by:    √  l1 = 0 1 0 , l2,3 = ∓ γ 0 0 1    1  ∓√ 0 0 1 γ  d1 =  1  , d2,3 =  . 0 2 0 1

(8a) (8b)

The model (2)-(4) can be also recast as a system of ODEs by looking for solutions of the form U(x, t) = U(z), with z = x − st, representing waves which, during the motion, travel with a constant migration speed s > 0 and preserve the spatial profile of the wavefront. Such a ODE system reads 4

ACCEPTED MANUSCRIPT

dU = N. (9) dz The spatially homogeneous equilibrium states U∗ ≡ (u∗ , w∗ , J ∗ ) admitted by both PDEs and ODEs systems depend upon the value of the rainfall A, namely

CR IP T

(M −s I)

for A < 2B ⇒ one equilibrium: U∗0 = (0, A, 0)

for A = 2B ⇒ two equilibria: U∗0 , U∗min = (1, B, 0)     B B ∗ ∗ ∗ for A > 2B ⇒ three equilibria: U0 , U1 = u1 , , 0 , U2 = u2 , , 0 u1 u2 where

√ A2 − 4B 2 A + A2 − 4B 2 u1 = < 1 < u2 = . (11) 2B 2B From the ecological viewpoint, the fixed point U∗0 corresponds to the desert state whereas the other equilibria represent uniformly vegetated ecosystems.

AN US

A−

√

(10)

3

System response to homogeneous and non-homogeneous perturbations

M

The local stability character of the equilibria (10) under small non-homogeneous perturbations is investigated by linearizing the governing system (2)-(4) through

which leads to



ED

˜ exp (ωt + i ξx) U = U∗ + U fu∗ − ω  gu∗ i γ 0∗ ξ

fw∗ ∗ gw − ω + i νξ 0

 −iξ  0 ω + γ 0∗

   u ˜ 0 w ˜  =  0 . 0 J˜

(12)

(13)

AC

with

CE

PT

In (13) the asterisk denotes that the functions are evaluated at U∗ while ω and ξ represent the growth factor and the wavenumber, respectively. The response of the system to spatially homogeneous perturbations can be obviously evaluated by setting ξ = 0. It is easy to verify that the growth factor ω satisfies the following characteristic equation

being

ω 3 + D1 ω 2 + D2 ω + D3 = 0

b 1 − i νξ, D1 = D b 2 + γ 0∗ ξ 2 + i νξ (fu∗ − γ 0∗ ) , D2 = D
∗ 2 b D3 = D3 − γ 0∗ gw ξ + i νξγ 0∗ fu∗ − ξ 2 , ∗ b 1 = γ 0∗ − (fu∗ + gw D ), ∗ ∗ ∗ ∗ ∗ b D2 = fu gw − fw gu − γ 0∗ (fu∗ + gw ), 0∗ ∗ ∗ ∗ ∗ b D3 = γ (fu gw − fw gu ) .

Let us now analyze in detail the stability character of each steady state.

5

(14)

(15)

(16)

ACCEPTED MANUSCRIPT

3.1

Equilibrium U∗0

By evaluating the partial derivatives of the source terms with respect to the field variables at U∗0 we obtain fu (U∗0 ) = −B < 0, gw (U∗0 ) = −1 < 0, fw (U∗0 ) = gu (U∗0 ) = 0 (17)

CR IP T

so that the equation (14) can be easily factorized and admits the explicit solutions given by   q 1 2 0∗ 0∗ 0∗ 2 ω1 = (−1 + i νξ) ω2,3 = − B + γ ∓ (B − γ ) − 4γ ξ . 2

(18)

Consequently, U∗0 is always a stable node or focus under homogeneous or nonhomogeneous perturbations, respectively. From the ecological point of view, this means that, when A < 2B, the rainfall level is too low with respect to the vegetation loss and whatever initial vegetation biomass will simply die out.

Equilibrium U∗min

AN US

3.2

In the particular situation A = 2B the system also admits the equilibrium U∗min characterized by fu (U∗min ) = B > 0, gw (U∗min ) = −2 < 0,

fw (U∗min ) = 1 > 0,

gu (U∗min ) = −2B < 0

(19)

The linear stability analysis for this steady state is carried out by distinguishing the behavior under homogeneous (ξ = 0) and nonhomogeneous (ξ 6= 0) perturbations. In the former case, the characteristic equation reduces to

M

ω (ω + γ 0∗ ) (ω − (B − 2)) = 0

0∗

(20)

U∗min

PT

ED

so that, being γ > 0 and B < 2, the equilibrium is always stable (not asymptotically) in response to homogeneous perturbations. Unfortunately, for non-homogeneous perturbations, the characteristic equation (14) is in the form of a third-order equation with complex coefficients which, except for the limiting case of a flat ground, cannot be easily managed. In fact, for ν = 0, (14) becomes an equation with real coefficients and the linear stability analysis can be carried out by using the Routh-Hurwitz criterion, namely Re ω < 0 ∀ω ⇐⇒ D1 > 0 D3 > 0 D1 D2 − D3 > 0 ∀ξ. (21)

3.3

CE

Evaluating these terms at U∗min , the conditions D1 > 0 and D3 > 0 are always satisfied whereas the last one is fulfilled for γ 0∗ > B. Therefore, in the case of a flat ground, U∗min is a stable focus under the restriction γ 0∗ > B. In the most general case ν 6= 0, significant information on the linear stability character of this steady state could be extracted through numerical investigations.

Equilibria U∗1,2

AC

For A > 2B the system admits, apart from the desert state, two additional equilibria U∗1,2 characterized by
fu (U∗1,2 ) = B > 0, gw (U∗1,2 ) = − 1 + u21,2 < 0, fw (U∗1,2 ) = u21,2 > 0, gu (U∗1,2 ) = −2B < 0 (22) In the case of homogeneous perturbations, the characteristic equation (14) reduces to h    2 i 2 (ω + γ 0∗ ) ω 2 − B − u1,2 − 1 ω + B u1,2 − 1 = 0

(23)

so that, taking into account that u1 < 1 < u2 , we may conclude that U∗1 is always unstable (a saddle point) whereas U∗2 is a stable node or focus depending on the values of A and B. In the case of non-homogeneous perturbations and a flat ground (ν = 0), the Routh-Hourwitz criterion applied to (14) establishes that the state U∗1 is always unstable whereas U∗2 is stable if γ 0∗ > B. In the general case ν 6= 0, since we are mostly interested in the occurrence of diffusion-advectiondriven instability giving rise to vegetation patterns, hereafter we will restrict our attention to the dynamics observed at U∗2 only. To this aim, by looking for roots of (14) with null real part, the following sixth-order equation is obtained 6

ACCEPTED MANUSCRIPT

a3 ξ 6 + a2 ξ 4 + a1 ξ 2 + a0 = 0

(24)

with

CR IP T

∗ b a3 = γ 0∗ (γ 0∗ − fu∗ )2 (γ 0∗ − ν 2 )2 gw D1 ,

AN US

a2 = γ 0∗ (γ 0∗ − fu∗ )×  ∗ ∗ b ∗ b 1 )(γ 0∗ D b 12 fu∗ − gw b 3 − ν2D b 1 fu∗ )(ν 2 gw b 1 )+ × 2ν 2 (gw −D D3 ) + (γ 0∗ − fu∗ )(D − γ 0∗ D h i ∗2 ∗ b b3 − D b 1D b 2 ) ν 2 (gw b 12 ) − 2γ 0∗ gw + (D +D D1 , h i 0∗ ∗ b ∗ b 0∗ b 2 ∗ b 2 f ∗ − g∗ D b 2 b b a1 = ν 2 (γ 0∗ D 1 u w 3 ) − (γ − fu )(D3 − D1 D2 )(gw D3 + γ D1 fu ) + h i ∗ b b 1 (D b3 − D b 1D b 2 ) 2D b 3 (γ 0∗ − fu∗ ) + gw b 1D b 2) , γ 0∗ D (D3 − D b 1D b 3 (D b3 − D b 1D b 2 )2 a0 = −D

(25)

and the bifurcation locus can be implicitly defined as

27a20 a23 − a21 a22 + 4a31 a3 + 4a0 a32 − 18a0 a1 a2 a3 = 0.

(26)

AC

CE

PT

ED

M

Owing to the algebraic complexity of such a locus, significative information are extracted through numerical investigations. In Fig.1 we report the bifurcation diagram obtained by varying the plant loss B and the rainfall A for different values of slope gradient ν and square characteristic velocity γ 0∗ . For simplicity, hereafter we consider γ(u) = γ0 u, being γ0 > 0 an arbitrary constant. More precisely, in Figs.1(a-c) we compare the results obtained by varying γ0 for fixed values of slope gradient (ν = 200, 30, 5) whereas in Figs.1(df) we present the dual situation obtained by varying ν and fixing the square characteristic velocity (γ0 = 2, 5, 10). In these figures, the dash-dot straight lines represent the critical condition A = 2B whereas all the other curves denote the bifurcation locus (26). These curves separate the B−A parametric plane into three domains for fixed values of ν or γ0 . In particular, the region below the line A = 2B is characterized by the absence of vegetation, being the desert state U∗0 the only equilibrium here admitted. The intermediate region delimited by A = 2B and the bifurcation locus corresponds to the scenario in which patterns of vegetation can be observed. Finally, in the region lying above such a bifurcation locus, the ecosystem develops a uniform vegetation. From a direct inspection of this figure, we can deduce that both the increase of the slope gradient ν and the decrease of the square characteristic velocity γ0 contribute in enlarging the domain in which patterns can be observed. Therefore, differently from the parabolic model, the hyperbolicity of the system provides an additional degree of freedom which can be used to tune the amplitude of the instability region. In addition, our results reveal that the behavior obtained for γ0 > 10 does not change quantitatively from the one drawn for γ0 = 10 which, in turn, mimic quite well that of the parabolic model [13]. In other words, by considering a linear dependence of γ on u, our results allow to estimate a critical value for the relaxation time which satisfactorily approximates the parabolic limit. Further confirmations on the value of such a threshold will be provided later on in the manuscript. The different dynamical behaviors depicted in Fig.1 can also be explained through the analysis of the wavenumber dependence of the real and imaginary parts of the three roots of (14) evaluated at U∗2 . Without loss of generality, in Fig.2 we carry out such an investigation for B = 1.75, γ0 = 5 and ν = 30, considering three different values of the rainfall parameter A corresponding to the points P1 ,P2 ,P3 indicated in Fig.1(b). As it can be noticed, two roots present negative real part independently of the wavenumber value ξ (see the dotted and dashed lines), so that the third root (solid line) determines the stability character of this steady state. In detail, at the point P3 (corresponding to A = 12.0), even the third root has negative real part (see Fig.2(e)). Therefore, for such a set of parameters, U∗2 is confirmed to be a stable focus and a homogeneous vegetation can be observed. At the point P2 (corresponding to A = 7.956), the third root is tangent to the axis Re {ω} = 0 at a given wavenumber ξ0 (see Fig.2(c)). In such a situation, a diffusion-advection-driven instability originates and vegetation patterns having 7

ACCEPTED MANUSCRIPT

35

35 0

30 0

A Rainfall

(a)

= 2

= 50

Homogeneous

= 100

Vegetation

10

(d)

= 30

30

Homogeneous

= 5

0

25

20

CR IP T

wavenumber ξ0 can develop. If the rainfall value is further decreased, the operating point lies below the bifurcation locus (as it happens for A = 5, which is represented by the point P1 ) so that there exists a range of wavenumbers [ξ1 , ξ2 ] corresponding to unstable modes (see Fig.2(a)). The ecological

A = 2B

Vegetation

= 200

25

A = 2B

= 200

20

0

= 2

Patterned

Vegetation 15

Patterned

10

Vegetation

5

AN US

15

10

5

No Vegetation

No Vegetation

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

25

15.0

12.5

0

= 5

Homogeneous

10

0

P

= 30 7.5

ED

P

5.0

P

= 30

3

2

20

= 100

Homogeneous

= 200

Vegetation

A = 2B

15

0

= 5 Patterned Vegetation

0

Patterned

(e)

= 50

Vegetation

A = 2B

10.0

A

(b)

= 2

M

0

Rainfall

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

10 P

1

Vegetation

2.5

5

PT

No Vegetation

No Vegetation

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

7

CE

0

6

5

25

(c)

= 2 = 5

= 50

10

0

4

3

(f)

= 30

Homogeneous

20

Vegetation

A = 2B

= 100

Homogeneous

= 200

Vegetation

A = 2B

= 5

A Rainfall

AC

0

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

15 0

Patterned

10

Patterned

Vegetation 10

Vegetation

2 5 1

No Vegetation No Vegetation 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

Plant Loss

B

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

Plant Loss

B

Figure 1: Bifurcation diagrams related to the equilibrium state U∗2 in the B-A parametric plane, for different values of slope gradient ν and square characteristic velocity γ0 .

8

ACCEPTED MANUSCRIPT

3000

1.5 (a)

(b)

wavenumbers of unstable modes

2000

} -3.0

Im {

Re {

}

-1.5

20

10

1500

0

1000

-4.5

-10 -3 10

500

-2

10

-6.0

10

-3

10

-2

10

-1

10

0

10

1

3

-1

0

10

10

AN US

0

CR IP T

2500

0.0

10

2

10

-3

10

-2

10

-1

10

0

10

1

10

2

3000

(c)

(d)

onset of instability

0

2500

20

-3

2000

-6

}

M

-9

Im {

Re {

}

10

ED

-12

-15

-18

-3

(e)

-2

10

-1

10

0

10

1

-10 -3 10

-2

-1

10

0

10

10

0

10

2

10

-3

10

-2

10

-1

10

0

10

1

10

2

3000 (f)

2500 20

-10

2000 10

}

-20

AC

Im {

Re {

}

0

1000

500

CE

0

10

PT

10

1500

-30

1500 0

1000

-10 -3 10

-2

-1

10

0

10

10

500

-40

0 -50 10

-3

10

-2

10

-1

10

0

10

1

10

2

10

Wavenumber

-3

10

-2

10

-1

10

0

10

1

10

2

Wavenumber

Figure 2: Wavenumber dependence of the real (left) and imaginary (right) part of the three roots of (14) evaluated for γ0 = 5 at the points P1 ,P2 ,P3 indicated in Fig.1(b). Figs.(a) and (b) correspond to P1 ; (c) and (d) to P2 ; (e) and (f) to P3

9

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

scenario corresponding to P1 is an illustrative example in which vegetation patterns with ξ ∈ [ξ1 , ξ2 ] can be observed depending on the choice of initial and boundary conditions. Remark. It should be pointed out that the results obtained through the present linear stability analysis are expected to hold when the patterns exhibit low amplitude, namely for parameters close to the bifurcation curves. For larger amplitude patterns, nonlinear effects could be relevant and induce dynamics significantly different from the ones here described. However, at this stage of investigation, we aim to remark the differences existing, in the linear regime, between the present hyperbolic model and the parabolic one [13]-[18]. Nevertheless, information on nonlinear dynamics are deduced from numerical simulations. Considering the complexity of the locus (26), in the following we provide an approximation through the use of power laws. In such a way, we achieve a twofold goal: first, we simplify the functional dependence of such a bifurcation locus on the main parameters here involved (plant loss, slope steepness and square characteristic velocity); second, we address a direct comparison with the results obtained in the parabolic case (see eq.(14) in [13] and eq.(5) in [14]). More precisely, we represent the surface A = κ1 ν α1 B β1 for α γ0 = 2 in Fig.3(a) and γ0 = 10 in Fig.3(b), whereas A = κ2 (γ0 ) 2 B β2 for ν = 5 in Fig.4(a) and for ν = 200 in Fig.4(b). The parameters κi , αi , βi (i = 1, 2) are obtained through a best fit procedure

Figure 3: Comparison between the exact bifurcation locus (26) (dotted lines) and its power law approximation A = κ1 ν α1 B β1 (surfaces) as a function of slope gradient ν and plant loss B, for fixed values of the square characteristic velocity: (a) γ0 = 2, (b) γ0 = 10. to the data extracted from the exact bifurcation locus (26). Results of this fitting procedure, shown in Tables 1 and 2, demonstrate a larger sensitivity of the rainfall level to reduced values of the square characteristic velocity and to more steep inclines, respectively. As it can be noticed from Table 1, while 10

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

PT

Figure 4: Comparison between the exact bifurcation locus (26) (dotted lines) and its power law approxiα mation A = κ2 (γ0 ) 2 B β2 (surfaces) as a function of square characteristic velocity γ0 and plant loss B, for fixed values of the slope gradient: (a) ν = 5, (b) ν = 200.

AC

the prefactor and the exponent of the plant loss are not so significantly affected by the reduction of γ0 , the dependence on ν becomes more relevant. In fact, if γ0 is reduced from 10 to 2, the parameter β1 increases of about 30% (from 1.239 to 1.584). We point out that the results obtained for γ0 = 10 are also in good agreement with those obtained by Sherratt for the corresponding parabolic system: κ1 = 0.60, α1 = 0.50, β1 = 1.25 [13]. Furthermore, from Table 2 we can also infer that the change in the slope gradient does not affect significantly the exponent of the plant loss, while it induces more significant modifications on the exponent of the square characteristic velocity, which even doubles when ν rises from 5 to 200. Finally, a direct inspection of Figs. 3 and 4 reveals that, when the plant loss approaches the limit value B = 2 and the square characteristic velocities is γ0 . 3, the power law approximations yield results which slightly deviate from the bifurcation locus (26). Let us now characterize more accurately the structure of the patterns arising from the diffusionadvection-driven instability. To this aim, we integrate numerically the system of PDEs (2)-(6) for γ0 = 10 and γ0 = 2 by considering the parameter set (B = 1.25, A = 5.4, ν = 30) corresponding to the point P0 indicated in Fig.1(b) which lies in the proximity of the bifurcation locus. The computation is performed for t ∈ [0, tend ] over the spatial domain x ∈ [0, xend ] by using periodic boundary conditions and with the following initial conditions

11

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 5: Spatial dependence of the field variables evaluated at tend = 200 with B = 1.25 and A = 5.4 [point P0 indicated in Fig.1(b)] for (a) γ0 = 10 and (b) γ0 = 2. The numerical integration of the system of PDEs uses the initial condition (27),(28) with xend = 1000.

u(x, 0) = 2 w(x, 0) = win (x) J(x, 0) = −0.001 h i 2 0.003 100 − |x − (xend /2)| 0

ED

win (x) =

(

M

where

κ1 0.815 0.810

α1 0.472 0.460

for |x − (xend /2)| ≤ 10 otherwise

(28)

β1 1.584 1.239

PT

γ0 2 10

(27)

ν 5 200

κ2 2.439 11.31

α2 −0.05 −0.112

AC

CE

Table 1: Parameters of the approximated bifurcation locus A = κ1 ν α1 B β1 β2 1.187 1.33

Table 2: Parameters of the approximated bifurcation locus A = κ2 (γ0 )

α2

B β2

The use of periodic boundary conditions is needed to allow a direct comparison with the results obtained in [14]. On the other hand, the initial conditions (27),(28) describe an uniformly vegetated habitat which is perturbed via a localized distribution of water in the middle of the hill. Despite analogous results have been obtained by using different (for instance, random) initial conditions, our choice aims to minimize possible spurious effects arising from the periodicity imposed at the boundaries. It should be also noticed that the wave features computed within a numerical framework which necessarily uses a finite domain are discrete. On the contrary, the analytical calculations performed over a theoretically unlimited domain imply continuous quantities. Nevertheless, our comparative analysis aims to show that, once the control parameters are fixed, the results obtained numerically are in qualitative agreement with the ones predicted analytically. In Fig.5 we show the spatial profiles of the field variables obtained at t = tend in the two investigated cases. Of course, owing to the downhill flow of water, only the region beneath the perturbation x ≤ xend /2 is meaningful and here represented. As shown in this figure, for γ0 = 10 the initial ecosystem 12

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

CR IP T

converges toward the homogeneous vegetation state U∗2 while, for γ0 = 2, it undergoes a dynamical regime characterized by spatial periodic oscillations for all the field variables. These behaviors support the results depicted in Fig.1(b), being the operating point P0 outside or inside the diffusion-driven instability region, respectively. Such an analysis also confirms that the reduction of the characteristic velocity can be used to enlarge the range of rainfall values at which vegetation patterns can originate. Moreover, the results depicted in Fig.5(b) give some insights into the mechanism of pattern migration. In fact, while the peaks associated to the water distribution are shifted upwards with respect to those of the vegetation biomass, very poor concentrations of water are present on the bottom edge of the patterns. These observations could indeed provide a clear evidence for the previously mentioned conjecture on the different moisture levels of the soil around the pattern which favor the advection-driven mechanism of uphill migration. Now, we choose the following set of parameters A = 1.25, B = 0.5 in such a way that the system falls within the instability region of pattern generation whatever the values of ν and γ0 . In this case, by integrating the system of PDEs with the abovementioned boundary and initial conditions, we analyze the spatio-temporal evolution of the vegetation biomass obtained by varying the slope gradient or the characteristic velocity. The corresponding results are shown in Figs.6 and 7, respectively.

AC

CE

Figure 6: Spatio-temporal dependence of the vegetation patterns observed for (a) ν = 200, (b) ν = 100 and (c) ν = 30. The parameters are: B = 0.5, A = 1.25, γ0 = 10, xend = 600, tend = 30.

Figure 7: Spatio-temporal dependence of the vegetation patterns observed for (a) γ0 = 10, (b) γ0 = 1 and (c) γ0 = 0.1. The parameters are: B = 0.5, A = 1.25, ν = 30, xend = 200, tend = 30. These figures clearly indicate the formation of non-stationary patterns which are thus associated with an instability which breaks both the temporal and the spatial symmetry. Such an instability is variously 13

ACCEPTED MANUSCRIPT

CR IP T

(m)

1.4

34

1.2

32

1.0

30

0.8

28

0.6

26 0.1

1

AN US

Wavelength

Migration speed s (m/year)

1.6

36

10

Square characteristic velocity

0.4

100

0

ED

M

Figure 8: Dependencies of the pattern wavelength Λ (filled squares) and uphill migration speed s (open circles) on the square characteristic velocity γ0 (depicted in log scale), computed by integrating the original system of PDEs. Dashed and dotted lines are representative of the logarithmic and constant trend lines which applies in the range γ0 < 3 and γ0 > 10, respectively.

AC

CE

PT

known as ”oscillatory Turing”, ”Turing-Hopf” or ”Wave” instability [30]. Moreover, our results show that, as time moves on, the creation of new patterns is a top-down process, according to the choice of the initial data, while the subsequent pattern migration is a bottom-up one. These dynamics are characterized by two physical quantities: the wavelength and the migration speed. More precisely, from Fig.6 we notice that the distance between two consecutive patterns increases with the slope steepness. Consequently, for a fixed time interval, a vegetation biomass is able to reach larger distances if it grows along steeper slopes. We also remark that the resulting monotonic behavior between the slope gradient and the wavelength is in full agreement with the results obtained in the parabolic model [see Eq.(15) in [13]]. On the other hand, Fig.7 also reveals that the reduction of the characteristic velocity affects significantly the uphill migration velocity while yields only a slight contraction on the pattern wavelength. In the following, we aim to quantify the functional dependence of the main pattern features on the square characteristic velocity. For this reason, we choose realistic parameter values that are valid for grass in semiarid regions, as reported by Klausmeier [9]. In detail, the mortality, the rainfall and the gradient slope are set to B = 0.45, A = 1.8 and ν = 182.5, respectively, whereas the characteristic velocity is swept in the range γ0 ∈ [0.1, 100]. According to the results shown in Fig.8, both wavelength and migration speed become independent of γ0 for γ0 ≥ 10, in agreement with our previous estimation of the parabolic limit. In the range 4 < γ0 < 10, these quantities exhibit a weak dependence on the square characteristic velocity whereas, for smaller values, they obey a logarithmic law. From this analysis we deduce that, by sweeping γ0 over three decades, the wavelength varies by about 15% only, whereas the migration speed can be tuned up to 70%. In other words, while the pattern wavelength is mostly determined by the slope gradient, the uphill migration speed exhibits a stronger dependence on the characteristic velocity. Therefore, we infer that it is possible to modulate, almost separately, the two main pattern features, wavelength and migration speed, by acting on two independent parameters, the slope gradient and characteristic velocity. The results depicted in Fig.8 allows us to address some additional comments from the ecological viewpoint. Indeed, the numerically computed values of wavelength Λ ∈ (29.5, 35) m and migration speed s ∈ (0.5, 1.6) m/year are in quantitative agreement with those reported in literature for grass (Λ ∈ (1, 40) 14

ACCEPTED MANUSCRIPT

4

CR IP T

m, s ∈ (0.3, 1.5) m/year, see [9]) and satisfactorily model the behavior observed in several semiarid regions of Mexico, Sudan, Mali and Niger (see, for instance, Table 12.2 in [19]). It is worthy noticing that, within the proposed hyperbolic model, the main pattern features are modulated by fixing the rainfall, the plant mortality and the slope gradient, i.e. without altering the equilibrium configurations. Such a possibility, precluded in the parabolic limit, reflects the assumption of a quasi-uniform distribution of the main environmental parameters within a bounded semiarid region. In this context, a proper identification of the characteristic velocity γ(u) would allow to model the different propagating character exhibited by each species of vegetation.

Periodic Travelling Waves

AN US

As shown in the previous Section, owing to the presence of the advection term in the evolution equation for the surface water, vegetation patterns move in the positive x direction. In view of this uphill migration, it is interesting to search for patterned solutions in the form of periodic travelling waves under the ansatz U(x, t) = U(z) where z = x − st, being s the wave speed. Substituting this solution into the governing PDE system (2) gives the ODE system (9). Therefore, in order to complete the analysis of formation and modulation of vegetation patterns, in this Section we investigate the occurrence of Hopf bifurcations in such a travelling waves system since, as known, patterned solutions of (2) correspond to limit cycles of (9) [14]-[18]. Linearizing the (9) around U∗ through the assumption

M

b az U = U∗ + Ue

ED

allows us to recast the original system as  as + fu∗ fw∗ ∗ ∗  gu a (s + ν) + gw aγ0 0

 −a  0 γ0 − as

 u b w b =0 Jb

(29)

(30)

so that the following characteristic equation for the eigenvalue a holds 3

2

PT

(as) + C1 (as) + C2 as + C3 = 0

C1 =

s (s+ν)(s2 −γ0 )

C2 =

s2 (s+ν)(s2 −γ0 )

AC

CE

where

C3 =

3

s (s+ν)(s2 −γ0 )

 2 ∗  ∗ ∗ − γ0 ) + νs (fu∗ − γ0 ) − gw γ0 , s (fu + gw

∗ ∗ {s [fu∗ gw − fw∗ gu∗ − γ0 (fu∗ + gw )] − νfu∗ γ0 } ,

(31)

(32)

∗ (fw∗ gu∗ − fu∗ gw ) γ0 .

It should be mentioned that the stability character of each equilibrium configuration has been addressed in [20] for the more general model involving a source term f1 (u) = up with p > 1. Here we specify the results therein obtained for p = 2. In particular, bearing in mind (17), we deduce that the desert state U∗0 is always unstable, namely √ 0 < s < γ0 ⇒ saddle point √ (33) s > γ0 ⇒ unstable node

so that no bifurcations can originate about such an equilibrium. At the equilibrium U∗min the coefficient C3 = 0 so that there exists a null eigenvalue which makes ∗ Umin stable (not asymptotically) iff C1 (U∗min ) > 0 and C2 (U∗min ) > 0, i.e. under the following conditions sb < s < min



 νB √ , γ0 , 2−B

B < γ0 15

(34) (35)

ACCEPTED MANUSCRIPT

0

0

= 2.0

CR IP T

12

= 100

= 0.75

10

8

6

P

4

P P

2

3

P

1

0 1.0

1.5

4

AN US

Travelling speed s

14

2.0

2

2.5

3.0

3.5

4.0

M

Rainfall A

ED

Figure 9: Hopf bifurcation diagram in the A − s parameter plane for γ0 = 0.75, 2, 100. The thick curves depict the condition (38), whereas the horizontal thin lines represent the upper limits provided by the √ right implication of (37), i.e. s = γ0 . The other parameters are: B = 0.45, ν = 182.5.

PT

where

sb =

ν (B − γ0 ) +

q

2

ν 2 (B − γ0 ) + 8γ0 (γ0 + 2 − B) 2 (γ0 + 2 − B)

.

(36)

AC

CE

For what concerns the equilibria U∗1,2 , by using the Routh-Hurwitz criterion and taking into account (22), we deduce that U∗1 is always unstable whereas U∗2 is asymptotically stable if the following restrictions hold se < s <

√

γ0

(C1 C2 − C3 )U∗ > 0

where

se =

ν (B − γ0 ) +

q

2

(37) (38)

2

ν 2 (B − γ0 ) + 4γ0 (1 + u22 ) (1 − B + u22 + γ0 ) 2 (1 − B + u22 + γ0 )

(39)

According to the local stability analysis carried out in the previous section, let us investigate the occurrence of a Hopf bifurcation at U∗2 which manifests itself as a small amplitude periodic solution of the ODEs system (9). As it can be verified, the Hopf bifurcation locus is defined by (C1 C2 − C3 )U∗ = 0 2

(40)

which is hereafter investigated numerically owing to its algebraic complexity. In particular, we depict in Fig.9 such a locus for different values of the square characteristic velocity by using the same set of parameters considered by Sherratt in [15]. We represent the locus (40) through √ the thick curves whereas the condition s = γ0 via horizontal thin lines. The left implication of (37) 16

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

PT

Figure 10: (Left) Dependence of the field variables (u, w) on the travelling wave coordinate (z). (Right) Trajectories in the u − w phase plane obtained by using A = 2, s = 0.7 (corresponding to the point P2 in Fig.9) and (a,b) γ0 = 0.75, (c,d) γ0 = 2, (e,f) γ0 = 100. The other parameters are the same as the ones used in Fig.9. The dotted curves are representative of the unstable limit cycles.

AC

CE

is, instead, not represented here since it is less restrictive than (38). Therefore, the resulting region in √ which patterns can be observed is delimited by (40) and s = γ0 . First of all, the bifurcation locus (40) obtained for γ0 > 10 does not differ significantly from the one obtained for γ0 = 10, confirming further our previous results. On the other hand, as γ0 is reduced, the pattern region extends towards larger values of the rainfall parameter A, consistently with the results shown in Fig.1. However, this is strictly valid for travelling speed values s & 0.8 because, for s . 0.8, the pattern region shrinks. From a √ direct inspection of Fig.9, we also argue that s = γ0 introduces an additional cutoff for the travelling ∗ speed s at which the equilibrium U2 loses its stability. In fact, as expected from the hyperbolic nature of the model, the existence of a non-null relaxation time determines an upper bound to the wave speed s. It is worth mentioning that such a restriction becomes particularly significant with the reduction of the characteristic velocity, namely far from the parabolic limit. Under such circumstances, the range of possible wave speeds is drastically reduced and, consequently, the possibility to determine what wave speed is selected by given initial and boundary conditions exhibits less uncertainties. From this viewpoint, our study could provide some hints for the solution of the open question reported in literature and known as ”pattern selection” issue [17]. In order to verify the results illustrated in Fig.9, let us now integrate numerically the ODE system (9) by using the set of parameters A = 2, s = 0.7 characterizing a state falling into the Hopf region whatever the value of γ0 . Then, the different dynamical behaviors obtained by varying the characteristic velocity and the initial conditions are shown in Fig.10. From Figs.10(b,d,f) it is possible to notice the existence of unstable limit cycles associated to subcritical Hopf bifurcations; indeed, any nearby trajectory spirals away from them. This result provides a further evidence of patterned solutions in the original hyperbolic

17

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

Figure 11: Real (left) and imaginary (right) part of the three roots of (31) evaluated for A = 2 and γ0 = 2 [figs.(a,b)] or γ0 = 100 [figs.(c,d)].

AC

CE

PT

PDE system. In addition, by analyzing such solutions, we can also confirm the weak dependence of the pattern wavelength on the square characteristic velocity. In fact, by reducing γ0 by almost two orders of magnitude (from 100 to 0.75), the resulting wavelength Λ varies from 16 to 12, namely it has undergone a reduction of 25% only, as shown in Figs.10(a,c,e). From these figures, we also get the further validation on the shift existing between the distributions of vegetation and water which could originate the pattern migration. Finally, we inspect the dependence of the roots of the characteristic equation (31) on the travelling speed s for fixed values of rainfall level A and square characteristic velocity γ0 . As an illustrative example, in Fig.11 we evaluate numerically the real and imaginary parts of such roots for A = 2 (vertical line in Fig.9). In particular, we represent the results obtained for γ0 = 2 in Figs.11(a),(b) whereas those corresponding to γ0 = 100 in Figs.11(c),(d). In detail, for s < s(P1 ) (see Fig.9 for the labels of the operating points), two roots exhibit positive real part, so that in this range the state U∗2 is unstable, as expected. Such roots change the sign of their real part in correspondence of the point P1 , which lies indeed on the bottom of the Hopf bifurcation locus. Then, the state U∗2 becomes stable in the region s(P1 ) < s < s(P3 ) for γ0 = 2 or in s(P1 ) < s < s(P4 ) for γ0 = 100, having all the roots negative real part. It should be noticed that the occurrence of different cutoff conditions can also lead to different √ dynamical behaviors. More precisely, for γ0 = 2, the state U∗2 loses its stability at s = γ0 (point P3 ) via a stationary bifurcation since the real part of one root change sign from negative to positive. On the contrary, for γ0 = 100, the instability occurs via a Hopf bifurcation since two conjugate roots exhibit null real part at the point P4 , which lies along the bifurcation locus (40).

5

Exact solutions via differential constraints method

In this Section we investigate a class of exact solutions for the more general model (2)-(4). As known, exact solutions of mathematical models, apart from their own theoretical value, are of relevant interest

18

ACCEPTED MANUSCRIPT

AN US

CR IP T

since they provide significant guidance in elucidating the properties of a system and can be also useful for testing numerical procedures for the system of PDEs under interest. Among the reduction methods to find exact solutions of nonlinear PDEs, we choose the differential constraints approach for a twofold reason: (i) the hyperbolicity of the model facilitates a systematic search of such solutions; (ii) the considered approach includes most of known methods for constructing solutions of PDEs: classical and non-classical symmetries [31, 32], weak symmetries [33], conditional symmetries [34]. In this context, the classical simple wave solutions to quasi-linear hyperbolic first-order homogeneous systems of PDEs can also be recovered [35, 36]. The differential constraint method can be briefly summarized as follows. Let us consider a quasilinear hyperbolic system in the general form (2), where U ∈ Rn is a column vector denoting the dependent field variables. As already mentioned, the system (2) is strictly hyperbolic, namely the n − th order matrix M admits n real eigenvalues λi 6= λj (i, j = 1, .., n) to which correspond n right eigenvectors di and n left eigenvectors li spanning the Euclidean space E n . Following the investigations carried out in [37]-[39] we append the following set of m differential constraints to the system (2) li · Ux = q (i) (x, t, U),

i = 1, ...m ≤ n

(41)

ED

M

with q (i) scalar functions to be determined. The first step of this procedure requires the compatibility between equations (2) and (41). This leads to an overdetermined set of conditions to be fulfilled by the functions q (i) and, in turn, restricts the functional form of the matrix coefficients M and N. Once these consistency conditions have been satisfied, then a class of exact solutions for the model under investigation is obtained by integrating the overdetermined system (2),(41). The solutions so deduced, which depend on n − m arbitrary functions, generalize the so-called multiple wave solutions and may be useful to solve suitable classes of initial and/or boundary value problems [40]. From a theoretical viewpoint, the most relevant case involves m = n − 1 differential constraints [38],[39]. In fact, in this case, it is straightforward to ascertain that (2) and (41) lead to Ut + λ n Ux = N +

n−1 X i=1

q (i) (λn − λi ) di

(42)

AC

CE

PT

so that the searched solutions to the governing equations can be obtained through integration along the characteristic curves associated to λn . Thus, since the solutions in point have inherent wave features, they may describe propagating signals. Moreover, the n − 1 constraint equations (41) also select the class of initial value problems since the initial data U(x, 0) associated to (42) must also satisfy the (41). Therefore, taking into account (8a), the possible first order differential constraints (41) compatible with (2)-(3) become wx = q (1) (u, w, J, x, t),

(43a)

Jx = q

(2)

(u, w, J, x, t) + λ2 ux ,

(43b)

Jx = q

(3)

(u, w, J, x, t) + λ3 ux .

(43c)

In this work we report the results obtained when the differential constraints (43b),(43c) are appended to (2)-(4). Details of such an analysis are given in the Appendix. Hereafter we focus the attention on the case I reported in the Appendix, whose results might provide meaningful source terms for the model described by (2)-(4). In particular, since the proposed approach offers the possibility to select the response functions that better mimic the reactions involved in the ecological scenario under investigation, in the following we try to characterize the exact solutions in order to describe the dynamics occurring in the Klausmeier model. To this aim, by requiring that g1 (u) and g2 (u) fulfil (5), from (A.8) the response functions f1 (u) and f2 (u) specialize to 1

f1 (u) =

(1 − νK0 )|νK0 u + c1 |1− νK0

1

(1 − νK0 )F1 − |νK0 u + c1 |1− νK0

f2 (u) = (νK0 u + c1 ) F2 + A(1 − νK0 )

R 19

, 1

|νK0 u + c1 |−1− νK0

1

(1 − νK0 )F1 − |νK0 u + c1 |1− νK0

!

du ,

(44)

ACCEPTED MANUSCRIPT

with

1   |νK0 u + c1 |2− νK0 |νK0 u + c1 | + F1 |c1 | + u < 0, νK0 − 1 νK0 u + c1

(45)

being c1 , K0 , F1 and F2 constants. Finally, integrating (A.7) along y = x + νt we have

CR IP T

 1  b (y) , c1 + exp (−νK0 t) U νK0 K1 b exp(−νK0 t), J(x, t) = exp(−γ0 t) + J(y) γ0 − νK0   1− νK1  1 0 b (y) exp(−νK0 t) U × w(x, t) = F1 exp(−t) − 1 − νK0   Z exp(t) c  × W (y) + A(1 − νK0 )  1− νK1 d t 0 b (y) (1 − νK0 )F1 − exp(t − νK0 t) U

(46)

AN US

u(x, t) = −

ED

M

b (y), J(y) b c (y) are specified by applying the constraints where K1 is a constant. The functions U and W (A.1a) on the initial data u(x, 0), J(x, 0) and w(x, 0). This leads to       b (0) exp νK0 − γ0 y − U b (y) + b = J(0) b exp νK0 − γ0 y − γ0 U J(y) ν ν 2 K0 ν Z y    γ0 (γ0 − νK0 ) νK0 − γ0 γ 0 − νK0 b − y σ d σ, exp U (σ) exp (47) ν 3 K0 ν ν 0 ! 2 b 1 1 c (y) = (U b (y)) νK0 −1 νK0 − γ0 J(y) b − ν − γ0 d U + (1 + F2 )(U b (y)) νK0 . W 2 ν ν K0 d y

AC

CE

PT

b , J, b W c is not subject to any restriction, so that a class of where, as expected, one of the functions U solution for the system (2)-(5) is selected. The results so obtained exhibit the advantage of offering several potential broad applications. Indeed, once the response functions have been determined, by means of a proper choice of the arbitrary function b (y) and of the parameters involved in (46),(47), it is possible to select that particular class of exact U solutions which better elucidates the properties of the investigated system. As an illustrative example, in order to enable a direct comparison with the response functions appearing in (6), we choose c1 = 0, F1 > 0 and νK0 = −1, so that we have  √   2u2 A 2F1 − u √ √ f1 (u) = , f (u) = u ln − F (48) 2 2 2F1 − u2 2F1 2F1 + u  √  with u ∈ 0, 2F1 . Of course, since the source terms appearing in the original Klausmeier model (5),(6) differ from those deduced via the differential constraints method (48), the comparison is addressed at qualitative level only. Bearing in mind such a goal, we analyze first the fixed points admitted by these models through the intersection of the corresponding nullclines. As shown in Fig.12, both models exhibit the same qualitative behavior and, indeed, when the rainfall value is varied in such a way A < 2B, A = 2B and A > 2B, we get one, two or three equilibria, respectively. It should be noticed that, in all cases, the desert state is a common steady state and that the condition u2 < 1 = umin < u2 is always fulfilled. These results are in full agreement with those reported in (10),(11). Finally, according to the previous choice of the arbitrary constants, the exact solution (46) specializes

20

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 12: Comparison between the nullclines of the Klausmeier model (left) and those obtained via the differential constraints method (right). Solid and dashed lines are representative of f (u, w) = 0 and g(u, w) = 0, respectively. The parameters are: B = 0.6, F1 = 6, F2 = 0.05. (a),(b) are obtained for A = 0.5; (c),(d) for A = 1.2; (e),(f) for A = 1.4.

21

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

b (x + νt), u(x, t) = exp(t)U

CE

to

PT

Figure 13: Spatio-temporal evolution of the field variables obtained through (49),(50). The parameters b = 0.02, K1 = 0, A1 = are the same as the ones used in Fig.12 together with A = 1.5, ν = 200, γ0 = 5, J(0) 3, A2 = 0.5, k2 = 0.1, k3 = 1.3k2 . (a) A3 = 0, (b) A3 = 0.75A2 .

AC

b + νt) + K1 exp(−γ0 t), J(x, t) = exp(t)J(x γ0 + 1   1 b 2 (x + νt) × w(x, t) = exp(t) 2F1 exp(−2t) − U 2 !! √ b (x + νt) A 2F + exp(t) U 1 c (x + νt) + √ × W ln √ b (x + νt) b (x + νt) 2F1 U 2F1 − exp(t)U

(49)

b is the arbitrary function, whereas Jb and W c are obtained from (47). where U The possibility to handle exact solutions allows us to explore the dynamical behavior of the system. In particular, we are interested in evaluating whether the solution (49) can be used to generate patterned configurations. To this aim, in Fig.13 we depict the spatio-temporal evolution of the three field variables obtained for a set of parameters which falls into the condition A > 2B and by using the following expression for the initial data b (y) = [A1 + A2 cos (k2 y) + A3 cos (k3 y)] exp (−y/ν) U

(50)

with A1 ,A2 ,A3 ,k2 and k3 arbitrary constants through which we also manage the harmonic character of the initial conditions. From the comparison of Fig.6, where the solutions are anharmonic in space, with Fig.13, where the 22

ACCEPTED MANUSCRIPT

Conclusions

AN US

6

CR IP T

solutions can be either harmonic (a) or anharmonic (b), we get a further confirmation on the possibility of the extended model to capture a larger variety of pattern dynamics. Moreover, by means of a proper choice of the arbitrary function (50), it is possible to reproduce an oscillatory behavior which qualitatively recalls the pattern formation observed in Section 3. Nevertheless, some quantitative differences have to be pointed out. In fact, the patterns predicted by the Klausmeier model are typically composed by vegetated regions separated by bare ground areas (see Figs. 6,7). On the contrary, the particular exact solution (49) reflects the functional dependence of the initial condition (50) and thus describes spatially-periodic aggregations of vegetation biomass not separated by bare ground areas and exhibiting an exponentially-decreasing trend. It is worthy noticing that such patterns undergo an analogous mechanism of uphill migration, although the corresponding migration speed values appear to be larger than the typical ones.

AC

CE

PT

ED

M

In this work the occurrence of patterned solutions in a hyperbolic reaction-diffusion-advection model of vegetation valid for semiarid environments is investigated analytically and numerically. Our results have shown that the hyperbolicity of the model does not bring the only advantage of removing the unphysical instantaneous diffusive effect typical of parabolic systems. More interestingly, it can be efficiently used to provide an additional degree of freedom to describe the complex vegetation dynamics. It is indeed shown that the main pattern features, wavelength and uphill migration speed, can be modulated almost independently by acting on slope steepness and characteristic velocity, without affecting the equilibrium states. In particular, the phenomenological response function γ(u) plays a fundamental role in this context since it allows to modify the bifurcation region in which patterns can be observed and to mimic the specific propagating character of a given species of vegetation. It should be emphasized that, by using a well established set of environmental parameters, it has been possible to reproduce quantitatively the values of the main pattern features observed in semiarid regions. Our analysis has also permitted to estimate a threshold value for the relaxation time (τ = 0.1) below which the behavior of the hyperbolic system approximates quite well that of the corresponding parabolic one. Moreover, by analyzing the spatial profiles of vegetation biomass and water distribution, a possible explanation for the mechanism which originates the uphill migration of the patterns has been provided. The hyperbolicity of the model has facilitated the search for a class of exact solutions via the differential constraint method. In particular, we have presented those solutions which are meaningful from the ecological viewpoint and are also in qualitative agreement with the results predicted by the Klausmeier model. Finally, we expect that a larger availability of long-term experimental data could support the results predicted by our hyperbolic model and may also contribute to a more accurate identification of the phenomenological response function in terms of the ecological quantities here involved.

7

Acknowledgments

This work was supported by INdAM-GNFM.

Appendix: Differential constraints analysis The requirement of consistency of the model outlined in (2)-(4) with two differential constraints leads to an overdetermined system which allows to characterize q (i) as well as to select classes of response functions f (u, w), g(u, w) involved in the model. We consider the set of equations (2)-(4) together with the differential constraints (43b),(43c) which can be recast in the form

23

ACCEPTED MANUSCRIPT

q (2) − q (3) q (2) + q (3) √ = r(u, w, J, x, t) , ux = = s(u, w, J, x, t) 2 2 γ0 ut = f (u, w) − r, wt = νwx + g(u, w), Jt = −γ 0 (u) (J + s)

Jx =

(A.1a) (A.1b)

CR IP T

whereupon, by cross differentiation, the resulting consistency conditions specialize to rw + νsw = fw ,

(A.2a)

νrw + γ 0 sw = 0,

(A.2b)

rx + st + (f − r)su + gsw − γ 0 (J + s)sJ + rrJ + s(r − f )u = 0, 0

00

0

(A.2c)

rt + γ sx + (f − r)ru + grw + γ s(s + J) + γ [r − (s + J)rJ + rsJ + ssu ] = 0.

From equations (A.2a,A.2b) we obtain νG − γ 0 (f + F ) , ν2 − γ0

s=

ν(f + F ) − G , ν2 − γ0

AN US

r=

(A.2d)

(A.3)

which, inserted into (A.2c,A.2d), give the following equations for F (u, J, x, t) and G(u, J, x, t) νGx − γ 0 Fx + νFt − Gt + (νf − G)Fu + F Gu − γ 0 (νJ + f + F )FJ + +(γ 0 J + G)GJ + νgfw − νF fu −

γ 00 F ν 2 −γ 0

[ν(f + F ) − G] = 0,

M

 νGt + γ 0 (νFx − Gx − Ft ) + (νf − G)Gu + γ 00 [ν(f + F ) − G] J +

νF

ν 2 −γ 0



(A.4) +

+γ 0 [F Fu − (νJ + f + F )GJ + (G + γ 0 J)FJ + F fu − gfw + νG − γ 0 (f + F )] = 0 and, in turn, differentiated with respect to w, lead to

f1 (u) (H(η, x, t) + Φ(u)) , ν2 − γ0

ED

F (u, J, x, t) =

PT

γ0 G(u, J, x, t) = − F (u, J, x, t) − γ 0 J + K(η, x, t) ν with

Φ(u) =

Z

γ0 − ν2


g1 (u) d u, f1 (u)

η = νJ + γ(u).

(A.5)

(A.6)

AC

CE

Further insertion of (A.5) into (A.4) allows to determine the functions H(η, x, t), K(η, x, t) involved in (A.5). Once H(η, x, t) and K(η, x, t) are determined, the solution u(x, t) , w(x, t) and J(x, t) of (2)-(4) is dx = −ν, i.e. obtained by integrating the set of equations (42) along the characteristic curve dt  ut − νux = −F (u, J, x, t)    wt − νwx = g(u, w) (A.7)    Jt − νJx = −γ 0 J − G(u, J, x, t)

The relations (A.4) induce also restrictions on the response functions γ(u), f (u, w) and g(u, w) leading to the following two cases. Case I)  γ(u) = γ0 u + γ1       f 0 (u)  g1 (u) = −νK0 + (νK0 u + c1 ) 1 f1 (u)     f2 (u) f 0 (u)   + (νK0 u + c1 ) 2  g2 (u) = −νK0 f1 (u) f1 (u) 24

(A.8)

ACCEPTED MANUSCRIPT

with f1 (u) and f2 (u) arbitrary functions, whereas K(η, t) = K0 η + K1 exp(−γ0 t) − K0 γ1 +

H = H0 ,

c1 γ0 ν

(A.9)

Case II)

(A.10)

AN US

 γ(u) = γ0 u + γ1 ,    f1 (u) = α1 , f2 (u) = α2 u + α3 ,    g1 (u) = −γ0 , g2 (u) = α4 ,

CR IP T

and γ0 , γ1 , c1 , H0 , K0 , K1 arbitrary constants. It is interesting to notice that the functional dependencies (A.8) could be meaningful from the ecological viewpoint. For this reason, we used them to provide the solution of the original system, as reported in the main text.

with γ0 , γ1 and αi (i = 1, ..4) constants, whereas H(η, x, t) and K(η, x, t) are solutions of
  α1 Ht − γν0 Hx − γ0να2 1 νK + α1 H + (ν 2 − γ0 )(η − γ1 ) Hη + h i γ2 +(ν 2 − γ0 ) Kx + K Kη − γν0 K + δ1 + ν02 (η − γ1 ) + α1 δ2 H = 0 2

γ2

(A.11)

M

0 1 2 α1 H + ν 2α−γ (−2γ0 Hx + ν +γ Ht ) νKx − Kt + νKKη − 2γ0 K + ννδ2 −γ ν 0 0 i h 2 2 γ γ0 γ0 −ν + να21−γ (η − γ1 ) − αν1 H − 2K Hη + ν0 (η − γ1 ) + νδ1 = 0 ν 0

References

PT

ED

where δ1 = α1 α4 + α3 γ0 , δ2 = ν02 − α2 . Unlike the previous case, the functional form (A.10) implies the independence of the vegetation density on the water redistribution. Owing to the lack of ecological meaning, such a scenario will not be further investigated.

CE

[1] J. von Hardenberg, E. Meron, M. Shachak, Y. Zarmi, Diversity of Vegetation Patterns and Desertification, Phys. Rev. Lett. 87 (2001) 198101

AC

[2] F. Borgogno, P. D’Odorico, F. Laio, L. Ridolfi, Mathematical models of vegetation pattern formation in Ecohydrology, Reviews of Geophysics 47 (2009) RG1005 [3] K. Gowda, H. Riecke, M. Silber, Transitions between patterned states in vegetation models for semiarid ecosystems, Phys. Rev. E 89 (2014) 022701 [4] K. Siteur, M.B. Eppinga, D. Karssenberg, M, Baudena, M.F.P. Bierkens, M. Rietkerk, How will increases in rainfall intensity affect semiarid ecosystems?, Water Resources Research 50 (2014) 59806001 [5] D. D’Onofrio, M. Baudena, F. D’Andrea, M. Rietkerk, A. Provenzale, Tree-grass competition for soil water in arid and semiarid savannas: The role of rainfall intermittency, Water Resources Research 51 (2015) 169-181 [6] Y. Mau, L. Haim, E. Meron, Reversing desertification as a spatial resonance problem, Phys. Rev. E 91 (2015) 012903 [7] D. Escaff, C. Fernandez-Oto, M. G. Clerc, M. Tlidi, Localized vegetation patterns, fairy circles, and localized patches in arid landscapes, Phys. Rev. E 91 (2015) 022924

25

ACCEPTED MANUSCRIPT

[8] S. K´efi et al., Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems, Nature 449 (2007) 213-217 [9] C.A. Klausmeier, Regular and irregular patterns in semiarid vegetation, Science 284 (1999) 1826-1828

CR IP T

[10] R. HilleRisLambers, M. Rietkerk, F. van de Bosch, H.H.T. Prins, H. de Kroon, Vegetation pattern formation in semi-arid grazing systems, Ecology 82 (2001) 50-61 [11] J. van de Koppel et al., Spatial heterogeneity and irreversible vegetation change in semiarid grazing systems, Americ. Nat. 159 (2002) 209-218 [12] K. Siteur et al., Beyond Turing: The response of patterned ecosystems to environmental change, Ecological Complexity 20 (2014) 81-96

AN US

[13] J.A. Sherratt, An analysis of vegetation stripe formation in semi-arid landscapes, J. Math. Biol. 51 (2005) 183-197 [14] J.A. Sherratt and G.J. Lord, Nonlinear dynamics and pattern bifurcations in a model for vegetation stripes in semi-arid environments, Theoretical Population Biology 71 (2007) 1-11 [15] J.A. Sherratt, Pattern solutions of the Klausmeier model for banded vegetation in semi-arid environments I, Nonlinearity 23 (2010) 2657-2675

M

[16] J.A. Sherratt, Pattern solutions of the Klausmeier model for banded vegetation in semi-arid environments II: patterns with the largest possible propagation speeds, Proc. R. Soc. A 467 (2011) 3272-3294

ED

[17] J.A. Sherratt and A.D. Synodinos, Vegetation patterns and desertification waves in semi-arid environments: mathematical models based on local facilitation in plants, Discrete and Cont. Dyn. Syst. Series B 17 (2012) 2815-2827 [18] J.A. Sherratt, Pattern solutions of the Klausmeier model for banded vegetation in semi-arid environments III: the transition between homoclinic solutions, Physica D 242 (2013) 30

PT

[19] D.J. Tongway, C.Valentin and J. Seghieri (editors), Banded Vegetation Patterning in Arid and Semiarid Environments Ecological, Springer, New York, 2001

CE

[20] E. Barbera, C. Curr` o, G.Valenti, On discontinuous travelling wave solutions for a class of hyperbolic reaction-diffusion models, Physica D 308 (2015) 116 [21] I. Muller, T. Ruggeri, Rational Extended Rational Thermodynamics. Springer, New York, 1998

AC

[22] S.R. Dumbar, H.G. Othmer, On a nonlinear hyperbolic equation describing transmission lines, cell movement, and branching random walks. In Nonlinear Oscillations in Biology and Chemistry, Othmer HG (ed.). Lecture Notes in Biomathematics, Springer, Berlin 66 (1986) 274 [23] J. Fort, V. M´endez, Wavefronts in time-delayed reaction-diffusion systems. Theory and comparison to experiment, Reports on Progress in Physics 65 (2002) 895-954 [24] E. Barbera, C. Curr` o, G. Valenti, A hyperbolic reactiondiffusion model for the hantavirus infection, Math. Methods Appl. Sci. 31 (2008) 481-499 [25] E. Barbera, C. Curr` o, G. Valenti, A hyperbolic model for the effects of urbanization on air pollution, Appl. Math. Model. 34 (2010) 2192-2202 [26] E. Barbera, C. Curr` o, G. Valenti, Wave features of a hyperbolic prey-predator model, Math. Methods Appl. Sci. 33 (2010) 1504-1515 [27] T. Ruggeri, Can Constitutive Relations be represented by non-local equations?, Q. Appl. Math. 70 (2012) 597-611

26

ACCEPTED MANUSCRIPT

[28] E. Barbera, G. Consolo, G. Valenti, Spread of infectious diseases in a hyperbolic reaction-diffusion susceptible-infected-removed model, Phys. Rev. E 88 (2013) 052719 [29] E. Barbera, G. Consolo, G. Valenti, A two or three compartments hyperbolic reaction-diffusion model for the aquatic food chain, Math. Biosciences and Engineering 12 (2015) 451-472

CR IP T

[30] E.P. Zemskov and W. Horsthemke, Diffusive instabilities in hyperbolic reaction-diffusion equations, Phys. Rev. E 93 (2016) 032211 [31] G.W. Bluman and S. Kumei, Symmetries and Differential Equations, Berlin, Springer, 1989

[32] G.W. Bluman and J.D. Cole, The general similarity solution of the heat equation, J. Math. Mech. 18 (1969) 1025-1042

AN US

[33] P.J. Olver and P. Rosenau, The construction of special solutions to partial differential equations, Phys. Lett. A 114 (1986) 107-112 [34] D. Levi and P. Winternitz, Non-classical symmetry reduction: example of the Boussinesq equation, J. Phys. A: Math. Gen. 22 (1989) 2915-24 [35] G. Boillat, Symtrisation des systmes d’quations aux drives partielles avec densit d’nergie convexe et contraintes, C. R. Acad. Sci. Paris 295 (1982) 551-554

M

[36] G. Boillat, On nonlinear plane waves, Rendiconti del Circolo Matematico di Palermo, Serie II Suppl. 45 (1996) 69-72 [37] S.V. Meleshko, V.P. Sapeev, N.N. Janenko, Soviet Math. Dokl. 22 (1980) 2

ED

[38] N. Manganaro, S.V. Meleshko, Reduction Procedure and Generalized Simple Waves for Systems Written in Riemann Variables, Nonlinear Dynamics 30 (2002) 87-102 [39] C. Curr`o, D. Fusco, N. Manganaro, A reduction procedure for generalized Riemann problems with application to nonlinear transmission lines, J. Phys. A: Math. Theor. 44 (2011) 335205

AC

CE

PT

[40] L.V. Ovsiannikov, Group Analysis of Differential Equations. Academic Press, New York, 1982

27