D Grasso

E Tassi

F L Waelbroeck

Nonlinear gyrofluid simulations of collisionless reconnection

The Hamiltonian gyrofluid model recently derived by Waelbroeck et al. [Phys. Plasmas 16, 032109 (2009)], is used to investigate nonlinear collisionless reconnection with a strong guide field by means of numerical simulations. Finite ion Larmor radius gives rise to a cascade of the electrostatic potential to scales below both the ion gyroradius and the electron skin depth. This cascade is similar to that observed previously for the density and current in models with cold ions. In addition to density cavities, the cascades create electron beams at scales below the ion gyroradius. The presence of finite ion temperature is seen to modify, inside the magnetic island, the distribution of the velocity fields that advect two Lagrangian invariants of the system. As a consequence, the fine structure in the electron density is confined to a layer surrounding the separatrix. Finite ion Larmor radius effects produce also a different partition between the electron thermal, potential, and kinetic energy, with respect to the cold ion case. Other aspects of the dynamics such as the reconnection rate and the stability against Kelvin-Helmholtz modes, are similar to simulations with finite electron compressibility but cold ions.

I. INTRODUCTION

The investigation of collisionless magnetic reconnection by means of Hamiltonian reduced fluid models has proved to be useful in various ways, for instance in the interpretation of nonlinear structures observed in simulations [1][2][3][4][5][6], the derivation of stability criteria [7], the extension of the model to include external fields while preserving a Hamiltonian structure [6], and the identification of negative energy modes [7]. The above mentioned results refer, however, to cold-ion Hamiltonian models. Energy-conserving, hot-ion (in particular, gyrofluid) versions of such cold-ion models have been derived, for instance in Refs. [8][9][10][11][12]. In particular, a gyrofluid two-field version of the model presented in Ref. [9] was investigated numerically in Ref. [13] and a dissipative version of it was studied in Ref. [14]. A separate application of this model to the study of Alfvénic turbulence appeared after completion of the present work. [15] Although different energy-conserving, hot-ion models are available, the nonlinear numerical investigation of Hamiltonian gyrofluid models for collisionless reconnection is relatively less developed with respect to that in the cold-ion limit. An example of results in this direction was presented in Ref. [13]. The analysis presented in that article, however, was limited to the linear and early nonlinear phase of the dynamics. The purpose of this article is to investigate the nonlinear dynamics of a Hamiltonian gyrofluid model and compare the results with those obtained in Ref. [13] and in previous investigations of cold-ion models. Particular emphasis will be given to the analysis of the structures that form nonlinearly in the fields. The structure of the electric field, for example, is of interest in theories of electron [16,17] and ion [18] energization during magnetic reconnection in the magnetopause, the solar corona, and in laboratory experiment [19,20]. Note, however, that energization is a kinetic process that cannot be described rigorously within the fluid model used here. We will investigate, with the help of the Hamiltonian structure of the model, the role played by the finite ion Larmor radius on the electron density distribution. The influence of finite ion temperature on the secondary fluid instabilities observed in recent years [4,5,21,22] will also be considered.

The paper is organized as follows: in Sec. II the model equations are introduced. Sec. III and Sec. IV are devoted to the analysis of the results of the numerical simulations, with focus on field structures and energy partition, respectively. Conclusions are drawn in Sec. V.

II. MODEL EQUATIONS

We consider the following Padé approximant, two-dimensional version of the Hamiltonian gyrofluid model of Ref. [12]:

∂n i ∂t + [Φ, n i ] = 0, (1)
∂n e ∂t + [φ, n e ] -[ψ, ∇ 2 ψ] = 0, (2)
∂ ∂t (ψ -d 2 e ∇ 2 ψ) + [φ, ψ -d 2 e ∇ 2 ψ] + ρ 2 s [ψ, n e ] = 0, (3)
n e = 1 1 - ρ 2 i 2 ∇ 2 n i + ∇ 2 1 -ρ 2 i ∇ 2 φ. (4)
Given a Cartesian coordinate system (x, y, z), we assume that all the fields are translationally invariant along z. The variables in Eqs.( 1)-( 4) are written in a dimensionless form and their relations with dimensional quantities are given by:

t = v A L t, x = x L , n i = L di ni n 0 , (5)
d e = de L , n e = L di ne n 0 , ψ = Âz BL , φ = ρ2 s L 2 L di e φ T e ,
where the carets denote dimensional variables. In ( 5) n 0 indicates a background density, L is a characteristic magnetic equilibrium scale length, v A is the Alfvèn speed based on a characteristic poloidal magnetic field intensity B, ni,e represent the ion guiding center and the electron density, respectively, Âz is the z component of the vector potential, ρs is the sonic Larmor radius, φ is the electrostatic potential, di,e the ion and electron skin depth, T i,e the ion and electron temperature, respectively, and e is the unit charge. In addition to this, we indicate with ρ i = T i /T e ρ s the ion Larmor radius, and with

Φ = 1 1 - ρ 2 i 2 ∇ 2 φ, (6)
the gyroaveraged electrostatic potential. Note that this normalization differs from that adopted in Ref. [12].

It is easy to see that in the limit ρ i → 0, the ions decouple and the above equations effectively reduce to the twofield, cold-ion model used in [2,23]. This model as well as the more complete model used here, describe the so-called "inertial-Alfvén" regime in the limit β e ≪ m e /m i (or equivalently, ρ s ≪ d e ) as well as the "kinetic Alfvén" regime in the opposite limit. In the previous expression, we indicated with m e,i the electron and ion mass, respectively, and with β e the ratio between the electron pressure and the magnetic pressure. These two regimes have been observed

in laboratory experiments [24] as well as in the auroral region with the Freja satellite. [25] The characteristics of magnetic reconnection in the two limits β e < m e /m i and β e > m e /m i have been compared and contrasted by Rogers et al. [26] III. ANALYSIS OF THE FIELD STRUCTURES

We perform numerical simulations of the model ( 1)-( 4) on the domain {(x, y) : -π ≤ x < π, -π ≤ y < π}, with a grid of 1024×128 points and imposing double periodic boundary conditions. The initial equilibrium is given by

n ieq (x) = n 0 , n eeq (x) = n 0 , ψ eq (x) = 11 n=-11
a n exp(inx), (7) where n 0 is a constant background density and the a n are the Fourier coefficients of the function f (x) = 1/ cosh 2 x. In

Ref. [27] it has been shown that such truncated Fourier series provides a very good representation of the equilibrium flux function. The equilibrium of the form 1/ cosh 2 x makes it possible to avoid the early cross-talking between magnetic islands that prevented the observation of a developed nonlinear phase in Ref. [13]. The equilibrium (7) allows to reach a nonlinear phase, in which the ratio between the island width and the electron skin depth is approximately twice as much as that permitted in Ref. [13], and eventually a saturation phase. The equilibrium is destabilized by perturbing the n i field with a four-cell pattern disturbance of the form ñi ∝ cos(x + y)cos(xy). The field φ is perturbed accordingly, in such a way that the initial perturbation on n e is zero. The choice (7) for the equilibrium implies that n i grows only very weakly compared to the other fields. Although this choice makes the role of n i on the dynamics essentially negligible, it permits to make a more direct comparison with previous results. An investigation in the presence of an ion guiding center nonuniform equilibrium will be the subject of a future publication.

Simulations have been run fixing d e = 0.2 and varying the values of ρ s and ρ i , in order to investigate the influence on the dynamics of finite ion temperature with respect to cold ion models or to the gyrofluid model of Ref. [13].

Linear growth rates observed in the simulations have been compared with the asymptotic formula γ ≈ 2(2d e ρ 2 τ /π) 1/3 , derived from linear theory [28], and a good agreement was found (here ρ 2 τ = ρ 2 i + ρ 2 s and the multiplicative factor 2 comes from the choice of the magnetic equilibrium adopted here).

We first carried out a comparison between simulations obtained by setting alternatively ρ s and ρ i equal to 0. From comparing the plots of the electron density one can see that the structures observed in the two simulations are qualitatively similar, in particular in the lobes around the magnetic separatrices. Such structures have already been observed in previous nonlinear simulations of different cold ion models (e.g. [2,4,6]). The inner regions look somewhat different, with a more accentuated quadrupolar structure in the ρ s = 0.4 case.

In order to understand the observed behavior we note the following: let us assume that the contribution coming from n i is negligible in (4). One can indeed verify that the amplitude of such term is typically much smaller than that of the other terms in the equation, namely due to the homogeneity of n i in the initial state. We then have the relation

φ = ∇ -2 n e -ρ 2 i n e . (8)
Through this relation we can eliminate φ in ( 2)-( 3) and obtain

∂n e ∂t + [∇ -2 n e , n e ] -[ψ, ∇ 2 ψ] = 0, (9)
∂ ∂t (ψ -d 2 e ∇ 2 ψ) + [∇ -2 n e , ψ -d 2 e ∇ 2 ψ] + (ρ 2 s + ρ 2 i )[ψ, n e ] + ρ 2 i d 2 e [n e , ∇ 2 ψ] = 0. ( 10
)
We can now see that in the evolution equation for n e (9) there is no explicit dependence on ρ i and ρ s . These parameters appear explicitly only in the electron momentum equation (10). Note, however, that the coefficient of [ψ, n e ] in (10), takes the same value in the two cases of our simulations. Therefore, apart from the term ρ 2 i d 2 e [n e , ∇ 2 ψ], the two systems we solve in the simulations are very similar. Thus, when it comes to determining the reconnection rate and qualitative field structures, interchanging the role of ρ i and ρ s can be expected to produce little differences, as already observed in Ref. [13]. This is reminiscent of the case with zero guide-field investigated in the "GEM challenge." In that case a comparison of kinetic, hybrid, and fluid (cold ion) codes showed that the reconnection rate is insensitive to the details of the plasma dynamics. [29] Subsequent work found that for all the models investigated in the GEM challenge, fast reconnection occurs when the phase velocity of the wave mediating the reconnection increases with the wavevector in a certain range. [26] Note, however, that more recent studies of reconnection in electron-positron plasmas have cast doubt on the role of the dispersion properties on the reconnection rate (see [30] and references therein).

In the case of the model considered here, when considering homogeneous equilibria n ieq = n 0 , n eeq = n 0 , ψ eq (x) = α ψ x, with constant α ψ , the linear dispersion relation resulting from Eqs. ( 1)-( 4) reads

