# Lotka–Volterra models and the global vegetation pattern

## Lotka–Volterra models and the global vegetation pattern

Ecological Modelling 135 (2000) 135 – 146 www.elsevier.com/locate/ecolmodel Lotka – Volterra models and the global vegetation pattern Yuri Svirezhev ...

Ecological Modelling 135 (2000) 135 – 146 www.elsevier.com/locate/ecolmodel

Lotka – Volterra models and the global vegetation pattern Yuri Svirezhev * Potsdam Institute for Climate Impact Research, PO Box 601203, Telegrafenberg, D-14412 Potsdam, Germany Received 26 November 1999; received in revised form 10 May 2000; accepted 6 June 2000

Abstract How can such a classic object of mathematical ecology as the Lotka – Volterra competition model for a community of n competing species be used for description of the global vegetation dynamics under the climate change? We assume that the present spatial distribution of global vegetation is a stable equilibrium solution of corresponding Lotka–Volterra equations. The problem is how to construct some discrete structure on a continuous set of parameters, in this context remaining in a framework of a continuous model description? The problem can be solved if a dynamical system with multiple equilibria is considered. In this case, a continuous quasi-stationary change of parameters induces a jump from one to another equilibrium, since they have different stability domain in the parametric space. On this conceptual base a special class of Lotka – Volterra equations is developed and a biologically interpreted procedure for estimating their coefficients is suggested. For this we must have maps of the annual production and biomass, a life span of each vegetation type, and a map of the current geographical distribution of different vegetation, i.e. the current global vegetation pattern (GVP). The latter is needed for the construction of the ‘elementary map’, which is an important part of the model. Another important part is the formula, which describes an annual production dependence on the temperature and precipitation (for instance, Lieth’s formula can be used). Considering a one-dimensional particular case for n =2 we have the analytical formulas describing a shift of borders between two vegetation zones under the climate change. It was shown that two different types of the transition zones, namely, ‘soft’ and ‘hard’, exist in this case. If the soft zone is characterised by the continuous and smooth replacement of one type of vegetation by another, then the hard zone is a typical ‘fractal’ structure with a mosaic of different vegetation. A special calibration procedure for the Lotka – Volterra model describing the dynamics of one-dimensional GVP with n types of vegetation is suggested. It is based on the stability theorem for a system with multiple equilibria and a special averaging method applying to real geographical distributions of the annual production and the biomasses. The results are used for estimating the shift of one transition zone between taiga and steppe in the Central Siberia under the climate change (CO2-doubling scenario). © 2000 Elsevier Science B.V. All rights reserved. Keywords: Global modelling; Global vegetation pattern; Climate change; Lotka – Volterra equations

1. Introduction

* Tel.: +49-331-2882671; fax: +49-331-2882695. E-mail address: [email protected] (Y. Svirezhev).

How can a classic object of mathematical ecology, the Lotka–Volterra equations, be used for the description of the dynamics of global vegeta-

0304-3800/00/\$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 4 - 3 8 0 0 ( 0 0 ) 0 0 3 5 5 - 0

136

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

tion under the climate change? Previously Svirezhev (1999) suggested the model of the global vegetation pattern (GVP) based on a probabilistic approach, which uses the so-named ‘urn scheme’. Svirezhev’s approach allows the creation of a family of different dynamic models for the GVP (‘evolutionary’ models). However, these models do not explicitly use such a concept as production, but the GVP dynamics very strongly depend on this value. Besides, these models operate only with two types of vegetation: forest and grass (although, there are no principal constraints to expand them for an arbitrary number of species). An attempt to apply Lotka – Volterra’s formalism to the description of the spatial dynamics of two competing plant species has been carried out by Svirezhev (1994) and Churkina and Svirezhev (1995). If in the first work a movement along one spatial co-ordinate is considered as some slow dynamics in a parametric space then in the second work a movement in the transition zone between two types of vegetation was considered as a wave solution of two reaction-diffusion equations. The similar approach was suggested by Jesse (1999) for modelling of the climate induced shifting of tundra and forest in North America. A simple Lotka – Volterra model is used in the simulation model TRIFFID for the description of competition between different biomes (Huntingford, 1999). One traditionally regards that in comparison with the ‘prey–predator’ models the competition models possess less rich behaviour. But the resent investigations (see, for instance, Tilman, 1982) shown that the competition is responsible for such effects as the appearance of dissipative structures and spatially organised patterns in plant communities. Ebenhoh (1994) and Ekschmitt and Breckling (1994) have carried out a sufficiently detailed survey of different models describing the competition and, as a consequence, coexistence. Vandermeer (1996), using the new concept of neutral stability, could explain the thinning process in rain forest. In this work an attempt is doing to create some special class of the Lotka – Volterra competition models, which directs on the solution of one very concrete problem: the description of the GVP dynamics. If the model contains n state variables,

