Alexandra Skripchenko

Serge Troubetzkoy
email: serge.troubetzkoy@univ-amu.fr

POLYGONAL BILLIARDS WITH ONE SIDED SCATTERING

We study the billiard on a square billiard table with a one-sided vertical mirror. We associate trajectories of these billiards with double rotations and study orbit behavior and questions of complexity.

Introduction

The table Π we consider consists of the square [0, 1/2] × [0, 1/2] with a vertical wall I connecting the points (a, 0) and (a, b). We play billiard on this table, with I acting as a one-sided mirror. That is, we consider a point particle in Π and a direction θ in S 1 , the particle moves at unit speed in the direction until it reaches the boundary of the table . If it arrives at the left side of I it passes through it Figure 1. The green side is transparent while the blue side is reflective.

unperturbed, while if it arrives at the right side of I or at the boundary of the square it is reflected with the usual law of geometric optics, the angle of incidence equals the angle of reflection (see Figure 1). Polygonal billiards with one-side straight mirrors were briefly described by M. Boshernitzan and I. Kornfeld [BK] in connection with a special kind of piecewise linear mapping of a semi-interval, called interval translation mappings (ITMs). Interval translation mappings are a natural generalization of interval exchange transformations.

In this article we prove that billiard flow on such table can be associated with special interval translation mappings called double rotations (Proposition 1). The term double rotations was introduced in [SIA], they have also been studied in [BT], [B], [BC]. We show that up to a natural involution there exists a bijection between double rotations and billiard map we work with, and therefore almost all of our billiard map are of finite type due to the corresponding results for double rotations (Theorem 2 part 2)) (see [SIA] and [BC]). In the other parts of Theorem 2 we collect various interesting implications of this result on unique ergodicity, non-unique ergodicity, and the Hausdorff dimension of the attractor (concretizing a suggestion of Boshernitzan and Kornfeld [BK]). Our main result is an exact linear formula for complexity of billiard trajectories in a given direction in the case a = 1 4 (Theorem 4), we also give a linear estimate in case of other rational values of a (Theorem 9). The main result also generalizes to certain other rational polygons with one sided scatterers located at an axis of symmetry (Theorem 7). The proof is based on an extension of combinatorial arguments introduced by J. Cassaigne ([C]) for languages with bispecial words (Theorem 11) and also uses a certain symmetry of orbits of the unfolded billiard.

The results

A double rotation is a map T : [0, 1) → [0, 1) of the form:

T y = y + α (mod 1) if y ∈ [0, c) y + β (mod 1) if y ∈ [c, 1).
Consider the billiard in the table Π described in the Introduction. There is a well known construction of unfolding the billiard in a rational polygon to a translation surface (see for example [MT]). The same construction applied to our setting yields a unit torus consisting of four copies of the billiard table Π with slits corresponding to I which are identified according to the billiard flow: when we hit a right copy of I (depicted in blue) we jump to the corresponding left copy of I (depicted in green), while the left (green) copies of I are transparent. With these identifications our billiard is equivalent to the linear flow on the unit slitted torus (see Figure 2). Consider the first return map to the vertical section x = a in the unfolded slit torus (see Figure 2, where this section is drawn by yellow dashed line). We denote this map by T (a,b,θ) : {(x, y) : x = a} → {(x, y) : x = a}. We identify its domain of definition with the circle [0, 1). It is easy to check that there are only two possibilities for our trajectories (here we assume that the the trajectory is not vertical: θ ∈ {± π 2 }): either the orbit of the point hits the blue wall, and immediately goes back to the section (this is depicted by the orange line in Figure 2); or the orbit of the point does not touch the blue wall (such an orbit is depicted by the purple line in Figure 2). In the first case T (a,b,θ) (y) = (1 -2a) tan θ + y (mod 1); in the second case T (a,b,θ) (y) = tan θ + y (mod 1). Therefore, we have the following:

Proposition 1. First return map on the vertical section x = a is a double rotation with the following values of parameters:

α = (1 -2a) tan θ (mod 1), β = tan θ (mod 1), c = b.
Let us consider the cubic polynomial P (x) = x 3x 2 -3x + 1 and let γ be the unique root of P (x) in [0, 1]. Let us consider the following values of parameters of our billiard:

a = (1 -γ) 2 , b = γ θ = arctan γ.
The map T (a,b,θ) is then the double rotation introduced by Boshernitzan and Kornfeld in [BK], the first example of an ITM whose dynamics differs from that of an interval exchange. To explain this difference, and to collect the most interesting implications of known results to our situation, we need to introduce some notation.

For any (α, β, c) the attractor of T is the set

Ω := J ∩ T J ∩ T 2 J ∩ • • • , where J = [0, 1). If there exists n ∈ N such that J ∩ T J ∩ T 2 J ∩ • • • ∩ T n J = J ∩ T J ∩ T 2 J ∩ • • • ∩ T n+1
J then we say that T is of finite type. Otherwise T is of infinite type. Informally, the interval translation map being of finite type means that it can be reduced to interval exchange transformations; if it is of infinite type, then the attractor is a Cantor set. The Boshernitzan-Kornfeld example given above is of infinite type. Fix the parameters a, b and a direction θ. Let T := T (a,b,θ) be the associated double rotation. Let X 0 := [0, c) and X 1 := [c, 1) be the two intervals of continuity of T . The code ω(y) ∈ {0, 1} N of the T orbit of y is the sequence of intervals it hits, i.e. w(y) n = i if and only if T n y is in the interval X i . The language L := L a,b,θ is the set of all infinite codes obtained as y varies, and L(n) is the set of different words of length n which appear in L. Let p(n) := p a,b,θ (n) := #L a,b,θ (n). Note that one could also consider p ∞ (n) ≤ p(n) the number of different words of length n which appear in L as x varies in the attractor.

We must exclude directions for which there is a billiard orbit from an end point of I to an end point of I. We call such directions exceptional. There are at most a countably many exceptional directions since for any positive constant N , there is a finite number of billiard orbits which start and end at end points of I and have length at most N .

The following theorem follows directly from several known results about double rotations, interval translation mappings and piecewise continious interval maps.

Theorem 2. 1) For all (a, b, θ) with θ non-exceptional the billiard/double rotation is minimal.

2) For almost all (a, b, θ) the double rotation is of finite type.

3) There exists an uncountable set of (a, b, θ) so that the Hausdorff dimension of the closure of the attractor is zero, in particular the double rotation is of infinite type. 4) For each (a, b, θ) with θ non-exceptional the billiard/double rotation has at most two ergodic invariant measures. 5) There exists an uncountable set of (a, b, θ) with θ non-exceptional such that the billiard/double rotation is not uniquely ergodic. 6) For all (a, b, θ) with θ non-exceptional the complexity p(n) grows at most polynomially with degree 3.

In [BK] it was suggested that there may exist configurations on a rational billiard table with mirrors which force the light to get concentrated in some arbitrarily small portions of the table. This suggestion is confirmed by Part 3) of the theorem.

Proof. It is easy to check that varying the parameters (a, b, θ) we can obtain all double rotations with only one restriction: α ≤ β. There exists an involution between two parts of parameter space (α ≤ β and β ≤ α) and the orbit behavior of ITM from these two parts are completely the same: the involution is measure preserving and does not change the dynamics. Thus all known results on double rotations hold in our setting. 1) We think of a double rotation as an ITM defined on an interval. First suppose that this interval is [0, 1) (in the coordinates of the definition of double rotation). Either this is an ITM on with 3 intervals of continuity, or it has 4 intervals of continuity. In the later case we choose the origin to be the point c, and the double rotation now has at most 3 intervals of continuity, thus we can always choose coordinates so that the ITM has at most 3 intervals and part 1) follows from Theorem 2.4 of [ST]. Part 2) follows from [SIA] Theorem 4.1 or from Theorem 1 of [BC], part 3) follows from Theorem 10 of [BT], part 4) follows from Theorem 3 of [BH], part 5) follows from Theorem 11 of [BT], and part 6) follows Theorem 1 of [Ba] (see also Corollary 8 of [BH]).

The billiard flow φ t is defined on the phase space Π := Π × S 1 , with proper identifications on the boundary, and the attractor of the billiard flow is

t≥0 φ t Π.
An immediate application of the Fubini Theorem to Theorem 2 part 2) yields Corollary 3. For almost every a, b the billiard attractor has full measure.

Question: For which polygons with one sided mirrors does the billiard attractor have full measure? Positive measure? Zero measure?

The main results of this article are improvements of part 6) of the theorem. We begin with the case a = 1 4 , i.e. the one sided mirror is in the middle of the square, where we get a complete description of the complexity. We define three sets of directions, in the below definitions all the parameters k i and l i are integers.