ω = ±|α ψ k y | 1 + (ρ 2 i + ρ 2 s )k 2 ⊥ 1 + d 2 e k 2 ⊥ . (11)
In Eq. ( 11) ω is the frequency of propagation of the wave, k y is the wave vector along the y direction and k ⊥ is the perpendicular wave vector. Note that this dispersion relation includes the so-called "inertial-Alfvén" regime in the limit β e ≪ m e /m i or equivalently, ρ s ≪ d e , as well as the "kinetic Alfvén" regime in the opposite limit. [24,25] The expression (11) indicates that the inclusion of finite ion temperature does not modify the dispersive properties of the wave, with respect to the ρ i = 0, ρ s = 0 case. It simply amounts to introducing an effective sonic radius

ρ τ = ρ 2 i + ρ 2 s .
In this respect, the cold plasma case is qualitatively different.

Another similarity between the effect of the ion and electron temperature is that the secondary Kelvin-Helmholtz instability, observed in cold plasma models, [4,5,21,22] is suppressed when ρ i = 0, even if ρ s is very small. Indeed, even in the cold electron case, if finite ion temperature is taken into account, the formation of thinning layers prone to the fluid instability is suppressed by the appearance of lobes enclosed in the separatrices. We can then conclude that finite temperature effects, of either the electron or of the ion species, inhibit the onset of a secondary Kelvin-Helmholtz instability.

The similarities in the growth rates, electron density structures around the separatrices and Kelvin-Helmholtz stability properties, however, should not be allowed to distract from an important qualitative difference between the cold-ion and hot-ion cases: for hot ions, the electrostatic potential develops very fine structure, whereas for cold ions, it varies no faster than on the scale of ρ s . This difference, which is clearly visible in Fig. 1, is confirmed by comparing profiles of |∇φ| in Fig. 2. Such plots indeed show that, as ρ i increases, strong electric fields develop in the plane perpendicular to the guide field, in regions surrounding the separatrices. In order to highlight this aspect, we chose to show the profiles at the chord y/π = -1/2, which crosses such regions. Formation of fine structures and steep gradients in the electric field accompany this growth in amplitude. Evidently, the same trend can be read also in terms of an increase in amplitude and in gradients of the E × B drift. The fine structure of the electrostatic potential influences also the parallel electric field. Since E = -∂ψ/∂t -[φ, ψ], the latter field too develops small scales as ρ i increases. As an example, the contour plot of E is drawn in fig. 3 for the case with ρ i = ρ s = 0.4. Fig. 3 shows clearly that all the area enclosed by the magnetic island is dominated by fine structures. This behavior of the parallel electric field, combined with the formation of density cavities, is suggestive [16] of a possible enhancement of the particle acceleration along magnetic field line and, consequently, of the plasma heating as ρ i is increased.

To understand the reason for fine-scale structure of φ, when ρ s is large, it is necessary to recall the conservation properties of the system. The study of these conservation properties [2,12] leads to the conclusion that the dynamical equations ( 1)-( 3) can be cast in the form of the following set of convection equations:

∂G j ∂t + v j • ∇G j = 0, (12)
where the convection velocities v j are prescribed in terms of stream-functions φ j by v j = ẑ × ∇φ j . Here the index j takes the values 0, +, and -. The fields G j are Lagrangian invariants given by G 0 = n i and G ± = ψ + d 2 e J ± d e ρ s n e , with J = -∇ 2 ψ, and the corresponding advecting stream-functions are φ 0 = Φ and φ ± = φ ± ρ s ψ/d e .

The Lagrangian representation of the dynamics leads to the "phase-mixing" picture of collisionless reconnection presented in Ref. [2]. According to this picture, the reconnection is enabled by the mixing of the invariants in a way analogous to the mixing of the distribution function during Landau damping. As a result of this mixing, the G ± fields develop increasingly fine structure, like cream stirred in a cup of coffee. In order to determine the smoothness of the original fields, such as n e and ψ, we need to solve for these fields in terms of the invariants.

ψ = 1 2 (1 -d 2 e ∇ 2 ) -1 (G + + G -); (13
)

n e = 1 2 (G + -G -)/d e ρ s . (14)
The first of these expressions shows that ψ is a smoothed version of the mixed-up G ± fields, the inverse Kelvin-Helmholtz operator (1d 2 e ∇ 2 ) -1 having the effect of suppressing all scales below d e . The second expression, Eq. ( 14), by contrast, shows that n e is not smoothed at all, and will contain all the fine scales produced by the mixing of the G ± . Applying the ∇ 2 operator to Eq. ( 13) shows that the axial current J, like n e , retains the fine structure of the G ± . The dark and light structures (six of each, see Fig. 1, case with ρ s = 0.4) in n e lying on and inside the separatrix can be recognized as being part of two spirals resulting from the roll-up of the invariants of the system in the counterclockwise and clockwise sense, respectively. [2,23] These spirals represent a historical record of the reconnection up to that point, much like geological strata preserve the record of tectonic subduction and buckling.

Returning to the question of the behavior of the electrostatic potential, inspection of Eq. ( 8) shows that φ will itself only be smoothed when ρ i = 0. For ρ i > 0, by contrast, φ contains an unfiltered contribution from n e that exhibits progressively finer structure as the reconnection progresses. That is, φ as well as n e and J, exhibit a cascade towards smaller scales. It is easy to see that the fine structure is retained by the electron perpendicular velocity, since v e = ẑ × ∇(φρ 2 s n e ). Thus, the gyrofluid model describes the formation of electron beams much like the ones observed in kinetic simulations. [16,17,31] The width of the beams is determined by the stretching and folding of the G ± , and their parallel velocity is determined by the electron momentum conservation in the ẑ direction.

Note also that if, on the other hand, one considers the limit ρ s → 0, then G ± ∼ ψ + d 2 e J and φ ± ∼ φ. Consequently, for small ρ s one expects not to see the stretching in opposite directions of G + and G -under the action of φ ± , because the magnetic contribution to the advection velocity fields v ± , which is what causes the opposite circulation of the two flows, gets suppressed. This can be seen in Fig. 1, when comparing the two plots of G + . In the case with small ρ s one does not see the spiral arms that are visible in the case with large ρ s and which are a result of the stretching in the clockwise sense caused by φ + (the case of G -would be identical but with a spiral winding up anticlockwise [1]). Therefore, in spite of the similar features in the electron density around the separatrices, the two cases differ considerably in the structure of the underlying invariants and of the corresponding mixing process.

The Fourier representation offers an alternative explanation for the behavior of φ and the gyro-averaged potential Φ. Let us write the periodic fields in Fourier series, so that n e = +∞ k=-∞ n e k (t) exp(ik • x) and analogously for the other fields. From (4) one obtains

φ k = 1 + ρ 2 i k 2 k 2 1 + ρ 2 i 2 k 2 n i k - 1 + ρ 2 i k 2 k 2 n e k , (15)
Φ k = 1 + ρ 2 i k 2 k 2 1 + ρ 2 i 2 k 2 2 n i k - 1 + ρ 2 i k 2 k 2 1 + ρ 2 i 2 k 2 n e k . ( 16
)
From these relations we deduce that for ρ i → 0,

φ k ∼ Φ k ∼ (n i k -n e k )/k 2 . ( 17
)
In this case, both φ and Φ are smoothed with respect to n e , which is proportional to the sum of the convected fields G + and G -, and with respect to n i which is itself a convected field.

For ρ i → +∞, by contrast,

φ k ∼ -ρ 2 i n e k ; ( 18
)
Φ k ∼ - 2 k 2 n e k . ( 19
)
Thus we can see that, as ρ i increases, φ tends to become proportional to n e whereas Φ remains screened, an effect of the gyroaveraging. The similarity in the structure of n e and φ is indeed what Fig. 1 shows. Finally note that the above argument, based on the quasi-neutrality condition, is independent of the value of ρ s and holds also in the presence of nonuniform background density, in which case the perturbations to n i , that were neglected in ( 9)-( 10), must be retained. Another effect related to the presence of ρ i is the flattening of the electron density as far as the central region of the island is considered. This effect can be explained in terms of the different behavior of the velocity fields v ± = ẑ × ∇φ ± which advect the Lagrangian invariants G ± . Indeed, as far as ρ i increases we observe progressively weaker velocity fields inside the island, as shown in fig. 4 (first column plots) for two different values of ρ i (0.01, 0.4) at the same value of ρ s = 0.4. Hence, the G ± are less strechted and filamented far from the separatrices and tend to coincide inside the island, as shown for a particular y-chord in the second column of fig. 4. Here we note that peaks, which occur where the velocity field is stronger, develop well inside the island for small ρ i , while they tend to be localized around the separatrices when ρ i increases. Since the electron density is recovered from the difference between the G ± fields through relation (14), it follows that its profile tends to become flatter and flatter as we approach the O-point of the magnetic island.

IV. ENERGY CONSIDERATIONS

The 3-field model ( 1)-( 4) possesses the following conserved energy integral (Hamiltonian) [12]:

H = 1 2 d 2 x |∇ψ| 2 + d 2 e (∇ 2 ψ) 2 + ρ 2 s n 2 e + Φn i -φn e . (20)
The various terms in (20) represent the magnetic energy (E B ), the parallel electron kinetic energy (E ke ), the electron thermal energy (E the ), the ion electrostatic energy (E eli ) and the electron electrostatic energy (E ele ), respectively.

We consider the temporal evolution of the different contributions to the total energy for ρ s = 0.4, ρ i = 0 and ρ s = 0, ρ i = 0.4. First we note that, with our choice of equilibria and initial perturbations, the total energy (corresponding to its value at t = 0) is given by

H = 1 2 d 2 x |∇ψ eq (x)| 2 + d 2 e (∇ 2 ψ eq (x)) 2 + Φ(x, y)ñ i (x, y) + 1 2 ρ 2 s 4πn 2 0 , (21)
where tildes denote the initial perturbations. The amount of energy turns out to be essentially the same for both cases, because the difference between the cases ρ i = 0 and ρ i = 0.4, present in the ion electrostatic energy, is less than 10 -5 and because the presence of the constant term in (21) depending on n 2 0 is irrelevant for the investigation of the energy deviations from their initial values. By comparing the two plots of Fig. 5 we can see that, although the total energy is the same in both cases, the distribution into the various forms is quite different. In both cases we decided to consider the simulations reliable only when the loss of total energy due to numerical dissipation is less than 1%. If we consider the plot for ρ i = 0 at t = 19 and the plot for ρ s = 0 at t = 16 (corresponding to magnetic islands of the same size) we see that essentially the same amount of magnetic energy has been lost in the two cases. This loss is, however, compensated in two different ways. For the cold ion case the reconnection process converts magnetic energy mainly into thermal electron energy (which prevails as ρ s is increased) and electron electrostatic energy. A smaller portion of magnetic energy also goes into parallel kinetic energy, although the latter starts to decrease when the total energy starts to be no longer conserved. We observe that in the cold ion limit the Hamiltonian can be written as