which are usually interpreted as number of individuals or biomass (density) of corresponding species, then here the competing species are interpreted as n different types of vegetation and the ith state variable is the density of living biomass of ith type at geographical point (x, y). Different stable equilibria of the model interpret as either the existence of only one tree species (or, a form) in a given territory or the coexistence of two species. The evolution of GVP under climate change is described by the change of stability domains under the change of model parameters. In turn, these parameters are considered as values depending on the geographical co-ordinates.

2. Special class of Lotka–Volterra models At the beginning we introduce the following definition:  The ‘GVP’ is a map, which represents a geographical distribution of different types of vegetation. This map is a space discrete structure, which changes continuously with time. Its mapping onto the set of continuous climate parameters (temperature, precipitation, soil moisture, etc.) is called the ‘Holdridge diagram’(see, for instance, Prentice et al., 1992). The mapping is not one-to-one: the same combination of climatic parameters can correspond to different types of vegetation. We take also into consideration the biomass densities of different vegetation types, their production and another elements of the carbon cycle. In fact, all these values are measured in carbon units. We assume that a continuous model describes all these values. But, on the other hand, the GVP is a discrete structure! Then the following mathematical problem arises: how to construct some discrete structure on a continuous set of parameters, remaining within a framework of a continuous model description? We shall prove that it can be done using the formalism of the Lotka–Volterra model for competing species (Svirezhev and Logofet, 1983). Firstly, we consider the ‘point’ dynamics of vegetation. The main part of any vegetation (and carbon cycle) model is the productivity function

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

or the net annual production. The general equation for living biomass, B, is the following: dB = P−mRB dt

(2.1)

where P is the net production, tR =1/mR is the real residence time for carbon contained in the biomass. In fact, P=P(T, WS, I, C . . .) where T is the temperature, WS is the soil moisture, which is a function of precipitation H, I is the light, C is the carbon concentration in the atmosphere, etc. All these variables might be either explicit or implicit functions of time. We assume that the real residence time differs from its maximal value, which can be attained in the absence of inter- and intraspecies competition. The latter reduces the residence time by means of increase of the speciesspecific mortality coefficient m, so that: mR = m +G(B, BS)

(2.2)

#Bi = Pi (Bi ; T, H)−m R i Bi, #t

137

i= 1, . . . , n;

(2.3)

where T and H are the annual temperature and precipitation at this point. Some additional assumptions must also be formulated: (a) Net primary production Pi (Bi ; T, H) is presented as: Pi = gi (T, H)Bi, where gi depends only on T and H. This means that the net production, depending on local climatic parameters, depends implicitly, through their geographical localisation, on geographical coordinates. For instance, T =T(x, y), H= H(x, y), etc. If these parameters change with time, then the production must also depend on time. (b) Local values of gi are so-called P/B-coefficients. They have known values at the equilibrium, which is identified with present state of the GVP. (c) Dependence of the coefficients m R i on competition is considered as a simple linear form: n

where Bs are the biomasses of other competing species. The value t =1/m can be interpreted as a mean lifetime for a specimen of given species without the negative influence of competition (so-called a ‘mean physiological life span’). We consider that as a species-specific parameter, and, therefore, assuming that this value is known. Let the GVP contain n different types of vegetation and Bi (x, y, t) be the density of living biomass of ith type at geographical point (x, y). By analogy with the previous consideration we assume that the dynamics of these variables at any point is described with the following equations:

mR i = mi + % gij Bj

(2.4)

j=1

The value mi = 1/ti, where ti is a mean residence time of carbon in living biomass of ith type of vegetation. As it was mentioned above, these values can be interpreted as a mean life span for specimen of ith type without competition. They are the known species- specific parameters. (d) Competition coefficients gij do not depend on spatial co-ordinates, i.e. gij = const.ij. Then the system (3) is rewritten as

n

n #Bi = Bi gi (T, H)−mi − % gij Bj , #t j=i

i= 1, . . . , n,

(2.5)

i.e. in the form of the standard Lotka–Volterra system.

3. Elementary map (elementary vegetation pattern)

Fig. 1. Elementary map (elementary vegetation pattern).

Let us consider the following elementary map (Fig. 1), which is the superposition of three domains: vi, vj and vk. That part of domain vi

138

& &

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

ni =

1 Si

nij = Fig. 2. The one-dimensional elementary map for two species.

