D De Ligny

P Loiseau

D R Neuville
email: neuville@ipgp.jussieu.fr

Control of nucleation and crystal growth of a silicate apatitic phase in a glassy matrix

Nucleation and growth of crystal in an oxide glass was studied in a Si B Al Zr Nd Ca Na O system. The nucleation and growth process were monitored by thermal analysis and isothermal experiments. The effect of the network modifier was studied. Therefore for a Ca rich sample the crystallization is homogeneous in the bulk showing a slow increase of crystallinity as temperature increases. On the other hand, a Na rich sample undergoes several crystallization processes in the bulk or from the surface, leading to bigger crystals. The activation energy of the viscous flow and the glass transition are of same magnitude when that of crystallization is a lot smaller. Early diffusion of element is done with a mechanism different than the configurational rearrangements of the liquid sate. The global density and small size of the crystals within the Ca rich matrix confirmed that it would be a profitable waste form for minor actinides.

INTRODUCTION

The containment of long-lived separated radionuclides, minor actinides, would be profitably realized in a silicate apatite glassceramic obtained by controlled crystallization of parent glass. Indeed, the advantages of this waste form are numerous : glasses are easily handled and processed, both the glass and apatite can be used to immobilize actinides and the chemical fluctuations of the waste are absorbed in the glassy matrix [1]. Moreover apatite type compounds have shown their ability to incorporate Pu, minor actinide as Np, Am and Cm or lanthanide as Gd in a large range of composition. [START_REF] Zhao | Material Research Society Symposium Proceedings[END_REF] However, to be of real interest such an apatite glass-ceramic needs to be internally stress less and to present a good partition coefficient of the actinides between crystals and residual glass. These two last criteria meet if the formed crystals are in the bulk, highly concentrated and small in size. Therefore the ability to control the crystal nucleation and growth from the parent glass is crucial. In this work, trivalent minor actinides were simulated introducing Nd 2 O 3 in the glass composition. Intense apatite crystallization is obtained in the bulk for an alumino borosilicate glass containing up to 15 wt.% of Nd 2 O 3 . Electron microscopy, X-ray diffraction (XRD) and differential thermal analysis are used to characterize the devitrification.

The description of the nucleation and growth mechanism is done in term of activation energy. An important physical property is the atomic diffusion inside the glass matrix; it is often related with the atomic reorganization controlling the viscous flow. From viscosity measurement activation energy was determined and compared to the crystallization one. As to evaluate the effect of local chemistry, two glasses are studied, a glass containing Na as network modifier and a second glass containing a mixture of Na and Ca (see

DIFFERENTIAL THERMAL ANALYSIS

Two main observations can be done from thermal analysis. Grain size effect on the crystallization temperature characterizes a surface nucleation. The evolution of a temperature with the heating rate can always be interpreted in term of activation energy [START_REF] Moynihan | [END_REF]. As seen in Fig. 1, the two glasses undergo a twostage crystallization above the glass transition. For Ca.Na.gl., a first exothermal peak is observed followed by a large and diffuse exothermal event. For Na.gl. the two events are closer and more difficult to identify.

Grain size effect

In the case of Ca.Na.gl., modifications of the size of the sample did not modify the position of Tc. The more diffuse effect at higher temperature however showed various shapes not very reproducible. It is deduced that Ca.Na.gl. undergoes a crystallization in the bulk, the sluggish variations at higher temperature would need further studies to determine if they are due to the sample or thermal diffusion effects as the powder melt and move inside the crucible.

In the case of Na.gl. at the opposite, the secondary peak was constant with grain size, when the first peak, notified Tc, increases and overlaps the constant exothermic effect. This last variation shows that a surface crystal growth mechanism coexists with growth in the bulk for Na.gl.

Heating rate effect

On a constant size fraction, 80-125 mm, the samples were heated up at rate from 1°C/mn to 55°C/mn. The temperature of crystallization was taken at the maximum of the first exothermic peak. The glass transition temperature of Na.gl. was noticeable at all rates, corresponding with a strong endothermic effect. It was not possible to trace Ca.Na.gl. Tg which is not surprising since heat capacities of Ca containing liquids are a lot smaller than Na ones. The equation ( 1) was used as well for the temperature of crystallization and the glass transition. q h represents the heating rate.

d d T H R h x x q ln . * 1 ( ) = - D (1)
The energy of activation DH * c is 345±7 kJ/mol for Ca.Na.gl. and 273±19 kJ/mol for Na.gl. For the glass transition only for Na.gl., DH * g is 630±70 kJ/mol. This strong contrast between activation energy of crystallization and of the glass transition shows that atomic reorganization around the growing crystal is a lot easier than the change of configurational states freed inside the liquid.

Fig. 3. Backscattered electron images obtained on samples heated at the notified temperature for 2 hours. Pictures on the left column are from Ca.Na.gl., those on the right column from Na.gl.

SCANNING ELECTRON MICROSCOPE

Isochronal preparation of 2h at different temperatures was done on glass pieces of 2g. The sample without initial nucleation stage was brought directly at high temperature in a tubular furnace. A thermocouple near the sample showed that the final temperature was reached within 4 minutes. At the end of the two hours, the sample was quenched at room temperature, part of it ground for X-ray diffraction and another mounted and polished for S.E.M. observations. Since the crystals we wanted to observe are rich in Nd which has a strong Z compared to the residual glass, backcscattered electron observation gives the best contrast.

The pictures, obtained and reported on Fig. 3, show a very different behavior for the two glasses. At low temperature, Ca.Na.gl. shows very fine granular structure under the resolution of our microscope. As temperature raises, the granular structure got sharper indicating a better repartition between rich and poor Nd phase. It is only above 1030°C that the silico apatite crystals could be seen still at a sub micron size. At higher temperature, as densification goes on, fine elongated shapes appear. Na.gl. shows a more complex behavior. Isothermal heating at low temperature (680°C) leads to devitrification in the bulk of dense round crystallites where just 50°C higher, only small surface crystals could be observed. Large dendritic apatite crystals develops through out the sample from the surface at 830°C. At the highest temperature well crystallized apatite could be seen in the bulk. Depending on temperature, distinct growth mechanism took place in Na.gl. with nucleation on the bulk or on the surface. This variability of behavior of the Na rich glass seems to demonstrate the Na, as a strong network modifier, homogenizes a lot more the matrix and allows large cation mobility. All the X-ray patterns showed the presence of a single silicate apatite phase. However Na or Ca Nd-apatite has mostly the same patterns and cannot be identified in a quick study. For Na.gl. the crystals were big enough to allow microprobe analysis that showed the stochiometric composition of NaNd 9 (SiO 4 ) 6 O 2 . The microprobe analysis of the residual glass for the sample prepared at 1180°C showed that one third of the Nd was jailed inside the crystal phase. No such measures was possible on Ca.Na.gl. due to the too small size of the crystals. However previous study showed the presence of a Ca apatite.

VISCOSITY MEASUREMENTS

Viscosity measurement was realized at low temperature below the crystallization region. The effect of temperature on viscosity is related only to the structural relaxation and is not perturbed by nucleation. The measurements were made with a creep apparatus, whereby the rate of deformation of the sample is measured as a function of an applied constant stress at a fixed temperature.

The value are 649±19 kJ/mol for Res.gl. and 615±19 kJ/mol for Ca.Na.gl. The nature of the internal rearrangement is similar for the two glasses. However as it can be seen the Res.gl. has a smaller viscosity. It can be foreseen that, as the silico apatite crystal appeared, atomic mobility increased in the residual glass, which help the growth mechanism.

If this activation energy is compared to those find earlier we can notice that the Na.gl. glass transition energetic (629±75 kJ/mol) is similar leading that structural relaxation and viscous flow have similar origin. The activation energies deduced from crystallization however are much lower, so they are underlined by an atomic diffusion that is not limited by the viscosity.

CONCLUSIONS AND PERSPECTIVES

If nucleation and growth mechanisms appear to be simple in the case of the Ca rich glass, the Na glass showed several behaviors. Calcium holds a highly clustered environment helping crystal nucleation. This Ca function is in agreement in the higher crystallization activation energy for Ca.Na.gl. The ratio Ca/Na is critical to control nucleation process. The poorly crystallized Na.gl. leaded to a segregation of one third of the Nd inside the crystal phase, it can be expected a better result for Ca.Na.gl. Since shape of the crystals, small and dense, and chemical separation are promising, this system is confirmed to be a profitable waste form for minor actinides.

Fig. 1 .

 1 Fig. 1. Thermogram of the two glasses for a 10°C/mn. heating rate on a 20-40 m fraction. Tg and Tc show respectively the glass transition and the maximum of crystallization.

Fig. 2 .

 2 Fig. 2. Dependence of the temperature of crystallization Tc and of the glass transition temperature Tg versus heating rate.

Fig. 4 .

 4 Fig. 4. Dependence of the viscosity with temperature. Using the simple and well-know Arrhenius equation (2), the energy activation of the structural relaxation, DH h * can be estimated from the viscosity h. d d T

ATALANTE 2004 Nîmes (France) June 21-25, 2004 2

  TABLE I).

	500°C for 5 days to relief internal strain and
	allow cutting for sample preparation. Being well
	below the Tg, around 600°C, this thermal
	treatment is taken as ineffective on the
	nucleation process. No further thermal treatment
	was done before experiment.
	TABLE I. Samples Nominal Composition in
	mol %.
	Oxydes	Ca.Na.gl.	Na.gl.	Res.gl.
	SiO 2	53.7	53.7	54.2
	B 2 O 3	9.3	9.3	10.5
	Al 2 O 3	9.6	9.6	10.8
	Na 2 O	15.1	21.7	17.1
	CaO	6.6	0.0	5.1
	ZrO 2	2.0	2.0	2.3
	Nd 2 O 3	3.7	3.7	0.0
				EXPERIMENTAL
				Three glasses were prepared from reagent grade
				oxide (SiO 2 , Al 2 O 3 , ZrO 2 , Nd 2 O 3 ), carbonate
				(Na 2 CO 3 , CaCO 3 ) and boric acid (H 3 BO 3 )
				powders. All glasses were melted in platinum
				crucibles at 1400°C. The chemical homogeneity
				was achieved by regrinding at least once each
				sample after a first quench. Microprobe analysis
				and wet chemistry analysis confirm the nominal
				composition given in Table I. The Res.gl. sample
				is a residual end member of Ca.Na.gl. assuming
				that all the Nd was captured in a silicate apatite,
				CaNd 8 (SiO 4 ) 6 O 2 as observed in a previous
				work [3]. The Na.gl. was obtained by
				substituting in a molar basis all the Ca by Na
				from the Ca.Na.gl glass.
				Samples prepared for isochronal recrystallization
				and viscosity measurement were annealed at


