Frédéric Mazenc
email: frederic.mazenc@lss.supelec.fr

Olivier Bernard
email: olivier.bernard@inria.fr

Interval observers for linear time-invariant systems with disturbances ⋆

Keywords: Interval observer, linear system, exponential stability

It is shown that, for any time-invariant exponentially stable linear system with additive disturbances, time-varying exponentially stable interval observers can be constructed. The technique of construction relies on the Jordan canonical form that any real matrix admits and on time-varying changes of coordinates for elementary Jordan blocks which lead to cooperative linear systems. The approach is applied to detectable linear systems.

Introduction

Interval observers, appeared in the last decade as an alternative technique for robust state estimation [START_REF] Gouzé | [END_REF], [Rapport et al., 2003]. It is based on differential inequalities that allow to derive bounds for the solutions of a differential system at any time instant. This approach is significantly different from the technique based on classical observers (see for instance the contributions [Luenberger, 1964], [Besançon and Hammouri, 1996], [Andrieu et al., 2009]): when only the initial condition is assumed to be unknown, both classical observers and interval observers converge asymptotically to the solutions of the considered system but, in addition, interval observers supply certain information at any instant, provided that bounds on the initial conditions are known.

Interval observers have been introduced to cope with uncertainties that are known to characterize some classes of systems. That is one of the reasons why this technique is becoming more and more popular [START_REF] Moisan | [END_REF], [Raissi et al., 2005], [Mazenc and Bernard, 2010], [Mazenc et al.], and has been especially successful in the domain of biotechnological processes [Bernard and Gouzé, 2004], [START_REF] Moisan | [END_REF], [Moisan, 2007] where the uncertainties are large.

Interval observers can be straightforwardly constructed for any asymptotically stable cooperative systems (i.e. for systems whose Jacobian matrices have non-negative off-diagonal entries [START_REF] Gouzé | [END_REF]). Indeed, these systems have the distinctive property to keep the partial order between two trajectories [Smith, 1995]. Thus, if the (unknown) initial condition of the real system can be bounded between two known values, the trajectories of the same system starting from these bounds will enclose the real trajectory. Observe that cooperative continuous-time linear systems are always positive systems but this is not the case for discrete-time linear systems. These systems are positive systems if and only if the system matrices are nonnegative (see [Haddad et al., 2010] and [START_REF] Germani | [END_REF]). However, cooperativity is a rather specific feature, and most of the systems are not cooperative. This is the main limitation of the interval observer theory. In the present paper, we overcome it for a fundamental family of systems: the linear time-invariant exponentially stable systems with additive disturbances. We succeed in doing this by showing how any of these systems can be transformed through a change of coordinates into a cooperative system. This question has been recently solved for linear systems of dimension two in the work [Mazenc and Bernard, 2010] where two key results are established. First, it is proved that for a specific class of systems no time-invariant changes of coordinates can turn these systems into cooperative systems. Second, it is shown that any planar exponentially stable time-invariant linear system can be transformed into a cooperative exponentially stable linear system through a time-varying change of coordinates.

The purpose of this paper is to extend the main construction proposed for two dimensional systems in [Mazenc and Bernard, 2010]: we present a technique of construction of interval observers which applies to any linear exponentially stable time-invariant system of arbitrary finite dimension.

Our design relies on the fact that any real matrix can be transformed into a matrix of the Jordan canonical form (see [Perko, 1991, Section 1.8]). This property enables us to find changes of coordinates which transform any time-invariant exponentially stable linear system into a cooperative and exponentially stable system. Surprisingly, although the changes of coordinates we use are time-varying, the transformed systems are timeinvariant. However, the interval observers we obtain are time-varying because they involve a time-varying change of coordinates, and thus they give lower and upper bounds for the state of the system studied which depend on the time. It is worth noticing that we focus our attention only on exponentially stable linear systems because a system which is not asymptotically stable never admits open-loop interval observers. Hence our main result offers a complete picture of what can be done or not for linear systems: it implies that a timeinvariant linear system admits an interval observer if and only if it is exponentially stable.

We illustrate the broad applicability of our approach, by combining it (for observable systems) with a Luenberger observer [Luenberger, 1964], to end-up with a real-time estimate of the Luenberger observer accuracy.

The paper is organized as follows. Section 2 gives definitions and a basic result. In Section 3.1, we state that any time-invariant exponentially stable linear system can be transformed into a block diagonal system with cooperative and exponentially stable subsystems. The result of Section 3.1 is proved in Section 3.2. The main result of the paper is given and illustrated through an example in Section 4. The result is then applied to systems endowed with a Luenberger observer in Section 5.1. An illustration of this interval Luenberger observer in large dimension is presented in Section 5.2 and its efficiency is assessed with simulations. Concluding remarks are drawn in Section 6.

Classical definitions and results

Notations and definitions

• Any k × n matrix, whose entries are all 0 is simply denoted 0.

• All the inequalities must be understood component by component (partial order of R r ) i.e. ℧ a = (℧ a1 , ..., ℧ ar ) ⊤ ∈ R r and ℧ b = (℧ b1 , ..., ℧ br ) ⊤ ∈ R r are such that ℧ a ≤ ℧ b if and only if, for all i ∈ {1, ..., r}, ℧ ai ≤ ℧ bi .

• max(A, B) for two matrices A = (a ij ) and B = (b ij ) of same dimension is the matrix where each entry is

m ij = max(a ij , b ij ).
• A square matrix is said to be cooperative or Metzler if its off-diagonal entries are nonnegative.

• A system ẋ = f (x, w) is said to be cooperative if the Jacobian matrix of f with respect to x, denoted ∂f ∂x (x, w), is a cooperative matrix for all (x, w).

• The notations will be simplified whenever no confusion can arise from the context.

Interval observer: definition and basic result

For the sake of generality, we first introduce a general definition of interval observer for time-varying nonlinear systems.

Definition 1 Consider a system ẋ = F (t, x, w(t)) ,

(1) (2)

with x ∈ R n , w ∈ R ℓ ,
Moreover, the initial condition x(0) = x 0 is assumed to be bounded by two known bounds:

x - 0 ≤ x 0 ≤ x + 0 . (3)
Then, the dynamical system

Ż = ϕ (t, Z, w(t)) , (4)
associated with the initial condition Z 0 = G(t 0 , x + 0 , x - 0 ) ∈ R nz and bounds for the solution x:

x + = H + (t, Z) , x -= H -(t, Z) (5)
(with ϕ, H + , H -, G Lipschitz continuous of appropriate dimension), is called an interval observer (resp. an exponentially stable interval observer) for ( 1) if (i) system ( 4) is globally uniformly asymptotically stable (resp. globally uniformly exponentially stable) when w is identically equal to zero, (ii) for all Lipschitz continuous function w(t), all the solutions of ( 4) are defined over [0, +∞), (iii) for any vectors x 0 , x - 0 and x + 0 in R n satisfying (3), the solutions of ( 1), ( 4) with respectively x 0 , Z 0 = G(t 0 , x 0 + , x 0 -) as initial condition at t = t 0 , denoted respectively x(t), Z(t), are defined for all t ≥ t 0 and satisfy, for all t ≥ t 0 , the inequalities

x -(t) = H -(t, Z(t)) ≤ x(t) ≤ H + (t, Z(t)) = x + (t) . (6)
In the sequel we focus on linear systems. We recall a classical result for cooperative linear systems. For these systems the design of interval observers is straightforward.