which does not contain any intersections with other domains, is occupied only by ith type of vegetation, the analogous parts of vj and vk are occupied by jth and kth types, respectively. We shall denote them by Vi, Vj and Vk. The domain vij =vi Svj may be occupied by both ith type and jth type, the domain vijk =vi S vj S vk may be occupied by three types of vegetation, etc., if intersections with other domains are considered. We assume that Sijk min{Si, Sj, Sk, Sij, Sik, Skj }, where Sijk, Si, Sj, Sk, Sij, Sik, Sjk are the areas of corresponding domains vijk and Vi, Vj, Vk, vij, vik, vkj. The similar assumptions are valid for the intersections of higher order. Note that the domain vij is a typical transition zone. By the same token, at any point (x, y), belonging to some transition zone, only one or two (not more) types of vegetation can exist. If two types co-exist at all points of the transition zone then it is ‘soft’. The definition is explained by the fact that the biomasse changes within the zone sufficiently smooth, without jumps and discontinuities. The biomass of one species increases during a movement from one to other borders, while the biomass of other species decreases along the same path. A spatial dimension of the soft transition zone is equal to 2. But the other, ‘hard’, type of transition zone could exist as well. This is when only one of two types of vegetation is presented in every point of zone, but the points with other type might exist within its arbitrary small vicinity. Such sort of transition zone may contain a lot of jumps and discontinuities. A spatial dimension of the hard transition zone can change from 1 to 2 (fractal dimensions). On the map the following operators (‘averaging operators’) . . .i and . . .ij may be defined: let be a variable n(x, y) then we have

Vi

1 Sij

n(x, y, t) dx dy and

vij

n(x, y, t) dx dy

(3.1)

It is noted that the GVP can be presented as a join of such elementary maps.

4. One dimensional and two species case Before the consideration of calibration process, i.e. the estimation of the coefficients in (Eq. (2.5)) in a general case, we consider one particular case. This helps us better to understand the general method. Let the vegetation pattern containing two types of vegetation (two species) be ordered along one spatial co-ordinate, x, for instance, along a meridian. Then from (Eq. (2.5)) we have #B1 =B1(o1(x)− g11B1 − g12B2) #t #B2 =B2(o2(x)− g21B1 − g22B2) #t

(4.1)

where oi (x)= oi (T(x), H(x))= gi (T(x), H(x))− mi = gi (x)− mi, i= 1, 2. This is the classic Lotka– Volterra model for two competing species. The corresponding elementary map is shown in Fig. 2. An elementary analysis shows that the system (7) has three non-trivial equilibria: 1.

(1) {B (1) 1 , 0} where B 1 = o1/g11,

2.

(2) (2) {B (2) 1 , 0} where B 1 = (o1/g21 − o2g22)/D, B 2

= (o2g11 − o1g21)/D, D= g11g22 − g12g21, 3.

(3) {0, B (3) 2 } where B 2 = o2/g22.

(4.2)

Their stability domains in the plane {o1, o2}are presented in Fig. 3a,b. We assume (and this is very important!) that F there is the mapping {o1, o2} {x} which maps the positive orthant of the plane {o1, o2} onto the interval of the axis x, so that: F

F

F

domain IV1, domain IIv12, domain III V2. This mapping transforms the axis o1(o2 = 0) onto the point x1 and the axis o2(o1 = 0) onto the

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