H = 1 2 d 2 x |∇ψ| 2 + d 2 e (∇ 2 ψ) 2 + ρ 2 s (∇ 2 φ) 2 + |∇φ| 2 . ( 22
)
Because in this limit n e = n i +∇ 2 φ the ion and electron electrostatic energy can be combined into a single electrostatic energy term which is also proportional to the perpendicular kinetic energy due to the E × B flow. The conversion can then of course be interpreted also as partial transformation into perpendicular kinetic energy.

In the opposite limit (ρ s = 0, ρ i → +∞), the Hamiltonian becomes

H ∼ 1 2 d 2 x |∇ψ| 2 + d 2 e (∇ 2 ψ) 2 + n i Φ + ρ 2 i 4 (∇ 2 Φ) 2 . ( 23
)
In this limit the electron thermal energy vanishes and almost all the magnetic energy is converted into electron electrostatic energy, which in this limit becomes proportional to the square of the laplacian of the gyroaveraged potential. In particular, if one assumes that ψ = O(1) and n e = O(1) as ρ i → +∞, then φ = O(ρ 2 i ), and the electron electrostatic will tend to grow and dominate the other terms in that limit. We note also that, compared to the cold-ion case, the presence of finite ion temperature determines a decrease in the parallel kinetic energy as shown in Fig. 6, where a comparison between two cases with ρ s = 0.2 and ρ i = 0, 0.2, respectively, is considered.

In the limit of negligible n i we can recover the following expressions:

d dt 1 2 d 2 x |∇ψ| 2 + d 2 e (∇ 2 ψ) 2 = d 2 x∇ 2 ψ[∇ -2 n e , ψ] + (ρ 2 i + ρ 2 s ) d 2 x∇ 2 ψ[ψ, n e ], (24)
d dt ρ 2 s 2 d 2 xn 2 e = -ρ 2 s d 2 x∇ 2 ψ[ψ, n e ], (25)
- d dt 1 2 d 2 xφn e = -d 2 x∇ 2 ψ[∇ -2 n e , ψ] -ρ 2 i d 2 x∇ 2 ψ[ψ, n e ]. (26)
These relations show how the sum of the magnetic and parallel kinetic energy, which is essentially the only form of energy available at the initial state, is transferred into different channels. Electron temperature terms are a source for electron thermal energy, whereas the corresponding ion temperature term and the convective term on the right-hand side of ( 24) provide a source for the electron electrostatic energy.

V. CONCLUSION

We have shown that, in some respects, ion gyration (parametrized by ρ i ) plays a role similar to electron parallel compressibility (parametrized by ρ s ). Both give rise to an acceleration of reconnection compared to the case of cold plasma (ρ s = ρ i = 0), and both have the effect of stabilizing the Kelvin-Helmholtz instabilities that have been observed in the collisionless reconnection of cold plasma. [4,5,21,22] In both cases, also, the acceleration of reconnection is consistent with changes in the dispersion relation, fast reconnection occurring when the dispersion curve has a convex segment. [26] Lastly, when ρ s is finite, also in the presence of finite ion temperature, the acceleration observed in hot plasma can be interpreted as resulting from the roll-up of the two Lagrangian invariants G ± by the convection cells associated with the reconnection. Such roll-up is energetically favorable since the magnetic energy is given by the integral of |∇ψ| 2 , and ψ is determined from the difference of the G ± by a smoothing operation. The smoothing filters the short scales, so that the magnetic energy in the rolled-up state is reduced from that in the initial state [2].

Beyond their superficial similarities, however, the cases of cold and hot ions exhibit important differences. Whereas the electrostatic potential is smooth in the case of cold ions, for hot ions it retains the fine structures created by the roll-up of the G ± . As a result, the finite-ρ i model exhibits electron beams similar to those observed in particle simulations. [16][17][18]31] This feature is important for the study of the particle energization occurring during magnetic reconnection. Moreover, another feature of finite ion temperature is a flattening of the electron density inside the magnetic island. This can be explained in terms of a decreasing of the speed of v ± going from the separatrices to the inner region of the island. Also the energy distribution is affected by the presence of hot ions. Indeed, when finite ion temperature is taken into account a relative increase in the electron electrostatic energy is compensated by a decrease in the parallel kinetic and electron thermal energies.

Lastly, we note that while particle codes are highly effective for small guide fields, their advantages recede for large guide fields and for current systems of macroscopic size. For such problems, the present paper shows that gyrofluid codes may prove to be valuable research tools. A good understanding of gyrofluid dynamics may also prove useful in the development of implicit moment algorithms [32] for gyrokinetic codes, since the gyrofluid equations are given by the moments of the gyrokinetic equation [33]. It is possible to observe that, as ρi increases, the amplitude of the peaks becomes larger and steeper gradients form in the perpendicular electric field. At the times which the three plots refer to, the corresponding magnetic islands have reached the same width. In both cases ρs = 0.4 and de = 0.2. By comparing the plots of the velocity fields one observes that for small ρi there exist high speed regions also inside the island (overplotted in red). This favours the formation of regions in which the difference between the local values of G+ and G-is enhanced. For large ρi, on the other hand, the velocity field is weaker inside the island, with respect to the separatrix regions. As a consequence, G+ and G-tend to coincide inside the island and therefore the electron density is locally flattened.

FIG. 5: Temporal evolution of energy deviations from the initial values, according to eqs. ( 20) and ( 21), for ρs = 0.4, ρi = 0 (left) and ρs = 0, ρi = 0.4 (right). In both cases de = 0.2.

  4 and 0.01. The case ρ i = 0.4, ρ s = 0.01, although unphysical, makes it possible to observe and isolate the effect of ion temperature on the cold plasma state. Results are shown in Fig.1. Because the lengths of the initial transient phases are different in the two cases, we did not compare the fields at the same time but rather when the magnetic island has reached approximately the same width, which indicates that the two dynamics have reached the same degree of advancement.

FIG. 1 :FIG. 2 :

 12 FIG.1:(Color online). Contour plots of electron density ne, electrostatic potential φ, gyroaveraged electrostatic potential Φ and the Lagrangian invariant G+ for ρs = 0.4, ρi = 0.01 (left column) and ρs = 0.01, ρi = 0.4 (right column). Contour plot of the magnetic island at the corresponding time have been superimposed onto each plot. The color scale goes from black (minimum value) to white (maximum value). In the plot of ne for ρs = 0.4, ρi = 0.01 dashed black and white lines are superimposed in order to show the presence of the two spirals resulting from the roll-up of the G± invariants. The black (white) spiral arm connect regions of high (low) electron density. The value of the electron skin depth is de = 0.2.

FIG. 3 :

 3 FIG.3:(Color online). Contour plot of the parallel electric field E = -∂ψ/∂t -[φ, ψ] for ρs = 0.4, ρi = 0.4 and de = 0.2. The plot exhibits fine structures which are determined by the filamentation of φ that occurs as ρi increases.

Acknowledgments

The authors acknowledge fruitful discussions with Dario Borgogno. This work was supported by the European Community under the contracts of Association between EURATOM and ENEA and between EURATOM, CEA, and the French Research Federation for fusion studies. FW was supported in part by a fellowship from the Katholieke Universiteit Leuven and by the US Department of Energy under contract DE-FG05-96ER-54346. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Financial support was also received from the Agence Nationale de la Recherche (ANR EGYPT).


F Bache

J Moreau

J L Rubino

C Gorini

B Van-Vliet

Lanoë

The Late Palaeozoic glaciation subsurface record, Chaco Basin (Bolivia)

Keywords: Late Palaeozoic Glaciation, Carboniferous, Chaco Basin, Bolivia, Palaeovalley

Late Palaeozoic glaciation is the longest of the Phanerozoic era. It is recorded in numerous Gondwanian basins, some having a high petroleum potential like the Chaco Basin.

In this basin, the quality of the available seismic, well and outcrop data permits to characterise the Late Palaeozoic glacial record. Palaeovalleys >500 m deep and ~7 km wide have here been analysed. Focusing on the glaciogenic Carboniferous deposits, the seismic data with well-ties and their outcrop analogues provide new sedimentological insights. The palaeovalley infill is imaged as a chaotic seismic facies overlain by an aggrading-prograding prism, interpreted as tillites covered by a fluvio-deltaic system respectively. Tillites form both under the ice and during rapid ice recession whereas fluvio-deltaic systems can only originate from a stable ice margin and last until the ice sheets withdraw inland. These two depositional

INTRODUCTION