Ã1 := θ : ∃k 1 , k 2 , l 1 , l 2 : ∀ k1 ≤ k 1 , k2 ≤ k 2 , l1 ≤ l 1 , l2 ≤ l 2 : tan θ ∈ 2l 1 1 + 2k 1 , 2l 1 + 4b 2k 1 + 1 ∩ k1, l1 l1 + 2b k1 + 1 , l1 + 1 k1 + 1 ∩ 2l 2 + 2 -4b 1 + 2k 2 , 2l 2 + 2 2k 2 + 1 ∩ k2, l2 l2 + 1 k2 + 1 , l2 + 2 -2b k2 + 1 , Ã2 := θ : ∃k 1 , k 2 , l 1 , l 2 : ∀ k1 ≤ k 1 , k2 ≤ k 2 , l1 ≤ l 1 , l2 ≤ l 2 : tan θ ∈ 2l 1 1 + 2k 1 , 2l 1 + 4b 2k 1 + 1 ∩ k1, l1 l1 + 2b k1 + 1 , l1 + 1 k1 + 1 ∩ 1 + l 2 -2b 1 + k 2 , 1 + l 2 k 2 + 1 ∩ k2, l2 2 l2 k2 + 1 , l2 + 2 -4b k2 + 1 } θ : ∃k ′ 1 , k ′ 2 , l ′ 1 , l ′ 2 : ∀ k′ 1 ≤ k ′ 1 , k′ 2 ≤ k ′ 2 , l′ 1 ≤ l ′ 1 , l′ 2 ≤ l ′ 2 : tan θ ∈ l ′ 1 + 1 1 + k ′ 1 , 1 + l ′ 1 + 2b k ′ 1 + 1 ∩ k′ 1 , l′ 1 2 l′ 1 + 2b 2 k′ 1 + 1 , 2 l′ 1 + 2 2 k′ 1 + 1 ∩ 2l ′ 2 + 2 -4b 1 + 2k ′ 2 , 2l ′ 2 + 2 2k ′ 2 + 1 ∩ k′ 2 , l′ 2 l′ 2 + 1 k′ 2 + 1 , l′ 2 + 2 -2b k′ 2 + 1 Ã3 := θ : ∃k 1 , k 2 , l 1 , l 2 : ∀ k1 ≤ k 1 , k2 ≤ k 2 , l1 ≤ l 1 , l2 ≤ l 2 : tan θ ∈ l 1 + 1 1 + k 1 , 1 + l 1 + 2b k 1 + 1 ∩ k1, l1 2 l1 + 2b 2 k1 + 1 , 2 l1 + 2 2 k1 + 1 ∩ 1 + l 2 -2b 1 + k 2 , 1 + l 2 k 2 + 1 ∩ k2, l2 2 l2 k2 + 1 , l2 + 2 -4b k2 + 1 , For i = 1, 2, 3 let A i := Ãi \ {θ exceptional}.
The set are pairwise disjoint and

3 i=1 A i = S 1 \ {θ exceptional}.
We give a geometric description of these sets. We work in the unfolded model (Figure 2) and consider the two singular points of the forward billiard map. We consider the backward trajectories of these two singular points, the sets are distinguished by the fact which part of the boundary of the semi-transparent mirror is first touched by the backwards trajectory, the transparent side or the reflecting side. This identifies how many preimages a corresponding point of the orbit has (none or two). In particular, A 3 corresponds to the case when both of the trajectories first touched the reflecting side, A 1 to the case when both of the trajectories first touched the transparent side and A 2 is the intermediate case, one trajectory first touches the reflecting side and the other tranjectory first touches the transparent side.

Theorem 4. Suppose a = 1 4 and that θ is non-exceptional. Then either, θ ∈ A 1 , and p(n) = n + 1 for all n ≥ 0, or there exists a positive integer constant C θ so that

(1) for θ ∈ A 2 , p(n) = 2n -C θ for all sufficiently large n, or (2) for θ ∈ A 3 , p(n) = 3n -C θ for all sufficiently large n.

Remark, in cases 1 and 2 the double rotation T θ is of infinite type where by T θ we mean T 1/4,b,θ for any arbitrary b. The behavior of the complexity for small n will also be described in the proof. Since the map T θ is minimal, we can apply a theorem of Boshernitzan, a minimal symbolic system with lim n→∞ sup p(n) n < 3 is uniquely ergodic ( [FM][Theorem 7.3.3]) to conclude that

Corollary 5. If θ ∈ A 1 ∪ A 2 then T θ is uniquely ergodic.
Note that Ã2 is open, and exceptional directions are countable, thus we have shown that for all but countably many θ in an open set of θ the map T θ is of infinite type, uniquely ergodic and of linear complexity. Corollary 6. In the case a = 1 4 , the billiard attractor has positive, but not full measure.

Except for the exact computation of the sets A i , Theorem 4 is a special case of a more general result. Take any rational polygon, reflect it in one side. Erase part of the side, make the other part a one sided mirror (which will be part of the line x = 0 to be concrete), to produce a table P . Consider the slitted flat surface M := M (P ) associated with P (see for example [MT]) with the slits identified as in the square case. For the moment consider P without the one-sided mirror, and M without the slits, let mi ni π be the angles between the sides of P , and let N be the least common multiple of the n i . Then a standard computation shows that M (P ) has R := N 1 ni vertices (not counting the endpoints of the slits) [MT]. Consider the section x = 0 as a subset of P . In M there are 2N copies of this section. If we fix θ non-exceptional (no orbit from any vertex of P to any vertex of P ) then there are N copies of the section for which the linear flow on M jumps via an identification (the linear flow goes through the other N copies of the section as if they where not there). We consider the first return map T θ to the N copies which produce a jump. The top and bottom of the section have already been counted as vertices of P . The map T θ is an interval translation map with R + N intervals of continuity, the R coming from the R singular points on M , and the N from the end point of the one-sided mirror (which we assume starts at the bottom of P ).

We assign a symbol to each of these intervals and code the billiard orbit by these R symbols. Let p(n) be the complexity of the associated language.

Theorem 7. If P as above and θ is non-exceptional then there exists a positive constant

C θ and k ∈ {0, 1, 2, • • • , R + N } so that p(n) = (R + N -1 + k)n -C θ for all sufficiently large n.
Example 1: In the double hexagon (Figure 3) all angles have n i = 3, so N = 3, R = 18. In this case all m i = 1, so we do not have any removable singularities. Example 2: In the U shaped figure (Figure 3) one can calculate N = 2, R = 12. However, on M only 4 of these singularities are actual singularities (the two points with angle 3 2 π and the two copies of the end point of the one sided mirror on M ), the others are removable singularities. We also remark that in the original case of the square N = 2, R = 4 and there are 4 removable singularities.

We remind the reader of two results: an aperiodic ITM with r intervals of continuity has at most {r/2} minimal components (Theorem 2.4 in [ST]) and if K ≥ 3 is an integer, then a minimal symbolic system satisfying limsup n→∞ p(n) n < K admits at most K-2 ergodic invariant measures (Theorem 7.3.4 in [FM]). Combining these two results with Theorem 7 yields Corollary 8. There are at most [(R + N )/2] minimal components and on each minimal component the number of ergodic invariant measures for the ITM/billiard is at most 2N + 2R -2.

We turn to the question of complexity for rational position of the mirror in the square case, i.e. a ∈ Q. Suppose a = p q . We consider the return map to the vertical sections x = k q with k ∈ {0, 1, . . . q -1}. The first return map S to these vertical sections is given by S(i, y) = (p, y + tan θ/q (mod 1)) if i = qp and y ∈ [-b, b], (i + 1 (mod q), y + tan θ/q (mod 1)) otherwise.

We code orbits of this map by a 3 letter alphabet, 2 letters for the section x = a (where the map is discontinuous) and a third letter for all the other sections.

Theorem 9. If a = p q , θ is non-exceptional then p θ (n) ≤ (2 + 2q)n for all n. Theorem 2.4 in [ST] and Theorem 7.3.4 in [FM] imply that Corollary 10. There are at most q minimal components and on each minimal component the number of ergodic invariant measures for the ITM/billiard is at most 1 + 2q.

Cassaigne's formula

The main technical tool will be a variant of Cassaigne's formula [C]. Let A be a finite alphabet, L ⊂ A N be a language, L(n) the set of words of length n which appear in L, and p(n) := #L(n). Note that p(0) = #{∅} = 1. For any n ≥ 0 let s(n) := p(n + 1)p(n), and thus We remark that while m r (u) ≥ 1 the other two quantities can be 0.

p(n) = 1 + n-1 i=0 s(i).
A word u ∈ L(n) is called left special if m l (u) > 1, right special if m r (u) > 1 and bispecial if it is left and right special. Let BL(n) := {u ∈ L(n) : u is bispecial}. Let L np (n) := {v ∈ L(n) : m l (v) = 0}.
In this section we show that Theorem 11.

s(n+1)-s(n) = v∈BL(n) m b (v)-m l (v)-m r (v)+1 - v∈Lnp(n):mr(v)>1
m r (v)-1 .

Remark: Cassaigne proved this theorem in the case of recurrent languages (i.e. L np (n) = ∅) [C] (see [CHT] for a English version of the proof). We use the same strategy of proof.

Proof. Since for every u ∈ L(n + 1) there exist b ∈ A and v ∈ L(n) such that u = vb we have

s(n) = u∈L(n) (m r (u) -1). Thus s(n + 1) -s(n) = u∈L(n+1) m r (u) -1 - v∈L(n) m r (v) -1 .
Now, for u ∈ L(n + 1) we can write u = av with a ∈ A and v ∈ L