B. BONNARD Toufik Bakir
email: toufik.bakir@u-bourgogne.fr

Bernard Bonnard
email: bernard.bonnard@u-bourgogne.fr

Jérémy Rouot
email: jeremy.rouot@epf.fr

A CASE STUDY OF OPTIMAL INPUT-OUTPUT SYSTEM WITH SAMPLED-DATA CONTROL: DING ET AL. FORCE AND FATIGUE MUSCULAR CONTROL MODEL

Keywords: Mathematics Subject Classification. 49K15, 93B07, 92B05 Muscular Force and Fatigue model, Optimization of sampled-data control, Non linear observer, Model predictive control

The objective of this article is to make the analysis of the muscular force response to optimize electrical pulses trains using Ding et al. force-fatigue model. A geometric analysis of the dynamics is provided and very preliminary results are presented in the frame of optimal control using a simplified inputoutput model. In parallel, to take into account the physical constraints of the problem, partial state observation and input restrictions, an optimized pulses train is computed with a model predictive control, where a nonlinear observer is used to estimate the state-variables.

(Communicated by the associate editor name) 1. Introduction. Functional Electrical Stimulation (FES) consists of applying an electrical stimulation to the muscle, in order to produce functional movements. It can be used for the muscular reinforcement, reeducation of the muscle and in the case of paralysis to activate the paralyzed muscles to produce movements. Mathematically FES leads to a sampled-data control problem which can be analyzed in this framework.

The simulations of muscular response to electrical stimulations are based on dynamics models. The origin comes from the Hill-Langmuir equation in the context of biochemistry and pharmacology, see [START_REF] Gesztelyi | The Hill equation and the origin of quantitative pharmacology[END_REF]. More recent models in the framework of model identification in non linear control are to due Bobet and Stein [START_REF] Bobet | A simple model of force generation by skeletal muscle during dynamic isometric contractions[END_REF] and Law and Shields [START_REF] Law | Mathematical models of human paralyzed muscle after long-term training[END_REF] and in this context a more sophisticated model was proposed by Ding et al. [START_REF] Ding | Two-step, predictive, isometric force model tested on data from human and rat muscles[END_REF][START_REF] Ding | Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains[END_REF][START_REF] Ding | A predictive model of fatigue in human skeletal muscles[END_REF][START_REF] Ding | Mathematical models for fatigue minimization during functional electrical stimulation[END_REF] where the force model is coupled to a fatigue model. This led to a set of five differential equations with a sampled-data control which can be used to stimulate and to control the force response to a pulses train of electrical stimulations and we shall refer to this model as the Ding et al. force-fatigue model in the sequel.

Only few works used this model to design an optimized train of pulses to control the force level [START_REF] Ding | Two-step, predictive, isometric force model tested on data from human and rat muscles[END_REF][START_REF] Ding | Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains[END_REF][START_REF] Ding | A predictive model of fatigue in human skeletal muscles[END_REF][START_REF] Ding | Mathematical models for fatigue minimization during functional electrical stimulation[END_REF]. Our aim being to make a more complete study of the problem in the framework of non linear optimal control, using sampled-data controls.

A geometric analysis of the model is provided and preliminary results in the frame of optimal control with sampled-data control [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF]. It is based on a simplified dynamics using the geometric properties of the force-fatigue model and control reduction to simplify the physical control constraints. The complete system is analyzed in details using a control predictive strategy (MPC), see [START_REF] Richalet | Algorithmic control of industrial processes[END_REF][START_REF] Wang | Model Predictive Control System Design and Implementation Using MATLAB[END_REF] coupled with a non linear observer technique based on [START_REF] Gauthier | A simple observer for non linear systems Application to bioreactors[END_REF] to estimate the state variable. This study completes [START_REF] Bakir | Predictive control based on nonlinear observer for muscu lar force and fatigue model[END_REF] and gives non trivial optimized pulses trains which can be implemented in practise up to specific devices to speed on line numerical computations.

The article is organized as follows. In Section 2, we make a brief presentation of Ding et al. force-fatigue model based on [START_REF] Wilson | Force response of locust skeletal muscle[END_REF]. In Section 3, the dynamics of the force model is briefly investigated to describe the input-output properties. In Section 4, the force-fatigue model is analyzed in the framework of geometric optimal sampleddata control systems and preliminary results are presented using a simplified model using a model reduction and an input transformation. Section 5 is devoted to the observer description. In Section 6, MPC method is presented using a further discretization of the dynamics to conclude by the algorithm used to compute in practise the optimized pulses trains. Numerical results are presented in the final Section 7.

2. Mathematical force-fatigue model. We refer to [START_REF] Wilson | Force response of locust skeletal muscle[END_REF] for a complete description and discussions of the model. The Ding et al. model studied in this article is presented next in the framework of model dynamics construction based on input-output observation and using the so-called tetania phenomenon in muscles responses. The first part of the model is the force (output) to pulses electrical stimulations (input). The pulses are normalized Dirac impulses at times 0 = t 1 < t 2 < . . . < t n :

v(t) = n i=1 δ(t -t i ) and I i = t i -t i-1
is the interpulse and convexifying leads to apply the input:

v(t) = n i=1 η i δ(t -t i ).
for some amplitudes η i ∈ [0, 1], i ∈ {1, . . . , n}. Such pulses train feeds a first order model to produce an output according to the dynamics

u(t) + u(t) τ c = n i=1 R i η i δ(t -t i ) (1)
where R i is a scaling factor associated to the phenomenon of tetania and corresponds to a accumulated effect of successive pulses and is modeled as

R i = 1 for i = 1 1 + (R 0 -1) exp -ti+1-ti τc for i > 1
where the magnitude is characterized by R 0 which is the limit case to high-frequency pulse. We denote E s (t) = u(t) the response to the pulses train and which is the (physical) control variable. It is of the form

E s (t) = 1 τ c n i=1 R i η i H(t -t i ) exp - t -t i τ c (2)
where

H(t -t i ) = 0 if t < t i 1 if t ≥ t i
is the Heaviside function.

The force response to such a train is modelled by the two equations of the socalled force model:

dC N dt + C N τ c = E s (t) (3)
which corresponds to a first order (resonant) linear dynamics which can be integrated with C N (0) = 0 and the force response F (t) is described by the equation

dF dt = A a -F b ( 4
)
where A is a fatigue parameter. Non-linear features of the model are described by the two mappings a, b:

a = C N K m + C N , b = 1 τ 1 + τ 2 a (5)
where K m , τ 1 are fatigue parameters and τ 2 an additional constant.

The complete model is obtained by describing the evolutions of the parameters associated to fatigue and corresponds to the linear dynamics

dA dt = - A -A rest τ f at + α A F (6) dK m dt = - K m -K m,rest τ f at + α Km F (7) dτ 1 dt = - τ 1 -τ 1,rest τ f at + α τ1 F. (8)
The full set of equations ( 3)-( 4)-( 6)-( 7)-( 8) are the force and fatigue model. We refer to Table 1 for the definitions and details of the symbols of the forcefatigue model.

For the purpose of our analysis the force fatigue model is written as the singleinput control system:

dx dt = F 0 (x) + u F 1 (x) (9)
with x = (x 1 , x 2 , x 3 , x 4 , x 5 ) = (C N , F, A, K m , τ 1 )
where u denotes the control u = E s (t) corresponding to the sampled physical control data:

(I 1 , I 2 , . . . , I n , η 1 , . . . , η n ) ( 10
)
with constraints I min ≤ I i ≤ I max , 0 ≤ η i ≤ 1 corresponding to interpulses bounds and amplitude convexification. This leads to a non linear model with sampled control data with prescribed convex control constraints (η, I) ∈ C; η = (η 1 , . . . , η n ), I = (I 1 , . . . , I n ).

3. The force model. The force model can be briefly investigated to describe preliminary results. Proof. Integrating (3), the equation ( 4) can be written as dF ds = c(s) -F (s) [START_REF] Gauthier | A simple observer for non linear systems Application to bioreactors[END_REF] with ds = b(t) dt, c = A a/b.

Thus this gives an explicit force response in the time parameter s

(I, η, s) → F (I, η, s). ( 12
)
Clearly we have Lemma 3.1. The above mapping is smooth with respect to I, η and piecewise smooth with respect to s.

Smoothing process.

For the sake of providing a smooth response for the observer it is sufficient to smooth the physical sampled control date as follows: use a bump function to smooth each Heaviside mapping H(t -t i ) at the sampling time t i .

3.3. Input-simplification. For the sake of the geometric analysis of the dynamics and to simplify the control constraints on (I, η) the FES signal E s (t) is taken as the input u(•) of the control system. Using (2) one can write

u(t) = E s (t) = e -t/τc v(t) with v(t) = 1 τ c n i=1 H(t -t i )R i η i e ti/τc . ( 13
)
Whence 0 = t 1 < t 2 < . . . < t n are given v(t) is a piecewise constant control depending upon the parameters η 1 , . . . , η n and the dynamics of the force model can be analyzed in the frame of geometric control. Fixing the interpulse led to a sample-data control model.

4.

Force-fatigue control model. First of all the control system (9) is analyzed in the framework of geometric control.

4.1.

The concepts of geometric control system. Consider a (smooth) control single-input control system.

dx dt = X(x) + u Y (x) y = h(x) (14)
where x ∈ R n , y = h(x) ∈ R corresponds to a (smooth) single-observation mapping.

The following concepts rely on seminal results of geometric control, see [START_REF] Sussmann | Controllability of nonlinear systems[END_REF][START_REF] Hermann | Nonlinear controllability and observability[END_REF].

We denote [U, V ] the Lie bracket of two (smooth) vector fields of R n :

[U, V ](x) = ∂U ∂x V (x) - ∂V ∂x U (x)
and a vector field U acts on (smooth) mappings f with the Lie derivative:

L U f = ∂f ∂x U (x).
Let D 1 = span{X, Y } and define recursively :

D k = span{D k ∪ [D, D k-1 ]}, k > 1.
Hence D k represents the Lie brackets of X, Y with lengths smaller than k + 1 and