The Late Palaeozoic Ice Age, developing from the Late Devonian (c. 390 Ma) to the Mid Permian (c. 270 Ma [START_REF] Hambrey | Earth's Pre-Pleistocene Glacial Record[END_REF]), generated the largest cyclic sedimentary deposits in Phanerozoic history [START_REF] Caputo | Late Devonian and Early Carboniferous glacial records of South America[END_REF]. It affected a large part of the supercontinent Gondwana, leaving evidences of ice activity in Australia, Antarctica, Arabic Peninsula, India, South Africa and South America [START_REF] Crowell | The late Paleozoic glaciation[END_REF]. The ice masses first developed in South America and South Africa during the Early Carboniferous.

During the Late Carboniferous, the ice sheets were situated in India, Australia and Antarctica.

The peak of the ice age occurred in Antarctica and Australia during the Permian [START_REF] Crowell | The late Paleozoic glaciation[END_REF][START_REF] Caputo | Migration of glacial centers across Gondwana during Paleozoic Era[END_REF]. The Late Palaeozoic glaciation coincides with the greatest episode of coal accumulation on Earth and consequent atmospheric CO 2 sequestration.

In South America, Late Palaeozoic glaciogenic rocks were first observed in the early 19 th Century [START_REF] Wegener | Die entstehung der kontinente und ozeane (The origin of continents and oceans)[END_REF][START_REF] Du Toit | Our Wandering Continents[END_REF]. Ice ages were usually examined by field studies, based on the occurrence of diamictites containing striated clasts attributed to icerafted debris or on glacial striated pavements recognition. Glacial valleys later became the main evidence of direct ice activity and, consequently, have been searched for within glaciogenic successions. Field mapping of the glacial surfaces led to the recognition of different types of glacial valleys associated with continental ice-sheets: (1) tunnel valleys, related to the incision of overpressurised meltwater; (2) cross-shelf troughs, related to ice-stream activity and abrasion processes; and (3) fjords which result from melt-water and abrasion erosion processes. Despite the development of subsurface seismic studies, palaeofjord incisions have rarely been described for ancient ice ages.

This paper aims at describing the seismic signature of the Late Palaeozoic succession in the Bolivian Chaco Basin which contains some very large-scale glacial incisions. The seismic stratigraphy provides new insights into the tectonics, the palaeoenvironment and the palaeoglaciology of the Chaco Basin.

REGIONAL SETTING

The Bolivian Chaco Basin is part of foreland basins associated with the Andean orogenic system and is identified as a backarc flexural intracratonic basin of Cenozoic age [START_REF] Sempere | Late Oligocene-Early Miocene major tectonic crisis and related basins in Bolivia[END_REF][START_REF] Decelles | Early to Middle Tertiary foreland basin development and the history of Andean crustal shortening in Bolivia[END_REF]. The basin is bounded to the north and north-east by the Brazilian craton, and to the west by the subandean fold and thrust belt zone which constitutes the eastern border of the Andes (Fig. 1). Since the late Oligocene, the deformation propagated toward the E and NE, uplifting the eastern part of the Andean depocentre [START_REF] Sempere | Phanerozoic evolution of Bolivia and adjacent regions[END_REF]. Consequently, the Phanerozoic strata are well exposed in the Bolivian Andes. A Precambrian-Tertiary clastic stratigraphic section is exposed in the Andean and the sub-Andean fold belts of southern Bolivia (Fig. 2), where shallow marine Ordovician-Devonian clastics are conformably overlain by a Carboniferous to Upper Permian sequence.

The Palaeozoic succession is unconformably overlain by Cretaceous fluvial sediments and Cenozoic progradational foreland basin deposits [START_REF] Dunn | Structural styles and hydrocarbon potential of the sub-andean thrust belt of southern Bolivia[END_REF].

In the Chaco Basin, the Izozog High (Fig. 3) is the result of an uplift that occured prior to the Cretaceous Period [START_REF] Gohrbandt | Paleozoic paleogeographic and depositional developments on the central proto-Pacific margin of Gondwana : their importance to hydrocarbon accumulation[END_REF][START_REF] Uba | Evolution of the late Cenozoic Chaco foreland basin, Southern Bolivia[END_REF]. Paleozoic units are tilted toward the W and NW and eroded at the top of the High (Fig. 4). This erosion is in part responsible for the variation in the thickness of the Paleozoic units (Fig. 5). The Carboniferous section is characterized by a continuous succession of well-developed massive sandstone units (Tupambi, Chorro, and Escarpment Fms.; Fig. 2) that alternate with thinner, muddy diamictite intervals (Itacuami, Tarija, Taiguati, and San Telmo Fms.; Fig. 2). This alternating succession is interpreted to result from episodic tectonic or climatic events [START_REF] Helwig | Stratigraphy, sedimentation, paleogeography, and paleoclimates of Carboniferous ("Gondwana") and Permian of Bolivia[END_REF]. Deeply incised palaeovalleys, typically 500 m deep and several kilometers wide, have been described in the Tarija and Escarpment Formations (Fig. 6; [START_REF] Helwig | Stratigraphy, sedimentation, paleogeography, and paleoclimates of Carboniferous ("Gondwana") and Permian of Bolivia[END_REF][START_REF] Salinas | Exploration del Sistema Carbonifero en la Cuenca Oriental de Bolivia[END_REF][START_REF] Tankard | Structural and tectonic controls of basin evolution in southwestern Gondwana during the Phanerozoic[END_REF]. The palaeovalley infills have been interpreted as continental tillites and subaerial meltwater channels, deposited by ice margins advancing and retreating, respectively [START_REF] Helwig | Stratigraphy, sedimentation, paleogeography, and paleoclimates of Carboniferous ("Gondwana") and Permian of Bolivia[END_REF][START_REF] Salinas | Exploration del Sistema Carbonifero en la Cuenca Oriental de Bolivia[END_REF]. The ice masses flowed across the Chaco Basin and were fed from ice centres to the S and E, situated within the Brazilian Craton [START_REF] Helwig | Stratigraphy, sedimentation, paleogeography, and paleoclimates of Carboniferous ("Gondwana") and Permian of Bolivia[END_REF][START_REF] Salinas | Exploration del Sistema Carbonifero en la Cuenca Oriental de Bolivia[END_REF]. To the contrary, [START_REF] Eyles | Hydrocarbon-bearing late Paleozoic glaciated basins of southern and central South America[END_REF] stress the importance of marine sedimentation and the paucity of evidence for any direct glacial influence on sedimentation. The Late Devonian-Early Carboniferous orogeny and Gondwana glaciations develop simultaneously and therefore, are usually considered to be associated with each other [START_REF] Amos | Las Cuencas Carbónicas y Pérmicas de Argentina[END_REF][START_REF] Eyles | Glacial geologic confirmation of an intraplate boundary in the Paraná basin of Brazil[END_REF][START_REF] Eyles | Glacioepochs and the supercontinent cycle after ~3.0 Ga: tectonic boundary conditions for global cooling: Palaeogeography[END_REF].

DATA AND METHOD

The dataset consists of 9000 km of 2D seismic lines and 4 wells tied to the seismic (Fig. 1). Seismic interpretations were carried out based on the principles of seismic stratigraphy [START_REF] Vail | Seismic stratigraphy and global changes of sea-level, part 3: Relative changes of sea level from coastal onlap[END_REF]. Seismic units are identified based on reflection terminations and the configurations of seismic reflections (Fig. 7). The units can be mapped in 3D throughout the basin thanks to the extensive coverage of the seismic data.

Additional information (datings, e-logs) were obtained from boreholes that sampled the Carboniferous units (Fig. 8). The gamma ray logging records the radioactivity of a formation. Shales (or clay-minerals) commonly have a relatively high gamma radioactive response and, consequently, gamma ray logs are considered to reflect the main grain size. Sonic logs measures the velocity of sound waves in rock. In complement to the gamma ray, sonic logs help to determine the lithology. The sonic logs have also been used to tie the well to the seismic lines.

The Paleozoic strata have been deformed after their deposition (Izozog deformation).

As the tectonic activity is supposed to have occurred at least after the Devonian, a Silurian horizon (Figs. 3 and4) has been used as a datum for the flattening of the seismic lines. At a regional scale, this method allows restoring the genuine morphology of the seismic stratigraphic surfaces.

This paper focuses on the upper part of the Carboniferous succession where two major stratigraphic surfaces have been picked and correlated throughout the basin (surfaces 4 and 5; Fig. 7). The geometric parameters of the observed incisions marked on these surfaces have been extracted from the studied (1) seismic strike sections for the width and depth, and (2) the flattened regional maps of the picked surfaces to define their longitudinal extension and drainage patterns. The described methodology allowed a detailed investigation of the seismic architecture of the Late Palaeozoic deposits of the Chaco Basin.

IZOZOG DEFORMATION

Glacial activity is commonly considered to be associated with a tectonic trigger [START_REF] Eyles | Glacioepochs and the supercontinent cycle after ~3.0 Ga: tectonic boundary conditions for global cooling: Palaeogeography[END_REF]. However, the presence of tectonic activity in the Chaco Basin during the deposition of the glaciogenic rocks of Late Palaeozoic is debatable. Previous mapping of the Carboniferous palaeovalleys (Fig. 6) suggested an influence of the Izozog High during their formation. [START_REF] Eyles | Hydrocarbon-bearing late Paleozoic glaciated basins of southern and central South America[END_REF] attributed the convergent direction of palaeovalleys to their formation in a confined basin.

Two main observations suggest, however, that the Izozog High formed after the Carboniferous, i.e. after the formation of the palaeovalleys. The first observation concerns the stratigraphic relationships between the Palaeozoic and the Mesozoic units. A major unconformity can be recognized on the seismic lines (yellow line, Fig. 3). This unconformity is of a regional extent and separates the Palaeozoic units from the Mesozoic ones (Fig. 8). The unconformity is deeply incised into the Palaeozoic rocks at the apex of the Izozog High (Fig. 5). The topography of the High is progressively onlapped and concealed by the Mesozoic deposits (Fig. 3). This erosion/deposition relationship is the evidence of a post-Palaeozoic tectonic uplift in the Izozog area. The second observation concerns the orientation of the Carboniferous palaeovalleys. The new maps of the palaeovalleys (Figs. 9 and 10) do not show any deflection over the Izozog High. On the contrary, the palaeovalleys are progressively more eroded toward the Izozog apex, indicating that the uplift is postdating the formation of the incisions.

PALAEOVALLEY MORPHOLOGIES

Ice masses are else flowing from highland areas (mountains) or lowland areas (quiet basins) and accordingly generate different types of valleys. In the Chaco Basin, the ice centres are situated onto the Brazilian craton [START_REF] Helwig | Stratigraphy, sedimentation, paleogeography, and paleoclimates of Carboniferous ("Gondwana") and Permian of Bolivia[END_REF][START_REF] Salinas | Exploration del Sistema Carbonifero en la Cuenca Oriental de Bolivia[END_REF]. Therefore, the here analysed glacial palaeovalleys are associated with lowland-ice sheets. Lowland-ice sheets usually generate three types of glaciogenic incisions visible on the seismic data (Table 1): (1) the tunnel valleys, (2) the ice stream cross-shelf troughs and (3) the fjords (not restricted to lowlands). (1) The tunnel valleys are frequently sinuous and are known to be organised in anastomosing to tributary branching patterns [START_REF] Huuse | Large-scale glaciotectonic thrust structures in the eastern Danish North Sea[END_REF]. They reach up to 4 km in width, 500 m in depth, 100 kilometres in length (e.g. [START_REF] O'cofaigh | Tunnel valley genesis: Progress in Physical Geography[END_REF][START_REF] Huuse | Large-scale glaciotectonic thrust structures in the eastern Danish North Sea[END_REF][START_REF] Ghienne | Palaeovalleys and faultcontrolled depocentres in the Late-Ordovician glacial record of the Murzuq Basin (central Libya[END_REF][START_REF] Praeg | Seismic imaging of mid-Pleistocene tunnel-valleys in the North Sea Basinhigh resolution from low frequencies[END_REF]). Tunnel valleys have been mainly observed in unconsolidated sediments but can also incise rocks. It has been assumed that tunnel valleys originate from subglacial pressurised meltwater incisions [START_REF] Wingfield | The origin of major incisions within the Pleistocene deposits of the North Sea[END_REF][START_REF] Brennand | Tunnel channels and associated landforms, south-central Ontario: their implication for ice-sheet hydrology[END_REF][START_REF] O'cofaigh | Tunnel valley genesis: Progress in Physical Geography[END_REF]. (2) The ice-stream cross-shelf troughs are U-shaped as well, they reach extreme widths (100' of km) and lengths (>600 km). Usually, their depths are subconstant along their flow lines at a regional scale, but they can vary significantly in tectonically active areas. On relatively passive margins like offshore Norway, they generate incisions between 200 and 500 metres deep and at least 10 km wide [START_REF] Ottesen | Ice-sheet dynamics and ice streaming along the coastal parts of northern Norway[END_REF]. These cross-shelf troughs are preferentially incised into soft unconsolidated sediments (e.g. the Norwegian channel northern boundary on hard rocks [START_REF] Ottesen | Submarine landforms and the reconstruction of fast-flowing ice streams within a large Quaternary ice sheet: The 2500-km-long Norwegian-Svalbard margin (57°-80[END_REF]). ( 3) As opposed to of (1) and ( 2), the fjords are deep glacial valleys incised into bedrock. Their depth reaches several kilometres while they stay relatively narrow in width (<10 km). Their longitudinal extent is 10's of kilometres but rarely exceeds 100 km (maps in Ottesen et al., 2008 for comparison). The drainage patterns of fjords are relatively straight, often guided by faults. They have a tendency to follow the structural grain onto the basement.

Fjords often correspond to the upstream part of ice-streams where they incise crystalline bedrock [START_REF] Ottesen | Submarine landforms and the reconstruction of fast-flowing ice streams within a large Quaternary ice sheet: The 2500-km-long Norwegian-Svalbard margin (57°-80[END_REF][START_REF] Ottesen | Ice-sheet dynamics and ice streaming along the coastal parts of northern Norway[END_REF].

In the Chaco Basin, the seismic analysis of the Late Palaeozoic sedimentary architecture shows an intricate succession of unconformities bounding the base of depositional sequences (Fig. 7). The peculiarity of these unconformities is the presence of large palaeovalleys, predominantly in the upper part of the Carboniferous succession. Both Uand V-shaped palaeovalleys are observed (Fig. 11). Their width ranges between 3 and 21 km with a mean width of 7.36 km (n: 24, standard deviation: 4.19 km). The incisions are deep with 500-700 m between the shoulders and the bottom of the thalweg. Two Upper Carboniferous surfaces have been mapped in detail (surfaces 4 and 5; Figs. 9, 10) and compared (Fig. 12). The mapped incisions are visible between the western limit of the dataset to the West and the Izozog High to the East, where they have been removed by erosion (Figs. 9, 10). They form extremely elongated depressions that were at least >100 km long. The palaeovalleys of surface 4 are subparallel with a predominent SE-NW direction. The palaeovalleys of surface 5 present a radial pattern with two different orientations in the North (SE-NW direction) and in the West (E-W direction, Fig. 13) of the Chaco Basin. The surface 5

shows two valleys merging and diverging twice to the South (Fig. 10).

The geometric parameters of the Chaco Basin incisions have been summarised and compared with the three known types of glaciogenic valleys (Table 1). Table 1 highlights that the studied glaciogenic incisions share characteristics with different types of glacial palaeovalleys under lowland-ice sheet conditions:

• The section shapes are comparable with fjords, indicating a mixed process abrasion (U-shape) and hydraulic (V-shape).

• The widths are clearly over the width of tunnel valleys but do not permit discriminating between fjords and cross-shelf troughs.

• The depths are similar to fjords and cross-shelf troughs but bigger than tunnel valleys.

• The lengths are comparable with ice stream troughs. Especially knowing that the lengths are underestimated because of the dataset extent and the Izozog High.

• The basement is assumed to be unlithified as in the tunnel valleys and cross-shelf troughs.

• They are straight to slightly sinuous like cross-shelf troughs.

• Except for two valleys on surface 5, they show no branching or anastomosing drainage patterns. This characteristic is also shared with cross-shelf troughs.

The size of the Chaco basin incisions rules out their possible interpretation as tunnel valleys. The anastomosed network to the South of Surface 5 (Fig. 10) is instead interpreted to result from the cross-cut between two generations of glacial valleys. This cross-cut is not resolvable with the available seismic data. However, it is difficult to discriminate whether the incisions represent cross-shelf troughs or fjords. The lengths of the Chaco incisions are an order of magnitude bigger than fjords. On the other hand, the section shapes and the depths are closer to fjords than to cross-shelf troughs. The Scandinavian ice sheet shows spatial transition between fjords on crystalline and lithified basements to cross-shelf troughs on sediments [START_REF] Ottesen | Ice-sheet dynamics and ice streaming along the coastal parts of northern Norway[END_REF].The shared characteristics of both cross-shelf troughs and fjords with the Chaco Basin valleys is therefore, interpreted to be the result of basement variations probably linked to varying lithification degrees of the subglacial substrate. These lithological variations are probably favouring abrasion-dominated processes (unconsolidated; cross-shelf trough type) or hydraulic-dominated incisions (lithified: fjord type).

The comparison of the drainage patterns on surfaces 4 and 5 shows a change in the orientation of the Chaco Basin incisions through time (Fig. 12). This indicates that the palaeoglaciological setting evolved between the two generations. The first generation (surface 4) is well organised with parallel troughs whereas the second generation of valleys (surface 5) radiates from a c. NW-SE axis (Fig. 12). This pattern change has been be interpreted to reflect one large ice mass forming surface 4 and individual ice sheets forming surface 5. The dispersive pattern may alternatively be the result of an ice margin closer during the formation of the surface 5 than the surface 4. This drainage pattern evolution at a regional scale is probably the result of a tectono-climatic evolution affecting the ice centres sourcing the Chaco Basin ice sheets during the Late Carboniferous.

PALAEOVALLEYS INFILL

The most detailed dataset concerns the sediments filling the palaeovalleys of surface 5.

The sediments have been attributed to the Escarpment and the San Telmo Formations (Fig. 2,8). Consequently, the analysis is focused on this stratigraphic interval based by surface 5 and topped by the next glaciogenic unconformity. The interval consists of three distinct sedimentary units with specific seismic facies, wireline log signatures, and sedimentary environments (Fig. 14). Facies 1 occurs at the base of the palaeovalleys (Fig. 14). It corresponds to a chaotic seismic facies suggesting an unsorted and unstratified accumulation of sediments. Gamma ray and sonic logs typically associated with diamictites can be observed (Fig. 14). Facies 2 is overlying Facies 1 (Fig. 14). Facies 2 is made of prograding clinoforms that can be seen downlapping on Facies 1. Gamma ray logs indicate coarsening-up shale to sandstone cycles (Fig. 14). On the seismic data, Facies 3 overlays unconformably facies 2.

Facies 3 has a channelised to slightly chaotic seismic facies. It is characterized by finingupward cycles of sandy units on gamma ray logs.

Facies 1, 2 and 3 have been interpreted as tills, delta progradations and fluvial deposits, respectively. This interpretation is supported by the field observations of the outcrop analogue (Fig. 14). Consequently, it is suggested that after the subglacial incision of palaeovalleys, the sedimentary environment passed through different phases: (1) subglacial to ice-marginal environment forming the till when ice is in the Chaco Basin, (2) the ice leaves the basin (3) the remnant glacial topography is progressively buried by an advancing fluviodeltaic system. This sedimentary system which concealed the 500-700m deep valleys was probably fed by ice margins further inland to the SE.

Analogues to the palaeovalleys infill patterns in the Chaco Basin can be found within Upper Ordovician palaeovalleys. The Ordovician valleys are filled by glacial sequences usually showing a motif with coarse-grained deposits (conglomerates, diamictites) at the base, passing upsection into prograding deltas and further into aggrading fluvial deposits [START_REF] Le Heron | Ancient glaciations and hydrocarbon accumulations in North Africa and the Middle East[END_REF]. The nature of the fill varies in style depending upon the palaeogeographic setting [START_REF] Le Heron | Sedimentary architecture of Upper Ordovician tunnel valleys, Gargaf Arch, Lybia: implications for the genesis of a hydrocarbon reservoir[END_REF]. In proximal settings, non-marine glaciofluvial sandstones overlie an initial fill of ice-proximal deposits. In deeper water settings, the initial, locally developed, iceproximal deposits are overlain by a transition from ice-distal diamictites to sand-dominated underflow fans. A deposition of this fill during deglaciation or marine transgression has been suggested [START_REF] Powell | Late Ordovician-Early Silurian glaciofluvial deposits preserved in palaeovalleys in South Jordan[END_REF][START_REF] Ghienne | Large-scale channel fill structures in Late Ordovician glacial deposits in Mauritania, western Sahara[END_REF][START_REF] Le Heron | Sedimentary architecture of Upper Ordovician tunnel valleys, Gargaf Arch, Lybia: implications for the genesis of a hydrocarbon reservoir[END_REF].

Considering the similarities between the successions in the Chaco Basin and the Ordovician, the palaeovalley infill of the Chaco Basin is considered to correspond to a glacial sequence.

As hundreds of metres of deltaic sediments overlay the basal till (Fig. 14), this glacial sequence is first developed in relatively deep water directly after ice retreat.

Terminal moraines are important indicators of former ice-front positions. Studies of exposed marine moraines have demonstrated that their internal facies architectures bear a high-resolution record of ice-front evolution, with direct implications for glacier dynamics and regional paleoclimatic conditions [START_REF] Lønne | Sedimentary facies and depositional architecture of ice-contact glaciomarine systems[END_REF](Lønne, , 1997)). [START_REF] Lønne | Sedimentary facies and depositional architecture of ice-contact glaciomarine systems[END_REF][START_REF] Lønne | Sedimentary Architecture and Dynamic Stratigraphy of a Marine Ice-Contact System[END_REF] and [START_REF] Lønne | Sedimentary Architecture and Dynamic Stratigraphy of a Marine Ice-Contact System[END_REF] have synthesized the development of marine ice-contact systems in the form of an allostratigraphic model. The model has implies a sedimentary architecture comparable to the sedimentary infill observed in the Chaco Basin palaeovalleys. In the southern part of the Chaco Basin, deltaic deposits (Facies 2) prograde in opposite directions (Fig. 15). In the Oslofjorden, the area between two clinoform orientations localise a former ice front [START_REF] Lønne | Sedimentary Architecture and Dynamic Stratigraphy of a Marine Ice-Contact System[END_REF]. Despite of the difference in the scale, this Quaternary example serves as an analogue; the pattern of clinoform orientations in the Chaco Basin is interpreted to mark the location of an ice front of a Late Carboniferous age (Fig. 16). Based on the Oslofjorden analogue, the glaciological and sedimentary evolution is deduced from the seismic analysis of surface 5 and its relative glacial sequence is summarised in a conceptual model presented in Fig. 16. The following multistage formation of the glacial sequence of surface 5 is suggested in the model:

(1) the ice advance and the creation of the valleys; (2) a still stand during ice recession in the basin and the generation of the southward-dipping proglacial fan-delta; (3) the recession of the ice from the basin associated with the isostatic uplift of the former ice front, leading to the reworking of the uplifted area and the creation of the northward-dipping clinoforms; (4) some minor readvance permitting the development of another fluvio-deltaic system prograding to the South and blanketing the northward-dipping clinoforms; and (5) the final ice recession associated with a last advance of the fluvio-deltaic system on top of the glacial sequence. This model honours the observed complexity of the glacial sequence and its potential palaeoglaciological record. In addition to ice-flow patterns deduced from the valley orientations, the localisation of ice-front positions and their consecutive depocentres are good constrains for understanding the evolution of ancient ice sheets.

CONCLUSIONS

The Late Palaeozoic sedimentary architecture of the Chaco Basin (Bolivia) is characterized by a succession of erosional surfaces and seismic sequences. This sedimentary architecture records the Late Palaeozoic glaciation events, and the modality of the ice age in the area. The main results of this study are:

• The Izozog High did not influence the Late Palaeozoic palaeovalley orientations.

Considering the progressive onlap of the Late Cenozoic deposits on the High, the Izozog Uplift is considered to be postdating the ice age and predating the Cretaceous fluvial deposits.

• Erosional surfaces present different morphologies from the base to the top of the succession. These surfaces can be flat (base of Carboniferous) whereas they form very large palaeovalleys (up to 20 km wide and 700 m deep) in the upper part of the succession.

• A lowland glacial environment is suggested, but the erosional processes responsible for the incisions remain debatable. A coupling between subglacial meltwater and mechanical abrasion processes is considered as a possible explaination for the observed valley geometries. In addition, the amount of meltwater at the base of the glacier (cold-/ warm-based or polythermal ice sheet), the nature and the lithification degree of the incised sediments probably induced the variability of the glacial erosional surfaces.

• The analysis of the glacial surface successions highlights the glacial history of the area, including ice flow rotations associated with changes in the glaciological parameters of the ice sheets.

• The palaeovalleys are covered by diamictites in the deepest part of the thalwegs overlain by deltas fed by fluvial depositional systems. As well as the Ordovician glacial sequences, the sand-prone Chaco Basin sediments present a good reservoir potential. This succession, corresponding to a glacial sequence, has been interpreted as the effect of glacier retreat. Depending of the distance to the ice margins, all or parts of the facies can be observed.

This paper aims at showing the intricate architecture of the glaciation record in the Late Palaeozoic sediments of the Chaco Basin. Although the origins of the glaciogenic features need to be investigated in more detail, this dataset highlights the importance of constructing a precise seismic stratigraphy for ancient ice ages. This work unravels a part of the glacial history in the area, built up of distinct events with different glacial settings and, thus, different climatic and tectonic triggers. This study localises ancient grounding lines of ice fronts which are palaeohighs containing proximal sandy facies. As they form isolated coarse-grained clastic depocentres, they represent good targets for hydrocarbon exploration.

Extending the analysis on each glacial depositional sequence is the next step to establish the context of the Late Palaeozoic ice age and its petroleum potential in the Chaco Basin. [START_REF] Helwig | Stratigraphy, sedimentation, paleogeography, and paleoclimates of Carboniferous ("Gondwana") and Permian of Bolivia[END_REF][START_REF] Salinas | Exploration del Sistema Carbonifero en la Cuenca Oriental de Bolivia[END_REF][START_REF] Tankard | Structural and tectonic controls of basin evolution in southwestern Gondwana during the Phanerozoic[END_REF]. Table . 1. Comparison between palaeovalleys observed in the Chaco Basin and different types of glaciogenic incisions associated with lowland-ice sheets.* From "pure" ice streams [START_REF] Bennett | Ice streams as the arteries of an ice sheet: their mechanics, stability and significance[END_REF]. M. =Mean.

FIGURE CAPTIONS

Fig. 1 .

 1 Fig. 1. Topographic map of Bolivia and subsurface dataset used for this study. Dotted lines

Fig. 2 .

 2 Fig. 2. Generalised stratigraphic column of the southern sub-Andean region with the source

Fig. 3 .

 3 Fig. 3. Time-depth of the top of the Silurian unit showing the location of the Izozog High

Fig. 4 .

 4 Fig. 4. Regional seismic lines showing the deformation of Paleozoic units by the Izozog High

Fig. 5 .

 5 Fig. 5. Time-thickness maps of Paleozoic and post-Paleozoic units in the Chaco Basin.

Fig. 6 .

 6 Fig. 6. Palaeovalleys described in the Carboniferous unit of the Chaco Basin Modified from

Fig. 7 .

 7 Fig. 7. Erosional surfaces identified in the Carboniferous. Palaeovalleys of surfaces 4 and 5

Fig. 8 .

 8 Fig. 8. Datation of seismic reflectors with boreholes. Location of boreholes on figure 1.

Fig. 9 .

 9 Fig. 9. Thickness map (in TWTT, left) between a continuous reference surface (top Devonian)

Fig. 10 .

 10 Fig. 10. Thickness map (in TWTT, left) between a continuous reference surface (top

Fig. 11 .

 11 Fig. 11. Seismic lines showing palaeovalleys of surface 5. TC: Top of the Carboniferous. BC:

Fig. 12 .

 12 Fig. 12. Comparison between palaeovalleys directions of surface 4 and surface 5.

Fig. 13 .

 13 Fig. 13. Detail of a palaeovalley of surface 5.

Fig. 14 .

 14 Fig. 14. Palaeovalleys infill and corresponding seismic facies, sedimentary environment and

Fig. 15 .

 15 Fig. 15. Seismic profile showing progradational facies (facies 2) in two different directions

Fig. 16 .

 16 Fig. 16. Conceptual model to explain the Carboniferous sedimentary depositional pattern of

ACKNOWLEDGMENTS

This study was funded by Total. We acknowledge Mads Huuse and journal reviewers for advice and comments on the manuscript. We would also like to thank Taija Torvela for its thoughtful comments.


S Kubler

W Derigent

A Thomas

E Rondeau

Problem definition methodology for the "Communicating Material" paradigm

Keywords: Communicating Material, System Engineering, Ubiquitous Computing

Today, more and more new technologies of information and wireless communication are useful for social, economic and industrial life. Auto-ID technologies such as RFID are a good example of that. The actual limits of these new systems are relative to their granularity, since items are the smallest parts of an ubiquitous system (Internet of things). In this paper we propose a "communicating material" new paradigm. Working on the credibility of this concept since some years, we highlight here its main technological and research perspectives.

INTRODUCTION

Currently, barcodes are progressively replaced by RFID tags (Radio-Frequency Identification) which enable a wireless communication with the "ambient environment". Products equipped with such devices are called communicating products and improve traceability, quality monitoring and logistic network performances because they have some information capacity offering to the system a better ability to react to unforeseen changes. An important aim of industrial management essential was to control the product material flow all along the supply chain. Right from the sixties, [START_REF] Plossl | La nouvelle donne de la gestion de la production[END_REF] showed that MRP2-based (Manufacturing Resource Planning) information systems could implement models reducing the gaps between the material and informational flows.This problem is still a hot topic today. An obvious way to suppress the synchronization problem between these two flows is to merge them : the material flow instrumentation via the use of technologies such as RFID is one solution, enabling products to be active by communicating information. A "communicating product" is thus an independent system, composed by a unique and single product connected to a communication device. As a result, it can wirelessly communicate with some computer facilities. This product can eventually carry information going from a simple reference (like a link to a database) to a more or less complex data set, accessible to the user when needed, by reading the tag attached to the product. Moreover, this product can eventually take simple decisions so that, under certain conditions, it could adapt its behavior to current circumstances, in spite of what was planned initially.

These last years, IFAC (International Federation of Automatic Control) [START_REF] Ollero | Milestone report of the manufacturing and instrumentation coordinating committee: From mems to enterprise systems[END_REF], [START_REF] Nof | From plant and logistics control to multi-enterprise collaboration[END_REF]] promoted the development of these concepts. More recently, the Holonic Manufacturing System (HMS) paradigm proposed the product should be given data storage and data processing capacities to interact with its environment [START_REF] Van Brussel | Reference architecture for holonic manufacturing systems[END_REF]]. The Intelligence in Manufacturing community [L Monostori . and Morel (2003), [START_REF] Morel | Editorial of special issue[END_REF]] demonstrated the pertinence of reconfigurable, modular and adaptable manufacturing systems. In order to implement the capacities exposed previously, Banaszak and Zaremba (2003) showed the product needs technologies not only to carry data but to carry information and knowledge. These products are either called intelligent or communicating products. Beyond these interpretations we will make no distinction between these two terms in this article.

Even so, this concept already showed a certain number of limits which are discussed later : discrete reading, risk of tag damage, problem of information transfer when cutting the product in several parts, . . . To overcome the previous problems, the concept of "communicating material" is proposed. It leads to an important paradigm change because the product does not communicate using a tag, but becomes intrinsically communicating, in his wholeness. In fact, in our vision, all the material which composing the product has this communicating property thanks to built-in specific micro-components. Besides, these components can give the product some capacities such as the data storage (the material would act as a computer hard disk and would update its data all along its life cycle), a decision making capacity (smart composite), or others.

The rest of article is organized as follows : A discussion about the paradigm change (passing from the communicating product to the communicating material) is introduced in section 2. Section 3 presents a state of the art on the concept of communicating products, by presenting works related to the intelligent product paradigm. Section 4 focuses on developing a rigorous innovation approach in order to organize hierarchically the objectives to achieve and to decline them in various research and technological problems. Finally, the last section uses this approach within the context of the communicating material.

OPPORTUNITY OF A NEW PARADIGM

Depending on studies that we have pursued for several years, we detail in this second section some reasons leading to the opportunity of this new paradigm (passing from the "communicating product" to the "communicating material").

Current research works about Intelligent Products within the framework of logistic suppose that each product is a small and non decomposable unit. In this context, each product has its own informational part, most of the time located on a single RFID tag which only stores the ID of the product, referring to its complete information saved in an external database. However, in a manufacturing process, products are rarely non decomposable and they are subjected to different transformations of several types (assembly/disassembly, welding, milling, sawing, . . . ). Depending on the transformation type, the information located on a certain area of the material might be moved to another location and perhaps modified. As a result, tools to manage data over the communicating product are needed. It is really difficult to implement these tools when considering the actual intelligent product architecture, because the informational part is stored in a fixed location on the product preventing it to be moved to another location. To solve this problem, the information will be continuously and intrinsically spread into the material. Moreover, accessing to the product information located in a distant database is not always possible and it could be interesting to store all the information into the product itself. But, currently, there are no commercial tag offering sufficient capacities coupled with small sizes. As explained before, the communicating material consists of a multitude of tag which offers an extensible and higher memory, according to the material dimensions and the tag density (formulated rows 1 & 2 Table . 1). The resistance to losses on a product integrating a single RFID tag is more critical than a object whose its material is entirely (and intrinsically) composed of tags. Indeed, if single tag is broken then all the information related to the product is lost (row 3 Table. 1). Similarly to communicating products, communicating material could provided several type of services (e.g. the monitoring, actuating notion or still information processing) [START_REF] Akhras | Smart composites -application and perspective[END_REF]].The main contribution of the communicating material is the capacity to provide services all over its surface (formulated row 4 Table . 1).

Thus, we aim to develop a new material which will be communicating from his natural state and in its entirety; it is the major distinction with the communicating product paradigm. Afterwards, it would be interesting to implement on the material the functions necessary for information management (to solve the previously mentioned problems). By the way, numerous projects interested in similar research topics are published as the SFIT axis (Smart fabrics and Interactive Textiles) from the next European FP7 call, which aims to develop new textile material with interactive functionalities. ISIS (a society hosted by the University of Brême) is currently working on sensory materials. The EMRS (European Materials Research Society) is also interested in this field and hosted a symposium around this topic during Spring 2010. Thus, one can easily see that giving enhanced functionalities to materials is a hot topic today, which seems to be promised to important research opportunities.

STATE OF THE ART

Today, works to make the material "communicating" are extremely rare. However, they find their sources in works aiming to give to products an active role in their life cycle (intelligent products) and those which goal is to functionalize the material, without making it intelligent (sensorial textiles). In what follows, we present works on these two research axes.

3.1

The "intelligent product" paradigm [START_REF] Meyer | Intelligent products: A survey[END_REF] offers a very complete state of the art on this theme. It seems that the term "intelligent product" has been mentioned for the first time in [START_REF] Ives | After the sale: leveraging maintenance with information technology[END_REF] in a context of after-sales services. In the IMS community, [START_REF] Mcfarlane | Auto id systems and intelligent manufacturing control[END_REF] is the first to show interest in intelligent products that he defines as being a physical product to which an informational representation is added. According to the intelligence type to implement within the product, various technologies are used. The "intelligence" concept, when it is interpreted as "information management", is often implemented through selfidentification technologies such as RFID. When this concept is used for decision making, more advanced technologies are used, like wireless sensors for instance.

Ubiquitous computer paradigm and functionalization of the material

The efforts of material functionalization have been numerous since the birth of the ubiquitous computing concept, formulated by [START_REF] Weiser | The computer for the twenty-first century[END_REF]. This concept proposes a world where computing systems, miniaturized to the extreme, are invisibly integrated in objects around us. These latest become then able to receive, store, process and transmit information.

Functionalizing the textile materials by adding electronic components takes place in the developed works within the topic of digital ubiquity. The smart textiles, or e-textiles, are materials combining microsystems, information technology and textiles. This e-textile participates in ambient intelligence and is related to the "wearable computer". Given our interpretation of "communicating material", e-textiles is to date the form of ambient intelligence that is the most similar to communicating material.

The ECE [Bradley Department of Electrical & Computer Engineering] from Virgina tech has worked for a long time on e-textiles. The issues concern computer architecture specifications, communication study (wired or wireless), software to implement in the e-textile, . . . [START_REF] Zeh | Softwear: A Flexible Design Framework For Electronic Textile Systems[END_REF]]. No instance of etextiles integrating solutions based on RFID tags spread into the material has been listed. Indeed, these tags are difficult to connect on the cloth. The works on e-textile are nevertheless numerous and concern a wide range of military [START_REF] Wilson | Electro-optic fabrics for the warrior of the 21st century[END_REF], [START_REF] Park | The wearable mother-board: An information infrastructure or sensate liner for medical applications[END_REF], [START_REF] Lind | A sensate liner for personnel monitoring applications[END_REF]], academic [START_REF] Stanley-Marbell | Modeling, analysis and self management of electronic textiles[END_REF]] or commercial applications.