point x2. The ray o2 =(g21/g11)o1 is mapped onto the point x12 (for the situation depicted in Fig. 3a and: onto the point x21 (for the situation depicted in Fig. 3b. The ray o2 =(g22/g12)o is mapped onto the point x21 and onto the point x12 for the situations depicted in Fig. 3a,b, respectively. Then we get immediately the following relations:

139

o1(x2)= 0, o2(x1)= 0, (4.3) o2(x12)= (g21/g11)o1(x12), o2(x21)= (g22/g12)o1(x21), for case Fig. 3a, (4.4a) o2(x12) = (g22/g12)o1(x12),o2(x21) = (g21/g11)o1(x21), for case Fig. 3b.

(4.4b)

It is easy to see that the domain II in Fig. 3a corresponds to a soft transition zone and the domain II in Fig. 3b corresponds to a hard one.

5. Soft transition zone

Fig. 3. Stability domains for different equilibria of (4.1). (a) Domain I: equilibrium (B (1) 1 , 0)is stable (globally, i.e., for any positive initial points). Domain II: non-trivial equilibrium (2) (B (2) existing and stable (globally), equilibria 1 , B 2 )is (1) (B 1 , 0)and (0, B2(3)) are unstable. Domain III: equilibrium (0, B2(3)) is stable (globally). (b) Domain I: equilibrium (B (1) 1 , 0)is stable (globally). Domain II: non-trivial equilibrium (2) (1) (B (2) 1 , B 2 )is existing and unstable, equilibria (B 1 , 0)and (3) (0, B2 ) are stable (locally, i.e., not for any initial points). Domain III: equilibrium (0, B2(3)) is stable (globally).

The intervals [x1, x12] and [x21, x2]in Fig. 2 are locations of two ‘biomes’ each of them is presented by only one type of vegetation. We assume (3) that the stable equilibria {B (1) 1 , 0} and {0, B 2 }are realised in the first and second intervals, respectively. The interval [x12, x21] is the transition zone. We assume that it is soft, i.e. the non-trivial (2) equilibrium {B (2) 1 , B 2 } is stable within the zone. The problem is how to calibrate the model and what kind of information is there for this? Really we have the information about the mean value of net annual production, Pi, the living biomass storage, Bi, and the species-specific life span, ti = 1/mi for each ‘biome’, i.e. for its determining species. We assume also that the spatial co-ordinates of x1, x12, x21 and x2 are known. If considering the present state as a steady one then we get from (Eq. (4.1)) under the stationary conditions #Bi /#t =0, i= 1,2: o1 − g11B (1) 1 = 0, o2 − g22B (3) where oi = [P*(x)/B *(x)]−m 2 = 0: i i i= gi (x)− mi ; for i= 1, x[x1, x12], and for i= 2, x[x21, x2], the asterisk means the observed values. In fact that both the production and the biomass may be changed along a biome (if it is large, too) but we consider the case when these values are constant inside each biome, gi (x)= g*= consti, i= 1,2. Then the coefficients of ini traspecific competition, gii, will be mi )/B*= o*/B *, gii = (g*− i i i i

i= 1, 2

(5.1)

Because of areas of the real transition zones are relatively small it is very problematic to get the

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

140

B1(x, t) B*, 1 B2(x, t) 0 for all x[x1, x12] B1(x, t) 0, B2(x, t) B*2 for all x[x21, x2] (2) (2) where B (2) B1(x, t) B (2) 1 , B2(x, t)B 2 1 , B2

" 0 for all x[x1, x12].

Fig. 4. Functions o1(x) and o2(x) (one of possible forms).

corresponding information about these values. Thus, another way will be used: we assume that the function o1(x) linearly decreases within the interval [x12, x2], and function o2(x) linearly increases within the interval [x1, x21] (see Fig. 4). Then o1(x21)=o*1

l2 l1 and o2(x12) = o*2 l2 +l12 l1 +l12

(5.2)

where l1, l2, and l12 are the lengths of intervals [x1, x12], [x21, x2] and [x12, x21], respectively. It is easy to see that conditions (4.3) are fulfilled. We require that the conditions (4.3a) are also fulfilled and use them for calculation of the coefficients of intraspecific competition g12 =

l2 o*1 , B*2 l2 + l12

g21 =

l1 o*2 B*1 l1 +l12

(5.3)

Let the system (Eq. (4.1)) be given on the interval [x1, x2]1), the functions o1(x) and o2(x) are as in Fig. 4 and the coefficients g11, g12, g21, and g22, are calculated by formulas (5.1) and (5.3). Then, in accordance to the construction, for any positive initial data B 01(x), B 02(x) we have for t  : 1

There is not a problem to extend the consideration on the whole axis x; for this we must define the functions o1(x) and o2(x) outside of the interval [x1, x2], for instance, setting them by zero. As far as only the dynamics of transition zone and, more correct, the shift of its border is interesting for us, we restrict ourselves by the consideration of the problem solely on the interval [x1, ,x2].

(5.4)

This means that the GVP is a limit set for the corresponding Lotka–Volterra system. If now local climate parameters are changed (as a result of the Global Change, in accordance with a climate scenario so that T= T(x, t),H = H(x, t))) are known functions of time at any spatial point) then, in the first place, the functions of net primary production will be changed for both species and oi (T(x, t),H(x, t))= gi (T(x, t),H(x, t))mi = gi (x, t)mi,i= 1,2. We assume that the change of these climate parameters does not change the characteristics of inter- and intraspecific competition, i.e. the coefficients g11, g12, g21, and g22, are retained the same. Then the evolution of the GVP will be described by the following non-autonomous system #B1 =B1(g1(x, t)−m1g11B1 − g12B2)+ j1(t), #t #B2 =B2(g2(x, t)−m2 − g21B1 − g22B2)+ j2(t) #t (5.5) where the random perturbations j1(t) and j2(t) (supposed sufficiently small) are used for description of random propagation of generative material. It is necessary since the system (Eq. (5.5)) without perturbations can not leave zero initial state (B 01 = 0 and B 02 = 0 are solutions of Eq. (5.5) under j1 = j2 0), i.e. a species can not spread to a new territory without such sort of mechanism. It is interesting to estimate a shift of borders between different zones under the climate change. We shall suppose that all these changes are sufficiently small, and the linearisation is correct.

5.1. Left border of the transition zone (x12) At the present equilibrium the equality o 02(x 012)= (g21/g11)o 01(x 012) has to be fulfilled. Using Eq. (5.2) and Fig. 4 the relation can be re-written in the form l1o*= k1o*, 2 1 where l1 = l1/l1 + l12 and

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