denotes D L.A. = span ∪ k≥0 D k is the Lie algebra generated by {X, Y }. The system is called weakly controllable if for each x ∈ R n , D L.A. (x) = R n .
The observation space is the set of mappings:

Θ = {L G h; G ∈ D L.A. } and the system is called observable if for each x 1 , x 2 ∈ R n , x 1 = x 2 there exists g ∈ Θ such that g(x 1 ) = g(x 2 ).
Taking x 0 ∈ R n , a frame at x 0 is a set of elements X 1 , . . . , X n of D L.A. such that: X 1 , . . . , X n are linearly independent at x 0 and

n i=1 length(X i ) is minimal. The system is called feedback linearizable in U if ẋ = X(x) + u Y (x) is feedback equivalent to the linear system ẋ = Ax + ub, that is there exists a diffeomorphism ϕ on U , ϕ(0) = 0 and a feedback u = α(x) + β(x) v, β(x) = 0 such that g • (X, Y ) = (A, b) with g = (ϕ, α, β)
acting by change of coordinates and (affine) feedback.

Fix x(0) = x 0 and T > 0 and consider the extremity mapping:

E x0,T : u ∈ L ∞ ([0, T ]) → x(T, x 0 , u) where x(•) is the response of ẋ = X(x) + u Y (x) to u defined on [0, T ].
The control u(•) is called singular on [0, T ] if the extremity mapping E x0,T is not of maximal rank n when evaluated at u(•).

Geometric analysis of an observation system of the form ( 14) amounts to compute D L.A. , the observation space, the singular controls and the feedback equivalence properties. Achievements of geometric optimal controls amounts to synthesize the optimal control in relation with the Lie brackets properties of a frame.

4.2.

Optimal control a force-fatigue model with sampled control data.

Concepts

. By the physical nature of the problem, the control in the force/fatigue model ( 1) is fixed at specific times and remains constant during some period of time. This falls into the framework of sampled-optimal control problem that we recall briefly and we refer to [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF] for more details.

Consider a system ẋ = f (x, u). When the state x(•) and the control u(•) evolves continuously in time, we speak of a continuous-time optimal control problem and the control is said permanent. When the state x(•) evolves continuously in time whereas the control u(•) evolves in a discrete way, we speak of an optimal sampled-data control problem. Denoting by t f the final time, on [0, t f ] the control u takes its values in a discrete set. That means, for T > 0 fixed and k ∈ N, once the value of u(kT ) is chosen, then u(t) = u(kT ) for t ∈ [kT, (k + 1)T ]. T is called the sampling period and [kT, (k + 1)T ] is the sampling time interval where the control "is freezing". 4.2.2. Application to the force-fatigue model. We present the geometric methods presented in [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF] and numerical results on a reduced model in dimension 3 related to the force-fatigue model to validate the method.

Consider the minimization problem

t f 0 u(s) 2 + (F -F ref ) 2 ds → min
subject to the dynamic in x = (C N , F, A) defined by the equations

dC N dt = - C N τ c + u dF dt = a A -b F dA dt = - A -A rest τ f at + α A F (15)
and to the initial conditions

F (0) = A(0) = C N (0) = 0.
The variables t f , K m , τ 1 , τ 2 , τ f at , α A , A rest are constant and fixed to some ad-hoc values.

The system (15) can be written into the form ẋ = F0 (x) + u F1 (x).

Remark 1. The model ( 15) is a simplification of the force-fatigue model. In particular, without loss of generality, we don't take into account the factor exp(-t/τc) τc appearing in [START_REF] Bobet | A simple model of force generation by skeletal muscle during dynamic isometric contractions[END_REF]. Control constraints |u|≤ 1 are not the physical constraints, see Section 3.3. Also the fatigue dynamics is reduced to a single equation, motivated by the controllability properties of the system ( 6)-( 7)-( 8).

The permanent control case. The problem considered can be summarized as a permanent optimal control problem as follows

           min t f 0 u(s) 2 + (F -F ref ) 2 ds ẋ = F0 (x) + u F1 (x) u(t) ∈ [-1, 1] (C N (0), F (0), A(0)) = (0, 0, 0). ( 16
)
The pseudo-Hamiltonian of the system is

H(x, p, p 0 , u) = p • ( F0 (x) + u F1 (x)) + p 0 (u 2 + (F -F ref ) 2 )
where (p, p 0 ) : [0, t f ] → R 4 is the adjoint vector.

We denote by H i , i = 0, 1 the Hamiltonian lifts of the vector fields Fi , i = 0, 1. Normal case: p 0 = -1/2. Applying the Pontryagin maximum principle (PMP), the optimal (permanent) control is given in the normal case by

u(t) =    1 if H 1 (x(t), p(t)) ≥ 1, -1 if H 1 (x(t), p(t)) < 1, H 1 (x(t), p(t)) otherwise. ( 17
)
Abnormal case: p 0 = 0. Abnormal controls are characterized by the equation H 1 = 0, which differentiated twice by t leads to

p • [ F0 , F1 ] = 0, p • [[ F0 , F1 ], F0 ] + u [[ F0 , F1 ], F1 ] = 0.
Also, we have

D = det( F1 , [ F1 , F0 ], [[ F1 , F0 ], F0 ]) = -α A (a (C N ) A -b (C N ) F ) 2 D = det( F1 , [ F1 , F0 ], [[ F1 , F0 ], F1 ]) = 0 and D = 0 ⇔ A(C N (τ 1 + τ 2 ) + K m τ 1 ) 2 + F τ 2 (C N + K m ) 2 = 0.
Differentiating one more time led to:

Lemma 4.1. There are no admissible singular trajectories for the problem [START_REF] Law | Mathematical models of human paralyzed muscle after long-term training[END_REF].

The sampled-data control case. The corresponding optimal sampled-data control problem is given by

           min t f 0 u(kT ) 2 + (F -F ref ) 2 dt, with k = E(t/T ) ẋ = F0 (x) + u F1 (x) u(k T ) ∈ [-1, 1] (C N (0), F (0), A(0)) = (0, 0, 0). ( 18
)
where T > 0 is a fixed sampling period such that t f = j T for some j ∈ N. Normal case: p 0 = -1/2. Following [START_REF] Bourdin | Optimal sampled-data control, and generalizations on time scales[END_REF], the optimal (sampled-data) control is

u(kT ) ∈ arg max y∈[-1,1] 1 T (k+1)T kT H(s, x(s), p(s), p 0 , y) ds ( 19
)
for all controlling times kT, k ∈ {0, . . . , d}. Write

H1 = 1 T (k+1)T kT p 1 (s) ds
then, the optimal sampled-data control is

u(kT ) =    1 if H1 > 1, -1 if H1 < -1, H1 otherwise. ( 20
)
Numerical results. In Fig. 1, we represent numerical results for several values of t f /T . The optimal permanent control is represented with thin continuous line. A Gaussian quadrature rule is used to compute H1 using approximations of terms of the form p 1 (kT + T ), p 1 (kT + 1/2) . . . . We observe the convergence of the optimal sampled-data control to the permanent control as T tends to 0.

5.

Observer. The model used for this study (five state variables) is based on force measurements collected from a set of subjects. Thus, the accuracy of the calculated parameters is directly related to these persons. In this study, we suppose that initial conditions are different following the subject under study. That is why, we need to estimate some of these five initial conditions. Indeed, from a first analysis we can deduce that the rest values C N (t 1 ) = 0 is well known and don't need estimation. The remaining parameters are A rest , K m,rest and τ 1,rest . Concerning K m,rest , a sensibility study is realized in order to determine its effect on the force variation.

5.1. Sensibility study of the force versus K m : The force evolution is compared for K m,rest and different values K m,rest (±30% of error) for I = 10ms, 50ms and 100ms (see Fig. 2 for I = 10ms). In the case of I = 10ms, the maximum force error is of -0.3%. Following Interpulse value, the maximum force error is obtained for I = 100ms (-1.3%) which means that a tolerance of ±30% gives force evolution with good accuracy. 5.2. High-gain observer synthesis for the estimation of A and τ 1 . In this section, we design a modified version of the standard high-gain observer given in [START_REF] Gauthier | A simple observer for non linear systems Application to bioreactors[END_REF] taking into account the specific structure of the problem.