Theorem 1 Consider the system

ẋ = Ax + w(t) , (7
)

with x ∈ R n , w ∈ R n , where A ∈ R n×n is a cooperative,
Hurwitz and constant matrix and let w be a Lipschitz continuous function bounded by two known Lipschitz continuous functions: for all t ≥ 0, w -(t) ≤ w(t) ≤ w + (t).

The system ( 7) admits the following exponentially stable interval observer

ż+ = Az + + w + (t) , ż-= Az -+ w -(t) , x + = H + (t, Z) = z + , x -= H -(t, Z) = z -, z + (t 0 ) = G + (t 0 , x + 0 , x - 0 ) = x + 0 , z -(t 0 ) = G -(t 0 , x + 0 , x - 0 ) = x - 0 , (8
)

with Z = (z + , z -) ⊤ and G = (G + , G -) ⊤ .
The proof of this result is in [START_REF] Gouzé | [END_REF]; it is directly related to the property of cooperative systems to keep the partial order between two trajectories.

Transformations of linear systems into cooperative systems

In this section, we state and prove that any exponentially stable time-invariant linear system can be transformed into a cooperative and exponentially stable timeinvariant system through a linear time-varying change of coordinates. We will use this result extensively in the next section to construct interval observers.

Change of coordinates: main result

The main result of the section is the following:

Theorem 2 Consider the system

ẋ = Ax , (9)
with x ∈ R n , where A ∈ R n×n is a constant Hurwitz matrix. Then there exists a time-varying change of coordinates ξ = P (t)x, where P is a matrix function whose entries are functions of class C ∞ bounded in norm, which transforms ( 9) into a cooperative and exponentially stable time-invariant linear system.

Proof of Theorem 2

The proof splits up into three steps. First we recall that any real matrix admits a real Jordan canonical form. In the second step, we transform, through a time-varying change of coordinates, systems with an elementary Jordan form into cooperative systems. In the last stage, we combine the results of the first two parts to end the proof.

Step 1: Jordan canonical forms.

From [Perko, 1991, Section 1.8], (see also [Hsu, 2005], p.73), we deduce that there exist two integers r ∈ {0, 1, ..., n}, s ∈ {0, 1, ..., n -1} and a linear timeinvariant change of coordinates

Y = Px , (10)
which transforms (9) into

Ẏ = JY , (11)
with

J =         J 1 0 . . . 0 0 J 2 . . . . . . . . . . . . . . . 0 0 . . . 0 J s         ∈ R n×n , (12)
where the matrices J i are partitioned into two groups: the first r matrices are associated with the r real eigenvalues of multiplicity n i of A and the others are associated with the imaginary eigenvalues of multiplicity m

i of A. Therefore n = r i=1 n i + s r+1 2m i and, for i = 1 to r, J i =         -µ i 1 . . . 0 0 -µ i . . . . . . . . . . . . . . . 1 0 . . . 0 -µ i         ∈ R ni×ni , (13)
where the µ i 's are positive real numbers and, for i = r+1 to s,

J i =            Λ i I 2 0 . . . 0 0 Λ i I 2 . . . . . . . . . . . . . . . . . . 0 . . . . . . . . . I 2 0 . . . . . . 0 Λ i            ∈ R 2mi×2mi , (14)
with

Λ i = -κ i ω i -ω i -κ i ∈ R 2×2 , (15)
and

I 2 = 1 0 0 1 , ( 16
)
where the κ i 's are positive real numbers and the ω i 's are non-zero real numbers. Notice that the fact that the constants µ i and κ i are positive is a consequence of the fact that A is Hurwitz. Notice also that if A has no real eigenvalues, then r = 0 and if all of the eigenvalues of A are real, then n = r i=1 n i . For j = 1 to r, the matrices J j are cooperative, but for j = r + 1 to s, they are not and one can prove (see [Mazenc and Bernard, 2010]) that, in the general case, the corresponding systems Ẏj = J j Y j cannot be transformed into cooperative systems through a linear time-invariant change of coordinates. On the basis of the result proposed in [Mazenc and Bernard, 2010] for planar systems, we show below that these systems can be transformed into cooperative systems through a linear time-varying change of coordinates.

Step 2: time-varying change of coordinates.

We consider the system ẋ = Γx , (17)

with x ∈ R 2p where Γ =            Λ I 2 0 . . . 0 0 Λ I 2 . . . . . . . . . . . . . . . . . . 0 . . . . . . . . . I 2 0 . . . . . . 0 Λ            ∈ R 2p×2p (18)
and

Λ = -κ ω -ω -κ ∈ R 2×2 , (19)
with κ > 0, ω = 0. We state and prove the following result:

Lemma 1 The time-varying change of coordinates

̺ = N (t)x (20)
with

N (t) =         ̟(t) 0 . . . 0 0 ̟(t) . . . . . . . . . . . . . . . 0 0 . . . 0 ̟(t)         ∈ R 2p×2p (21)
with

̟(t) = cos(ωt) -sin(ωt) sin(ωt) cos(ωt) ∈ R 2×2 , ( 22
)
transforms system ( 17) into

̺ = M̺ (23) with M =            -κI 2 I 2 0 . . . 0 0 -
           ∈ R 2p×2p . ( 24
)
Proof. Let us rewrite system (17) as follows

                               α1 = -κα 1 + ωβ 1 + α 2 , β1 = -ωα 1 -κβ 1 + β 2 , . . . αp-1 = -κα p-1 + ωβ p-1 + α p , βp-1 = -ωα p-1 -κβ p-1 + β p , αp = -κα p + ωβ p , βp = -ωα p -κβ p . ( 25
)
For each (α i , β i )-subsystem, we perform the change of coordinates

a i b i = ̟(t) α i β i . ( 26
) It gives ȧi = -ω[sin(ωt)α i + cos(ωt)β i ] + cos(ωt) αi -sin(ωt) βi , ḃi = ω[cos(ωt)α i -sin(ωt)β i ] + sin(ωt) αi + cos(ωt) βi . (27) It follows that, for i ∈ {1, ..., p -1}, ȧi = -ω[sin(ωt)α i + cos(ωt)β i ] + cos(ωt)[-κα i + ωβ i + α i+1 ] -sin(ωt)[-ωα i -κβ i + β i+1 ] , ḃi = ω[cos(ωt)α i -sin(ωt)β i ] + sin(ωt)[-κα i + ωβ i + α i+1 ] + cos(ωt)[-ωα i -κβ i + β i+1 ] . (28) Therefore, for i ∈ {1, ..., p -1}, ȧi = -κ cos(ωt)α i + cos(ωt)α i+1 +κ sin(ωt)β i -sin(ωt)β i+1 , ḃi = -κ sin(ωt)α i + sin(ωt)α i+1 -κ cos(ωt)β i + cos(ωt)β i+1 .
(29) Thus, we obtain, for i ∈ {1, ..., p -1}, ȧi = -κa i + a i+1 , ḃi = -κb i + b i+1 .

(30)

Similarly, we obtain ȧp = -κa p , ḃp = -κb p .

We deduce that the change of coordinates (20) transforms system (25) into system (23). This completes the proof of the Lemma 1.

Step 3: end of the proof of Theorem 2.