k1 = (g21/g11). We assume that under some climate change (P/B)-coefficients for these two species are also changed, so that o L1 =o*1 +Do L1 and o L2 = o*2 + Do L2 . The upper index ‘L’ is shown that the change is considered on the left border. In order that the equality l L1 o L2 =k1o L1 will be again fulfilled the borders have to be shifted, so that x12 =x 012 + Dx12 and x21 =x 021 + Dx21. We also assume that the external borders of biomes (points x1 and x2) are not changed and only the borders of transition zone are shifted. Then l L1 = and

l1 + Dx12 Dx12 Dx21 : l1 1 + − l1 +l12 +Dx21 l1 l1 +l12

Dx12 Dx21 1 k1 L − = Do 1 −Do L2 l1 l1 + l12 o*2 l1

We neglected here terms with exponents (Dx12)2, (Dx12Dx21), (Dx12Do L2 ),etc. Since k1 = (g21/g11) then from Eqs. (5.1) and (5.3) o* l1 k1 = 2 o*1 l1 +l12

Dx12 Dx21 Do Do − = − =L l1 l1 + l12 o*1 o*2

(5.7)

Substituting Eq. (5.8) into Eq. (5.6) we get

L 2

(5.8)

5.2. Right border of the transition zone (x21) At this point we have o 02(x 021)= (g22/g12)o 01(x 012) under present condition. It is re-written in the form o*= k2l2o*1 where l1 = l2/l2 + l12 and k2 = 2 (g22/g12). After some climate change o R 1 = o*+ 1 R R Do R , o = o *+ Do (the upper index ‘R’ is shown 1 2 2 2 that the change is considered on the right border) R R and o R 2 = k2l 2 o 1 , where lR 2 =

(5.6)

141

L 1

l2 − Dx21 Dx12 Dx21 : l2 1+ − l2 + l12 − Dx12 l2 + l12 l2

And finally, in the analogous way, we get Dx12 Dx21 Do R Do R 2 1 − = − =R l2 + l12 l2 o*2 o*1

(5.9)

It is easy to see that Eqs. (5.8) and (5.9) form the system of two linear algebraic equations in respect to Dx12 and Dx21, and there is not a problem to solve it. Thus, solving this system we can calculate the shift of the both borders of transition zone. It is interesting to estimate the shifts for some concrete case, for instance, for the transition zone between taiga and steppe.

6. Case study: transition zone between taiga and steppe in central Siberia Let us consider a very idealised presentation of real structure of vegetation along the Yenisei meridian (90° E): the biome of the Southern taiga is located from 59 to 56° N, the transition zone (forest-steppe)-from 56 to 54° N, and the steppe biome-from 54 to 51° N (Fig. 5). Then l1 = 3°, l12 = 2°, l2 = 3°, (all the distances are expressed in degrees: 1°: 111 km). Solving Eqs. (5.8) and (5.9) for these values we get: Dx12 = 4.7L − 2.8R, Fig. 5. Idealised structure of vegetation along the Yenisei meridian (90° E): the biome of the Southern taiga is located from 59 to 56° N, the transition zone (forest-steppe)-from 56 to 54°N and the steppe biome-from 54 to 51°N.

Dx21 = 2.8L − 4.7R

(6.1)

The problem is how to estimate the deviations Do L,R and Do L,R when the temperature and the 1 1 precipitation are changed? There are different ways for this. Let us consider one of them using

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

142

Lieth’s formula for the annual production or for the productivity in Lieth’s formulation (Lieth, 1975). Since oi = (P/B)i −mi =gi −mi, i =1, 2 then Doi g*i : g*− o*i mi i +

n

# ln gi (T, H) #T

# ln gi (T, H) #H

DT

where f(T0)= 1.5× 103 ln[1+ exp(0.119T0 − 1.315)], (6.5)

0

DH

(6.2)

0

where the derivatives are calculated for present temperatures and precipitation on the both borders. There is not a problem to calculate the factors *−m bi = g*/g i i i in Eq. (6.2). In accordance with Bazilevich’s data (Bazilevich et al., 1986) the net primary productions (NPP) are equal to P*1 = 6.35 tons/ha year for the Southern taiga (in basis, Larix czekano6skii ) and P*2 =13 tons/ha year for the steppe (the dominant species are Stipa capillata and Caragana arborescens); the biomasses are equal to B*=125 tons and B*2 = 9 tons, respectively. Natu1 rally assuming thatt1 =180 years for trees (free stand) and t2 = 2 year for perennial grass, then m1 =0.005 and m2 =0.5 (1/year). Thus we have the full information for the calculation of b1 =1.11 and b2 = 1.53. Since there is not any explicit dependence of the biomass on the temperature and precipitation (this dependence is realised by dynamic equations) then we can write that