(n). Let L p (n) := {v ∈ L(n) : m l (v) ≥ 1}. Thus s(n + 1) = v∈Lp(n)   av∈L(n+1) m r (av) -1   . Let L np (n) := L(n) \ L p (n), and thus s(n + 1) -s(n) equals v∈Lp(n)   av∈L(n+1) m r (av) -1 -m r (v) -1   - v∈Lnp(n) m r (v) -1 .
For any word v ∈ L p (n) with av ∈ L(n + 1) any legal prolongation to the right of av is a legal prolongation to the right of v as well thus if m r (v) = 1 then m r (av) = 1. Thus words with m r (v) = 1 do not contribute to any of the above sums. Thus s(n + 1)s(n) is equal to the above sum restricted to those v such that m r (v) > 1. For the left sum, if furthermore m l (v) = 1 then there is only a single a such that av ∈ L(n + 1). For this a we have m r (av) = m r (v) thus such words do not contribute to the left sum. Thus the only terms which contribute to the left sum are the bispecial words, and to the right the words for which m r (v) > 1; in other words s

(n + 1) -s(n) equals v∈BL(n)   av∈L(n+1) m r (av) -1 -m r (v) -1   - v∈Lnp(n):mr(v)>1 m r (v) -1 . For any v ∈ BL(n) we have m b (v) = av∈L(n+1) m r (av) and m l (v) = av∈L(n+1) 1, thus s(n + 1) -s(n) equals v∈BL(n) m b (v) -m l (v) -m r (v) + 1 - v∈Lnp(n):mr(v)>1 m r (v) -1 .

The proofs

Proof of Theorem 4. We use Theorem 11. In our case m r (v) = 1 or 2, so the second term of the equation reduces to

#{v ∈ L np (n) : m r (v) > 1}, thus s(n + 1) -s(n) equals v∈BL(n) m b (v) -m l (v) -m r (v) + 1 -#{v ∈ L np (n) : m r (v) > 1}.
We have p(0) = 1 and p(1) = 2. Suppose that y has two preimages, and that T n θ (y) is in the boundary of one of the two intervals of continuity of T θ , i.e. the billiard orbit of y arrives at the top of the reflecting side of I after n steps. Consider the code v of length n. Clearly m l (v) = m r (v) = 2. Since θ is non-exceptional then m b (v) = 4. (Note that if θ is exceptional then m b (v) = 3 and thus the orbits would not contribute to the sum.) Thus for non-exceptional directions s(n + 1)s(n) = #BL(n) -#{v ∈ L np (n) : m r (v) > 1}.

We switch back and forth between the language of double rotations and that of the billiard. Let {e, f } be the two (common) endpoints of X 0 and X 1 (in the original definition of double rotations these points are called 0 and c). Let {e

(n) i } := T -n
θ (e) be the collection of nth-preimages of e which we will denote by T e (n) for short. Then we will define the tree T e of preimages of e to be the set T e = ∪ n≥0 T e (n), with a directed arrow from e n i to e n-1 j if T θ e n i = e n-1 j . We define similarly the preimage tree T f . Any right special word corresponds to a billiard orbit which hits e or f , thus we can decompose s(n + 1)s(n) into two parts, those words contributing to this difference corresponding to a billiard orbit arrives at e, and those which arrive at f ; we note these two contributions by (s(n + 1)s(n)) e = #BL e (n) -#{v ∈ L np (n) : m r (v) > 1} e for the point e, and similarly for the point f . This formula can be seen as counting the number of leaves at level n of a weighted tree, a vertex e n i has weight +1 if it has two preimages (contributing one element to BL e (n)), weight 0 if it has one preimage, weight -1 if it has no preimages (contributing one element in L np (n)) (Figure 4).

Proof of Proposition 12. By assumption there are no right special words, thus no bispecial words. Thus s(n) = const. We have s(0) = p(1)p(0) = 3 -1 = 2, thus s(n) = 1 and p(n + 1) = p(n) + 1. It follows that p ∞ (n) = n + 1.

5.

Acknowledgements.

Figure 2 .

 2 Figure 2. First return map T θ is a double rotation

Figure 3 .

 3 Figure 3. Polygons with linear complexity.

For

  u ∈ L(n) let m l (u) := #{a ∈ A : au ∈ L(n + 1)}, m r (u) := #{b ∈ A : ub ∈ L(n + 1)}, m b (u) := #{(a, b) ∈ A 2 : aub ∈ L(n + 2)}.

We thank Pascal Hubert for fruitful discussions and Henk Bruin for a critical reading of an earlier version of this manuscript. We gratefully acknowledge the support of ANR Perturbations and ERC Starting Grant "Quasiperiodic". We also thank the anonymous referee for several useful suggestions which improved our presentation.

For any n such that each of the e n j has a unique preimage, we have #BL e (n) = #{v ∈ L np : m r (v) > 1} e = 0 and thus (s(n 0 + 1)s(n 0 )) e = 1.

Now consider an n so that at least one point e n-1 j has two preimages. The first time this happens, n 0 := n e 0 , if this never happens then set n 0 := +∞. Thus #BL e (n 0 ) = 1 and #{v ∈ L np (n 0 ) : m r (v) > 1} e = 0. This implies (s(n 0 + 1)s(n 0 )) e = 1. By the symmetry, the next time, n e 1 one of these backwards orbits hits the mirror, the other one will also hit the mirror from the other side. Thus one of the orbits will have two preimages and the other will be non-left-extendable. Thus #BL e (n e 1 ) = #{v ∈ L np (n a 1 ) : m r (v) > 1} e = 1. The same holds for all times n e k for which #BL e (n e k ) > 0. We have thus shown that n e 0 is the unique time for which (s(n + 1)s(n)) e = 0. The same holds for (s(n + 1)s(n)) f = 0, but n e 0 and n f 0 are not necessarily equal.

If n e 0 = n f 0 =: n 0 then