Methods of requirement specifications

It seems that proposing a paradigm change of this type requires a system engineering approach rigorous enough to correctly setup the bases of a research field which could be broad. In order to do so, we wished to adopt a methodology to specify our Fig. 2. Schematization of the proposed method requirements and define the technological or scientific research fields inherent to this paradigm. As a result, in the following we expose some methodologies allowing to achieve these objectives by describing and presenting their limits, which justify our proposal. Thus, various approaches and communities have emerged, include :

• INCOSE (International Council on Systems Engineering) : Across all the system engineering tools offered by this community, there are methods and tools used by project executives, adequate to design, develop and check a system which offers an economic and highly efficient solution to the customer requirements.

• AFAV (French Association of Value Analysis) :

This organism proposes several methods and integrated tools, such as the functional and value analysis. The functional analysis allows to establish specifications, but it does not allow to raise the problems to process from the objectives. • 6 Sigmas : Originally, it is a quality self-assurance method relying on statistics and that which found an evolution as a design methodology. Since this method includes the requirement specifications and that our objective is to specify the research fields from product functionalities, this method is not appropriate in our case.

Consequently, we focus on an approach allowing to organize the objectives and to decline them in various technological and scientific problems. This method is illustrated in section. 4.

METHODOLOGY OF PROBLEM DEFINITION