# ln gi #T

=

0

# ln Pi #T

# ln Pi = #H

=ai and

0

# ln gi #H

(6.3)

=bi

0

!

NPP = 3 min [1− exp( −0.664 × 10 − 3H)], 1 1+ exp(1.35−0.119T)

n"

0.644× 10 − 3 Á , if H0 B f(T0), Ã −3 b1,2 = Íexp(0.644×10 H0)− 1 Ã0, if H0 ] f(T0). Ä (6.6) These coefficients are calculated in the points x 012 and x 021 where T0(x 012)= 0.2°C, H0(x 012)= 600 mm and T0(x 021)= 1.25°C, H0(x 021)= 400 mm, then a1,2(x 012)= 0.093, b1,2(x 012)= 0, a1,2(x 021)= 0, −3 0 b1,2(x 21)= 2.192×10 . Then L= (b1 − b2)[a1,2(x 012)DT(x 012)+ b1,2(x 012)DH(x 012)] = − 0.039DT(x 012), R=(b2 − b1)[a1,2(x 021)DT(x 021)+ b1,2(x 021)DH(x 021)] = 0.92× 10 − 3DH(x 021).

(6.7)

And finally, using Eq. (6.1) we get Dx12 = − 0.183DT(x 012)− 2.58× 10 − 3DH(x 021), Dx21 = − 0.11DT(x 012)− 4.32× 10 − 3DH(x 021). (6.8)

0

For estimation of these derivatives we use Lieth’s formula, which connect the NPP with annual temperature T (in °C) and precipitation H (in mm). We will use it in some different (not principally) form:

0.119 Á , if H0 \ f(T0), Ã a1,2 = Í1+exp(0.119T0 − 1.315) Ã0, if H0 5 f(T0) Ä

[kg dry matter] m2 ×year (6.4)

Differentiating (Eq. (6.4)) in respect with T and H we get

In accordance with the standard scenario of climate change for CO2-doubling (Petoukhov et al., 1999) in this region the annual temperature is expected to rise by 5°C and the precipitation increases by 10% (so that DT(x 012)= 5°C and DH(x 021)= 40 mm). Then both borders will be shifted northwards, where the border between taiga and ‘forest-steppe’ zones will be shifted by 1.018° : 113 km, and the border between ‘forest-steppe’ and steppe zones will be shifted by 0.723 : 80.25 km; as a result, the width of transition zone between taiga an steppe increases by 33 km under the climate change. 7. Hard transition zone We consider again that the intervals [x1, x12] and [x21, x2] in Fig. 2 are locations of two

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

143

B1(x,t) B*, 1 B2(x,t) 0 for all x[x1, x12], B1(x, t)0, B2(x, t)B*2 for all x[x21, x2]. (7.3) Relating to the situation in the transition zone [x12, x21] it depends on the values of initial data, and one of the equlibria {B*, 1 0} or {0, B*}can 2 be realised. This means that domains representing different ‘biomes’ (but not transition zones between them) of the GVP are limit sets for corresponding Lotka–Volterra system. It is interesting, what kind of structure will have the transition zone [x12, x21]? In order to answer the question we carry on the following numerical experiment. Let B*= B*= 1, o*= o*= 1, l*= 1 2 1 2 1 l*= l * = 1. Then, calculating g , g , g , and g22, 2 12 11 12 21 by Eqs. (7.1) and (7.2) we get g11 = g22 = 1, g12 = g21 = 2. The basic equation will be

Fig. 6. ‘Hard’ transition zone for two different amplitudes of random perturbations (ra = 0.006, rb = 0.008). The ‘fractal’ structures of zones are shown. Solid: 1st species, Dotted line: 2nd species. (1) 1

‘biomes’, and that the stable equilibria {B , 0} and {0, B (3) 2 }are realised in the first and second intervals, respectively. The interval [x12, x21] is the transition zone, but, on the contrary to the previous section, the zone is ‘hard’, i.e. the non-trivial (2) equilibrium {B (2) 1 , B 2 } is unstable within the zone. Further all the considerations are as in the previous section, so that gii = (g*−m o*/B *, i i )/B*= i i i

i =1, 2

(7.1)

The functions o1(x) and o2(x) are as in Fig. 4. Only one difference is between this and the previous section: conditions (4.3) must be fulfilled (not 4.3a). From this we get g12 =

o*1 l1 + l12 , B*2 l1

g21 =

o*2 l2 +l12 B*1 l2

(7.2)

Let the system (4.1) be given on the interval [x1, x2], the functions o1(x) and o2(x) are as in Fig. 4 and the coefficients g11, g12, g21 and g22 are calculated by Eqs. (7.1) and (7.2). Then, in accordance to the construction, for any positive initial data B 01(x), B 02(x) we have for t  :