The system defined by the force equation ( 4) and the fatigue model ( 6) and ( 8 be rewritten as the single input-output system:

ẋ(t) = β m (t, E s (t))f (x(t), E s (t)) y(t) = h(x(t)) = F (t), (21)
with x = (F, A, τ 1 ) ∈ R 3 , y ∈ R, E s ∈ R, β = C N /(C N + K m ), 0 < β < 1,
and m is a positive integer. Note that in ( 21), K m is not a state variable thanks to the weak sensibility of the solution with respect to this variable (see sensibility study of the force versus K m ). We introduce the change of variables φ:

φ : R 3 → R 3 x → φ(x) = [h(x), L f (h(x)), L 2 f (h(x))]. (22)
We have:

∂φ ∂x =    1 0 0 -1 τ1+τ2β β F (τ1+τ2β) 2 ∂L 2 f h(x) ∂x1 ∂L 2 f h(x) ∂x2 ∂L 2 f h(x) ∂x3   .
where

∂L 2 f h(x) ∂x1 = βα A + 1 (x3+τ2β) 2 (1 - x3-τ1,rest τ f at + 2α τ1 x 1 ) ∂L 2 f h(x) ∂x2 = -β τ f at +x3+τ2β τ f at (x3+τ2β) ∂L 2 f h(x) ∂x3 = 1 (x3+τ2β) 2 (x 2 β -x1 τ f at ) - 2x1 (x3+τ2β) 3 (1 - x3-τ1,rest τ f at + α τ1 x 1 ) and det( ∂φ ∂x ) = β(x 2 βτ f at (x 3 + τ 2 β) -τ f at x 1 + 2x 1 (x 3 -τ 1,rest -α τ1 τ f at x 1 ) τ f at (x 3 + τ 2 β) 3 (23)
A sensibility study in ( 6) and ( 8) concerning respectively -A-Arest τ f at and -τ1-τ1,rest τ f at , shows that neglecting these two terms induces a maximum error of 7% for A and τ 1 , and 6% for the force value. Using the simplified model:

dA dt = α A F (24) dτ 1 dt = α τ1 F. (25)
with (4) gives:

det( ∂φ ∂x ) = β(x 2 β(x 3 + τ 2 β) -x 1 -2α τ1 x 2 1 ) (x 3 + τ 2 β) 3 . ( 26
)
Considering the input E S (t) ≥ 0, we have

∀t ∈ [0, T ] : (C N , x 1 ) ≥ 0, (x 2 , x 3 , K m ) > 0 ⇒ (x 3 + τ 2 β) 3 > 0 (see 26). det( ∂φ ∂x ) = 0 ⇒ β = 0 ⇒ C N = 0 (C N = 0 ⇒ x 1 = 0) (rest time) x 2 β(x 3 + τ 2 β) -x 1 -2α τ1 x 2 1 = 0 ⇒ x 1 = 1∓ √ 1+8ατ 1 x2β(x3+τ2β) 4ατ 1 (27)
α τ1 ∼ 10 -5 , 0 ≤ β ≤ 1 and (x 3 + τ 2 β) ∼ 10 2 . In the second case of (27), we have

8α τ1 x 2 β(x 3 + τ 2 β) << 1 ⇒ x 1 = 0 (C N = 0), x 1 = 1 2α τ1 .
Thus ( ∂φ ∂x (x(t)) = 0 for: • (C N , x 1 ) = 0 which correponds to the rest time (no control and then no need of state variable estimation).

• x 1 = 1 2ατ 1
, x 1 ∼ 10 5 which is outside of force maximum value.

In the simplified model ( 4),( 24) and (25), there is no singularity of ∂φ ∂x (x(t)) -1 during stimulation period (C N , x 1 > 0), and ill-conditioning of the matrix is reduced.

Based on the simplified model, the modified high-gain observer is defined as:

ẋ(t) =β m f 1 (x(t), E s (t)) -β m ( ∂φ ∂x (x(t)) -1 S -1 θ C T (C x(t) -y(t)). (28)
Recall that m is a positive integer. It depends on the used pulses frequency for the stimulation of the muscle. S θ is a symmetric positive definite matrix given by the following Lyapunov equation:

θS θ (t) + A T S θ (t) + S θ (t)A = C T C ( 29
)
where θ is a tuning parameter,

A =   0 1 0 0 0 1 0 0 0   ,C = 1 0 0 .
The terms of this matrix S θ = [S θ (l, k)] 1≤l,k≤3 have the following form:

S θ (l, k) = (-1) l+k l + k -2 k -1 θ -(l+k-1) , n k = n! (n -k)! k! .
5.3. High-gain observer proof.

   ẋ(t) = β m (t)f 1 (x(t), E s (t)) = f (x(t), E s (t)) y(t) = h(x(t)) (30) z 1 = h(x), ż1 = L f h(x) = β m L f1 h(x) = β m z 2 , ż2 = L f (L f1 h(x)) = β m L f1 (L f1 h(x)) = β m z 3 , . . . żn-1 = β m z n żn = β m L n f1 h(x) = ϕ n (u, z) ⇒      ż1 ż2 . . . żn      = β m       0 1 • • • 0 0 . . . . . . . . . . . . . . . . . . 1 0 • • • • • • 0            z 1 z 2 . . . z n      +      0 . . . 0 ϕ n (u, z)      (31) ⇒ ż = β m Az + ϕ(u, z) y = Cz. ( 32
)
Take C = 1 0 • • • 0 and consider the system (32). We make the following assumptions 5.4. Assumptions.

• H1) ϕ n (u, z) is globally Lipschitz with respect to z and uniformly with respect to u.

• H2) ∃ U ∈ R of admissible controls, a compact K ∈ R n and two positive constants β min and β max such that: for all u ∈ U and all y(t) associated to u and initial conditions z(0) ∈ K, we have

β min ≤ β(t) ≤ β max . • H3) β(t) is C 1 • H4) m β(t) β(t) < 1
The observer has the following expression:

ż(t) = β(t) m Aẑ(t) -β(t) m S -1 θ C T (C ẑ(t) -y(t)) ( 33
)
where S θ is the solution of the Lyapunov equation:

θS θ + A T S θ + S θ A -C T C = 0. ( 34
)
The solution of (34) can be rewritten as

S θ = 1 θ ∆ θ S 1 ∆ θ . ( 35
)
and

∆ θ = diag(1, 1 θ • • • 1 θ n ) (36)
with S 1 is the solution of (34

) for θ = 1. Let ρ = β m and e = ẑ -z ⇒ ė = ż -ż, then ė = ρAẑ + φ(u, ẑ) -S -1 θ C T (C ẑ(t) -y(t)) -ρAz -ϕ(u, z) = ρ(A -S -1 θ C T C)e + φ -ϕ (37) Write ē = ρ∆ θ e, then ė = ρ∆ θ e + ρ∆ θ ė = ρ∆ θ e + ρ∆ θ (ρ(A -S -1 θ C T C)e + φ -ϕ) = ρ∆ θ e + ρ∆ θ (ρ(A -θ∆ -1 θ S -1 1 ∆ -1 θ ∆ θ C T C∆ θ )e + φ -ϕ), ∆ θ C T C∆ θ = C T C = ρ ρ ē + ρθAē -θS -1 1 C T C ē + ρ∆ θ ( φ -ϕ), ∆ θ A∆ -1 θ = θA (38)
Consider the Lyapunov function:

V = ēT S 1 ē. Then V = 2ē T S 1 ė = 2ē T S 1 ( ρ ρ ē + ρθAē -θS -1 1 C T C ē + ρ∆ θ ( φ -ϕ)) (39) 2ē T S 1 Aē = ēT (S 1 A + A T S 1 )ē and S 1 A + A T S 1 = -S 1 + C T C.
Then, we have

V = -θρV + θρē T C T C ē -2ē T C T C ē + 2ē T S 1 ρ ρ ē + 2ē T S 1 ρ∆ θ ( φ -ϕ) = -θρV + (θρ -2) || C ē || 2 +2 ρ ρ V + 2ē T S 1 ρ∆ θ ( φ -ϕ) = -(θρ -2 ρ ρ )V + (θρ -2) || C ē || 2 +2ē T S 1 ρ∆ θ ( φ -ϕ) (40)
We deduce that

θρ -2 ρ ρ > 0 ⇒ θ > 2 ρ ρ 2 , θρ -2 < 0 ⇒ θ < 2 ρ ⇒ 2 ρ ρ 2 < θ < 2 ρ ⇒ ( ρ ρ = m β β ) < 1. (41)
Using (41) and assymption H1:

V ≤ -γV + 2ē T S 1 ρ∆ θ || φ -ϕ ||, γ > 0 ≤ -γV + 2ē T S 1 ρ θ n | φn -ϕ n | ≤ -γV + 2ē T S 1 ρ θ n K ē ≤ -(γ - β m θ n K 1 )V. (42)
For m and θ sufficiently large:

γ > β m θ n K 1 ⇒ V ≤ -γ 1 V, γ 1 > 0.
In the particular case of the force fatigue model, β(t) is piecewise smooth, the lack of regularity is numerically bypassed by the choice of the integer m. For example, for I = 10ms, m = 3 is sufficient to estimate the whole variables. However, for I = 25ms, m must be at least equal to 5 (see observer simulations in Section 7).

Model predictive control (MPC).

MPC computes a sequence of future controls to optimize "future" plant behavior (see Fig. 4). This computation is based on the dynamic model of the system, and must respect a certain cost and associated constraints (optimal control). In this framework, Model Predictive Heuristic Control (MPHC) was introduced by Richalet et al. [START_REF] Richalet | Algorithmic control of industrial processes[END_REF] using impulse response type dynamic model. Dynamic Matrix Control (DMC) followed in 1980 (Cutler et al. [START_REF] Cutler | Dynamic matrix control: A computer control algorithm[END_REF]) using step response type dynamic model. State space formulation of MPC was introduced by Li et al. [START_REF] Li | A state space formulation for model predictive control[END_REF].