Thanks to the structure of matrix J defined in (12), we can decompose system (11) as follows:

Ẏa = GY a , (32)
Ẏb,j = J j Y b,j , j = r + 1, ..., s ,

with

G =         J 1 0 . . . 0 0 J 2 . . . . . . . . . . . . . . . 0 0 . . . 0 J r         ∈ R q×q (34)
with q = r i=1 n i and J j defined in ( 14), ( 15), for j = r +1 to s. Lemma 1 applies to each Y b,j -subsystem because ω j = 0. This leads us to perform the following timevarying change of coordinates

ξ = η(t)        Y a Y b,r+1 . . . Y b,s        (35)
with

η(t) =         I q 0 . . . 0 0 N r+1 (t) . . . . . . . . . . . . . . . 0 0 . . . 0 N s (t)         ∈ R n×n , (36)
where I q is the identity matrix of R q×q and, for j = r + 1 to s,

N j (t) =         ̟ j (t) 0 . . . 0 0 ̟ j (t) . . . . . . . . . . . . . . . 0 0 . . . 0 ̟ j (t)         ∈ R 2mj ×2mj , (37) with ̟ j (t) = cos(ω j t) -sin(ω j t) sin(ω j t) cos(ω j t) ∈ R 2×2 . ( 38
)
From Lemma 1, we infer that (35) transforms system (32), (33) into the system ξ = ℵξ ,

with

ℵ =         G 0 . . . 0 0 M r+1 . . . . . . . . . . . . . . . 0 0 . . . 0 M s         ∈ R n×n , (40)
with, for ℓ = r + 1 to s,

M ℓ =            -κ ℓ I 2 I 2 0 . . . 0 0 -κ ℓ I 2 I 2 . . . . . . . . . . . . . . . . . . 0 . . . . . . . . . I 2 0 . . . . . . 0 -κ ℓ I 2            ∈ R 2m ℓ ×2m ℓ .
(41) Since the matrices M ℓ and G are cooperative and Hurwitz (see ( 13) and ( 41)), we deduce that system (39) is cooperative and exponentially stable. We can thus apply Theorem 1 to system (39). We can conclude by observing that system (39) is obtained from (9) by performing the change of coordinates

ξ = P (t)x , (42)
with P (t) = η(t)P, where P is the matrix in (10) and η(t) is defined in (36).

Interval observers for linear exponentially stable systems

We can state and prove the main result of our work.

Statement of the main result

Theorem 3 Consider the system

ẋ = Ax + w(t) , ( 43
)
with x ∈ R n , w(t) ∈ R n ,
ż+ = ℵz + + P + (t)w + (t) -P -(t)w -(t) , ż-= ℵz -+ P + (t)w -(t) -P -(t)w + (t) , x + = H + (t, Z) = M + (t)z + -M -(t)z -, x -= H -(t, Z) = M + (t)z --M -(t)z + , z + (t 0 ) = G + (t 0 , x + 0 , x - 0 ) = P + (t 0 )x + 0 -P -(t 0 )x - 0 , z -(t 0 ) = G -(t 0 , x + 0 , x - 0 ) = P + (t 0 )x - 0 -P -(t 0 )x + 0 , Z = (z + , z -) , ( 44
)
where P + (t) = max(0, P (t)), P -(t) = P + (t) -P (t) and the matrix M (t) is the inverse of P (t), M + (t) = max(0, M (t)) and M -(t) = M + (t) -M (t), is an exponentially stable interval observer for system (43).

Remark. Theorem 2 applies to system (43) in the absence of the disturbance w and leads to a system to which Theorem 1 applies. However, since there is a disturbance in system (43), we cannot directly derive an interval observer for this system from these two theorems.

Proof of Theorem 3

Using the remark above, we observe that Theorem 2 provides a constant, Hurwitz and cooperative matrix ℵ and a function P of class C ∞ , bounded in norm, such that, for all t ∈ R, Ṗ (t) = ℵP (t) -P (t)A and P (t) is invertible. Next, we introduce the time-varying change of coordinates z = P (t)x, where all the entries of P are functions of class C ∞ bounded in norm, which transforms system (43) into ż = ℵz + P (t)w(t) .

(45)

We now prove that ( 44) is an exponentially stable interval observer for (43). The exponential stability of (44) when both w + and w -are identically equal to zero is a consequence of the fact that ℵ is Hurwitz. Moreover the solutions of ( 44) are defined over [0, +∞) for Lipschitz continuous functions w + , w -, because P + and P -are Lipschitz continuous.

Let t 0 be a nonnegative real number. Consider a solution (x(t), Z(t)) of ( 43) and ( 44) with initial conditions x(t 0 ),

Z(t 0 ) = (z + (t 0 ), z -(t 0 )) with z + (t 0 ) = P + (t 0 )x + 0 -P -(t 0 )x - 0 , z -(t 0 ) = P + (t 0 )x - 0 -P -(t 0 )x + 0 , where x + 0 , x - 0 are vectors such that x - 0 ≤ x(t 0 ) ≤ x + 0 .
We show that z(t 0 ) = P (t 0 )x(t 0 ) can be bounded, using both the equality z(t 0 ) = (P + (t 0 ) -P -(t 0 ))x(t 0 ), and the following inequalities:

P + (t 0 )x - 0 ≤ P + (t 0 )x(t 0 ) ≤ P + (t 0 )x + 0 , P -(t 0 )x - 0 ≤ P -(t 0 )x(t 0 ) ≤ P -(t 0 )x + 0 , (46)
which hold because the entries of P + (t) and P -(t) are nonnegative. We get thus:

P + (t 0 )x - 0 -P -(t 0 )x + 0 ≤ z(t 0 ) z(t 0 ) ≤ P + (t 0 )x + 0 -P -(t 0 )x - 0 , (47)
or, equivalently,

z -(t 0 ) ≤ z(t 0 ) ≤ z + (t 0 ) . (48)
Moreover, for all t ≥ 0,

P + (t)w -(t) ≤ P + (t)w(t) ≤ P + (t)w + (t) , P -(t)w -(t) ≤ P -(t)w(t) ≤ P -(t)w + (t) . (49)
It follows that, for all t ≥ 0,

P + (t)w -(t) -P -(t)w + (t) ≤ P (t)w(t) P (t)w(t) ≤ P + (t)w + (t) -P -(t)w -(t) . (50)
Then z(t) = P (t)x(t) is solution of ( 45) and since ℵ is cooperative and the inequalities (48) are satisfied, we deduce from Theorem 1 that for all t ≥ t 0 , the inequalities

z -(t) ≤ P (t)x(t) ≤ z + (t) (51)
hold. Since the entries of M + (t) and M -(t) are nonnegative, it follows that, for all t ≥ t 0 ,

M + (t)z -(t) ≤ M + (t)P (t)x(t) ≤ M + (t)z + (t) , (52)
M -(t)z -(t) ≤ M -(t)P (t)x(t) ≤ M -(t)z + (t) . (53) Since, for all t ≥ t 0 , M + (t) -M -(t) = P (t) -1 , it follows that, for all t ≥ t 0 , M + (t)z -(t) -M -(t)z + (t) ≤ x(t) x(t) ≤ M + (t)z + (t) -M -(t)z -(t) . (54)
This concludes the proof.

Three dimensional example

We illustrate our technique by constructing step by step an interval observer for the system

ẋ = Ax + w(t) (55)
with