Introduction

As previously stated, the paradigm change must be studied. In this paper, we develop an approach relying on some tools of previously mentioned methods. On the basis of a functional analysis while allows to list the functions that the product can realize, we implement a matrix approach being inspired by the QFD (Quality Function Deployment, [START_REF] Zaïdi | QFD. Quality Function Deployment[END_REF]). On the one hand, this approach allows us to identify the technological and scientific problems relating to the product and on the other hand, the correlations between these problems. Finally, using an approach based on spectral algorithms [Fiedler (1975a), [START_REF] Fiedler | A property of eigenvectors of nonnegative symmetric matrices and its applications to graph theory[END_REF]], these problems are gathered by problem families with strong correlations, that must be handled together. This approach is synthesized Fig. 1.

Quality Function Deployment (QFD)

The QFD approach takes into account the market needs and/or the future user wishes from the product design phase (or service), to develop the best design process in accordance with fixed quality imperatives. However, this approach initially needs to list the customer requirements. In order to do so, the product functional analysis is realized to extract all the requirements. These requirements are placed in the rows of the first level QFD matrix (Fig. 2 : Matrix level 1). They are also called What in the QFD method and form the set W 1 .

A value P w relative to the weight of the requirement w (with w ∈ W 1 ) is also specified. Next, we search to list all the product specifications (the How), composing the set H 1 enumerated by the columns of this matrix. Finally, we develop a relationship r w,h between What-How, based on 4 levels : 0, 1, 3, 9, where 0 means there is no relation between these elements and 9 a very strong one. Two influence indicators are calculated :