In the state space framework, the system takes the form:

   ẋ = f (x, u) x(0) = x 0 t ∈ [0, +∞[ . (43)
Consider the cost function to be minimized

J ∞ (u(.), x 0 ) = ∞ 0 q(x u (τ, x 0 ), u(τ ))dτ (44)
with

x u (t, x 0 ) = x 0 + t 0 f (x u (τ, x 0 ), u(τ ))dτ (45)
and

   q(0, 0) = 0, q ∈ C 2 q(x, u) ≥ c q (||x|| 2 +||u|| 2 ), c q > 0 u → q(x, u) is convex for all x. (46)
Without constraints, Bellman's principle of optimality (1950) gives the solution. In the case of constrained problem, instead of J ∞ , we minimize the following cost:

J(u(.), x 0 , T ) = T 0 q(x u (τ, x 0 ), u(τ ))dτ + v(x u (T, x 0 ) (47) with J * (x 0 , T ) = min u(.) J(u(.), x 0 , T ) (48)
and u * (t, x 0 , T ) = arg min u(.)

J(u(.), x 0 , T ).

This optimization over a finite horizon follows the algorithm: 1. Solve min J(u(.), x 0 , T ) and find u * (., x 0 , T ). 2. Apply u * (., x 0 , T ) for τ ∈ [0, T s [, 0 ≤ T s ≤ T , (T s : time sampling period). 3. Repeat using x(T s ) instead of x 0 .

6.0.1. Discrete linear system (basic case). To explain the method, we consider the discrete linear system:

x(k + 1) = Ax(k) + Bu(k) y(k) = Cx(k) (50)
where x(t) ∈ R n , u(t) ∈ R m , y(t) ∈ R p are respectively the state, input and output vectors, and t = kT s . This system (64) results from a discrete system modeling or continuous-discrete transformation. Practical considerations of implementation make this form more appropriate in the framework of MPC. We suppose that the system is both controllable and observable, and we call x(k + i/k) with i ≥ 0, the future state vector x(k + i) starting from t = kT s .

x(k + i/k) is assumed to be available through using linear state variables estimator. In the case of unconstrained problem, a quadratic cost is minimized in order to obtain optimal control sequence. For instance, choose the following cost:

J(k) = ∞ j=k ((y(j + 1/k) -y ref ) T Q(y(j + 1/k) -y ref ) + (u(j/k) -u ref ) T R(u(j/k) -u ref ) + ∆u(j/k) T S∆u(j/k)) (51) with Q = Q T 0, R =
x(k + i/k) = A i x(k/k) + i-1 j=0 A i-j-1 Bu(k + j/k). (52)
The input-output relation along the receding horizon is then:

     y(k + 1/k) y(k + 2/k) . . . y(k + N r /k)      =      CA CA 2 . . . CA Nr      x(k/k) +       CB 0 • • • 0 CAB CB . . . . . . . . . . . . . . . 0 CA Nr-1 B • • • CAB CB            u(k/k) u(k + 1/k) . . . u(k + N r -1/k)      . (53)
Write (53) as ȳ = Ψx(k/k) + Γū (54) where ȳk ∈ R pNr , ūk ∈ R mNr , Ψ ∈ R pNr×n , Γ ∈ R pNr×mNr .

Most frequently used constraints concern the minimum and maximum input, state and output variables:

• The control variables:

u min ≤ u(k) ≤ u max . (55)
• The change rate of the control variables:

∆u min ≤ ∆u(k) ≤ ∆u max (56) with ∆u(k) = u(k) -u(k -1).
(57) • Soft output variables constraints (relaxed output constraints using large slack variable s v to avoid constraints conflicts when solving control problem):

y min -s v1 ≤ y(k) ≤ y max + s v1 . (58)
• Soft state variables constraints (for the same raison as output constraints):

x min -s v2 ≤ x(k) ≤ x max + s v2 (59)
In the case of constrained problem, the cost function (51) becomes:

J(k) = k+Nr-1 j=k ((y(j + 1/k) -y ref ) T Q(y(j + 1/k) -y ref ) + (u(j/k) -u ref ) T R(u(j/k) -u ref ) + ∆u(j/k) T S∆u(j/k)) (60)
subject to (54) and a set of inequality constraints among (55) -(59). (left half plane) E s and force profile for applied amplitude and interpulse stimulation, (right half plane) Predicted E s and force using single move strategy to be optimized.

Compute

C N (t k ), F (t k ), Â(t k ), τ1 (t k ), K m (t k ) 3. Data F = N r * u k (1)
step int with: u i (1) = u k (1) for i = k, k + 1, ..., k + (N r -1), and step int is a submultiple of u i (1) = I i . 4.

F mean = 1 Data F Data F j=1 F (t k + j Step int , E s ) with: E s = 1 τ c k+(Nr-1) i=1 u i (2)H(t Nr -t i ) R i exp(- t Nr -t i τ c ) . (65)
We recall u i (1) = u k (1) and u i (2) = u k (2) for i = k, k + 1, ..., k + (N r -1), (single move profile) and t i = t i-1 + u i-1 (1), i ≥ 2, t Nr = k+(Nr-1) i=1

I i . 5. u * (k/k) = I * k η * k = arg min I k ∈[I min ,I max ] η k ∈[η min ,η max ] J 6. if t k+1 ≥ Final t ⇒ stop, else, k = k + 1, back to 2.
7. Numerical results.

7.1. High-gain observer. To exhibit the interest of the use of the power m in the nonlinear observer, let's consider MPC based nonlinear observer using only amplitude as control variable. We consider also the worst case of +30% of error of K m .

The following simulation results follow the stimulation protocol used in practise (a set of periods with stimulation and rest time slots in each period). During the stimulation time slot, the control is calculated to bring the force to F ref and in the rest time slot, the stimulation amplitude is set to 0. In this section, only two stimulation periods are considered.

In the cases of interpulse I = 10ms (100Hz) and I = 25ms (40Hz), the power values are m = 3 and m = 5, respectively. Fig. 6 and Fig. 7 represent the A estimates for I = 10ms and I = 25ms, respectively. Fig. 8 and Fig. 9 are the τ 1 estimates for I = 10ms and I = 25ms. Â converges after 50ms when I = 10ms and 100ms when I = 25ms. Concerning τ1 , it converges after 75ms when I = 10ms and 200ms when I = 25ms. Large I seems to delay the convergence of the observer. that maximum value of amplitude control (see Fig. 15) and computed interpulse (see Fig. 14) cannot maintain F at F ref . Maximum interpulse frequency is not used at the beginning of the forth and the fifth periods. Higher value of N r could correct this problem (supposing that global minimizer of MPC algorithm is reached at each iteration). However, increasing N r will render MPC algorithm time consuming which causes a problem for real time implementation.

Figure 1 .

 1 Figure 1. Time evolution of the permanent control (thin continuous line) and sampled-data control for several values of the sampling period T ∈ {t f /20, t f /40, t f /200}.

Figure 2 .

 2 Figure 2. Evolution of K m for different initial conditions (case of I = 10ms).

Figure 3 .

 3 Figure 3. Relative error of the force for a well known and erroneous K m initial condition (case of I = 10ms).

  R T 0 and S = S T 0. The couple (y ref , u ref ) corresponds to a set point and ∆u(k) = u(k) -u(k -1). The solution of this problem is obtained using the LQ controller. In the case of constrained problem,

Figure 4 .

 4 Figure 4. General MPC strategy diagram

Figure 5 .

 5 Figure 5. (left half plane) E s and force profile for applied amplitude and interpulse stimulation, (right half plane) Predicted E s and force using single move strategy to be optimized.

Figure 7 .Figure 8 .Figure 9 .Figure 10 .

 78910 Figure 7. Evolution of A and  for I = 25, 30% error of K m .

Figure 11 .Figure 12 . 8 .AmplitudeFigure 13 .

 1112813 Figure 11. Evolution of the force for a reference force of 425N and different receding horizons(3, 5 and 10)

Figure 14 .

 14 Figure 14. Evolution of the interpulse (control) for a reference force of 425N and a preditive horizon of 10.

Figure 15 .

 15 Figure 15. Evolution of the amplitude (control) for a reference force of 425N and a preditive horizon of 3.

Table 1 .

 1 Margin settings

	Symbol Unit	Value	description
	C N	-	-	Normalized amount of
				Ca 2+ -troponin complex
	F	N	-	Force generated by muscle
	t i	ms	-	Time of the i th pulse
	n	-	-	Total number of
				the pulses before time t
	i	-	-	Stimulation pulse index
	τ c	ms	20	Time constant that commands
				the rise and the decay of C N
	R 0	-	1.143	Term of the enhancement
				in C N from successive stimuli
	A	N ms	-	Scaling factor for the force and
				the shortening velocity
				of muscle
	τ 1	ms	-	Force decline time constant
				when strongly bound
				cross-bridges absent
	τ 2	ms	124.4	Force decline time constant
				due to friction between actin
				and myosin
	K m	-	-	Sensitivity of strongly bound
				cross-bridges to C N
	A rest	N ms	3.009	Value of the parameter A
				when muscle is not fatigued
	K m,rest	-	0.103	Value of the parameter K m
				when muscle is not fatigued
	τ 1,rest	ms	50.95	The value of the parameter τ 1
				when muscle is not fatigued
	α A	1 ms 2	-4.0 10 -7 Coefficient for the force-model
				parameter A in the fatigue
				model
	α Km	1 msN	1.9 10 -8 Coefficient for the force-model
				parameter K m in the fatigue
				model
	α τ1	1 N	2.1 10 -5 Coefficient for force-model
				parameter τ 1 in the fatigue
				model
	τ f at	s	127	Time constant controlling the
				recovery of (A, K m , τ 1 )
	3.1. Parameterization. Observe that in the force model the static non-linearities
	described by the mappings a and b can be absorbed by time reparameterization to
	provide an explicit form for the force responses.

supported by NSF grant xx-xxxx.

This minimization problem could be expressed as a quadratic cost function subject to linear equality and inequality constraints:

and ȳ = C x with X = [x, ū] T . Clearly, optimization problem (61) is a convex Quadratic Programming (QP) problem; X * = [x * , ū * ] T is than a global minimizer at each iteration. ū * being calculated, only u * (k/k) is applied at time t = kT s and optimization algorithm is repeated at t = (k+1)T s . In the case of NMPC (Nonlinear Model Predictive Controller), the dynamic is defined using a discretization of the general form (43). Using the same above cost function, optimization problem could be defined as follow: [START_REF] Wang | Model Predictive Control System Design and Implementation Using MATLAB[END_REF], Fletcher [START_REF] Fletcher | Practical Methods of Optimization[END_REF] and Boyd and Vandenberghe [START_REF] Boyd | Convex Optimization[END_REF]).

In the case of force or force-fatigue model, we consider the more general case of a varying T s (T s → I(i)), then u(i) = [I(i) η(i)] T , where I(i) = t(i) -t(i -1) is the interpulse between two successive pulses and η(i) is the pulse amplitude applied at time t i . In this case, ū(k) = [u(k/k) . . . u((k + N r -1)/k)] is the discrete receding horizon control vector (of length N r ) to be calculated by the NMPC at time t corresponding to iteration k. Then:

t(N r ) being the final time of the optimization horizon which is unknown a priori. For instance, using single move strategy u(k/k) = u(k + 1/k) = . . . = u((k + N r -1)/k), (see Fig. 5).

The force-fatigue model being nonlinear, an NMPC is used to solve the problem with Nonlinear Programming Algorithm (NPA). In this case, criterion to be minimized is:

subject to: 0 ≤ η(i) ≤ 1 and 0.01ms ≤ I(i) ≤ 0.1.

The estimation of the state variables vector (by the high-gain observer) is used as an initial variables vector to perform the NMPC over the horizon N r : ALGORITHM 1. Give F inal t , k = 1 of stimulations is dedicated to the control of the force level using both stimulation interpulse and amplitude controls. These simulations show the effect of the receding horizon on the control efficiency. Reasonable value of N r is however suitable to guarantee a short computation time for real time implementation.


Thierry Woignier
email: thierry.woignier@imbe.fr

Florence Clostre

Luc Rangon

Magalie Lesueur Jannoyer

Thierry

Natural Nano Materials Trapp Pollutants in Soils

Keywords: natural nanomaterials, pesticide, pollution, chlordecone, allophanic clay, andosols

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

HE fate of pesticides is an important environmental issue because large quantities of pesticides are applied to agricultural soils to control pest and insects. Organochlorine pesticides are persistent organic pollutants (POP) [START_REF] Madaj | Lindane, kepone and pentachlorobenzene: chloropesticides banned by Stockholm convention[END_REF], which are not only toxic but also are biomagnified through the food chain [START_REF] Dallaire | Cognitive, visual, and motor development of 7-month-old Guadeloupean infants exposed to chlordecone[END_REF]. Even if the use of such pesticides has been forbidden for decades, soils still are a reservoir of pollution and continue to contaminate rivers, groundwater, ecosystems and finally human beings [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF].

A wide range of pesticides has been used on tropical volcanic soils and different authors showed that andosol displayed higher sorption capacity than other soil types [START_REF] Paton | Distribution, bioavailability and behavior of persistent organic pollutants in Andosols: with specific reference to Iceland[END_REF]- [START_REF] Fernandez Bayo | Sorption characteristics of chlordecone and cadusafos in tropical agricultural soils[END_REF]. Paton et al. [START_REF] Paton | Distribution, bioavailability and behavior of persistent organic pollutants in Andosols: with specific reference to Iceland[END_REF], found that dichlorophenol and pentachlorophenol were tightly bound to the organic fractions within the andosols, which could have binding affinity several orders of magnitude greater than other soils. In Mexico, for two organophosphorous insecticides, parathion and cadusafos, the amounts of bound residues were greater in the andosols compared with the vertisols [START_REF] Olvera-Velona | Sorption and desorption of organophosphate pesticides, parathion and cadusafos, on tropical agricultural soils[END_REF]. Prado et al. [6] also demonstrated that andosol showed a higher atrazine sorption capacity than vertisol with Kd values roughly one order of magnitude higher.

In French West Indies, sorption of cadusafos was also higher in andosol than in nitisol [START_REF] Fernandez Bayo | Sorption characteristics of chlordecone and cadusafos in tropical agricultural soils[END_REF]. Similarly, for the organochlorine insecticide chlordecone, andosol showed higher contents than the vertisols and nitisols [START_REF] Brunet | Determination of soil content in chlordecone [organochlorine pesticide] using near infrared reflectance spectroscopy [NIRS[END_REF][START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF]. Chlordecone (CLD) was used years ago in order to control pests [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF][START_REF] Brunet | Determination of soil content in chlordecone [organochlorine pesticide] using near infrared reflectance spectroscopy [NIRS[END_REF]. This neurotoxic and carcinogenic molecule [START_REF] Dallaire | Cognitive, visual, and motor development of 7-month-old Guadeloupean infants exposed to chlordecone[END_REF] still persists in the soil and contaminates crops, water resources and ecosystems thus leading to human exposition through food. Andosols are generally more polluted than the other kinds of tropical soils but recent data have shown that andosols release less CLD to percolating water and crops than other volcanic soils [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF], [START_REF] Woignier | Chlordecone retention in the fractal structure of volcanic clay[END_REF]. One explanation of these retention effect proposed in the literature is the high organic content of andosols because of the high affinity between CLD and soil organic carbon (SOC). However, andosols contain clay (allophane) that presents unique structures and physical properties compared to crystalline clays.

In the case of parathion, cadusafos and 2,4dicchlorophenoxyacetic acid [START_REF] Olvera-Velona | Sorption and desorption of organophosphate pesticides, parathion and cadusafos, on tropical agricultural soils[END_REF], [8], the presence of allophane clay (the constitutive clay of andosols), seems to increase pesticide sorption. Allophane is also known to be an adsorbent for polluting chemical species [9] including phosphates [10], [11], arsenates [12], [13], sulphates, molybdates, chromates and selentes [START_REF] Dallaire | Cognitive, visual, and motor development of 7-month-old Guadeloupean infants exposed to chlordecone[END_REF]. The AlOH and SiOH groups on the allophane surface have interesting properties with respect to boron adsorption [9], [10]. Allophane has an adsorption ability for benzene derivatives [START_REF] Nishikiori | Adsorption of benzene derivatives on allophane[END_REF] but also for the removal of fluoride from water [START_REF] Reinert | Duclaux Characterization and boron adsorption of hydrothermally synthesised allophanes[END_REF]. Recently, allophane nanoparticles were chemically synthesized as an adsorbent of organic pollutants [14]- [16] and synthetic allophane were also studied for boron removal from water [START_REF] Kaufhold | Allophane compared with other sorbent minerals for the removal of fluoride from water with particular focus on a mineable Ecuadorian allophane[END_REF]. Allophane has been found to be effective in the stabilization of C in soil and composts, offering perspectives for carbon sequestration technology [START_REF] Calabi-Floody | Nanoclays from an Andisol: Extraction, properties and carbon stabilization[END_REF].

Allophanes are naturally occurring nanoparticles formed through the alteration of volcanic parent materials. Allophane is an amorphous alumino silicate, the unit cell of which appears as spheroids with diameter close to 5 nm [START_REF] Wadak | The distinctive properties of Andosol[END_REF] forming aggregates with a fractal structure [START_REF] Adachi | Application of a scaling law to the analysis of allophane aggregates[END_REF]. Some researchers described allophane aggregates as "gel like" [START_REF] Wallace | Structural behaviour of residual soils of the continually wet Highlands of Papua New Guinea[END_REF]. The clay in its natural state consists of aggregations of spherical allophane Thierry Woignier a,b , Florence Clostre c , Luc Rangon a,b , Magalie Lesueur Jannoyer particles [START_REF] Wadak | The distinctive properties of Andosol[END_REF], [START_REF] Lindner | Hollow nanospheres, allophanes 'All-organic' synthesis and characterization[END_REF]. Allophane spherules tend to form porous nano size aggregates [START_REF] Calabi-Floody | Nanoclays from an Andisol: Extraction, properties and carbon stabilization[END_REF], [START_REF] Garrido-Ramirez | Catalytic wet peroxide oxidation of phenol over iron or copper oxide-supported allophane clay materials: Influence of catalyst SiO2/Al2O3 ratio[END_REF], the average diameter of the aggregates is around 100 nm and it exists a clear analogy between allophane aggregates and synthetic silica gels [START_REF] Clark | Cation and anion retention by natural and synthetic allophane and imogolite[END_REF], [START_REF] Wadak | The distinctive properties of Andosol[END_REF], [START_REF] Adachi | Application of a scaling law to the analysis of allophane aggregates[END_REF]. Transmission electron microscopy image shows the occurrence of allophane nano spheres with a diameter of about 5 nm, the allophanic particles aggregate and form clusters [START_REF] Wadak | The distinctive properties of Andosol[END_REF], [START_REF] Calabi-Floody | Nanoclays from an Andisol: Extraction, properties and carbon stabilization[END_REF]. Small Angle X-rays Scattering (SAXS) results shows that the allophane aggregates are fractal clusters with fractal dimension close to 2.5-2.7 [START_REF] Chevallier | Fractal structure in natural gels: effect on carbon sequestration in volcanic soils[END_REF]. The tortuous structure of the allophane aggregates is typically in the mesopore range 5-100 nm.

Our study hypothesizes that, besides the well-known effect of soil organic content on CLD retention, the features of this natural gel "allophane" favor pollutants accumulation in soils. The purpose of this work is to: 1) describe some physical properties of the allophane clay, 2) put forth the importance of the natural gel properties of allophane clay (microstructure, large water content, pore collapse during drying) for pollutant trapping, 3) discuss the synergy "soil organic contentallophane microstructure" on the pesticides retention in tropical volcanic soils.

II. EXPERIMENTAL

A. Soil samples

We sampled several andosols (allophane soils) and for comparison nitisols (containing halloysite clays) in Martinique (14°40 N, 61°00 W) in the vicinity of the "Montagne Pelée" volcano. These different soils are representative of the main polluted soils found in French West Indies [START_REF] Dallaire | Cognitive, visual, and motor development of 7-month-old Guadeloupean infants exposed to chlordecone[END_REF][START_REF] Brunet | Determination of soil content in chlordecone [organochlorine pesticide] using near infrared reflectance spectroscopy [NIRS[END_REF]. The soils we selected were known to be historically polluted by CLD. The soil samples were retained in closed containers to avoid evaporation.

Organic carbon contents were measured by dry combustion using a CHN chromatograph analyzer (Thermo Finnigan Flash EA 1112, thermo Finnigan Italia, Rodano, Italy).

The allophanic content was measured by the method of Mizota and van Reewijk [START_REF] Mizota | Clay mineralogy and chemistry of soils formed in volcanic material in diverse climatic regions[END_REF] using Al and Si content extracted by oxalate and pyrophosphate and measured by ICP-AES (Inductively Coupled Plasma Atomic Emission Spectroscopy ICAP 6500, Thermo Scientific, France). The soils samples covered the allophane range 0 to 26 % (w/w basis).

The shrinkage curve during drying was measured according to the following procedure. The water saturated soil samples were placed on a balance to measure loss of water during drying. In parallel, the linear shrinkage in the height and diameter of each sample was measured by three spot lasers. The data were recorded and converted to specific volume and water content.

The crystalline and amorphous structure was studied by Xrays diffraction (Cu Kα) with a Philips PW 1830 and the presence of allophane was characterized by Infrared Spectroscopy. The apparatus used was an IR-FT Nicolet 510P spectrometer and the samples were diluted in KBr pellets with a 5 10 -3 mass ratio.

Physical fractionation by sieving in water was done after soil dispersion in water and hexametaphosphate. We isolated by size three different organo mineral fractions: < 50 µm size, 50-200 µm size and > 200 µm size. Three replicates of each fraction were prepared.

B. Transfer Experiments

Soil to plant transfer coefficients were determined on andosol and nitisol for two crop models: sweet potato [START_REF] Cabidoche | Contamination of Harvested Organs in Root Crops Grown on Chlordecone-Polluted Soils[END_REF] and carrot. Tubers and corresponding soil were sampled and the CLD content analyzed. Six replicates were prepared on nitisol and 15 on andosol. Soil to plant transfer coefficient is the mean value of the ratio of the pesticide concentration in plant to the pesticide concentration in soil expressed in mg.kg-1 of fresh matter / mg.kg-1 of dry soil. For the experiments the soil contamination ranges from 0.7 to 7.2 mg/kg of dry soil in andosols and from 0.1 to 2.7 mg/kg of dry soil in nitisols. CLD concentration in soils and crops samples were analyzed by the LDA26 (Valence, France) as described in [START_REF] Cabidoche | Contamination of Harvested Organs in Root Crops Grown on Chlordecone-Polluted Soils[END_REF]. All statistical analyses were performed using XL STAT 2012.6.08 (Addinsoft 1995-2012). ANOVA and Tukey Honestly Significant Difference tests (p<0.05) were used to compare means of transfer for the different soil types and to assess the difference of CLD content in the different soil fractions for each soil type. The linear regressions and ANOVA were used to assess the equation of the linear models and the significance of the explanatory variable.

C. Porous properties

Porous and structural measurements generally require dried solid samples. A preliminary study [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF] has shown that the CO 2 supercritical drying is necessary to preserve the mesoporosity of soils. A critical point dryer (CPD 010, Balzers, Liechtenstein) was used in a procedure such as previously described [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF]. The pore size distribution of the sample was calculated from the N2 adsorption-desorption curves (Micromeritics ASAP 2010) BJH (Barret-Joyner-Halenda) analysis [START_REF] Brunauer | Adsorption of gases in multimolecular layers[END_REF]. The soil powder was placed in a test tube to degas for 24 hours at 50°C in a vacuum of 2-4 μm Hg. The scanning electron micrographs were obtained with a Cambridge stereo scan 360 scanning electron microscope.

III. RESULTS AND DISCUSSION

A. Soils characterization

Allophane are normally referred to as X-rays amorphous materials but for some authors it could give weak and broad Xrays traces at a lattice spacing d close to 3.3 and 2.25 Å [START_REF] Kaufhold | Allophane compared with other sorbent minerals for the removal of fluoride from water with particular focus on a mineable Ecuadorian allophane[END_REF], [START_REF] Lindner | Hollow nanospheres, allophanes 'All-organic' synthesis and characterization[END_REF], [START_REF] Garrido-Ramirez | Catalytic wet peroxide oxidation of phenol over iron or copper oxide-supported allophane clay materials: Influence of catalyst SiO2/Al2O3 ratio[END_REF]. Figure 1 shows a typical X-rays spectrum that we obtained on the studied andosols and nitisols, the general trend of the andosol spectrum is a clear evidence of bulges of the base line between 2 and 7 Å, attributed to noncrystalline solids, allophane and/or imogolite [START_REF] Barois | Andosol-forming process linked with soil fauna under the perennial grass Mulhembergia macroura[END_REF]. The Xray diffraction spectrum for nitisols showed peaks at 4.45 and 7.4 Å that were attributed to the clay halloysite [START_REF] Dubroeucq | Weathering and soil forming processes under semi-arid conditions in two Mexican volcanic ash soils[END_REF], [START_REF] Onodera | materials prepared from allophane soil: I. Preparation and characterization of silver/phosphorus-silver loaded allophanic specimens[END_REF]. Volcanic soils contain different types of poorly crystalline aluminosilicate minerals like imogolite (Al/Si = 2), proto imogolite and allophane. Imogolite forms 2 nm hollow tubes, while allophane forms spheres [START_REF] Nanzyo | Volcanic Ash Soils. Developments in Soil Science[END_REF]. Different allophane structures have been identified based on their chemical composition: Al-rich allophane (Al/Si ~ 2) also called "proto imogolite allophane" since the local structure is similar to that of imogolite [START_REF] Parfit | [END_REF]. Allophane with a lower Al/Si ratio is called Si-rich allophane [START_REF] Brigatti | Developments in Clay Science[END_REF][START_REF] Montargès-Pelletier | [END_REF]. Transmission electron microscopy [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF] showed no imogolite structure features, such as hollow tubes of 2 nm in diameter. The X-Ray's data showed no diffraction peaks, 8 Å and 12 Å, related to imogolite or to proto imogolite respectively [START_REF] Cradwick | Imogolite, a Hydrated Aluminium Silicate of Tubular Structure[END_REF], [START_REF] Levard | Structure and distribution of allophanes, imogolite and proto-imogolite in volcanic soils[END_REF].

A method to reveal the occurrence of allophane and the absence of imogolite is the infrared spectroscopy (FTIR). In the literature [START_REF] Barois | Andosol-forming process linked with soil fauna under the perennial grass Mulhembergia macroura[END_REF], [START_REF] Gustafsson | Podzolisation mechanisms and the synthesis of imogolite in northern Scandinavia[END_REF], [START_REF] Denaix | Structure and affinity towards Cd 2+ , Cu 2+ Pb 2+ of synthetic colloidal amorphous aluminosilicates and their precursor Colloids and Surfaces A[END_REF], it has been shown that characteristic bands like 506, 577 and a shoulder at 970 cm -1 are the signature of allophane. These bands are also accompanied by a broad band in the silanol stretching region (around 3500-3700 cm -1 ). Figure 2 shows the typical FTIR spectrum recorded for the studied andosols. The spectrum reveals the shoulder at 970 cm -1 (Si-O stretching), the bands in the range 500-600 cm -1 and broad band in the silanol region, characteristics of allophane. Moreover the IR spectra of the different studied andosols do not present doublets of the band at 577 and 967 cm -1 characteristic of the imogolite structure [START_REF] Levard | Structure and distribution of allophanes, imogolite and proto-imogolite in volcanic soils[END_REF]- [START_REF] Denaix | Structure and affinity towards Cd 2+ , Cu 2+ Pb 2+ of synthetic colloidal amorphous aluminosilicates and their precursor Colloids and Surfaces A[END_REF]. Thus, we can deduce that the major part of the poorly crystalline content measured by the method by oxalate and pyrophosphate extraction [START_REF] Mizota | Clay mineralogy and chemistry of soils formed in volcanic material in diverse climatic regions[END_REF] is allophane and not imogolite or proto imogolite. Figure 3 shows the scanning electron microscopy of typical andosol and nitisol samples. Wada [START_REF] Wadak | The distinctive properties of Andosol[END_REF] describes the allophane particles as follows: the unit cell appears as spheroids with diameter between 3 and 5 nm clearly different of the platelet-like microstructure of phyllosilicate clays like halloysite clays which are 100-300 nm in size. Allophane particles are thus one to two orders of magnitude smaller than halloysite particles and can be considered as nanomaterials. These alumino silicates have a high pore volume (~ 70%) and develop large specific surface. From a structural point of view, light and X-rays scattering data demonstrate that allophane aggregates have fractal geometry [START_REF] Woignier | Chlordecone retention in the fractal structure of volcanic clay[END_REF]- [START_REF] Adachi | Application of a scaling law to the analysis of allophane aggregates[END_REF], [START_REF] Chevallier | Fractal structure in natural gels: effect on carbon sequestration in volcanic soils[END_REF]. Thanks to these data we proposed an analogy between these natural nano porous materials and synthetic mineral gels [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF].

B. Allophane: a natural mineral "gel "

A synthetic gel of silica, is a two phase solid-liquid material, amorphous and extremely porous with 70 to 99 % of porosity. Its bulk density is typically in the range 0.1 and 0.5 g.cm -3 and specific surface are in the range of several hundreds of m 2 g -1 [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF]. The gel network can be described as an assembly of fractal clusters (~50 nm) built by the aggregation of small particles (~1-2 nm). The porosity is totally open and spans over the range of the micro to the macroporosity [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF]. Allophane is derived primarily from the in situ weathering of volcanic materials. Allophane particles are obtained by the leaching of amorphous volcanic materials [START_REF] Levard | Structure and distribution of allophanes, imogolite and proto-imogolite in volcanic soils[END_REF] leading to oversaturated solutions which form sols [46]. Karube et al. [START_REF] Adachi | Application of a scaling law to the analysis of allophane aggregates[END_REF], [47]- [48] studied allophane aggregation using light scattering. These authors reported that allophane particles occurred in primary clusters of unit particles like a string of beads several nanometers in size. They also pointed out that secondary clusters could be formed when flocculation occurred and estimated the diameter of the cluster to be 150-200 nm.

Allophane aggregates have physical features very close to that of synthetic gels. As explained above allophane aggregates have fractal geometry and can be described as clusters formed by aggregation of small particles. These "natural gels" have a bulk density close to 0.5 g.cm -1 [START_REF] Wadak | The distinctive properties of Andosol[END_REF] and develop a large specific surface area as high as 700 m 2 g -1 [START_REF] Woignier | Allophane: a natural gel in volcanic soils with interesting environmental properties[END_REF]. The formation of the two kinds of network seems different: the polymerization of silicic acid for the silica gels and a corrosion mechanism for the allophane. In the case of the silica sol-gel process, the organosilicon is hydrolyzed and the silicic acid polymerizes into discrete particles that aggregate into chains and networks while for allophane, the weathering leads to over-saturated solutions which precipitates, the fine particle size and the predominance of glass favor formation of amorphous phases [START_REF] Shoji | Chemical kinetics of weathering in young andisosl from northeastern Japan using 10°C normalized soil age[END_REF]. Nevertheless, in both cases the porous network is obtained by the aggregation of colloidal particles whatever the way of formation of these particles. Moreover, the literature showed that it is possible to synthesize allophane by the sol-gel process [START_REF] Denaix | Structure and affinity towards Cd 2+ , Cu 2+ Pb 2+ of synthetic colloidal amorphous aluminosilicates and their precursor Colloids and Surfaces A[END_REF].

In conclusion, allophane aggregates have morphological features very close to that of synthetic gels: large pore volume and water content, a broad pore size distribution, a high specific surface area and a fractal structure. We will see later that these peculiar "gel" properties are likely involved in the pesticides trapping process.

C. Influence of allophane clay on soil pore properties

Recent work [50] have shown that the microstructure of allophane clay has an influence on the soil porous properties. The pore volume Vp and the specific surface area S determined after supercritical drying were clearly correlated to the allophane contents. The specific surface area could be as high as 180 m 2 .g -1 and pore volume close to 2 cm 3 .g -1 .

Figure 4: Pore size distributions of a set of andosols with different allophane contents %: 0.5 (▲), 5 (•), 12 (-), 23 (■).

The allophane aggregates are highly porous in the range of the mesoporosity [START_REF] Calabi-Floody | Nanoclays from an Andisol: Extraction, properties and carbon stabilization[END_REF]. Figure 4 shows that the pore size distribution of a set of andosols ranges from 5 to 70 nm. Pore size distribution and mesopore volume increase with an increase in allophane content. The pore size distribution also shifts to larger pore sizes with an increase in allophane content. These results are in agreement with hierarchical aggregation, where larger and larger allophane clusters lead to an increasing volume of pores and to a broadening of pore size distribution. These observations are compatible with fractal organization of the allophane aggregates. This kind of aggregation process leads to a microstructure similar to a labyrinth with a high tortuosity and the extent of the labyrinth structure (i.e. the extent of the fractal clusters) increases with the allophane content [50].

D. Results of transfer experiments

Soils derived from volcanic materials have a high potential for agricultural production because of their properties (high exchange capacity, high phosphate absorption and organic matter content). Pest control technology, which was originally developed from temperate climates, was introduced into the tropics without considering the fragility of the ecosystems. Some pesticides whose pollution potential has long been underestimated are now recognized to pose serious human health and environmental problems. For example, chlordecone was extensively used to control pests mainly in banana plantations. CLD persists in the soil and contaminates crops and water resources [START_REF] Dallaire | Cognitive, visual, and motor development of 7-month-old Guadeloupean infants exposed to chlordecone[END_REF], [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF], [START_REF] Brunet | Determination of soil content in chlordecone [organochlorine pesticide] using near infrared reflectance spectroscopy [NIRS[END_REF]- [START_REF] Woignier | Chlordecone retention in the fractal structure of volcanic clay[END_REF], thus leading to human exposure through food. Research has been focused on chlordecone transfer from soils to crops and data shows that the capacity of soils to contaminate plants depends not only on the chlordecone content of the soil, but also on the concerned type of soil [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF], [START_REF] Woignier | Chlordecone retention in the fractal structure of volcanic clay[END_REF], [START_REF] Cabidoche | Contamination of Harvested Organs in Root Crops Grown on Chlordecone-Polluted Soils[END_REF]. Table 1. Mean CLD transfer from soil to crop (mg.kg -1 dry soil) according to the soil type (andosol or nitisol), dispersion index in brackets

The literature [START_REF] Woignier | Chlordecone retention in the fractal structure of volcanic clay[END_REF], [START_REF] Cabidoche | Contamination of Harvested Organs in Root Crops Grown on Chlordecone-Polluted Soils[END_REF], [START_REF] Clostre | Organochlorine [chlordecone] uptake by root vegetables[END_REF], [START_REF] Clostre | Soil Type and Growing Conditions Influence Uptake and Translocation of Organochlorine [Chlordecone] by Cucurbitaceae Species[END_REF] shows that crops contamination shows different CLD contents according to the considered crop, organ and soil type. The most contaminated crops are roots and tubers then cucurbits and leaf vegetables (Table 1). Observing CLD transfer in soil and plant, a general assessment is that andosol tends to retain CLD higher than other soils. For the two crops of this studies, we observed a highly significant difference between the two types of soil (P < 0.0001): transfer from andosol is lesser than from nitisol (table 1) with a ratio close to 2 for sweet potato and 3-5 for carrot.

These results show that: 1) Andosols are generally more contaminated than nitisols and 2) andosols release less pollutant to plants than halloysite soils. In the literature, one explanation for this retention effect is the high organic content of andosols combined with the high affinity of the pesticide for soil organic matter [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF][START_REF] Cabidoche | Contamination of Harvested Organs in Root Crops Grown on Chlordecone-Polluted Soils[END_REF]. The soil-water partition coefficient relative to organic C content (Koc) partly accounts for this process [START_REF] Cabidoche | Long-term pollution by chlordecone of tropical volcanic soils in the French West Indies: A simple leaching model accounts for current residue[END_REF]. However, Koc is a global concept behind which lies a series of explanations (hysteresis; organic matter content; composition of organic matter, etc…), which may make it possible to compare soils, but do not explain the origin of the differences. So, the question remains: what is the direct influence of allophane properties on the retention of pesticides in the soil?

E. Influence of allophane properties on the pesticides retention

To put forth the importance of the allophanic microstructure we fractionated different nitisols and andosols soils in 3 classes: fraction size higher than 200 μm (class a), between 50 and 200 μm (class b) and lower than 50 μm (class c). In figure 5 we followed the CLD concentration in the different a, b and c fractions with increased soil allophane contents. For nitisol samples (0% allophane content on the figure 5), CLD concentrations are showing small differences between b and c fractions, these two fractions are not significantly different from the a fraction. The figure 5 shows that the CLD increases in the c fraction but is quite constant in the a and b fractions; 80% of the CLD content is retained in the c fraction ("allophane fraction") [START_REF] Basile-Doelsch | Mineralogical control of organic carbon dynamics in a volcanic ash soil on La Reunion[END_REF]. These differences are due to the fact that the c size class is rich in allophane aggregates. These results confirm the effect of the allophane content on the CLD retention [50]. As already explained, CLD and similar pesticides tend to bind with soil organic matter due to of their hydrophobicity. Previous study [START_REF] Basile-Doelsch | Mineralogical control of organic carbon dynamics in a volcanic ash soil on La Reunion[END_REF] have shown than the fraction c, the "allophane fraction" is rich in organic matter and we could assume that the large CLD concentration in fraction c is the consequence of a high soil organic content. As expected, figure 6 shows a positive correlation between CLD and SOC in andosols (P<0.01and r 2 =0.68). These results could indicate that the SOC is the explanation for the higher CLD content in andosols.

However, in figure 6 the CLD content of the andosols class c is also plotted versus SOC. These samples have clearly a higher CLD content at a given SOC. The CLD concentration in the andosols c fraction is systematically higher than the general trend. The results stress the fact that the high organic matter content in andosols is not the only effect and we propose that the allophane tortuous microstructure plays also a role in the CLD retention. The structural properties of the mineral matrix with which pollutants are associated may influence pollutant bioavailability and retention in soils [START_REF] Peters | Differences in p,p′-DDE bioaccumulation from compost and soil by the plants Cucurbita pepo and Cucurbita maxima and the earthworms Eisenia fetida and Lumbricus terrestris[END_REF]- [START_REF] Liu | Mechanisms Associated with the High Adsorption of Dibenzo-p-dioxin from Water by Smectite Clays[END_REF]. The decrease in pesticide bioavailability has often been explained by slow diffusion in small pores and can be physically explained as the result of the pollutant entrapment [57]- [59]. Data revealed greater physical protection and adsorption of persistent pesticides in nanosized structures. We have seen in section B that the allophane aggregate microstructure resembles a labyrinth at the nm scale. This description suggests that the accessibility in the allophane microstructure decreases because of this nano materials microstructure. The combination of pore size and network tortuosity may be a key feature in the inaccessibility of pesticides to microorganisms.

F. High water content and collapse during drying

Another distinctive property of allophane clay is involved in the high trapping of pesticides: the high water content and the associated pore collapse during drying. Figure 7 shows the large increase of the water content versus the allophane content, the results expresses that the allophane clay is like a sponge with water content as high as 180%.

Figure7: Water content in volcanic soils versus allophane content Moreover, figure 8 compares the drying curve between an andosol and a nitisol and shows clearly a specific volume and water content loss 4-5 times larger than for the nitisol. Figure 8 shows that during a classical drying, and because of the large water content, the andosol exhibits an important shrinkage which modifies the soil physical properties. This shrinkage comes from capillary stresses during the drying, in the pores and from the large compliance of the porous solids. This behaviour is typical of the mineral gels like silica gels [START_REF] Brinker | Sol-gel science: the physics and chemistry of sol-gel processing[END_REF]. It must be noted that if the shrinkage is quite reversible for nitisols, in the case of andosols, the shrinkage is irreversible and a large part of the pore volume is lost. To quantify the amplitude of the volume irreversibly lost after drying we have rewetted the soils samples at room temperature during two days and measured the final water content. The water content after drying and rewetting was compared with the initial water content and allowed to calculate the water content lost and consequently the lost pore volume (Fig. 9).

Figure 9 shows that the irreversible volume shrinkage is related to the allophane content in volcanic soil and can be higher than 50%. These results show that a large part of the volume collapse should be attributed to the shrinkage of the allophane aggregates. The loss of this pore volume affects the pesticide bioavailabilty. One would expect the pore structure to have an influence on the mechanism of pesticide diffusion and sorption. In this case, changes in the soil pore volume should affect the flow and diffusion processes and pesticide transport thereby changing pesticide bioavailability. The main structural consequence of the shrinkage is the progressive contraction of the allophane aggregates associated with closure of the mesopores. Consequently, the capillary strain collapses the clay microstructure and hence affect the transport properties within the clay aggregates. As a results, possible exchange and chemical reactions able to degrade the pesticides should be poor. This reduced accessibility might partially account for the containment process. The chemical affinity for organic matter associated with low transport properties inside allophane aggregates and the collapse of mesopores, play a role in the entrapment of pesticides in andosols.

IV. CONCLUSION

A wide range of pesticides have been used for tropical volcanic soils and the general agreement is that the pesticide sorption behavior is most pronounced in andosols. In this work we have confirmed that andosols retain and trap more pesticides than other kinds of soils. The study confirms that there are clear correlations between the soil chlordecone content, soil organic content and allophane microstructure (natural gel features). The bioavailability of each pesticide is related to the intrinsic physical-chemical properties of the pollutant (water solubility, chemical affinity with soil organic matter, and poor biodegradability). We show that the structural properties of the mineral matrix (nano size sponge) with which the pollutants are associated may also influence pollutant bioavailability and retention in soils. As mentioned above the more accepted explanation for the retention effect reported in the literature is the high affinity of the pesticide for soil organic carbon. However, physical parameters like poor transport properties inside clay aggregates and porosity closing could participate to the overall pesticide sequestration. The nano porous structure of the porosity slows down access to trapped pesticides by enzymes and microorganisms. Moreover, after a strong drying, the andosols porosity irreversibly collapses. So, during the cycle of dry and wet seasons the allophane aggregates could collapse because of the drying and rewetting process. This shrinkage, which increases with the allophane content, can participate to the whole pesticides trapping in andosols by the closing of the porosity. Future research for pollutants degrading microorganisms and other bioremediation tools to clean up polluted soils should take into account these different aspects of volcanic soils, large soil organic content, poor CLD accessibility and collapse of the meso porosity.

ACKNOWLEDGMENT

Funding was provided by the French Chlordecone National Plan ("JAFA" project), the French National Research Agency ("Chlordexco" project) and the French Ministry for Overseas development (MOM).

.

Figure1

  Figure1: X-rays typical diffractogram of volcanic soils samples. Counts versus lattice spacing d (Å)

Figure 2 :

 2 Figure 2: Infra-red typical spectra of volcanic soils. Transmittance in arbitrary units versus wavelength in cm -1 .

Figure 3

 3 Figure 3 Scanning electron Microscopy of volcanic soils: nitisols (a) and andosols (b)

Figure 5 :

 5 Figure 5: Evolution of the chlordécone content in the different classes of size versus allophane content. a (-), b (◊), and c (■).

Figure 6 :

 6 Figure 6: CLD concentration versus SOC in volcanic soil samples (▲) and in the c size class of andosols (■).

Figure 8 :

 8 Figure 8: Drying curves of typical andosol and nitisol

Figure 9 :

 9 Figure 9: Irreversible volume shrinkage versus allophane content