A = 1 3     -1 0 -4 6 -9 0 10 0 -5     . ( 56
)
This system is not cooperative. One can check that A is Hurwitz and therefore Theorem 3 applies to (55).

Step 1: Transformation of A into a matrix of Jordan form.

For the matrix

P = 1 3     1 -3 1 1 0 1 2 0 -1     (57)
the equality

PAP -1 = J (58)
with

J =     -3 0 0 0 -1 2 0 -2 -1     (59)
is satisfied. The Jordan matrix J is not cooperative.

Step 2: Construction of a diagonalizing time-varying change of coordinates.

For the particular case we consider, the matrix η(t) in (36) is

η(t) =     1 0 0 0 C(t) -S(t) 0 S(t) C(t)     (60)
with S(t) = sin(2t), C(t) = cos(2t) and the matrix P (t) = η(t)P is

P (t) = 1 3     1 -3 1 C(t) -2S(t) 0 C(t) + S(t) 2C(t) + S(t) 0 S(t) -C(t)     . (61)
It is solution of

Ṗ (t) = ℵP (t) -P (t)A (62)
with

ℵ = 1 3     -3 0 0 0 -1 0 0 0 -1     . (63)
Observe for later use that

M (t) = P -1 (t) =     0 C(t) -S(t) C(t) + S(t) -1 C(t) S(t) 0 2C(t) + S(t) 2S(t) -C(t)     .
(64)

Step 3: Interval observer.

According to Theorem 3, system (55) admits the following interval observer

ż+ = ℵz + + P + (t)w + (t) -P -(t)w -(t) , ż-= ℵz -+ P + (t)w -(t) -P -(t)w + (t) , x + = M + (t)z + -M -(t)z -, x -= M + (t)z --M -(t)z + , z + (t 0 ) = P + (t 0 )x + 0 -P -(t 0 )x - 0 , z -(t 0 ) = P + (t 0 )x - 0 -P -(t 0 )x + 0 , Z = (z + , z -) , (65)
with ℵ defined in (63),

P + (t) = 1 3     1 0 1 L(C(t) -2S(t)) 0 L(C(t) + S(t)) L(2C(t) + S(t)) 0 L(S(t) -C(t))     , (66)
P -(t) = P + (t) -P (t) , (67)
M + (t) =     0 L(C(t) -S(t)) L(C(t) + S(t)) 0 L(C(t)) L(S(t)) 0 L2C(t) + S(t)) L(2S(t) -C(t))     , (68)
with L defined by L(s) = max{s, 0} and

M -(t) = M + (t) -M (t) . ( 69
)
5 Application to systems endowed with Luenberger observers

Interval observers for detectable systems

In this section, we show how interval observers and classical observers can be used simultaneously to obtain complementary information on the solutions of a linear detectable system. We consider a system

Ẋ = αX + βu + φ 1 (t) , (70)
with X ∈ R n , where u ∈ R p is an input, φ 1 is an unknown locally Lipschitz function and α and β are constant matrices. We assume that system (70) is endowed with an output

y(t) = γX(t) + φ 2 (t) , (71)
where γ ∈ R q×n is a constant matrix such that the pair (α, γ) is detectable and φ 2 is an unknown locally Lipschitz function. It is worth observing that we do not assume that α is Hurwitz.

The detectability property of the pair (α, γ) ensures that there exists a constant matrix L ∈ R n×q such that the matrix α -Lγ is Hurwitz [Brogan, 1990]. Then we consider the observer

Ẋ = α X + βu + L(y -γ X) . ( 72
)
The error variable X = X -X gives the error equation

Ẋ = (α -Lγ) X + φ 1 (t) -Lφ 2 (t) . ( 73
)
Since α -Lγ is Hurwitz, Theorem 3 applies to (73) (with the notations of Theorem 3, we have A = α -Lγ, w = φ 1 -Lφ 2 ). Therefore an interval observer can be constructed for this system. Thanks to it, an estimate of the distance between a solution of ( 70) and a solution of (72) can be obtained. Indeed, if (73) admits the interval observer Ż = ϕ (t, Z, w(t)) , (74)

associated with Z 0 = G(t 0 , X+ 0 , X- 0 ) ∈ R nz and X+ = H + (t, Z), X-= H -(t, Z) then, if for t 0 ≥ 0, the solutions X(t), X(t) satisfy X- 0 ≤ X(t 0 ) -X(t 0 ) ≤ X+ 0 (75)
and, for all t ≥ t 0 , w -(t) ≤ w(t) ≤ w + (t), then, the solution of (74) Z(t) such that Z(t 0 ) = Z 0 is such that, for all t ≥ t 0 ,

H -(t, Z(t)) ≤ X(t) -X(t) ≤ H + (t, Z(t)) . (76)

Interval observer design and simulation

In this simulation example we illustrate the result of Section 5.1. In order to demonstrate the applicability of our approach to high dimensional systems, we have considered a simulation example in dimension six. Our purpose is to illustrate how the proposed design can be applied to a standard Luenberger observer to estimate, in real time, its convergence despite uncertainties (here, noise in the measurements). The example we consider is the following linear system: ẋ = αx + βu(t) with u(t) = sin(t), β = (-18, -13, -5, -4, -10, 22) ⊤ and

α =             -1 1 0 0 -1 0 -1 -2 0 -1 0 1 -2 0 -3 -2 0 0 -1 0 -2 -3 0 1 -1 0 2 0 -4 0 -1 -1 0 1 0 -1             . (77)
Matrix α has stable poles. Four of them have an imaginary part different from zero. We denote them ω 1 ± jκ 1 and ω 2 ± jκ 2 and µ 5 and µ 6 the two real eigenvalues.

We assume that the output y = x 1 + φ 2 (t) is available, at discrete times (every 0.05 hours) i.e. x 1 is measured, with a bounded noise : -2 ≤ φ 2 (t) ≤ 1, for all t ≥ 0 (see Figure 1). The considered noise for simulation is based on a uniformly distributed noise in the interval [-2, 1] at the sampling instants. Function φ 2 (t) thus results from linear interpolation of φ 2 (t i ) between two sampling time instants.

Through routine calculations, one can prove that the pair (α, γ) with γ = (1 0 . . . 0) is observable. With a pole placement algorithm, we select a matrix K so that the poles of α -Kγ give a desired observer closed-loop dynamics. More specially, the slowest pole of α being µ 6 = -0.26, we select K in order to move this pole to -1, keeping the five other poles unchanged [Brogan, 1990]. Thus, the observer dynamics is faster than the original system and the gains are, roughly speaking, not too large in order to limit observer over reactivity in the presence of large noises.

Then, the interval observer presented in Theorem 3 can be derived. Due to the dimension of the system, this in- (20,10,6,20,30,40) ⊤ as initial condition is simulated. We also draw the solutions of the interval observer associated with unknown initial conditions but with components belonging to the interval [0,40].

Figure 2 shows the estimates for this interval observer. Due to the uncertainty in the knowledge of the initial condition and to the initial change of variables, the interval is initially quite large. However its length becomes rapidly reasonably small despite the uncertainties on the measurements y. At the end, this approach provides an evaluation of the Luenberger observer estimation accuracy, i.e. an error interval associated to the Luenberger observer predictions. It is worth noting that we have considered a case of biased noise (the average value of the noise is -0.5) to highlight the generality of the situations where this approach applies.