Inf h that gives information about the relative influence of a single specification h (h ∈ H 1 ) on all requirements :

Inf h = w∈W1 (P w • r w,h )
Inf w that gives information about the relative influence of all specifications on a single requirement w (w ∈ W 1 ) :

Inf w = P w • h∈H1 r w,h
Secondly, we reiterate this approach to create the second level QFD matrix (Fig. 2 : Matrix level 2). The specifications (How) of the first level matrix become the What to solve of the second matrix, which compose the set H 2 . By this way, the "Quality" deployment is assured by constraints spreading, since the specifications (H 1 ) become the What of the second matrix (W 2 ). The underlined problems (H 2 ), once solved, will then help to answer to the customer needs. In addition to the previous study, we here add the correlation half-matrix, being situated above the matrix (which exists for the matrix of level 1 in the original method, also known as the "roof" of the House of Quality).

The aim is to seek the correlations between various problems or

Px

Competitive purchase price 1 9 9 1 090 28 Cheap maintenance price 010 3 090 12 0.8 0.5 0.6 1.6 0.8 1.7 1.4 0.9 0.5 0.5 0.2 0.7 1.4 1 0.6 1.5 0.9 1.2 0.7 0.4 0.9 0.6 0.1 0.1 Fig. 3 How, which may be positive or negative and of variable values. Indeed, some process can damage themselves, others reinforce each other. In our case, they are expressed according to 5 levels : ++, +, 0, -, --. Indicators of the Matrix level 2 (Inf h , Inf w , r w,h ) are calculated like previously, taking into account not the sets

Q 1 , C 1 but the sets Q 2 , C 2 .
Note that we do not use the QFD tools in a classic way. Indeed the matrix of level 2 highlighting the problems, it is important for us to show the positive and negative interactions for organize hierarchically the problems to solve.

Organization of the problems to solve

The challenge is to organize the problems to solve (H 2 ) by considering their interdependence. The purpose is to gather these problems according to their negative or positive correlations. In order to achieve this step, we propose to define clusters of problem by relying on some spectral algorithms and especially on the serialization algorithm developed by [START_REF] Atkins | A spectral algorithm for seriation and the consecutive ones problem[END_REF]. The principle of this algorithm is the following : Given a set of elements n to order, and a correlation function f (i,j) , which corresponds to an attraction level between two elements i and j, in order to find them side by side in the sequence, the aim is to find optimized sequence between elements according to their consistency specified by the correlation functions. So if π is the permutation of elements and π