#B1 =B1(o1(x)− B1 − 2B2)+ j1(t), #t #B2 = B2(o2(x)− 2B1 − B2)+ j2(t) #t

(7.4)

where Á Ã1, o1(x)= Í1 Ã (3− x), Ä2 Á Ã 1, o2(x)= Í1 Ã x, Ä2

if x[0,1] , if x[1,3] if x[2, 3] if x[0, 2]

and random perturbations, j1(t) and j2(t) are sufficiently small. The initial distributions of B 01(x) and B 02(x) are similar to o1(x) and o2(x). The stable distributions B1(x) and B2(x) (the solution of Eq. (7.4) for t ) are shown in Fig. 6a,b for two different amplitudes of random perturbations. It is interesting that the borders becomes ‘vague’ when perturbations increase (the ‘patches’ with other species arises within the area of basic species). Based on the reasoning of the previous section we can prove that the shift in zones does not depend on the type of borders between them: soft or hard ones. Therefore, we can use the formulas

144

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

of the previous section in order to estimate the shift values for vegetation zones caused by the climate change. However, it should be noted that the dynamics of the GVP depends on what kind of transition zone, ‘soft’ or ‘hard’, there is. It is natural since the coefficients g12 and g21 will be distinguished for these cases. This is a typical non-uniqueness, and in order to avoid it we have to use a priori information about the type of transition zone. Note that in case of a ‘hard’ transition zone, the term ‘borders’ looses its original sense as well as their shift does. Now there are not well outlined borders between two species: among an areal, which is completed by individuals of the first species, an individual (or, a small group) of the second species may occurs, and vice versa. You can see this very visually in Fig. 6b. Apparently, in this case, speaking about the shift of transition zone, we must imply the shift (as the whole) of some fractal pattern, the shift of the whole picture.

Applying to the both parts of (8.1) the averaging operator (3.1) we get g* i i − mi − gii B* i i= 0 and gii =

g* i i − mi B* i i

(8.2)

Secondly, for intermediate domains vij we assume that there is only one stable equilibrium {0, . . . , 0, B*, i 0, . . . , 0, B*, j 0, . . . , 0}. Then for any points (x, y)vij we have g*(x, y)− mi − gii B*(x, y)− gij B*(x, y)= 0, i i j g*(x, y)− mj − gji B*(x, y)− gjj B*(x, y)= 0. j i j (8.3) Applying to the both parts of (8.2) the averaging operator (3.1) we get g* i ij − mi − gii B* i ij − gij B* j ij = 0, g* j ij − mj − gji B* i ij − gjj B* j ij = 0, and 1 {g* i ij − mi − gii B* i ij }, B* j ij 1 {g* gji = j ij − mj − gjj B* j ij }. B* i ij gij =

8. Calibration of a general dynamic model If coming back to the system (2.5) then immediately we collide with a problem: how to determine n(n +1) coefficients gij, mi ; i, j =1,...,n and n concrete functions gi (T, H). How to decide the problem? Firstly, we assume that there is only one stable equilibrium {0,...,0,B*,0,...,0}in Vi (see Section 3 i and Fig. 1). Then we get from Eq. (2.5) for any points x, yVi : g*(x, y)− mi −gii B*(x, y) =0 i i

We formulate the following Theorem. Let gij, i, j= 1, . . . , n be calculated accordingly with Eqs. (8.2) and (8.4) and B. (x, y, t)= {B. 1(x, y, t), . . . , B. n (x, y, t)} be the solution of the system

n

n #B. i = B. i g*− mi − % gijB. j , i #t j=1

i= 1, . . . , n (8.5)

(8.1)

where g*(x, y) = g*(T*(x, y), H*(x, y)). Since the i i values of net primary production P*(x, y) and i equilibrium biomass B*(x, y) inside V are i i known, by the same token we know P/B-coefficients, i.e. the function g*(x, y) = P*(x, y)/ i i B*(x, y) is also known. In addition, we assume i that the value of mi =1/ti, where ti is a mean life span for specimen of ith type of vegetation, is the known species-specific parameter and it is constant inside Vi.

(8.4)

1. Then for any positive initial values B. (x, y, t) B. *(x, y) when t . B. *=

!

{0, . . . , 0, B. *, if x, yVi ; i 0, . . . , 0} . *, {0, . . . , 0, B. *, if x, yvij. i 0, . . . , 0, B j 0, . . . , 0} (8.6)

. * . * 2. B. * i i = B* i i, B i ij = B* i ij, B j ij = B* j ij,

i, j= 1, . . . , n.

(8.7)

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