Conclusions

We have presented a construction of exponentially stable interval observers for the fundamental family of the time-invariant exponentially stable linear systems with additive disturbances. This approach has several possible applications. One of them is the possibility to on-line assess the accuracy of classical observers such as Luenberger observers. These techniques can be improved by considering simultaneously several admissible changes of variables and by taking the intersection of all the predicted intervals [START_REF] Moisan | [END_REF], [Bernard and Gouzé, 2004]. This may be important when one wants to improve the performances through bundles of interval observers in the spirit of what is done in [START_REF] Moisan | [END_REF]. This will be the subject of further investigations. Besides, extensions to time-varying, nonlinear and time-delay systems are expected.

Fig. 1 .

 1 Fig. 1. Variable x1 is supposed to be measured with a biased noise of average -0.5.

Fig. 2 .

 2 Fig. 2. Estimation of the bounds associated with the Luenberger observer (dotted lines) for the variables x1 and x6, and comparison with the real state (continuous line). Bounds given by the interval observer are represented by dashed lines terval observer is of dimension twelve and given by complicated expressions. A trajectory with t 0 = 0, x 0 =(20, 10, 6, 20, 30, 40) ⊤ as initial condition is simulated. We also draw the solutions of the interval observer associated with unknown initial conditions but with components belonging to the interval[0, 40].


Morad Lakhssassi

Bernard Guy
email: guy@emse.fr

Eric Touboul

Jean-Yves Cottin

Keywords: Cristallisation fractionnée, solidification, échange chimique, sédimentation, répartition bimodale, discontinuité de composition, chambre magmatique Fractional crystallization, solidification, chemical exchange, sedimentation, bimodal distribution, discontinuity of composition, magmatic chamber

Bimodal distribution of the solid

Introduction

Associated plutonic intrusions or volcanic emissions within the same province frequently show composition discontinuities (all the compositions are not equally seen), on local or larger scale (see for example [START_REF] Halmlyn | Origin of chromite compositional variation in the Panton sill, Western Australia[END_REF][START_REF] Leblanc | Chromite cristallization in a multicellular magma flow : evidence from a chromitite dike in the Oman ophiolite[END_REF]). One also very often observes bimodal distributions; for instance for the volcanic rocks, on one side the basaltic terms, and the other side the felsic terms (trachytes, phonolites) whereas the intermediate compositions (andesites) are very little represented (for the equivalent plutonic rocks: gabbroic rocks on one side, granitic rocks on the other, with the relative scarcity of intermediate rocks such as quartz bearing diorites or syenites) [START_REF] Baker | Intermediate oceanic rocks and the "Daly Gap[END_REF][START_REF] Bonnefoi | The « Daly gap » as a magmatic catastrophe[END_REF][START_REF] Chayes | Relative abundance of intermediate members of the oceanic basalt-trachyte association[END_REF][START_REF] Clague | The oceanic basalt trachyte association: an explanation of the Daly gap[END_REF][START_REF] Daly | The geology of Ascension Island[END_REF][START_REF] Gagnevin | Open-system processes in the genesis of silica-oversaturated alkaline rocks of the Rallier-du-Baty peninsula, Kerguelen archipelago (Indian Ocean)[END_REF][START_REF] Marsh | On bimodal differentiation by solidification front instability in basaltic magmas, part I: basic mechanics[END_REF]). One speaks of Daly gap for the absence of intermediate terms in the volcanic series of oceanic islands, but the situation is similar in the continental series. In the present paper, we wish to explain these features, and particularly the bimodality of the rock compositions, by way of simple models.

Many models have already been proposed for the first understanding of the variety of the plutonic and/or volcanic rock compositions of the same province. The fractional crystallization model is the best example. It suffers important limitations: -in the models that are derived from it, nonlinear (composition dependent) chemical partition coefficients between the solid and the liquid are not taken into account (the models are adapted to trace rather than major elements); -the relative displacement between the solid and the liquid, and the possible chemical exchanges along this displacement are not taken into account: there is no space variable. This model can be developed within larger models. For convenience, we will distinguish three classes of processes: -processes of chemical exchanges between liquids and solids; -processes of thermal transfer and solidification; -processes of relative movement between the solid and the liquid (in relation to tectonics, or to the density differences between the solid and the liquid). Each one of these classes deserves more thorough developments (see [START_REF] Baker | Liquid fractionation. Part III: Geochemistry of zoned magmas and the compositional effects of liquid fractionation[END_REF][START_REF] Bons | Melt extraction and accumulation from partially molten rocks[END_REF][START_REF] Godard | Effects of mineralogical reactions on trace element redistributions in mantle rocks during percolation processes: a chromatographic approach[END_REF][START_REF] Jaupart | Dynamics of differentiation in magma reservoirs[END_REF][START_REF] Jellinek | Magma dynamics, crystallization and chemical differentiation of the 1959 Kilaulea Iki lava lake, Hawaii, revisited[END_REF][START_REF] Kuritani | Magmatic differentiation examined with a numerical model considering multicomponent thermodynamics and momentum, energy and species transport[END_REF] for example).

However, the models remain generally specialized in each of the three classes of processes.

Coupled models exist, but they only relate to partial aspects of the phenomena, their variety being so large. And, in the end, the models do not allow a satisfactory answer to the above questions (appearance of bimodality and of composition discontinuities). Many explanations were proposed (see previous papers for example) and we will not discuss them in detail. They generally take us away from the elementary phenomena of fusion/solidification, relative displacement and chemical exchange between the solid and the liquid, and utilize some various additional phenomena (multiple injections, contamination in an open system etc.). In this work, we would like to go back to the elementary phenomena, deciding to consider them together. This can only be done at the price of a great simplification of each of them, but the result is a fertile overall picture. In particular, these three phenomena will reduce to two independent dimensionless parameters. These will enable us to discuss the important types of behaviour of the systems. We will take as an example that of the evolution of a closed magmatic chamber. This will be a vertical column containing the magma (with or without some solids) at time zero, along one space dimension, the vertical axis (Fig. 1). The other case of an open system is simpler, in that it corresponds to the subsystem of a larger system, and will also be discussed briefly. The three classes of phenomena are simplified as follows:

solidification is imposed by the thermal situation of the chamber, in particular the temperature gradient at the top. Without modelling all the heat transfers, we will consider that we are able to define a solidification rate q: this is the solid/solid + liquid proportion production per unit of time. This rate can be variable from one point to another of the chamber (one can restrict the solidification zone to the top of the chamber, or give any other conceivable evolution in space and time).

relative movement between the solid and the liquid is described by its velocity: once the solid is formed, we will consider that it can go down and settle at the bottom because of the density differences. The velocity is a function of temperature (viscosity) and of the solid abundance. Other modes of relative displacement can be generated by the action of internal or external tectonics (filter press phenomenon for example), and/or when the solid crystallizes on the walls of the chamber. On the whole, we will consider that we are able to write a law for the relative movement velocity between the solid and the melt.

-chemical equilibrium between the solid and the melt may finally more or less achieve; we will consider that the kinetics law is governed by departure from equilibrium. The equilibrium partition law between the solid and liquid will be supposedly known in an analytical form.

Importantly, this law can be nonlinear. By simplification, we will consider first the case of one independent chemical component.