(i) < π (j) < π (k) then f (i,j) ≥ f (i,k) and f (j,k) ≥ f (i,k) .
The correlation half-matrix represents the relationship graph between problems. We have seen that these relationships can be positive or negative, so it is necessary to transform the half-matrix in a matrix suitable for the Atkins algorithm. In concrete terms, we substitute the symbols ++, +, 0, -, -by values 20, 10, 1, 50, 100 respectively. Only the matrix diagonal keeps the value 0. The minimal value 1 indicates that there is no interaction between problems but ensures that the graph is connected (essential condition to apply the spectral algorithm). The arbitrary choice of values 20, 10, 50, 100 highlights the correlation by overwhelmingly favoring the negative correlations, aiming to underline the problems which are the most difficult to solve. Once the problem sequence is determined via the Atkins algorithm, the purpose is then to form several problem clusters that will be simultaneously processed, these clusters composing the set C. To do so, the values 50, 100 are first gathered in at least one cluster. To facilitate the study, the clusters are limited to a reasonable size. Finally, we calculate for each cluster the contribution rate rt c to the final project (given below), where c ∈ C and H 2c is a subset of H 2 only composed of problems belonging to cluster c.

rt c = 100 • h∈H2 c Inf h h∈H2 Inf h
For instance, processing all problems belonging to cluster c whose contribution rate is 20 means answering to 20% of the whole problematic, i.e. the set of problems related to the product (illustrated before : H 2 ).

APPLICATION TO "COMMUNICATING MATERIAL"

In this section, we apply the approach elaborated before within the framework of the communicating material. We established functional specifications from a generic functional analysis resulting of various studies relating to user cases representative of the life cycle (scenario of traceability in factory to detect product defaults like holes, scenario of a custom officer checking a shipment of communicating materials in a truck, . . . ). Thus, we have formalized the customer requirements in various categories : Quality, Security and Price. The QFD matrix of first level is presented Fig. 3. In what follows, we will present the requirements, specifications and problems which are the most important. The aim of this paper is not to be exhaustive but to facilitate the understanding with some examples. We figure out from the indicator Inf h that the choice of the tag type, the volume of information or still the resistance against data losses are specifications having the greatest impact on the customer requirements. In the same way, the indicator Inf w shows that the requirements "must be able to be read/written", "must be able to contain the information", and "resist to undergone treatments by the material throughout its life cycle" are respectively the most important.

Next, we fill the second level QFD matrix and the correlation half-matrix, illustrated Fig. 4. Note that the specifications (H 1 ) related to requirements (W 1 ) concerning the "price" will be developed in the matrix of level 2, because the main objective is to raise the scientific and technological problems, without taking into account the financial aspect (for the present).

Then, we make the serialization of problems (via the Atkins algorithm), presented by the matrix Fig. 5. We define 4 clusters :

• Instrumented the material : it is mainly composed of issues related to the achievement of the communicating material in its natural state. It concerns the choice of tags and RFID reader to implement, or still the integration/the fixation of "tag systems" in a fibrous composite or in a textile. The encountered problems in this cluster are principally problems of choice of solution with respect to the technological constraints. • Accessing the material information : the aim is to develop a middleware allowing to exploit the physical layer (composed of tags and readers), to explore and optimize various architectures (tag-reader-middleware) for several applications. More specifically, it addresses issues like the accessibility of material, the tags order to read/write, . . . • Processing information : the previous clusters are responsible for design the communicating material as well as methods to access it. In this group, we wish to develop an embedded information system within the material, able to offer some services expected by users involved in the life cycle, from the design phase to the exploitation phase : for instance, the material could be used to traceability in a factory, to struggling against the brand piracy (custom officer), or to backing up data of final users. In order to implement this information system, it is necessary to define the conceptual specifications of the data model supporting the exchanged information throughout its life cycle, to study the mechanisms needed for allocating data on the communicating material and to specify the protocols and applicative services answering to actors needs. Furthermore, these specifications will be followed by a data model validation regarding to the reliance level offered by the information system (consistency of distributed data) and a validation of protocols and services by emulation on the one hand and on a software platform interfaced with the middleware developed earlier on the other hand.

• Ensuring security & completeness of the information :

like the cluster "Processing information", it studies the adaptability of data and focuses specifically on issues concerning information completeness. During its life cycle the material will be modified because of physical transformations (assembly/disassembly, . . . ), informational transformations (data access, data addition -deletion by several users) and of accessibility constraints (partial reading/writing of the material, . . . ). It raises issues like the choice of data to keep/modify within the material (at a given moment of its life cycle), or the reduction of the data volume via some methods of aggregation -desegregation.

It is also in this cluster that will be processed the data security (encryption, user identification, . . . ).

Finally, the contribution rate rt c of each cluster is computed (Fig. 5). We figure out that the cluster Processing information is therefore the heart of the project. Indeed, it has the greatest contribution rate with 58%, but it is also the only cluster in interaction with three others. Every cluster can be associated with a specific view of the system. Indeed, we evolve from a physical view of the system, intrinsically technological (Instrumented the material) towards an informational view of the system (Processing information).

CONCLUSION

In this study, we showed the relevance of this new paradigm. Thus, we have established an innovative approach in order to organize hierarchically the objectives to achieve and to decline them into various scientific and technological issues. This generic approach can be used in contexts other than the one exposed in this paper.

The different clusters require specific scientific skills. In order to implement the communicating material, it is necessary to identify the best scientific and industrial partners to answer problems highlighted by every clusters. As a result, it is important to implement a project structure between the partner teams in order to process these problems through mutual cooperation.

To validate the concept pertinence, we currently develop a specific demonstrator where the material is intrinsically communicating (i.e the material itself become communicating). Taking into account the technologies limits, we are building a prototype that integrate Hitachi µchip (dust of tags) spread into the material. As said previously, our perspectives is to give at the

  Fig. 4. Matrix QFD : Second level

Table 1 .

 1 Communicating product vs communicating material

	Characteristics	Communicating Product	Communicating Material
	Information location	localized	continuous
	Embedded memory	limited	extensible
	Resistance to the losses of tags	critical	less critical
	Services location	localized	continuous

  . Matrix QFD : First level

																						+	+
																						++
					++ ++				++					++	++ ++ ++	++ ++ ++ ++ ++ ++ + ++ ++	++ ++
			Problems	Area covered	Offered throughput	Frequency range	Type of reading	Power system	Tag memory	Insulate type	Controlled environment	Existing norm	Tag integration	Data allocating software	Data encapsulation	Tag density	Material memory	Adaptability to terminal devices	Access to useful info	Info to conserve/modify	Data integrity	Associated data model	Data classification	Data aggregation methods	Updates/propagation	Bit error correction code	Inter-tag redundancy	Areas accessible in read./writ.	Order of tags to read	Order of tags to write	Type of encryption key	Identification key
	Specifications C 1 ⇒ Q 2	P			Technological												Scientific
	Technological	Type of reader Reading/Writing short delays0.5 0.8090090090090030 090 090 Protocol of communication 0.6 090090 ...	090090 090030030030 010 010030090			090030030010 090030010090010010010000000000010010090090090 090010010 030	63 46 20

 20 0 50 50 50 1 20 50 1 50 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Type of reading 1 50 0 1 1 1 1 1 100 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Tag integration 1 50 1 0 50 50 50 50 50 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Power system 1 50 1 50 0 50 50 50 1 1 1 1 1 1 100 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Insulate type 1 1 1 50 50 0 50 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Area covered 1 20 1 50 50 50 0 100100 1 1 100 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Controlled environment 1 50 1 50 50 20 100 0 50 50 100100 1 100 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Tag density 1 1 100 50 1 1 100 50 0 1 1 50 100 1 50 1 20 1 20 1 20 1 1 1 1 1 1 1 1 Existing norm 1 50 1 1 1 1 1 50 1 0 1 50 1 1 1 50 1 1 100 50 50 1 1 1 1 1 1 1 1Order of tags to read 1 010 1 1 1 1 1 100 1 1 0 50 1 1 1 10 1 1 50 1 1 1 1 1 1 1 1 1 1 Areas accessible in read./writ. 1 1 1 1 1 1 100100 50 50 50 0 1 50 1 1 50 1 50 50 20 1 20 1 1 20 1 1 1 Data allocating software 1 1 1 1 1 1 1 1 100 1 1 1 0 1 100 50 1 1 50 1 1 1 1 1 1 1 1 1 1 Order of tags to write 1 1 1 1 1 1 1 100 1 1 1 50 1 0 1 10 1 100 50 1 1 1 1 1 1 1 1 1 1 Tag memory 1 1 1 1 100 1 1 1 50 1 1 1 100 1 0 100 1 50 1 1 100 50 1 1 1 1 1 1 1 Data encapsulation 1 1 1 1 1 1 1 1 1 50 10 1 50 10 100 0 1 1 20 1 100 1 1 1 1 1 1 1 1 Updates/propagation 1 1 1 1 1 1 1 1 20 1 1 50 1 1 1 1 0 1 50 1 1 1 10 1 1 1 1 1 1 Inter-tag redundancy 1 1 1 1 1 1 1 1 1 1 1 1 1 100 50 1 1 0 1 1 100 1 1 1 1 1 1 1 1 Associated data model 1 1 1 1 1 1 1 1 20 100 50 50 50 50 1 20 50 1 0 100100 1 20 100 50 20 1 1 20 Adaptability to terminal devices 1 1 1 1 1 1 1 1 1 50 1 50 1 1 1 1 1 1 100 0 1 1 1 1 1 1 1 1 20 Material memory 1 1 1 1 1 1 1 1 20 50 1 20 1 1 100100 1 100100 1 0 50 100 1 20 1 50 50 1 Bit error correction code 1 1 1 1 1 1 1 1 1 1 1 1 1 1 50 1 1 1 1 1 50 0 1 1 1 1 1 1 1 Info to conserve/modify 1 1 1 1 1 1 1 1 1 1 1 20 1 1 1 1 10 1 20 1 100 1 0 1 20 20 1 1 1 Data classification

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 50 1 1 1 1 1 1 0 1 Access to useful info 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 20 20 1 1 1 1 1 1 1 1 0 Fig. 5. Fiedler serialization matrix material itself functionalities (such as carrying a data model, capacities to react to external Event, . . . ).