By definition p(0) = 1 and thus s(0

In particular if n 0 = ∞, i.e. if the preimages of both points e and f disappear before being doubled, then

It is easy to check that n 0 = ∞ happens exactly when θ ∈ A 1 .

On the other hand if n e 0 = n f 0 , if follows that

If N + 0 = ∞ then exactly one of the two points e or f disappears before being doubled, thus θ ∈ A 2 , the ITM is of infinite type and

If N + 0 < ∞ then neither preimage disappers before being doubled, thus θ ∈ A 3 .

Proof of Theorem 7. The main difference with Theorem 4 is that instead of two points {e, f } which produce right special words, there are now R + N points. Call these point {g 1 , g 2 , . . . , g R+N }. The other difference is that there are R + N intervals of continuity, thus p(1) = R + N and thus s(0

Otherwise the proof is identical, for each of the g i we construct its preimage tree. We consider the first time n gi when the preimage is doubled. The symmetry argument is identical to the square case, once a g i has two preimages at some time, it has two preimages for larger times. The k in the statement of the theorem then corresponds to the number of g i which have two-preimages at a certain time, while for the other R + Nk points the preimages disappear before being doubled. As before, C θ is responsible for the events that happened before the moment when the first preimage was doubled. We conclude that s(n + 1)s(n) = k for n sufficiently large and thus p

Proof of Theorem 9. Remember that in this case the map S is the first return map to the union of vertical sections x = k q described before the Theorem. The proof follows the same line as the the case a = 1 4 , in which we argued that by symmetry that once there are two preimages of e, each time one disappears a new one appears. This is no longer true in the general rational case. As in the case a = 1 4 we need to consider the preimage tree under the map S of the points e and f (which as in case a = 1 4 are the points of discontinuity on the circle x = a = p q ), the map being continuous on the other circles.

As mentioned in the proof of Theorem 4, the difference (s ns 0 ) e counts the number of leaves of level n, but since all points in the set {S -n (i, y)} have the same second coordinate, we have #{S -n (i, y)} ≤ q, for any n ≥ 0, for any point (i, y).


Mangiarotti
email: sylvain.mangiarotti@ird.fr

P & Mougin

Predictability of vegetation cycles over the semi-arid region of Gourma (Mali) from forecasts of AVHRR-NDVI signals

Keywords: vegetation cycle, semi-arid region, horizon of predictability, spatial scale, NDVI satellite data, nonlinear prediction

The NOAA-AVHRR Normalised Difference Vegetation Index (NDVI) dataset is used to investigate the predictability of the vegetation cycle in an area centred on the Gourma region in Sahelian Mali at scales varying from 8 km 2 to 1024 km 2 over a period spanning from 1982 to 2004. The predictability of the vegetation cycle is analysed with a model based on a reconstruction approach that fully relies on the dataset. Two parameters deduced from the growth of the forecast error are considered: the horizon of effective predictability, H E , which is the horizon at which a satisfying prediction can be effectively forecasted at a given level of error, and the level of noise.

Predictability is therefore analysed with regard to the horizon of prediction and the spatial scale; the influence of the model's dimensions is also discussed. The analysis clearly indicates that the signal predictability increases, and the level of noise decreases with an expanding area. However, even though the signal is more regular, its complexity increases within the narrowing entangled trajectory, setting the level of error of any prediction at a minimum of 15%, which matches the level of noise characteristic of the AVHRR-NDVI data series.

The forecasting error quickly increases with the horizon of prediction, setting the optimum horizon of predictability in the range of 2 to 4 decades, with high intraannual variability. At the short horizon of 1 decade, a resolution of 16 km 2 is reasonable to achieve an accuracy of approximately 20%. At the longer horizon of 3 decades, only low resolutions (256 km 2 or lower) give an accuracy equal to or better than 35%.

Introduction

The annual cycle of Sahelian vegetation is split in two seasons: the rainy season, which is controlled by the African monsoon regime, with the earliest rains beginning in June and lasting until the end of September, and the dry season, consisting of the remaining 8 months. The economy and, therefore, the society of this region rely strongly on cropping and pastoral activities, which in turn, depend greatly on the rainfall regime driving primary production. This dependence gives practical importance to the predictability of the vegetation cycle.

This analysis focuses on the Gourma region (in Mali), located south of the large loop in the course of the Niger River between Timbuctu and Labezanga at the border between Mali and Niger, downstream of Gao. Although the climate in the region is characterised by a strictly monomodal pattern of rainfall events, the rainfall distribution within the rainy season is irregular, and its spatial distribution is patchy [START_REF] Lebel | Rainfall Estimation in the Sahel: what is the Ground truth?[END_REF][START_REF] Frappart | Rainfall regime over the Sahelian climate gradient in the Gourma, Mali[END_REF]. As a consequence, primary production is strongly variable in time and space, at least between years and locally, progressively smoothing as the area considered expands.

The Normalised Difference Vegetation Index (NDVI) is an efficient parameter for monitoring the spatial patterns of vegetation and their cycle timing from space [START_REF] Tucker | Red and photographic infrared linear combinations for monitoring vegetation[END_REF]. The NDVI takes advantage of the differential reflectance of green vegetation between the infrared and near-infrared bands. The NDVI provides information about the green leaf area index (LAI) and the fraction of photosynthetically Active Radiation (FAPAR) absorbed by vegetation (e.g., [START_REF] Myneni | On the relationships between FAPAR and NDVI[END_REF], which are used as main inputs in primary production models (e.g., [START_REF] Running | A general model of forest ecosystem processes for regional applications. I. Hydrological balance, canopy gaz exchange and primary production processes[END_REF]. Thus, NDVI data are used to estimate primary production (or yields) in space and over time through empirical regressions (e.g., [START_REF] Maselli | Use of NOAA-AVHRR NDVI data for environmental monitoring and crop forecasting in the Sahel. Preliminary results[END_REF][START_REF] Balaghi | Empirical regression models using NDVI, rainfall and temperature data for the early prediction of wheat grain yields in Morocco[END_REF]. Alternatively, these data are introduced into predictive models of vegetation growth to better constrain these models using assimilation and/or inversion approaches [START_REF] Jarlan | Assimilation of SPOT/VEGETATION NDVI data into a Sahelian vegetation growth model[END_REF][START_REF] Mangiarotti | Evolutionary Bi-Objective Optimization of a semi-arid vegetation dynamics model with NDVI and σ0 Satellite data[END_REF].

The predictability of natural systems depends first on the dynamics of the system itself: an unstable system will naturally present a lower predictability than a stable one. Dynamical instability can be characterised by the first positive Lyapunov exponents [START_REF] Wolf | Determining Lyapunov exponents from a time series[END_REF] that provide information about the rate of separation of infinitesimally closed trajectories. Dynamical instability is a typical characteristic of chaotic systems, which although deterministic, are characterised by low predictability. Of course, the ability of a model to simulate or to forecast behaviour will rely on the stability of the system: the more unstable the dynamics, the more difficult the forecasting. Stability is therefore an important indicator of the dynamic properties of a system and should be given particular attention with respect to forecasting. The predictability of the Sahelian vegetation cycle was considered in a previous study [START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF] using the theory of nonlinear dynamics as a background. This former study was purely statistical and did not require any forecasting: predictability was assessed through statistical analysis of time series providing estimates of geometric invariants from the correlation integral [START_REF] Grassberger | Characterization of strange attractors[END_REF][START_REF] Diks | Nonlinear time series analysis[END_REF]. The Gaussian kernel algorithm developed by C. Diks and his colleagues [START_REF] Yu | An efficient implementation of the Gaussian kernel algorithm in estimating invariants and noise level from noisy time series data[END_REF] was used for this purpose. Three variables were estimated and analysed: the correlation dimension, D c , which is an indicator of dynamical complexity (this parameter will not be considered in the present paper); the correlation entropy, K, which is an indicator of dynamical stability (to make this parameter more directly understandable, this entropy was converted to a horizon of predictability, H P , defined as the time for error to double); and an estimate of Gaussian noise, σ, was also provided by the algorithm. This previous analysis provided us with estimates whose order of magnitude should be confirmed and refined.

Although the analysis presented in this paper is also based on the theory of dynamical systems, the approach used to investigate the predictability of the vegetation cycle is quite different. Indeed, the NDVI signal is taken here as a proxy of the system's dynamics and used as a forecasting model; for this reason, the corresponding models will be referred to as proxy models hereafter. Forecasts are, thus, obtained from these models, and predictability is assessed from the growth of the forecast error. Because proxy models are directly derived from observational time series, the results are only slightly model dependent.

Two main parameters will be estimated in this analysis: the horizon of effective predictability, H E , which is estimated from effective forecasts; and the level of noise. These parameters will be compared with the results obtained in the previous study for the horizon of predictability, H P , and the level of Gaussian noise, σ, respectively.

The main aims of this work are to introduce a prediction approach based on a proxy model built directly on an NDVI series and to demonstrate its forecasting ability and its usefulness for studying the predictability of vegetation cycles in semiarid regions. Because spatial non-stationarity is expected to arise from dynamics when addressing varying scales, the balance between prediction accuracy and spatial scales is considered carefully. As this work represents a first step, it will not consider the entire West Africa region but will instead focus on a smaller region centred on Gourma.

After a brief description of the data and data pre-processing are presented in section 2, the theoretical background of the forecasting approach is described in section 3, as is a description of the methodology used to estimate the horizon of effective predictability, H E , and the level of noise. The results are presented and discussed in section 4, and finally, conclusions are drawn in section 5.

Data

The NDVI product used in this study c was produced by the Global Inventory Modelling and Mapping Study (GIMMS) of the Global Land Cover Facility (GLCF, www.landcover.org) [START_REF] Tucker | Global Inventory Modeling and Mapping Studies[END_REF]2005;[START_REF] Pinzon | Satellite time series correction of orbital drift artifacts using empirical mode decomposition[END_REF]. These NDVI data have been calculated using radiometric data from the AVHRR (Advanced Very High Resolution Radiometer) sensor onboard NOAA (National Oceanic Atmospheric Administration) satellites 7, 9, 11, 14, 16 and 17. The data used cover the period spanning from 1982 to 2004 at a spatial resolution of 8×8 km 2 resulting from subsampling of initial 1.2 km² pixel data [START_REF] Justice | Spatial degradation of satellite data[END_REF][START_REF] James | The pathfinder AVHRR land data set: an improved coarse resolution data set for terrestrial monitoring[END_REF]. The data are 10-day composites selecting the maximum NDVI value for each pixel over ten days to reduce atmospheric interference and observe angle effects on the signal [START_REF] Holben | Characteristics of maximum-value composite images for temporal AVHRR data[END_REF]. Corrections applied to the dataset account for sensor degradation over time as well as the geometric effects of the view and atmospheric aerosols resulting from eruptions of El Chichon and Mt Pinatubo. To fill the gaps encountered in the time series (less than 4% in the region of study), data pre-processing is applied, including temporal interpolation when 1 or 2 successive data points are missing in the series and spatial interpolation in other cases (see Mangiarotti et al., 2010 for details). Note that such treatments are no longer required with the last version of the GIMMS product. However, to achieve a rigorous comparison of the results of this study with those of the previous study, the same time series had to be used.

The effective predictability of the signal is considered at scales ranging from 8×8 km 2 pixels (aggregation index, I A = 1) to 1024×1024 km 2 pixels (aggregation index, I A = 8). Series representing scales larger than the original resolution are constructed via simple spatial aggregation of the original time series. To reduce the number of basic assumptions and to maintain coherency with the former study, no direction is privileged, and no land surface classification has been considered in the aggregation process. The number of aggregated series and the corresponding resolution are given in Table 1.

This study is centred on the Agoufou site of the AMMA program (15.3°N, 1.5°W, see [START_REF] Mougin | The AMMA Gourma observatory site in Mali: relating climatic variations to changes in vegetation surface, hydrology, fluxes and natural resources[END_REF] located in the centre of the Gourma region (Figure 1). The surface characterised by the NDVI time series varies greatly when the signal is spatially aggregated. A qualitative description of the aggregated surfaces is given in Figure 2 with regard to the scale. At I A = 1, the vegetation is mostly characterised by open grassland and shrubby steppes. As the aggregation scale increases, the area considered expands in both longitude and latitude, including an increasing diversity of landscapes. As a consequence, the percentage of grassland regularly decreases as the aggregation scale increases, reaching ~15% at I A = 8, whereas the contribution of shrubby steppes, together with that of woody and shrubby savannahs from the Sudanian and Guinean zones, increases up to ~15%. The percentage of bare soils linked to rock outcrops and to the bare sands of Saharan landscapes becomes significant beginning at I A = 3. Other contributions become clearly significant (> 5%) starting at I A = 7, especially open and park croplands (annual crops), further associated with perennial crops at I A = 8.

Theoretical background and methodology

The predictability of the vegetation cycle is investigated using the theory of nonlinear dynamical systems [START_REF] Bergé | Order within chaos: toward a deterministic approach to turbulence[END_REF] for support. This theory associates a geometric trajectory with the successive states experienced by a system. Practically, each state is characterised by a vector of values corresponding to the variables of the system. This trajectory can be reconstructed easily when all of the variables of the system are assessed. Most often, only a few variables are measured or usable, impeding the completion of the task. However, using the theorem of Floris Takens [START_REF] Takens | Detecting strange attractor in turbulence[END_REF], under ergodic assumptions, information regarding the whole system's dynamics can be retrieved from a single time series. In this application, the trajectory is reconstructed from the NDVI series and used as a proxy for the original dynamics. A forecasting algorithm is built from this reconstructed trajectory, with the behaviour of one current state being estimated by considering the neighbouring trajectories experienced by the system (see [START_REF] Farmer | Predicting chaotic time series[END_REF] for a development based on similar background; see [START_REF] Miksŏvský | Testing for nonlinearity in European climatic time series by the method of surrogate data[END_REF] for an application to climatic time series).

Reconstruction approach

Among the methods for reconstructing a trajectory in multi-dimensional space from a single time series, the delayed variable method consists of representing the NDVI time series   t x as a function as itself with time delays:

      1 ... , ,     e d j j j j x x x y ( 1
)
where j is the j th decade since January 1982; τ is the delay; and d e is the embedding dimension. An illustration of a 3-dimensional reconstruction is presented in Figure 3. Reconstructed trajectories of NDVI time series corresponding to spatial scales I A = 1, I A = 5 and I A = 8 are presented in Figure 4 at dimensions of d r = 2 and d r = 3. For this purpose, it is necessary to determine appropriate values of the time delay, τ, and the embedding dimension, d e . The appropriate time delay is estimated from the first minimum of the average mutual information function [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF], which is a generalisation of the correlation dimension. The resulting delay corresponds to the average time for information to be lost. The embedding dimension is the number of dimensions required for appropriate reconstruction of the trajectory. The embedding dimension is classically estimated using the Global False Neighbours method [START_REF] Kennel | False neighbors and false strands: A reliable minimum embedding dimension algorithm[END_REF]. In practice, this method consists of checking that two trajectories cannot cross each other in space. The lowest dimension value for which no crossing is observed is the embedding dimension d e . When applied to the present case, this method provides a dimension range of 3 to 4. However, this approach has limitations (see [START_REF] Letellier | Comparison of tests for embeddings[END_REF], and as an important aim of this study was to assess the optimal predictability of the vegetation cycle, it was decided not to fix the embedding dimension but to retain and explore a large range of reconstruction dimensions from d r = 2 to d r = 8 to obtain the optimal forecasting dimension.

Simpler approaches

For the sake of comparison, two alternative approaches, a seasonal model and a plain climatology model, were developed. The seasonal model is a 2-dimensional model in which one dimension is the seasonal time. One state of the system is therefore defined as follows:

  j j j x doy y ,  (2)
where j is the j th decade since January 1982; doy j is the corresponding decade (each month is split in three decades) of the year; and x j is the corresponding NDVI value. An illustration of the corresponding trajectory is provided in the top panels of Figure 4 for I A = 1, I A = 5 and I A = 8. The trajectory crossings of this representation clearly indicate that this approach does not allow embedding of the NDVI signal.

The climatology model is based on the intra-annual average of the full series as follows: for each decade of the year, an averaged value of the NDVI is computed from the NDVI values for the same decade of all other years available. The state of this model is also 2-dimensional but is characterised by a single cycle. Its trajectory is defined as follows:

  climato , j j j x doy y  (3)
where j is the j th decade since January 1982; doy j is the corresponding decade of the year; and climato j x is the averaged value. Because of its static definition, this last model does not allow for any inter-annual variations in the prediction. The corresponding trajectory is given by the individual dark lines in the top panels of Figure 4.

Forecasting

The trajectory reconstructed in the phase space is a geometric representation of the vegetation dynamics constructed from the successive dynamical states experienced by the system. This trajectory is used as a surrogate of the dynamical equations underlying the system. The general principle of the forecasting algorithm is as follows: starting from a current state of the system corresponding to one point in the reconstructed space, the neighbouring states already experienced by the system are identified in the vicinity of this state; the evolution of these neighbouring states is monitored over time (by simply following the trajectories) and used to forecast the future states of the system. An illustration of the method is presented in Figure 5 and is based on a local representation of the geometric trajectory reconstructed in a 3dimensional space. Corresponding states of the NDVI values are also shown in Figure 6 in a temporal representation. More precisely, forecasts are achieved through the following steps. (1) The neighbouring states of the currently experienced state are searched within a circle whose radius, r n , has a moderate length. This circle is represented in black in Figure 5. (2) The geometric distances of these neighbouring states from the current state are calculated. (3) Only the N n closest states are retained (if the number of closest neighbours is less than 40, the search radius is increased, and steps (2) and (3) are repeated accordingly). In Figure 5, these closest states are represented as small (blue) squares inside the black circle. (4) A weight is attributed to each of these N n closest states as a function of its distance to the current state (an L 1 norm is used for this purpose). ( 5) The evolution of these small squares over time is investigated at a chosen horizon of prediction, h, which is the term for which a forecast is targeted. The resulting states are obtained after h decades along the trajectory (the direction of the trajectory is illustrated with an arrow). In Figure 5, the resulting states are plotted after 11 decades as small dark (blue) circles. A weighted average of the resulting states provides the forecast of the new state, which is represented in Figure 5 as a thick black circle. The 10 preceding black dots represent the forecast obtained in a similar manner for h = 1 to 10 decades. This forecasted trajectory can be compared a posteriori to the real (temporarily unknown) trajectory, which is represented in the figure as a (red) dotted line (corresponding to h = 1 to 10 decades) ended with a (red) cross (corresponding to h = 11 decades).

The forecasting approach used based on the seasonal model is similar.

(1) The distances of the state currently experienced by the system from the states experienced in other years (in the same decade of the year) are estimated. (2) A relative weight is deduced from these distances. (3) The evolution of these neighbouring states is considered at the horizon h, basically following their trajectory. (4) The forecast is obtained from a simple weighted average of the resulting states.

Finally, the predictions of the climatology model do not depend on the current state of the system but only on the decade of the year for which the forecast is run. The forecast is therefore directly provided by the static value of the climatology in a given decade, doy j .

Effective predictability

Predictability is estimated through a series of forecasting experiments upon which a statistical analysis is compiled. One experiment is constructed as follows: any time series is divided into 2 subsets (see Figure 6). The first subset (the part of the NDVI curve that has been removed) covers the time interval over which the NDVI will be forecasted. These data are removed from the time series and set aside to assess the forecasting skill. The forecasting skill is estimated over a horizon ranging from 1 decade to 1 year. To avoid boundary effects, the length of the forecasting window is fixed at two years. The second subset (light green line) is used as the informative part of the forecasting model. To obtain sufficiently long data archives of a constant length to achieve a robust estimate of the statistical error, data posterior to the prediction zone are also included in this subset.

Statistics for the predictive skill are based on the forecasting error,   h e j , which is the forecasting error at time j+h when forecasting from time j for the horizon h. The forecasting error is defined as follows:

    h j j j x h x h e    ˆ (4)
where h j

x  represents the data at time j+h, and  

h x j ˆ
is the forecast from time j at horizon h, corresponding to time j+h.

The % of error   h p is defined as follows:

    100   x e h h p   (5)
where x  is the standard deviation of the original time series x j , and  

h e  is defined as         N j j e h e N h 2 1  (6) Note that in practice, the systematic error       N j j h e N h e 1 is very low (   3 10   h e
).

Results and discussion

Inter-annual variability and scale dependence of the forecasting error

Both the forecasting error   h e j and its absolute value   h e j are presented in Figure 7 together with the corresponding time series   t x with an original resolution of I A = 1. In these patterns, j is the date from which the forecast is estimated. The diagonal orientation of the pattern is the most salient feature and represents the error associated with temporally approaching events. This diagonal pattern exhibits strong intra-annual heterogeneity with a low predictability (larger error) during rainy seasons and a high predictability (lower errors due to less vegetation) in the dry season (see Figure 7, bottom panel). Large errors may also occur during dry seasons (e.g., at the end of 2003). Large inter-annual variability is also observed with maximum magnitudes of errors generally corresponding to the peak of the rainy season (Figure 7 middle panel). Positive errors correspond to low rainfall years while negative errors correspond to high rainfall years. These results highlight the limits of the method with respect to estimating extreme events (maximum or minimum) when based on a representation of the dynamics restricted by the length of the time series (23 years of data available).

Another important feature of the pattern is the low error values at short horizons. However, the error level quickly increases with the magnitude of the horizon. Indeed, at a horizon of 1 decade, the error does not exceed ±0.05 in magnitude, while at 3 decades, it remains in the range of ±0.15, and at 4 decades, it surpasses ±0.25 in magnitude.

Although less salient, some vertical patterns also appear in the error diagram (e.g., at the end of 1985). Vertical patterns are associated with a poorly experienced current state from which the prediction is achieved. These effects can be observed at any horizon but are generally characterised by a limited range of action of ±5 decades. These patterns can be interpreted to result from temporary meteorological conditions (sudden growth of vegetation associated with large rainfall events or, in contrast, sudden drying up of vegetation due to long dry spells during the rainy season). These conditions act as a stochastic perturbation and diminish (generally temporarily) the predictability of the dynamics. Vertical patterns may also be due to artefacts associated with meteorological conditions (e.g., aerosols, dust) that affect the NDVI observed from satellites. Although they are less probable at these latitudes, perturbations caused by clouds cannot be excluded.

The forecasting error

 

h e h j is plotted in Figure 8 for I A = 1, I A = 5 and I A = 8. In this representation, the time corresponds to the target date of the prediction j+h (instead of the date of the initial condition j). As a consequence, the diagonal (approaching events) patterns encountered in Figure 7 become vertical in Figure 8, whereas the vertical (current events) patterns become diagonal. At I A = 1 and I A = 5, the error patterns are dominated by the peak annual error corresponding to the rainy season and vegetation growth. The annual variation pattern of NDVI predictability is low during vegetation growth, which depends strongly on rain (date of germination, peak yield), and high during vegetation senescence, which depends only marginally on rainfall. Between 1985 and[START_REF] Running | A general model of forest ecosystem processes for regional applications. I. Hydrological balance, canopy gaz exchange and primary production processes[END_REF], the error patterns associated with current events (i.e., the diagonal patterns in Figure 8) are more salient and exhibit a higher level of error overall, which may arise from the lower quality of the sensor used in earlier NOAA missions. Note that the range of the error progressively decreases with scale from ±0.25 at I A = 1 to ±0.1 at I A = 8.

Although the error level is lower at I A = 8 (see Figure 8, bottom panel), the trajectory appears to be more complex (its predictability varies greatly over time). Indeed, the annual error patterns are less salient. Intuitively, such complexity of the scene may have led to a greater complexity of the dynamics and therefore to lower predictability [START_REF] Justice | Spatial degradation of satellite data[END_REF][START_REF] Lambin | Change detection at multiple temporal scales: seasonal and annual variations in landscape variables[END_REF]. This is not what is observed. In practice, the aggregation of the scene tends to concentrate the signal in a more delimited zone, appearing to be visually simpler (compare the top panels in Figure 4) and becoming more predictable. However, the complexity of the dynamics has not actually disappeared but has just been condensed in a delimited zone. The complex pattern shown in Figure 8 reflects the intricate mix of trajectories in a narrower channel that ensues from the large number of aggregated series. At this scale, a wide range of vegetation landscapes are contributing to the signal, including both annual and perennial plants, such as those found in croplands, woody and shrubby savannahs, steppes and forests in the south and open grasslands and bare sands in the north (see Figure 2). Practically, one consequence of this complexity is the necessity of reconstructing the signal into a higher dimensional space to obtain more precise forecasts.

The error level also varies with the forecasting horizon. The % error   h p values obtained with the applied models (the climatology model, the seasonal model and the proxy models for d r = 2 through 8) are presented in Figure 9 for I A = 1, I A = 5 and I A = 8. As expected, the errors obtained with the climatology model are independent of the horizon of prediction h (straight line), highlighting the low forecasting skill of this model. The other approaches exhibit a clear dependence with regard to the target horizon, showing low errors at short horizons and a rapid increase in error as the horizon expands, followed by stabilisation affected by some minor oscillations.

The seasonal model exhibits poor scores, with a reasonable error being obtained only at horizons of less than 3 to 4 decades. The errors associated with the proxy models exhibit a more complex dependence. At a low horizon (h = 1 decade), the best results are obtained with d r = 2 through d r = 4. This behaviour reflects the greater robustness of low-dimensional models due to their better ability to filter contamination caused by noise (although the resulting dynamics will consequently be simplified). At longer horizons, higher dimensional models clearly become more efficient, showing compensation for the sensitivity to noise through the use of a larger window of memory, leading to a better representation of the dynamics. The convergence of the error with regard to the reconstruction dimension thus provides some information about the dynamical dimension of the underlying dynamics (see section 4.2).

When passing from I A = 1 to I A = 5 (Figure 9), the behaviour of the error   h p is mostly a matter of magnitude. At I A = 5, errors become lower than at I A = 1 for both the seasonal model and proxy models. At intermediate horizons (h = 4 to h = 8), the best improvements are seen for the models with higher dimensions (d r = 5 to d r = 8), revealing that deterministic and nontrivial dynamics are accessible at longer time scales. This intermediate scale appears to be optimal for achieving a reasonable forecast of a nontrivial behaviour. This scale is in agreement with the 'optimal' nontrivial prediction scale (I A = 5) found in the previous study [START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF]. The lower contamination by stochastic and local perturbations may also contribute to explaining the higher deterministic portion of the signal. At I A = 8, the level of relative error is much lower, regardless of the modelling approach applied. Predictability is therefore much higher, but the dynamics are more trivial and thus, less interesting (a simple linear forecasting approach would be sufficient but would provide a forecast with a poor resolution).

4.2

Model dimension and forecast accuracy Convergence to optimal models is appreciated through the convergence of the error level. This latter type of convergence is achieved at a lower dimension when expanding the area studied (see Figure 10 top panel). However, depending on the targeted horizon of prediction, the convergence of the system dimension clearly differs. The dimension for which the convergence of the error is obtained is shown in Figure 10 (top panel) for h = 1, h = 10 and h = 15 decades. The closer the horizon of prediction, the fewer reconstruction dimension are required. This result contributes to explaining the greater robustness of the low-dimensional approach at short time scales. Another salient feature of this result is the decrease in the reconstruction dimension required with an increasing area of application. Indeed, at h = 15 decades and more, the dimension requirement decreases from d r = 8 at I A = 1 to d r = 4 at I A = 8, which shows that the complexity of the trajectory at a local scale is higher than the complexity arising from the data aggregation process. This result also suggests that despite the heterogeneous composition of scenes, the regional scale is dominated by the global behaviour of mixed vegetation dynamics controlled by regional gradients in rainfall, which are locally patchy but result from a coherent monsoon regime. More sophisticated methods of data aggregation may be interesting to test in further analyses to account for directional effects resulting from the climatic gradients. Alternative strategies of aggregation by eco-zone, land cover type or rainfall isohyets may also lead to a better or at least complementary description of the dynamics of scenes. Indeed, the changes in the vegetation phenology analysed from AVHRR data over the 1982-2005 series by [START_REF] Heumann | AVHRR derived phenological change in the Sahel and Soudan, Africa, 1982-2005[END_REF] presented a zonal pattern with an increased amplitude of the NDVI signal in northern Sahel while the duration of the green signal expands in southern Sahel. A zonal distribution was also observed when attempting to disentangle the effects of climate and population in the AVHRR-NDVI trend over time in West Africa [START_REF] Herrmann | Recent trends in vegetation dynamics in the African Sahel and their relationship to climate[END_REF][START_REF] Seaquist | Disentangling the effects of climate and people on Sahel vegetation dynamics[END_REF][START_REF] Fensholt | Analysis of trends in the Sahelian 'rain-use efficiency' using GIMMS NDVI, RFE and GPCP rainfall data[END_REF].

Another way of analysing the forecasting error is to plot the lowest error min p obtained among the model forecasts as a function of the aggregation scale for horizons of prediction ranging from h = 1 to h = 15 decades (Figure 10

bottom panel).

A deeper horizon of prediction leads to a larger error. However, the error clearly reaches saturation at the horizons of 5 to 6 decades, which means that, statistically, averaged forecasts only improve marginally at farther horizons. At shorter horizons (h ≤ 4 decades), the error in the forecast strongly depends on the aggregation scale. Moreover, the horizon of predictability also varies with the expected accuracy. For example, at a maximum acceptable error of 20%, (statistically) satisfying forecasts at the horizon of 1 decade are obtained for most aggregation scales (slightly less than 1 decade at I A = 1, slightly more at I A = 8). Fixing the maximum accepted error at 30% instead extends the horizon at 2 decades for aggregation scales ranging from I A = 2 to I A = 7 (slightly less than 2 at I A = 1 and much more (4 decades) at I A = 8).

One of the objectives pursued in the first part of the previous study [START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF] was to determine the aggregation scale that would optimise the predictability of a non-trivial behaviour. The optimal aggregation scale found for the site of Agoufou was I A =3, corresponding to an area of 32 km 2 . The results of the present analysis clearly indicate that the aggregation scale is not the only factor that should be considered when searching for an optimal forecast of the vegetation cycle; the horizon of prediction is another factor to consider. Indeed, at low horizons, the results are quite similar whatever scale is considered, and full resolution (I A = 1) or the first level of aggregation (I A = 2) may be preferred. At larger horizons such as h = 5 decades, the errors observed at scales I A = 1 to I A = 6 are larger than 40%, and thus, higher aggregation scales should be chosen to keep the level of errors below the threshold.

Horizon of effective predictability

The horizon of predictability H P aims to assess the term for which a forecast can reasonably be targeted. [START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF] defined this parameter as the horizon at which the perturbation doubles. In this previous study, the horizon of predictability H P was derived from the correlation entropy K, which was computed from the correlation integral with a Gaussian kernel developed by C. Diks and his colleagues [START_REF] Diks | Nonlinear time series analysis[END_REF][START_REF] Yu | An efficient implementation of the Gaussian kernel algorithm in estimating invariants and noise level from noisy time series data[END_REF]. The horizon of predictability was statistically deduced from the structure of the trajectories in the phase space. In contrast in the present study, the horizon of predictability is not based on a direct analysis of the trajectories but on forecasts. Here, the predictability is assessed from the growth of the error of the effective forecasts provided by the proxy models. The horizon of predictability is therefore referred to as the horizon of effective predictability, denoted H E , and defined as the horizon for which the forecasting error doubles, assuming an exponential error growth. In practice, this parameter is computed from the ratio of the increasing error between decades 1 and 2 (as the level of noise corresponding to h = 0 is unknown) as follows:

        1 2 log 2 log p p H E  (7)
The resulting values of H E obtained at I A = 1 through 8 are shown in Figure 11 (top panel) together with the values obtained for H P via the correlation integral method [START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF]. The averaged values of H P and H E obtained for d r = 3 through 5 are also given in Table 1. The horizons of effective predictability H E are deeper than H P when H E ≥ 2 decades, while H P < 1 decade at most of the aggregation scales. At high aggregation scales, especially at the scale I A = 8, the horizons of predictability expand, regardless of the approach applied, and the differences between H E and H P tend to decrease (H E ranges from 2 to 3 and H P from 2 to 4 decades). This expansion indicates the regularity of the overall behaviour of the vegetation cycle at the scale of a large region.

The residual error   0 p is also estimated by assuming exponential growth for the error and by extrapolating this growth in a backward direction. The percentage of error at h = 0 decades is estimated as follows:

      2 1 0 2 p p p  (8)
When considering the optimal models at each scale (i.e., the proxy models exhibiting the lowest prediction error 1  h p min at h = 1), a clear scale dependence of the error appears, with errors ranging from 14.5% at I A = 1 to 10.4% at I A = 8 (see 0  h p min in Table 1). This decrease in error is due to compensation for the Gaussian noise when aggregating the signal. This decrease was only marginally detectable when the error was assessed via a nonlinear invariant approach [START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF]. This poor detectability could be arise from the uncertainty associated with the invariant approach under noisy conditions and from the difficulty of confirming the embedding dimension (additionally, estimates obtained with reconstruction dimensions of d r = 3 through d r = 5 were averaged). To allow comparison with the error estimates σ obtained via the nonlinear method, the errors observed for d r = 3 through 5 were averaged in both methods. Averaging increases the corresponding average errors 1) and tends to attenuate the decreasing pattern associated with an increasing aggregation scale (Figure 11 bottom panel). However, the average errors σ obtained with the invariant approach are much higher than 1  h p and 0  h p , suggesting that the invariant approach tends to overestimate the error.

1  h p and 0  h p (Table
As noted in the introduction, one of the critical issues in the Sahelian region is the forecasting of seasonal primary production. The absolute values of the forecasting error are plotted separately for aggregation scales I A = 1, I A = 5 and I A = 8 and for both approaching events   h e j and current events   h e h j (Figure 12). These graphs also show the local horizon (white curve) of effective predictability

  dOY H l E
, which was deduced from the seasonal error as follows:

          1 2 2 dOY dOY l E p p dOY H log log  (9)
The representation utilising the forecast date j+h as the abscise (Figure 12, top panels) clearly highlights two seasons: one with a short horizon of predictability ( l E H ranges from 2 to 6 decades), stretching from the earliest rain event of the rainy season until the end of the growing season of annual vegetation, and one with a shorter length and a farther horizon of predictability ( l E H >> 10 decades), corresponding to the dry season, during which annual plants are dead and form straw and litter. However, the horizons of predictability, l E H , being defined by the doubling of the error (equation 9) are highly sensitive to the initial error (here,

  1 dOY p
). The particularly long horizon of predictability values reached at the end of the dry season is attributed to the low dispersion of the signal due to the absence of vegetation. As expected, this behaviour is less marked at the scale of a large region (I A = 8) because of the greater contribution of evergreen perennial vegetation in the south. The farther horizon of predictability associated with the end of the dry season does not considerably contribute to the forecasting of primary production, which requires a farther horizon during the rainy season. However, the decrease in the horizon of predictability at the onset of the rainy season is not abrupt, as intermediate low values are reached, followed by a much slower decrease (with possible fluctuations) before the wet season lower threshold is reached, and these intermediate values may be helpful in predicting germination and early growth. If predictability is defined by the maximum threshold error, rather than by the growth of the error, quite different horizon magnitudes are obtained (see the dashed white lines plotted in Figure 12). This definition highlights the differences arising from variations in scales, with large increases of predictability being observed when the area studied is expanded. In particular, at an average NDVI error of greater than ±0.07, the horizon of predictability becomes very far for I A = 8 (even disappearing from the figure). In contrast, below this limit of precision, the horizon of predictability collapses. Indeed, at this level of precision, the trajectory becomes highly entangled because of the aggregation of a large scene of tremendous complexity. Scale non-stationarity therefore contributes to the origin of the increasing complexity of the patterns. Alternate aggregation strategies will likely lead to different but complementary results. An aggregation based on directional (latitudinal and longitudinal) averages or on rainfall isohyets would help to understand the gradient effects. Similarly, an aggregation based on eco-zones or land cover types would represent a better approach for addressing the various components of the scene.

The representation utilising the current date j as the abscise (Figure 12 bottom panels) provides complementary information regarding the horizon of predictability arising from the approaching event at any point in the course of the year. Two seasons are also clearly identifiable in this representation: at the onset of the dry season, the trajectories are converging, and the horizon of predictability expands abruptly, then regresses as the next rainy season approaches, characterised by a building up of divergence in the trajectories.

Conclusions

The predictability of the vegetation cycle in the semi-arid region of Gourma in Mali is analysed using a forecasting approach based on a reconstruction method. The analysis of predictability is estimated from an analysis of the forecasting errors. The forecasting error is first analysed over time with regard to the horizon of prediction. The forecasting error reveals a strong annual variability in the error, which is closely linked to the alternating rainy/dry seasons, with higher errors corresponding to the rainy season. The analysis indicates a rapid increase in error with regard to the horizon of prediction and shows large interannual variability, revealing the weak ability of the approach to forecast extreme situations. Forecasting error due to the error associated with the current state of the system is also assessed but is found to be comparatively low.

The degree of forecasting error clearly decreases as the aggregation scale increases, revealing the higher predictability of the behaviour of vegetation at the scale of large regions. However, the horizon of predictability also depends on the expected level of precision. On average, at a threshold level of error of 15%, no satisfactory forecasting based on a decadal NDVI dataset can be performed, even at the scale of very large regions [1024 km] 2 .

The reconstruction approach is investigated for dimensions ranging from d r = 2 to d r = 8 and is also compared with a seasonal model and to a plain climatology model. Models based on a reconstruction approach are systematically more efficient, regardless of the horizon of prediction targeted. However, the optimal dimension of such models varies with respect to the horizon of prediction. At a very short horizon (1 decade), low-dimensional models (especially d r = 3 and 4) are more robust, whereas at farther horizons (2 decades and more), 'optimal' predictions are only obtained with models of higher dimensions. From a functional point of view, this contrast of behaviour reflects the more homogeneous response of vegetation to meteorological variations at short time scales, as captured by simple models, and more spatially differentiated responses at longer time scales, requiring more complex models. Spatially differentiated responses may result from the different phenological behaviours of plant species at small aggregation scales or from the heterogeneous landscape (e.g., regarding the topography, soil texture, redistribution of water at the soil surface, soil fertility, and land use) and patchy climate forcing that occur at larger scales.

The convergence of the forecasting error at long horizons of prediction is used as a proxy to investigate the model dimensions required for optimal forecasts, suggesting a decrease from 8-dimensional dynamics at an aggregation scale of 8 km 2 to 4-dimensional dynamics at 1024 km 2 . These values are higher than the dimensions estimated with the global false neighbours approach (from 4 to 3). These results show that from a forecasting perspective, a fully empirical estimate of dimensions should be preferred.

Finally, the horizon of predictability and level of noise is estimated and compared with the results obtained in a former study based on a statistical analysis. The horizons of predictability estimated in the present study are longer at all aggregation scales (especially at small and intermediate scales), and the levels of noise estimated are significantly lower. The results obtained with the present approach are more consistent and robust than in the previous study.

The semi-arid Sahel region is not the only area where a better forecasting of the vegetation cycle can be expected from satellite data. Vegetation prediction can be considered to be an important task in most areas where cropping and pastoral activities take place. It would therefore be interesting to extend this analysis of effective predictability over the entire West Africa window. Moreover, as a forecasting tool, this approach may also prove to be an effective and robust method. Therefore, this approach should also be tested in other semi-arid regions of the world. Finally, among other developments related to potential applications to the real world, this approach could be coupled to biophysically based models, or these forecasts could be used as an input in simple statistical models used to estimate primary production or yields. is the correlation dimension; H P is the horizon of predictability estimated from the correlation entropy K; and the % noise represents the percentage of additive Gaussian noise, in which Diks c D , K and the % noise σ were simultaneously estimated using an algorithm from [START_REF] Diks | Nonlinear time series analysis[END_REF]. Figure 5: The geometric trajectory of the NDVI time series (see Figure 6) reconstructed in a 3-dimensional space using the time delayed method. The black arrow shows the direction of the flow (increasing time). The last states experienced by the system are represented by the (red) thick line; the large (red) square is the ultimate state from which the forecasting will be launched and is referred to as the current state. The dotted (red) line continuing the trajectory represents measures that have been removed from the trajectory and retained to analyse the forecasting skill. The states experienced by the system in the close vicinity of the ultimate state are identified by small squares (blue). The resulting states at the horizon h are represented by small circles (blue circles). The observed value at horizon h is indicated by the thick (red) cross. The forecasting trajectory is delineated by the large (black) dots. Current points corresponding to these various states are also represented in the temporal representation (Figure 6); the same symbols are used.

I A # of pixels Surface (km 2 ) x x  1  h p min (%) 0  h p min (%) 1  h p (%) 0  h p (%) 1  h r d 10  h r d 15  h
Figure 6: The green line represents the informative signal. The last states experienced by the signal are shown by the thick solid (red) line, until the ultimate state from which the forecasting is launched is reached. The small (blue) circles represent the N p states closest to the ultimate state experienced by the signal (see Figure 5). The small (blue) circles represent the time evolution after a given horizon of h decades. The black circle is the model forecast obtained from the weighted averaged values (computed in the embedded space, see Figure 5) of the evolution of the neighbouring states. This forecast can be compared with the real data represented with a (red) cross. Black dots represent intermediate forecasts, whereas the dashed (red) line represents the intermediate signal that has been temporarily removed for testing of the forecasting skill.

Figure 7: Relative (middle panel) and absolute (bottom panel) forecasting errors obtained from NDVI time series (top panel) at an aggregation scale of I A = 1 (full resolution). The abscise gives the date j from which the forecasting is launched, and the ordinate gives the horizon h of the prediction (in decade) at which the forecast is estimated. The error is given in NDVI units (no unit). Diagonal patterns represent the error associated with approaching events. Vertical patterns represent the error associated with currently experienced states. 7, the abscise gives the date j+h at which the forecasting is performed. The ordinate gives the horizon h (in decade) at which the forecast is estimated. The error is given in NDVI values (no unit). The orientations of the patterns are therefore contrary: vertical patterns represent the error associated with approaching events, whereas diagonal patterns represent the error associated with currently experienced states.

Figure 9: The % forecasting error is given as a function of the horizon of prediction h (in decade) for the original non-aggregated time series I A = 1 (left panel), for the I A = 5 aggregated time series (middle panel) and for the I A = 8 aggregated time series (right panel). A logarithmic representation of the horizon of prediction was selected to make the initial development of error more visible and to retain the behaviour observed at long horizons. The (red) solid thick line is obtained from the climatology model and the dashed (red) line from the seasonal model. The black lines correspond to the lowdimension proxy models: d r = 2 (dashed line), d r = 3 (solid line) and d r = 4 (dasheddotted line). The thin (green) solid lines correspond to proxy models for higher dimensions, d r = 5 through 8. 9) as a function of the aggregation scale I A , for h = 1 decade (green dots and line), for h = 10 decades (red dots and line) and for h = 15 decades (black circles and line). The forecasting error corresponding to the optimal dimension (i.e., the reconstruction dimension for which the lowest forecasting error is obtained) is given as a function of the aggregation index I A (bottom panel) for h = 1 to 6 decades (alternating solid (green) lines and dashed (red) lines from bottom to top) as well as for h = 15 decades (thick black line). (in decades) is superimposed for each of the six patterns (white lines). Seasonal evolution deduced from an NDVI error threshold of ±0.07 is also overlaid (white dashed lines); note that these last lines do not appear at aggregation scale I A = 8 because the averaged error is always lower than 0.07.

Figures

  is reached are given for horizons of h = 1, h = 10 and h = 15 decades, respectively. the increase in the error of prediction between 1 decade and 2 decades. The other estimates come from[START_REF] Mangiarotti | The Vegetation Cycle in West Africa from AVHRR-NDVI data: Horizons of predictability versus Spatial Scales[END_REF]:GFN e dis the embedding dimension estimated using the Global False Neighbours method; Diks c D

Figure 3 :

 3 Figure 3: An example of the reconstruction approach using the delayed vector technique. The reconstructed trajectory plotted in the right figure is built from the time series shown in the left figure. The initial state of the time series is based on three values (dark, red) separated by a time delay, τ, and represented by a cross, a circle and a square. These values can be represented in a 3-dimensional space, one on each axis, designated   t x ,     t x

Figure 4 :

 4 Figure 4: seasonal NDVI time series (top panels) and its geometric reconstruction (at d r = 2, middle panels, and d r = 3, bottom panels) at a non-aggregated scale of I A = 1 (left panels), an aggregated scale of I A = 5 (vertical middle panels) and an aggregated scale of I A = 8 (right panels). In the top panel, the thicker line represents the intraannually averaged signal (climatology model), whereas the thinner lines represent the individual years from 1982 to 2004 on which the seasonal model is based.

Figure 8 :

 8 Figure 8: Relative forecasting error obtained from NDVI time series at an initial resolution of I A = 1 (top panel), at an aggregated scale of I A = 5 (middle panel) and at an aggregated scale of I A = 8 (bottom panel). In contrast to what is shown in Figure7, the abscise gives the date j+h at which the forecasting is performed. The ordinate gives the horizon h (in decade) at which the forecast is estimated. The error is given in NDVI values (no unit). The orientations of the patterns are therefore contrary: vertical patterns represent the error associated with approaching events, whereas diagonal patterns represent the error associated with currently experienced states.

Figure 10 :

 10 Figure 10: Reconstruction dimension d r (top panel) corresponding to the convergence of the level of forecast error (see Figure9) as a function of the aggregation scale I A , for h = 1 decade (green dots and line), for h = 10 decades (red dots and line) and for h = 15 decades (black circles and line). The forecasting error corresponding to the optimal dimension (i.e., the reconstruction dimension for which the lowest forecasting error is obtained) is given as a function of the aggregation index I A (bottom panel) for

Figure 11 :

 11 Figure 11: Top graph: the horizons of predictability H P (green, light lines) and effective predictability H E (red, dark lines) are given as a function of the aggregation scale I A . Bottom graph: the % noise σ (green lines) and % error 0  h p extrapolated at h = 0 (red lines) are plotted as a function of the aggregation scale I A . In both graphs, results obtained at three reconstruction dimensions are presented: d r = 3 (solid line), d r = 4 (dashed-dotted line) and d r = 5 (dashed line).

Figure 12 :

 12 Figure 12: The patterns represent the seasonal variations (x axis in decade of year -dOY) in the error as a function of the horizon of prediction h (y axis in decades) for aggregation scales I A = 1 (left patterns), I A = 5 (middle patterns) and I A = 8 (right patterns). The top graphs show the forecast error at the decade j+h for which the forecast is performed whereas the bottom graphs show the forecast error at the decade j from which the forecast is launched. The seasonal evolution of the horizon of effective predictability   dOY H l E

  Figure 1:

Figure 8 (

 8 Figure 8 (black and white version for print edition):

Figure

  Figure 12:

1 2 Table 1 :

 21 scale dependence of the time series average x , standard deviation

  Maps of the study area. The larger scale map (I A = 8, 1028×1028 km 2 , full left map) includes a climatic gradient ranging from the Soudanian region in the south to the arid region in the north with a large panel of land surface (see working Table1for details). Decreasing areas corresponding to scales of I A = 7 to I A = 4 are also plotted. Note that the I A = 6 window is entirely included within Sahel and Gourma. Smaller areas corresponding to I A = 5 to I A = 2 are shown on the right map. The cross represents the Agoufou site of the AMMA program, on which all of the square areas are centred.

	Figure Captions
	1 2 3 4 5 6 7 8 Figure 1: Figure 2: Qualitative description of the aggregated surfaces with regard to the I A scale 1 8×8 0.21 0.080 22.9 14.5 24.4 17.0 3 4 4 16×16 0.20 0.070 20.0 12.9 21.9 15.5 2-3 4 16 32×32 0.19 0.067 19.6 12.6 21.6 15.3 2-3 4 64 64×64 0.19 0.068 18.7 12.0 20.9 14.6 2-3 4 256 128×128 0.19 0.064 18.8 11.5 21.3 15.0 2 4 1024 256×256 0.20 0.063 19.5 12.3 21.3 15.0 2 4 4096 512×512 0.21 0.062 17.8 11.5 20.5 15.0 2 4 16384 1024×1024 0.23 0.060 15.9 10.4 18.1 14.5 2 3 (see Figure 1). The distribution of the land cover is separated into natural vegetation,	d 8 6 6 6 5 5 5 4	r	2 (dec.) 1 E H 2.0 2.0 2.0 1.9 2.0 2.1 2.3 3.3	d 4 4 4 4 4 4 4 3	GFN e	Diks c D 1.59 1.55 1.67 1.68 1.73 1.68 1.59 1.33	H P (dec.) 0.5 0.8 0.7 0.6 0.6 0.7 1.0 2.4	Σ (%) 26.0 23.7 21.7 22.5 22.6 23.4 23.0 24.5
	crops and bare soils. The classes of natural vegetation include (a) thickets + riverine
	forests + open woodlands (>> 50% woody plants); (b) woody savannahs + shrubby
	savannahs + shrubs steppes (>> 15% woody plants, north contribution); (c) open
	grassland (90% annual herbaceous plants, < 1% woody plants); and (d) wetland
	aquatic herbaceous, including perennial grasses + open water (<10% woody plants).
	Crops are separated into two classes: (e) cropland parks (>> 3% woody plants) and (f)
	open cropland (annual crops, < 3% woody plants). Bare soils (0% vegetation) are
	separated into two classes: (g) bare sands and (h) rock outcrops. Note the jump in
	composition that occurs between I

A = 6 and I A = 7 when the area expands beyond Sahel.

Figure 7 (colour version for Web edition):

Figure 8 (colour version for Web edition):

Acknowledgments:

This work was performed within the framework of the AMMA and ECLIS projects. Based on a French initiative, AMMA has been constructed by an international group and is currently funded by large number of agencies, especially from France, the UK, the US and Africa. It has been the beneficiary of a major financial contribution from the European Communities Sixth Framework Research Program. Detailed information on scientific coordination and funding is available on the AMMA international web site (https://www.amma-eu.org/).