The theorem is intuitively clear, but this is not a mathematically rigorous proof. Nevertheless, computer experiments give arguments in favour of the theorem.In other words, when t  , the solution of (Eq. (8.5)) converges on average to the real map of vegetation, i.e. to the map, which was used for estimation of coefficients. This theorem gives us the effective method for the constructing of the dynamic GVP; i.e. sequences of maps, which show the evolution of GVP, caused by the global change. How to do this? In order to remain in the frameworks of our model we have to assume that both the coefficients mi and the competition coefficients gij must be constant, only the functions gi must be changed with the climate change (which is described in our model by the scenario dynamics of annual temperature and precipitation). Since gi are simply P/B-coefficients then only the numerator in this ratio depends on T and H explicitly. The biomass depends on T and H implicitly, by means of dynamic equations. Therefore, we can present the dependence of gi on T and H in the form:

!     "

gi (T, H)

# ln P * DH(x,y,t) #H

climate change. It does not take into consideration such characteristics as soil properties, different physiological parameters, etc., which naturally influences on the GVP dynamics. By the same token, the model is very distinguished from other simulation models of the global vegetation (see, for instance, Cramer and Field, 1999). However, in spite of its external simplicity, it possesses a very complex behaviour with multiple equililibria, hysteresis and fractal topology. We can expect that the model can describe such sufficiently complex spatial structures as the GVP. There is an additional attractive property: the simplicity of requirements to basic information, which is needed for its calibration. At last, I would like to note that in the model an influence of different factors on vegetation dynamics is concentrated in the function, which describes a dependence of production on them. We used Lieth’ formula depending on the annual temperature and precipitation. However, nothing hinders us to use another formula, which contains some dependence of production, for instance, on soil properties at given loci.

Acknowledgements

# ln P * :g*i 1+ DT(x, y, t) #T +

145

(8.8)

For calculation of the logarithmic derivatives (# ln P/#T) and (# ln P/#H) the Lieth formula (see Eq. (6.4)) can be used. The values of DT(x, y, t) and DT(x, y, t), i.e. deviations of the future temperatures and precipitation are given by one of possible scenarios of the climate change. Finally, the dynamics of GVP will be described by the non-autonomous dynamical system (2.5) where the function gi (T, H) is given by Eq. (8.8). 9. Discussion At first glance, the suggested model seems too abstract and ‘theoretical’ in order to apply it to simulate a real dynamics of the GVP under the

I am grateful to S. Venevsky and S. Sitch for their helpful comments and criticism. I am indebted also to W. von Bloh for his help in the processing of the manuscript.

References Bazilevich, N.I., Grebentshikov, O.C., Tishkov, A.A., 1986. Geographical Regularities of the Ecosystems Structure and Functioning. Nauka, Moscow. Churkina, G., Svirezhev, Y., 1995. Dynamics and forms of ecotone under the impact of climate change: mathematical approach. J. Biogeog. 22, 565 – 569. Cramer, W., Field, C.B., (Eds), 1999. The Potsdam NPP Model Intercomparison. Global Change Biology, 5, Suppl. 1: 1 – 76. Ebenhoh, W., 1994. Competition and coexistence: modelling approach. Ecol. Model. 75-76, 83 – 95. Ekschmitt, K., Breckling, B., 1994. Competition and coexistence: the contribution of modelling to the formation of ecological concept. Ecol. Model. 75-76, 71 – 82.

146

Y. S6irezhe6 / Ecological Modelling 135 (2000) 135–146

Huntingford, C., 1999. Personal communication. Jesse, K.J., 1999. Modelling of a diffusive Lotka–Volterra system: the climate induced shifting of tundra and forest realms in North-America. Ecol. Model. 123, 53–64. Lieth, H., 1975. Modelling the primary productivity of the world. In: Lieth, H., Whittaker, R.H. (Eds.), Primary Productivity of the Biosphere. Springer-Verlag, New York, pp. 237 – 263. Petoukhov, V., Ganopolsky, A., Brovkin, V., et al., 1999. CLIMBER – 2: a climate system model of intermediate complexity. Part I: model description and performance for present climate. Climate Dyn. 15, 950–966. Prentice, I.C., Cramer, W., Harrison, S.P., Leemans, R., Monserud, R.A., Solomon, A.M., 1992. A global biome model

based on plant physiology and dominance, soil properties and climate. J. Biogeog. 19, 35 – 57. Svirezhev, Y.M., 1994. Non-linear Problems in Mathematical Ecology. WP-94-71, IIASA, Laxenburg. Svirezhev, Y.M., 1999. Simplest dynamic models of the Global Vegetation Pattern. Ecol. Model. 124, 131 – 144. Svirezhev, Y.M., Logofet, D.O., 1983. Stability of Biological Communities. Mir, Moscow. Tilman, D., 1982. Resource Competition and Community Structure. Monograph in Population Biology. Princeton University Press, Princeton NJ. Vandermeer, J., 1996. Disturbance and neutral competition theory in rain forest dynamics. Ecol. Model. 85, 99 – 111.

.