In this approach, we took as a starting point the chromatographic model, as it is used for the understanding of metasomatic rocks; it calls upon similar simplifications for what relates to relative movement and chemical exchange phenomena. A major characteristic of this model is to account for the appearance of composition discontinuities, i.e. of plurimodal distributions of the crystallized rock compositions. We already proposed in a preliminary way its application to magmatic rocks, in particular in open systems, as also did other authors (for example [START_REF] Baker | Liquid fractionation. Part III: Geochemistry of zoned magmas and the compositional effects of liquid fractionation[END_REF][START_REF] Boudreau | Iridium, a program to model reaction of silicate liquid infiltrating a porous solid assemblage[END_REF][START_REF] Godard | Effects of mineralogical reactions on trace element redistributions in mantle rocks during percolation processes: a chromatographic approach[END_REF][START_REF] Guy | Mathematical revision of Korzhinskii's theory of infiltration metasomatic zoning[END_REF][START_REF] Guy | The behavior of solid solutions in geological transport processes : the quantization of rock compositions by fluid-rock interaction, in: Complex inorganic solids, structural, stability and magnetic properties of alloys[END_REF][START_REF] Hauri | Melt migration and mantle chromatography, 1: simplified theory and conditions for chemical and isotopic decoupling[END_REF][START_REF] Navon | Geochemical consequences of melt percolation : the upper mantle as a chromatographic column[END_REF][START_REF] Sedqui | Echange chromatographique de deux constituants indépendants entre un fluide aqueux et une solution solide à trois pôles; application à la substitution Ca-Fe-Mn des grenats de skarn[END_REF]). We want to give a more general approach (use o dimensionless parameters etc.) and extend it to the study of a closed magmatic chamber. We will show that, under relatively simple physical assumptions, the crystallization products can present a non-homogeneous distribution, in particular bimodal, starting from a homogeneous liquid. We present an abstract model and give the indications for applying it to some particular systems. The chamber, considered as a vertical column, is modelled along one space dimension. The axis is directed to the downward direction. The abscissa is noted x, the origin is taken at the higher point of the chamber. The total height is noted h; p S (x, t) or p(x, t) represents the solid volumic proportion (also called compactness) at position x and time t. The liquid or melt volumic proportion is thus p m =1-p. Solids and liquids are incompressible. The hydrodynamic aspect of sedimentation is taken into account through the Richardson-Zaki law [START_REF] Richardson | Sedimentation and fluidization, part 1[END_REF]:

0 (1 ) U U p   
, where U 0 is the Stokes velocity, and =5. For p values higher than 0.7, the sedimentation velocity is almost zero. This accounts for what is observed in experimental work and nature. The U m liquid velocity is obtained by writing the total volume conservation.

At any point, solid phase and liquid phase fluxes are equal and have opposite signs:

(1 ) m p U pU   
The solidification rate is represented by a function q S (x) or q(x) for the creation of solid volume per unit of total volume and unit of time. The x-dependency accounts for the fact that the solidification speed depends on temperature, and so on depth. In order to separate the intensity of the solidification process and its space variation, we will write q in the form of the product of a constant and a space function. So q(x) =q 0 .s(x); s(x) will be a simple function of

x. In certain cases, s might also depend on p and t. The unknown function is p(x,t) and the conservation equation of the solid phase is written as:

5 0 ( (1 ) ) ( p ) pU p q t x        x
The equation can be transformed by considering dimensionless time and space coordinates t * and x * .

* 0 * * , x [0,1] U t t h x x h          
It is then written in the form:

(1)

5 * * * 0 * * 0 ( (1 ) ) ( ) ( ) q h p * p p s x A s t x U         x imension With s*(x*) = s(hx*). A = q 0 h
/U 0 is the analogue of a Damkhöler number, measuring the ratio of solidification rate to sedimentation rate. Equation ( 1) is a non-linear hyperbolic equation, likely to give place to "shocks", i.e. discontinuities of the solution, even for very regular initial conditions (constant p).

Initial and boundary conditions: the initial condition will be p(x,0) = 0: at t= 0, the chamber is filled with magma, which then solidifies gradually, giving place to sedimentation. The conditions at top and bottom of the chamber are zero matter flux. It is considered that the bed of the settled solid has a null porosity, p=1 (and thus U=0). Other choices would be possible and lead to change space and time scales.

Chemical exchanges

One independent chemical component is taken into account, its concentration is noted C in the solid phase and C m in the melt (kg/m 3 of solid or liquid). The rate of variation of the chemical component is supposed to be proportional to the surface of the solid-liquid contact (reaction surface), and to a function G(C m , C) describing the departure from equilibrium, the d of which is the same as for C. For a solid volume V with reaction surface S r :

. . ( , )

r m VC K S G C C t   
K is the kinetics constant (ms-1). In the magmatic chamber, the solid volume is controlled by p variable. In order to model the contact surface, we consider that it depends only on compactness p. To estimate it, we suppose that the solid is composed of similar spheres with radius r, the surface of the grains (geometrical surface) is then obtained by a simple calculation according to compactness:

3 g p S a r p  
With a = 3/r. We then consider that the effective reaction surface S r is lower than the surface of the grains (when p approaches 1, S r goes to 0), which leads us to adopt a qualitative law:

(1 )

r g S p S    or (1 ) r S ap p   
Where one will take β = 0.3. The function G is written:

( , ) ( ( ) m m G C C C g C )   g(C)
, the isotherm, is able to have various, linear or not linear shape.

Solidification by itself also generates a chemical transfer from the liquid phase to the solid phase, without any chemical reaction, the rate of corresponding change is equal to qC m . The conservation equations of the chemical component are finally written as:

( ) ( ( )) (1 ) ((1 ) ) ( ( )) r m m m m m r m m pC pUC KS C g C qC t x p C p U C KS C g C qC t x                         
Time, space, the sedimentation velocity are transformed to dimensionless parameters in the same way as previously. The dimensionless concentrations are defined by ratios to the

densities * * , m m s m C C C C    
The conservation equations are then written as:

(

( ) ( 1 ) ( ( )) ( ) (1 )
(( 1)

) (1 ) ( ( )) ( ) m m m m m s m m pC ahK pU C p p C g C As x C t x U p C ahK * * * m p U C p p C g C As x C t x U                                                              With: * * * 1 ( ) ( ) s m g C g C   
This makes a second dimensionless coefficient appear:

0 ahK B U 
, ratio of the reaction velocity to the settling velocity (Damköhler number). The system is finally governed by two dimensionless parameters, A and B, which reflect the ratios of preponderance between the three phenomena that come into play: sedimentation, solidification and chemical reaction.

There remain however "form factors" or adjustable functions which are the space function s * (x * ) and the isotherm g * (C * ). One will note that, in order the preceding equations are consistent, the equality of  m and  s is implicit: in a first step we have assumed that there is no volume effect associated to solidification. Dimensionless velocities U* and U m * are obtained by dividing corresponding velocities by U 0 .

Ranges of values for the useful physical parameters, form of the laws used

We need to define ranges and typical values for the principal physical parameters in the equations. Those are synthesized in Table 1. We used various works (see already quoted works and [START_REF] Jellinek | Magma dynamics, crystallization and chemical differentiation of the 1959 Kilaulea Iki lava lake, Hawaii, revisited[END_REF][START_REF] Kirkpatrick | Kinetics of crystal growth from silicate melts: anothite and diopside[END_REF][START_REF] Schwindinger | Particle dynamics and aggregation of crystals in a magma chamber with application to Kilaulea Iki olivines[END_REF][START_REF] White | Geochemistry[END_REF]) together with some semi-quantitative reasoning. The variety of th natural situations is great (the preferred parameters chosen match the expertise of the author and other values of the parameters could be considered. For instance, the densities chosen refer mostly to mafic systems while felsic melts can be lighter than 2.6; the bubble effect could also be taken into account. The law chosen for the spatial distribution of solidification i given in Fig. 2a and the isotherm chosen for the numerical experiments is given in Fig. 2b.

Numerical simulation and results

The sedimentation equation (1) and the system of chemical exchanges (2) are both solved numerically by a "finite volume" type method. A dissipative numerical scheme of Lax-Friedrich type was adopted in order to avoid as much as possible the numerical oscillations in the vicinity of the shocks. Non-linearities are treated implicitly for the sedimentation equation, explicitly for the reaction equations. The complete solidification of any zone located above the sedimentation front was prohibited, so that the front is never blocked. The two

Visualization of the bimodality zones

In order to visualize the appearance of bimodality according to parameters A and B, we defined a "bimodality criterion". A classification is done between two classes by a traditional method minimizing the intraclass variance. One can then propose as a criterion a real number involving the intraclass variances V 1 and V 2 and the distance between the centers of each class, C 1 and C 2 . One does not take into account the ratio of the populations. Thus for example: criterion

1 2 1 2 *(exp( ) exp( )) C C V V     
. This gives the map in Fig. 5 (Log to the base 10), where the bimodality is very marked when the criterion goes close to or exceeds 1, in agreement with the first look at the preceding graphs. We see that the bimodality is obtained for A = solidification/sedimentation < 1 (the log is negative) and B = reaction/ sedimentation >1, or in other words, solidification < sedimentation < reaction.

Discussion

Our main result is as follows: starting from very uniform initial conditions, the solidified rocks within the magmatic chambers may present a bimodal distribution of their chemical composition. This is understood within a modelling which takes into account the simultaneous action of the sedimentation, solidification and chemical reaction phenomena.

The hyperbolic character of the equations plays an important part: the dynamic process of reaction is modelled in the same time as the relative movement of the solid with respect to the liquid. Even if the model remains elementary for the three phenomena, the mere coupling made it possible to clarify the appearance of bimodality. In order to illustrate the results, a particular choice of isotherms and of the space variation of the solidification rate was taken, but bimodality is generally obtained for nonlinear isotherms, due to the role of major elements in crystallization. The appearance of composition spatial discontinuities or shocks (that give rise to bimodality) can be qualitatively understood in the simplest cases of non-linear convex or concave isotherms. A convex isotherm depicts an enrichment of the chemical component in the solid as compared to the liquid. The rapid solidification (quenching) of such a liquid will give an "initial" composition of the solid that will be impoverished as compared to the "final" composition, i.e. that in equilibrium with the liquid. The respective positions of the initial point and of the final point on the isotherm can thus make a shock appear, according to the standard chromatography rules. A concave isotherm would equally lead to the appearance of shocks.

What are the geological factors that can account for the bimodality? The systems must have sufficient size, the evolutions must be long enough at elevated temperatures as to permit the chemical exchanges. These must be efficient (high parameter B value). There is an optimum range of values for parameter A: a too rapid solidification as compared to transport does not allow the bimodality: the system solidifies without chemical transfer. When solidification is too slow, there is not enough disequilibrium for shock appearance. Our model relates to any situation where a magma crystallizes (or equivalently where a rock melts) and where a subsequent reworking of chemical element distribution is made thanks to a differential movement or a separation between the solid and the liquid whatever the mechanism. If there is a tendency to local equilibrium, the process may give rise to products with differing proportions, even for homogeneous starting conditions. If, at a given time, the liquid escapes from such systems, its composition will be "quantified" in space. As a consequence, the emitted compositions within a volcanic system will vary in time along more or less discrete steps. We must then stress that our model thus explains possible bimodality both for the liquids and the solids: from the knowledge of the solid and liquid profiles in the chamber in figures 3, one could draw histograms for the solids as well as for the liquid: both would show possible (correlated) bimodalities. Generally, the standard interpretation for the above geological systems, with varied products separated by discontinuities, calls upon several distinct events. We stress that one phenomenon may be enough. If the initial compositions are not homogeneous, the gradients can degenerate into discontinuities. In open system, one already knows that the arrival of a magma in a rock system in disequilibrium can give rise to the propagation of discontinuity waves [START_REF] Guy | Mathematical revision of Korzhinskii's theory of infiltration metasomatic zoning[END_REF]. Natural systems undoubtedly show intermediates between purely closed behaviours and purely open behaviours.

The simple model with one independent chemical component can be applied for systems with one solid solution phase, or a set of solid solution phases described by a single chemical parameter (e.g. the Fe/Fe+Mg ratio of olivines, the Al ratio of chromites, the non-linear fractioning of trace elements etc.), all other factors being equal. This may apply to small mono-mineralic magmatic chambers.

As a whole, the proposed model provides an intellectual framework to discuss the natural cases and the basis for future more complicated models. The development of our approach is related to our simplifications of the phenomena:

texture, crystal size, absence of size variation, absence of chemical zoning etc. Absence of size variation may not be a problem insofar as the increase of size during crystallization is compensated by faster descent. Generally speaking, one can say that the various effects related to the size of the crystals are hidden in the kinetics and the movement laws.

-Choice of chemical equilibration models: natural evidence of solid diffusion re-equilibration of crystals is generally poor. The more likely is that chemical disequilibrium induces solid melting (because of solidus variation with composition) and recrystallization makes the stable composition appear. This is similar to what occurs during reactive transfer of aqueous solution in porous media (dissolution recrystallization operates rather than true exchange of components by solid diffusion). Whatever the physical phenomenon taking place however, the modelling of chemical exchange by macroscopic laws involving surface transfer and thermodynamic disequilibrium gives a good agreement with reality. The heating effects are not taken into account insofar as the total balance of the transformation is practically null.

Modelling by diffusion well accounts for the mass transfers necessary for equilibration (diffusion coefficient may be chosen intermediary between that of the liquid and the solid).

The magnitude obtained for K parameter is similar to that for the water/rock interactions and can vary on several orders depending on the systems. This homogenizes the aspects of texture, size and mode of matter displacement.

-isotherm concept: taking into account more complicated systems of equations does not present an obstacle. The isotherm concept as a whole gathers all the mineralogical effects within a single formalism; the nature of the mineral which precipitates is simply a consequence of the value of the solution parameter for a more or less complex isotherm (with a plateau for example).

-taking into account of only one chemical component : we consider here that the phenomena involving that particular chemical component may be discussed independently from the other components of the magma. In reality, other chemical elements, like SiO 2 , play the role of solvent, which, at a given time in the chemical evolution, will give place to new mineral crystallizations and interfere with the other phenomena. The presence of several simultaneous minerals may also involve several simultaneous movements with possibly different speeds and complex interferences.

-simplification of the overall heating/cooling effects : the heating effects interfere with the other effects, in particular chemical. Part of the temperature pattern is controlled from the interior of the system and not from outside. The isotherm is independent of temperature (one could easily add a dependence). Similarly the movement law will depend on temperature.

These dependencies may be hidden within a dependency on p, the evolution of which globally corresponds to a decrease in temperature.

Caption of figures

Schéma d'une chambre magmatique. L'axe vertical représente la profondeur sans dimension figurée par la variable x, variant entre 0 et 1 (correspondant à des épaisseurs de quelques centaines de mètres pour des chambres naturelles). La solidification peut se faire dans la partie supérieure de la chambre ; les cristaux formés peuvent ensuite se sédimenter et remplir le bas de la colonne (partie solide).

Figure 1

Schematic view of a magmatic chamber: the vertical axis is the dimensionless x axis used in the text (variation between 0 and 1; it corresponds to thicknesses of several hundreds of meters for natural chambers). The solidification may take place in the upper part of the chamber; the crystals formed subsequently settle and fill the bottom of the chamber (solid).

Figure 2a

Fonction d'espace s*(x*) pour la solidification, en fonction de la variable sans dimension x*.

Cette fonction intervient en facteur de la vitesse de solidification q 0 exprimée en termes de production de proportion de solide p. La solidification est ici maximale au sommet de la chambre (s* = 1 pour x* = 0) et nulle en bas (s* = 0 pour x* = 1).

Figure 2a

Space function for dimensionless solidification s*(x*); this function intervenes as a factor to the speed of solidification q 0 expressed in production of proportion of solid p. Solidification is maximum here at the top of the chamber (s* = 1 for x* = 0) and null at bottom (s* = 0 for x* = 1).

Figure 2b

Isotherme sans dimension C m * = g(C S *). L'isotherme a deux parties, représentant chacune une phase solide, séparées par un plateau horizontal représentant l'équilibre entre les deux solides.

Figure 2b

Dimensionless isotherm C m * = g(C S *). The isotherm has two parts, representing each a solid phase, separated by a horizontal plateau corresponding to equilibrium between the two solids.

Figure 4a

Différents types d'évolution spatiale pour la composition du solide dans la chambre à la fin de l'évolution, pour différentes valeurs des paramètres sans dimension A et B; pour le profil plat, une composition est sélectionnée ; pour le profil avec une pente constante, une gamme de compositions équiprobables est sélectionnée. Pour le profil à deux plateaux, les compositions sont séparées en deux modes.

Figure 4a

Various types of spatial trends for the composition of the solid within the chamber at the end of the evolution, for various values of dimensionless parameters A and B: for the flat profile, one composition is selected; for the profile with constant slope, a range of compositions with same probability is selected; for the two-plateau profiles, the compositions are separated into two modes.

Figure 4b

Spectre de fréquences donnant les proportions des différentes compositions du solide dans la chambre (pour l'expérience précédente avec A = 1, B = 10). Il y a ici une forte bimodalité marquée par deux pics, l'un de composition proche de C = 0.9 et l'autre de composition proche de C = 0.3.

Figure 4b

Frequency spectrum giving the proportions of the various compositions of the solid within the chamber (preceding experiment, A = 1, B = 10). There is a strong bimodality marked by two peaks, one close to C = 0.9 and the other close to C = 0.3.

Figure 4c

Exemple de spectre unimodal (une seule composition de la roche, correspondant à celle du magma initial, est sélectionnée) pour les valeurs A = 0.1 et B = 0.1 des paramètres sans dimension.

Figure 4c

Example of uni-modal spectrum (one rock composition, corresponding to that of the initial magma, is selected) for the values A= 0.1 and B = 0.1 of the dimensionless parameters.

Figure 5

Carte des valeurs des paramètres sans dimension A (solidification / sedimentation) et B (reaction / sedimentation) permettant une bimodalité de composition des roches cristallisées. q or q s s -1 rate of solid proportion production q = q 0 s(x) 10 -7 to 10 -14 10 -12 q 0 s -1 average rate of solid proportion production 10 -7 to 10 -

f 2 .

 2 Modelling and defining dimensionless parameters2.1. Variation of the solid /liquid proportion by solidification and sedimentation

  the whole, the two dimensionless parameters A and B can vary between 0.1 and 10 for the most current geological situations, but the extreme values are in the bracket 10 -4 , 10 4 or even 10

  adjustable functions s * et g * were chosen as shown in Figs. 2a and 2b. The two dimensionless numbers A and B are chosen within the interval [0.1 -10]. The initial condition C m * (t=0), was fixed at 0.8 (an extreme value of 0 or 1 would not generate any chemical differentiation effect). A typical evolution of the compactness and composition profiles at three times are given in Fig. 3a, b, c, for A=1 and B=10. The distribution functions of the concentration in the solid at the end of the process are then drawn (obtained by kernel smoothing of the histograms) and we observe, in this case, a conspicuous bimodality in the distribution of the concentrations in the solid, whereas the initial condition in the magma was perfectly uniform. The profiles obtained, for various values of A and B (for the same s* and g* functions), give an idea of the various distributions of concentration: uni or bimodality according to the values of A and B, for the selected functions s* and g* (Fig. 4a). The corresponding bimodality is shown in Fig. 4b. In Fig. 4c, an example of rather uniform distribution of concentration for A=0.1, B=1 is shown. Going from dimensionless time to physical time, one sees that the end of evolution (for values of t* of several tens of units) corresponds to times of thousands to hundreds of thousands years for typical values of the chambers.

Figures 3a, b, c

 c Figures 3a, b, c Profils de différentes variables en fonction de la profondeur à différents temps. La loi de vitesse de solidification et l'isotherme sont donnés dans les figures 2a et 2b. Les paramètres sans dimension ont pour valeurs A = 1 et B = 10 respectivement. Le trait noir représente la proportion de solide p ; la sédimentation s'effectue au bas de la chambre, le front du solide remonte en fonction du temps. La condition initiale représente une concentration homogène avec C m = 0.8. La forme d'ensemble de l'isotherme est concave, le solide est enrichi par rapport au liquide. En remontant, on observe un appauvrissement du solide (courbe bleue),avec une inflexion représentant une discontinuité de composition entre deux plateaux. La courbe verte représente la composition du liquide et la courbe rouge celle du total liquide + solide. Si le système est figé, cette dernière courbe donnera la composition du solide final.

  Figure 1

Figure 3a

  Figure 3a

Figure 3b

  Figure 3b

Figure 3c

  Figure 3c

Figure 4a

  Figure 4a

Figure 4b

  Figure 4b

Figure 4c

  Figure 4c

Figure 5

 5 Figure 5

Table 1

 1

Acknowledgements : the authors thank Daniel Garcia, Jacques Moutte, Pierre Boivin, Marie-Christine Gerbe, Ariel Provost, the geologists of the UMR 6524 "magmas and volcanos" for the stimulative discussions, and two anonymous reviewers for helpful comments.

suivantes : vitesse de réaction chimique > vitesse de sédimentation > vitesse de solidification.

Figure 5

Map of the values of dimensionless parameters A (solidification/sedimentation) and B (reaction/sedimentation) allowing the chemical composition bimodality of crystallized rocks.

The zones of marked bimodality are in red. Bimodality is acquired in the condition: chemical reaction rates > sedimentation rates > solidification rates.

Caption of table Tableau 1

Liste des principaux symboles utilises et valeurs typiques des paramètres correspondant.

Table 1

List of principal symbols used and typical values of the corresponding parameters.